Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

draft implementation of how to adden leung dust emissions #43

Draft
wants to merge 1 commit into
base: noresm_develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions src/aero_model.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down
187 changes: 58 additions & 129 deletions src/oslo_aero_dust.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
!===============================================================================
Expand Down Expand Up @@ -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

!===============================================================================
Expand All @@ -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
Expand All @@ -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:
Expand All @@ -123,122 +129,45 @@ 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
! same formulation as MAM now and tune later
! 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