From 830ae29247047ace5f127b995ed433684fde458d Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 27 Nov 2024 19:50:33 +0100 Subject: [PATCH] draft implementation of how to adden leung dust emissions --- src/aero_model.F90 | 6 +- src/oslo_aero_dust.F90 | 187 +++++++++++++---------------------------- 2 files changed, 60 insertions(+), 133 deletions(-) diff --git a/src/aero_model.F90 b/src/aero_model.F90 index d757bb4..25dafda 100644 --- a/src/aero_model.F90 +++ b/src/aero_model.F90 @@ -50,7 +50,7 @@ module aero_model use oslo_aero_condtend, only: N_COND_VAP, COND_VAP_ORG_SV, COND_VAP_ORG_LV, COND_VAP_H2SO4 use oslo_aero_condtend, only: initializeCondensation, condtend use oslo_aero_seasalt, only: oslo_aero_seasalt_init, oslo_aero_seasalt_emis, seasalt_active - use oslo_aero_dust, only: oslo_aero_dust_init, oslo_aero_dust_emis, dust_active + use oslo_aero_dust, only: oslo_aero_dust_init, oslo_aero_dust_emis use oslo_aero_ocean, only: oslo_aero_ocean_init, oslo_aero_dms_emis use oslo_aero_share, only: getNumberofTracersInMode, getCloudTracerIndexDirect, getCloudTracerName use oslo_aero_share, only: getCloudTracerName, getTracerIndex, aero_register @@ -592,9 +592,7 @@ subroutine aero_model_emissions( state, cam_in ) integer :: icol integer :: pndx_fdms ! DMS surface flux physics index - if (dust_active) then - call oslo_aero_dust_emis( state%lchnk, state%ncol, cam_in%dstflx, cam_in%cflx) - endif + call oslo_aero_dust_emis( state%lchnk, state%ncol, cam_in%dstflx, cam_in%cflx) if (seasalt_active) then call oslo_aero_seasalt_emis(state%ncol, state%lchnk, & diff --git a/src/oslo_aero_dust.F90 b/src/oslo_aero_dust.F90 index 44d09a7..edfb603 100644 --- a/src/oslo_aero_dust.F90 +++ b/src/oslo_aero_dust.F90 @@ -8,16 +8,12 @@ module oslo_aero_dust use spmd_utils, only: mpicom, mstrid=>masterprocid, masterproc use spmd_utils, only: mpi_logical, mpi_real8, mpi_character, mpi_integer, mpi_success use namelist_utils, only: find_group_name - use ppgrid, only: pcols, begchunk, endchunk - use phys_grid, only: get_ncols_p, get_rlat_all_p, get_rlon_all_p + use ppgrid, only: pcols use constituents, only: cnst_name, pcnst - use interpolate_data, only: lininterp_init, lininterp, lininterp_finish, interp_type - use mo_constants, only: pi, d2r use cam_logfile, only: iulog use cam_abortutils, only: endrun - use cam_pio_utils, only: cam_pio_openfile - use ioFileMod, only: getfil - use pio, only: file_desc_t,pio_inq_dimid,pio_inq_dimlen,pio_get_var,pio_inq_varid, PIO_NOWRITE + use shr_dust_emis_mod,only: is_dust_emis_zender, is_zender_soil_erod_from_atm, shr_dust_emis_readnl + use soil_erod_mod, only: soil_erod_init ! use oslo_aero_share, only: l_dst_a2, l_dst_a3 @@ -29,26 +25,16 @@ module oslo_aero_dust public :: oslo_aero_dust_init public :: oslo_aero_dust_emis - ! private routines (previously in soil_erod_mod in CAM) - private :: soil_erod_init - character(len=6), public :: dust_names(10) - integer , parameter :: numberOfDustModes = 2 !define in oslo_aero_share? + integer , parameter :: numberOfDustModes = 2 real(r8), parameter :: emis_fraction_in_mode(numberOfDustModes) = (/0.13_r8, 0.87_r8 /) integer :: tracerMap(numberOfDustModes) = (/-99, -99/) !index of dust tracers in the modes - integer , parameter, public :: dust_nbin = numberOfDustModes - !Related to soil erodibility - real(r8) :: dust_emis_fact = -1.e36_r8 ! tuning parameter for dust emissions + real(r8) :: dust_emis_fact = 0._r8 ! tuning parameter for dust emissions character(len=cl) :: soil_erod_file = 'soil_erod_file' ! full pathname for soil erodibility dataset - logical, parameter, public :: dust_active = .TRUE. - - real(r8), allocatable :: soil_erodibility(:,:) ! soil erodibility factor - real(r8) :: soil_erod_fact ! tuning parameter for dust emissions - !=============================================================================== contains !=============================================================================== @@ -76,12 +62,30 @@ subroutine oslo_aero_dust_readnl(nlfile) end if close(unitn) end if + ! Broadcast namelist variables call mpi_bcast(dust_emis_fact, 1, mpi_real8, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(subname//" mpi_bcast: dust_emis_fact") call mpi_bcast(soil_erod_file, len(soil_erod_file), mpi_character, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(subname//" mpi_bcast: soil_erod_file") + call shr_dust_emis_readnl(mpicom, 'drv_flds_in') + + if ((soil_erod_file /= 'none') .and. (.not.is_zender_soil_erod_from_atm())) then + call endrun(subname//': should not specify soil_erod_file if Zender soil erosion is not in CAM') + end if + + if (masterproc) then + if (is_dust_emis_zender()) then + write(iulog,*) subname,': Zender_2003 dust emission method is being used.' + end if + if (is_zender_soil_erod_from_atm()) then + write(iulog,*) subname,': Zender soil erod file is handled in atm' + write(iulog,*) subname,': soil_erod_file = ',trim(soil_erod_file) + write(iulog,*) subname,': dust_emis_fact = ',dust_emis_fact + end if + end if + end subroutine oslo_aero_dust_readnl !=============================================================================== @@ -90,7 +94,9 @@ subroutine oslo_aero_dust_init() ! local variables integer :: imode - call soil_erod_init( dust_emis_fact, soil_erod_file ) + if (is_zender_soil_erod_from_atm()) then + call soil_erod_init( dust_emis_fact, soil_erod_file ) + end if ! Set module variables tracerMap(1) = l_dst_a2 @@ -109,7 +115,7 @@ subroutine oslo_aero_dust_emis(lchnk, ncol, dstflx, cflx) !----------------------------------------------------------------------- ! Purpose: Interface to emission of all dusts. ! Notice that the mobilization is calculated in the land model and - ! the soil erodibility factor is applied here. + ! the soil erodibility factor is applied here if the zender scheme is used !----------------------------------------------------------------------- ! Arguments: @@ -123,19 +129,6 @@ subroutine oslo_aero_dust_emis(lchnk, ncol, dstflx, cflx) real(r8) :: soil_erod_tmp(pcols) real(r8) :: totalEmissionFlux(pcols) - ! Filter away unreasonable values for soil erodibility - ! (using low values e.g. gives emissions in greenland..) - where(soil_erodibility(:,lchnk) < 0.1_r8) - soil_erod_tmp(:)=0.0_r8 - elsewhere - soil_erod_tmp(:)=soil_erodibility(:,lchnk) - end where - - totalEmissionFlux(:) = 0.0_r8 - do icol=1,ncol - totalEmissionFlux(icol) = totalEmissionFlux(icol) + sum(dstflx(icol,:)) - end do - ! Note that following CESM use of "dust_emis_fact", the emissions are ! scaled by the INVERSE of the factor!! ! There is another random scale factor of 1.15 there. Adapting the exact @@ -143,102 +136,38 @@ subroutine oslo_aero_dust_emis(lchnk, ncol, dstflx, cflx) ! As of NE-380: Oslo dust emissions are 2/3 of CAM emissions ! gives better AOD close to dust sources - do imode = 1,numberOfDustModes - cflx(:ncol, tracerMap(imode)) = -1.0_r8*emis_fraction_in_mode(imode) & - *totalEmissionFlux(:ncol)*soil_erod_tmp(:ncol)/(dust_emis_fact)*1.15_r8 - end do - - end subroutine oslo_aero_dust_emis - - !============================================================================= - subroutine soil_erod_init( dust_emis_fact, soil_erod_file ) - - ! arguments - real(r8), intent(in) :: dust_emis_fact - character(len=*), intent(in) :: soil_erod_file - - ! localvaraibles - real(r8), allocatable :: soil_erodibility_in(:,:) - real(r8), allocatable :: dst_lons(:) - real(r8), allocatable :: dst_lats(:) - character(len=cl) :: infile - integer :: did, vid, nlat, nlon - type(file_desc_t) :: ncid - type(interp_type) :: lon_wgts, lat_wgts - real(r8) :: to_lats(pcols), to_lons(pcols) - integer :: c, ncols, ierr - real(r8), parameter :: zero=0._r8 - real(r8), parameter :: twopi=2._r8*pi - - soil_erod_fact = dust_emis_fact - - ! Summary to log file - if (masterproc) then - write(iulog,*) 'soil_erod_mod: soil erodibility dataset: ', trim(soil_erod_file) - write(iulog,*) 'soil_erod_mod: soil_erod_fact = ', soil_erod_fact + if (is_zender_soil_erod_from_atm()) then + ! Filter away unreasonable values for soil erodibility + ! (using low values e.g. gives emissions in greenland..) + where(soil_erodibility(:,lchnk) < 0.1_r8) + soil_erod_tmp(:)=0.0_r8 + elsewhere + soil_erod_tmp(:)=soil_erodibility(:,lchnk) + end where + + totalEmissionFlux(:) = 0.0_r8 + do icol=1,ncol + totalEmissionFlux(icol) = totalEmissionFlux(icol) + sum(dstflx(icol,:)) + end do + + do imode = 1,numberOfDustModes + cflx(:ncol, tracerMap(imode)) = -1.0_r8*emis_fraction_in_mode(imode) & + *totalEmissionFlux(:ncol)*soil_erod_tmp(:ncol)/(dust_emis_fact)*1.15_r8 + end do + + else ! Leung emissions + + totalEmissionFlux(:) = 0.0_r8 + do icol=1,ncol + totalEmissionFlux(icol) = totalEmissionFlux(icol) + sum(dstflx(icol,:)) + end do + + do imode = 1,numberOfDustModes + cflx(:ncol, tracerMap(imode)) = -1.0_r8*emis_fraction_in_mode(imode) & + *totalEmissionFlux(:ncol)/(dust_emis_fact) + end do end if - ! read in soil erodibility factors, similar to Zender's boundary conditions - - ! Get file name. - call getfil(soil_erod_file, infile, 0) - call cam_pio_openfile (ncid, trim(infile), PIO_NOWRITE) - - ! Get input data resolution. - ierr = pio_inq_dimid( ncid, 'lon', did ) - ierr = pio_inq_dimlen( ncid, did, nlon ) - - ierr = pio_inq_dimid( ncid, 'lat', did ) - ierr = pio_inq_dimlen( ncid, did, nlat ) - - allocate(dst_lons(nlon)) - allocate(dst_lats(nlat)) - allocate(soil_erodibility_in(nlon,nlat)) - - ierr = pio_inq_varid( ncid, 'lon', vid ) - ierr = pio_get_var( ncid, vid, dst_lons ) - - ierr = pio_inq_varid( ncid, 'lat', vid ) - ierr = pio_get_var( ncid, vid, dst_lats ) - - ierr = pio_inq_varid( ncid, 'mbl_bsn_fct_geo', vid ) - ierr = pio_get_var( ncid, vid, soil_erodibility_in ) - - ! convert to radians and setup regridding - dst_lats(:) = d2r * dst_lats(:) - dst_lons(:) = d2r * dst_lons(:) - - allocate( soil_erodibility(pcols,begchunk:endchunk), stat=ierr ) - if( ierr /= 0 ) then - write(iulog,*) 'soil_erod_init: failed to allocate soil_erodibility_in, ierr = ',ierr - call endrun('soil_erod_init: failed to allocate soil_erodibility_in') - end if - - soil_erodibility(:,:)=0._r8 - - ! regrid - do c=begchunk,endchunk - ncols = get_ncols_p(c) - call get_rlat_all_p(c, pcols, to_lats) - call get_rlon_all_p(c, pcols, to_lons) - - call lininterp_init(dst_lons, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi) - call lininterp_init(dst_lats, nlat, to_lats, ncols, 1, lat_wgts) - - call lininterp(soil_erodibility_in(:,:), nlon, nlat, soil_erodibility(:,c), ncols, lon_wgts, lat_wgts) - - call lininterp_finish(lat_wgts) - call lininterp_finish(lon_wgts) - end do - deallocate( soil_erodibility_in, stat=ierr ) - if( ierr /= 0 ) then - write(iulog,*) 'soil_erod_init: failed to deallocate soil_erodibility_in, ierr = ',ierr - call endrun('soil_erod_init: failed to deallocate soil_erodibility_in') - end if - - deallocate( dst_lats ) - deallocate( dst_lons ) - - end subroutine soil_erod_init + end subroutine oslo_aero_dust_emis end module oslo_aero_dust