From 6273be20fb50ae819c5399bec038bf1a8d97d70e Mon Sep 17 00:00:00 2001 From: Alexander Date: Wed, 24 Mar 2021 17:51:50 +0300 Subject: [PATCH] removing dxy from GOCART dust emission schemes (#1431) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TYPE: enhancement KEYWORDS: GOCART, dust emission, dust_opt, WRF-Chem SOURCE: Alexander Ukhov (KAUST) DESCRIPTION OF CHANGES: This PR affects dust emission schemes (dust_opt=1,3,4) coupled with GOCART aerosol module. Current implementation of dust_opt=1,3 calculates air mass in the grid box using cell area and height, and air density, however a more elegant solution is now implemented (by analogy with dust_opt=4). This solution only requires height of the lowest model layer. Another issue have been also resolved when using dust_opt=1. Due to incorrect indexing in the tile loop dust emissions diagnostics (variables EDUST1…5) were zero near the boundaries of the domain. Each emission scheme (dust_opt=1,3,4) reported instantaneous dust emissions diagnostics (variables EDUST1…5) in different units: dust_opt=1 [kg/cell/dt], dust_opt=3 [gram/m^2/s], dust_opt=4 [ug/m^2/s]. For user convenience this discrepancy also has been fixed. now EDUST1..5 variables store accumulated dust emissions [kg/m^2] for dust_opt=1,3,4. For example, this can be useful when is necessary to estimate dust mass balance. In this case accumulated dust emission needs to be multiplied by cell area (m^2), see Appendix of [1]. Code clean-up has been also done in dust_opt=1,3. LIST OF MODIFIED FILES: M Registry/registry.chem M chem/emissions_driver.F M chem/module_gocart_dust.F M chem/module_gocart_dust_afwa.F M chem/module_qf03.F M chem/module_uoc_dust.F TESTS CONDUCTED: To demonstrate correctness of the implemented changes a numerical simulations of the dust over the Middle East (chem_opt=301 with dust_opt=1,3,4, 100km resolution) with and without changes has been run. --- Registry/registry.chem | 14 +- chem/emissions_driver.F | 19 +-- chem/module_gocart_dust.F | 238 ++++++++++----------------------- chem/module_gocart_dust_afwa.F | 106 +++++++-------- chem/module_qf03.F | 8 +- chem/module_uoc_dust.F | 25 ++-- 6 files changed, 160 insertions(+), 250 deletions(-) diff --git a/Registry/registry.chem b/Registry/registry.chem index 9b4abd6141..1089598679 100644 --- a/Registry/registry.chem +++ b/Registry/registry.chem @@ -314,11 +314,11 @@ state real e_hum i+jf emis_ant 1 Z i5 "E_H # dust and seas emission arrays state real - i{dust}jf emis_dust - - - - "Dust Emissions" "" -state real edust1 i{dust}jf emis_dust 1 Z - "EDUST1" "DUST emissions bin1" "" -state real edust2 i{dust}jf emis_dust 1 Z - "EDUST2" "DUST emissions bin2" "" -state real edust3 i{dust}jf emis_dust 1 Z - "EDUST3" "DUST emissions bin3" "" -state real edust4 i{dust}jf emis_dust 1 Z - "EDUST4" "DUST emissions bin4" "" -state real edust5 i{dust}jf emis_dust 1 Z - "EDUST5" "DUST emissions bin5" "" +state real edust1 i{dust}jf emis_dust 1 Z - "EDUST1" "Accumulated DUST emissions bin1" "kg/m2" +state real edust2 i{dust}jf emis_dust 1 Z - "EDUST2" "Accumulated DUST emissions bin2" "kg/m2" +state real edust3 i{dust}jf emis_dust 1 Z - "EDUST3" "Accumulated DUST emissions bin3" "kg/m2" +state real edust4 i{dust}jf emis_dust 1 Z - "EDUST4" "Accumulated DUST emissions bin4" "kg/m2" +state real edust5 i{dust}jf emis_dust 1 Z - "EDUST5" "Accumulated DUST emissions bin5" "kg/m2" state real - i{dust}jf emis_seas - - - - "Sea-Salt Emissions" "" state real eseas1 i{dust}jf emis_seas 1 Z - "ESEAS1" "Sea-Salt emissions bin1 " "" state real eseas2 i{dust}jf emis_seas 1 Z - "ESEAS2" "Sea-Salt emissions bin2 " "" @@ -684,8 +684,8 @@ state real sandfrac ij misc 1 - i01r "SAND state real clayfrac_nga ij misc 1 - i01r "CLAYFRAC_NGA" "Clay fraction in each grid cell (0-1)" "none" state real sandfrac_nga ij misc 1 - i01r "SANDFRAC_NGA" "Sand fraction in each grid cell (0-1)" "none" state real afwa_dustloft ij misc 1 - h02 "AFWA_DUSTLOFT" "AFWA Diagnostic dust lofting potential (U10-U10t)" "m s^-1" -state real tot_dust ikj misc 1 - h02 "TOT_DUST" "Total dust concentration (0.1-20 um)" "ug m^-3" -state real tot_edust ij misc 1 - h02 "TOT_EDUST" "Total dust emission flux (0.1-20 um)" "g m^-2 s^-1" +state real tot_dust ikj misc 1 - h02 "TOT_DUST" "Total dust concentration (0.2-20 um)" "ug m^-3" +state real tot_edust ij misc 1 - h02 "TOT_EDUST" "Total accumulated dust emission (0.2-20 um)" "kg m^-2" state real vis_dust ikj misc 1 - h02 "VIS_DUST" "Visibility due to dust only" "m" # These 3D arrays are output from the SOA module for different purposes diff --git a/chem/emissions_driver.F b/chem/emissions_driver.F index e2b9431eda..6648e01c50 100644 --- a/chem/emissions_driver.F +++ b/chem/emissions_driver.F @@ -5,6 +5,9 @@ ! Saulo Freitas (CPTEC), and Georg Grell ! ! +! A. Ukhov, 11 March 2021, remove unused parameters in gocart_dust_driver(), +! gocart_dust_afwa_driver(), and uoc_dust_driver() subroutines. + MODULE module_emissions_driver IMPLICIT NONE CONTAINS @@ -721,18 +724,18 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & dust_select: SELECT CASE(config_flags%dust_opt) CASE (DUSTGOCART) CALL wrf_debug(15,'Gocart dust emissions') - call gocart_dust_driver(ktau,dtstep,config_flags,julday,alt,t_phy,moist,u_phy, & - v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,dustin, & - ivgtyp,isltyp,vegfra,xland,xlat,xlong,gsw,dx,g,emis_dust, & + call gocart_dust_driver(dtstep,config_flags,alt,t_phy,u_phy, & + v_phy,chem,rho_phy,dz8w,smois,u10,v10,erod,dustin, & + isltyp,xland,g,emis_dust, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) CASE (DUSTGOCARTAFWA) CALL wrf_debug(15,'AFWA modified Gocart dust emissions') - call gocart_dust_afwa_driver(ktau,dtstep,config_flags,julday,alt,t_phy,moist,u_phy, & - v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,erod_dri,dustin,snowh,zs, & - ivgtyp,isltyp,vegfra,lai_vegmask,xland,xlat,xlong,gsw,dx,g,emis_dust, & - ust,znt,clayfrac,sandfrac,clayfrac_nga,sandfrac_nga,afwa_dustloft, &!EDH + call gocart_dust_afwa_driver(dtstep,config_flags,alt, & + chem,rho_phy,smois,u10,v10,p8w,dz8w,erod,erod_dri,dustin,snowh, & + isltyp,vegfra,lai_vegmask,xland,dx,g,emis_dust, & + ust,znt,clayfrac,sandfrac,clayfrac_nga,sandfrac_nga,afwa_dustloft,&!EDH tot_dust,tot_edust,vis_dust,dust_alpha,dust_gamma,dust_smtune,dust_ustune, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -750,7 +753,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & CASE DEFAULT imod = 2 END SELECT scheme_select - call uoc_dust_driver (ktau,dtstep,config_flags, & + call uoc_dust_driver (dtstep,config_flags, & chem,rho_phy,dz8w,smois,ust, isltyp,vegfra,g,emis_dust, & ust_t, imod, rough_cor, smois_cor, soil_top_cat, erod, & ids,ide, jds,jde, kds,kde, & diff --git a/chem/module_gocart_dust.F b/chem/module_gocart_dust.F index 5138600275..1049227c67 100644 --- a/chem/module_gocart_dust.F +++ b/chem/module_gocart_dust.F @@ -1,77 +1,59 @@ MODULE GOCART_DUST - + +! A. Ukhov, 11 March 2021, Now "emis_dust" is accumulated dust +! emission (kg/m2). Before was instantenious flux (kg/cell). +! Bug fix in the loop over cells: Cells near domain boundaries +! were not processed. Code cleanup, remove unused variables. USE module_data_gocart_dust CONTAINS - subroutine gocart_dust_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy, & - v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,dustin, & - ivgtyp,isltyp,vegfra,xland,xlat,xlong,gsw,dx,g,emis_dust, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + subroutine gocart_dust_driver(dt,config_flags,alt,t_phy,u_phy, & + v_phy,chem,rho_phy,dz8w,smois,u10,v10,erod,dustin, & + isltyp,xland,g,emis_dust, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) USE module_configure USE module_state_description - USE module_model_constants, ONLY: mwdry IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags - INTEGER, INTENT(IN ) :: julday, ktau, & - ids,ide, jds,jde, kds,kde, & + INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte - INTEGER,DIMENSION( ims:ime , jms:jme ) , & - INTENT(IN ) :: & - ivgtyp, & - isltyp - REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & - INTENT(IN ) :: moist + + INTEGER,DIMENSION( ims:ime , jms:jme ), & + INTENT(IN ) :: isltyp REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & - INTENT(INOUT ) :: chem + INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,& - INTENT(INOUT ) :: & - emis_dust + INTENT(INOUT ) :: emis_dust REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , & INTENT(INOUT) :: smois - REAL, DIMENSION( ims:ime , jms:jme, 3 ) , & + REAL, DIMENSION( ims:ime , jms:jme, 3 ) , & INTENT(IN ) :: erod - REAL, DIMENSION( ims:ime , jms:jme, 5 ) , & + REAL, DIMENSION( ims:ime , jms:jme, 5 ) , & INTENT(INout ) :: dustin REAL, DIMENSION( ims:ime , jms:jme ) , & - INTENT(IN ) :: & - u10, & - v10, & - gsw, & - vegfra, & - xland, & - xlat, & - xlong + INTENT(IN ) :: u10,v10,xland REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & - INTENT(IN ) :: & - alt, & + INTENT(IN ) :: alt, & t_phy, & - dz8w,p8w, & + dz8w, & u_phy,v_phy,rho_phy - REAL, INTENT(IN ) :: dt,dx,g + REAL, INTENT(IN ) :: dt,g ! ! local variables ! - integer :: nmx,i,j,k,ndt,imx,jmx,lmx + integer :: nmx,i,j,k,imx,jmx,lmx integer,dimension (1,1) :: ilwi real*8, DIMENSION (1,1,3,1) :: erodin real*8, DIMENSION (5) :: tc,bems - real*8, dimension (1,1) :: w10m,gwet,airden,airmas - real*8, dimension (1) :: dxy + real*8, dimension (1,1) :: w10m,gwet,airden + real :: dz_lowest real*8 conver,converi - real dttt - integer ibegc,jbegc,iendc,jendc - jbegc=max(jts,jds+3) - jendc=min(jte,jde-4) - ibegc=max(its,ids+3) - iendc=min(ite,ide-4) -! conver=1.e-9*mwdry -! converi=1.e9/mwdry conver=1.e-9 converi=1.e9 ! @@ -82,9 +64,9 @@ subroutine gocart_dust_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy, lmx=1 nmx=5 k=kts - do j=jbegc,jendc - do i=ibegc,iendc -! + DO j=jts,jte + DO i=its,ite + ! no dust over water!!! ! if(xland(i,j).lt.1.5)then @@ -99,62 +81,43 @@ subroutine gocart_dust_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy, tc(5)=chem(i,kts,j,p_dust_5)*conver endif w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) - airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g ! ! don't trust the u10,v10 values, is model layers are very thin near surface ! if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j)) - erodin(1,1,1,1)=erod(i,j,1)!/dx/dx - erodin(1,1,2,1)=erod(i,j,2)!/dx/dx - erodin(1,1,3,1)=erod(i,j,3)!/dx/dx + erodin(1,1,1,1)=erod(i,j,1) + erodin(1,1,2,1)=erod(i,j,2) + erodin(1,1,3,1)=erod(i,j,3) ! ! volumetric soil moisture over porosity ! gwet(1,1)=smois(i,1,j)/porosity(isltyp(i,j)) - ndt=ifix(dt) airden(1,1)=rho_phy(i,kts,j) - dxy(1)=dx*dx -! if(erod(i,j,1).gt.0.)write(0,*)'er1=',p_dust_1,num_chem,erod(i,j,1),tc(2) -! if(erod(i,j,1).gt.0.)write(0,*)'er1=',dt,dxy(1),u10(i,j),w10m(1,1) -! erodin(1,1,1,1)= 0.149748762553862 -! erodin(1,1,2,1)= 7.487438878070708E-002 -! erodin(1,1,3,1)= 7.487438878070708E-002 -! ilwi(1,1)= 1 -! dxy(1)= 54585850453.7552 -! w10m(1,1)= 10.6305338763678 -! gwet(1,1)= 9.136307984590530E-002 -! airden(1,1)= 1.16423276395132 -! airmas(1,1)= 8114017750938.79 -! tc (1) = 1.000000000000000D-030 -! tc (2) = 1.000000000000000d-030 -! tc (3) = 1.000000000000000d-030 -! tc(4) = 1.000000000000000d-030 -! tc(5) = 1.000000000000000d-030 -! dttt=3600. + dz_lowest = dz8w(i,1,j) + call source_du( imx,jmx,lmx,nmx, dt, tc, & - erodin, ilwi, dxy, w10m, gwet, airden, airmas, & - bems,config_flags%start_month,g) -! write(0,*)tc(1) -! write(0,*)tc(2) -! write(0,*)tc(3) -! write(0,*)tc(4) -! write(0,*)tc(5) -! if(erod(i,j,1).gt.0.)write(0,*)'er2=',i,j,erod(i,j,1),tc(2) - if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then - dustin(i,j,1:5)=tc(1:5)*converi - else - chem(i,kts,j,p_dust_1)=tc(1)*converi - chem(i,kts,j,p_dust_2)=tc(2)*converi - chem(i,kts,j,p_dust_3)=tc(3)*converi - chem(i,kts,j,p_dust_4)=tc(4)*converi - chem(i,kts,j,p_dust_5)=tc(5)*converi - endif -! for output diagnostics - emis_dust(i,1,j,p_edust1)=bems(1) - emis_dust(i,1,j,p_edust2)=bems(2) - emis_dust(i,1,j,p_edust3)=bems(3) - emis_dust(i,1,j,p_edust4)=bems(4) - emis_dust(i,1,j,p_edust5)=bems(5) + erodin, ilwi, w10m, gwet, airden, & + dz_lowest,bems,config_flags%start_month,g) + + if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then + dustin(i,j,1:5)=tc(1:5)*converi + else + chem(i,kts,j,p_dust_1)=tc(1)*converi ! tc(1...5) is (kg/kg), p_dust_1...5 (ug/kg) + chem(i,kts,j,p_dust_2)=tc(2)*converi + chem(i,kts,j,p_dust_3)=tc(3)*converi + chem(i,kts,j,p_dust_4)=tc(4)*converi + chem(i,kts,j,p_dust_5)=tc(5)*converi + endif + + ! A. Ukhov + ! for output diagnostics + ! bems (kg/m2) per dt + ! p_edust1...5 is accumulated dust emission (kg/m2) + emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1) + emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2) + emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3) + emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4) + emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5) endif enddo enddo @@ -164,8 +127,8 @@ end subroutine gocart_dust_driver SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, & - erod, ilwi, dxy, w10m, gwet, airden, airmas, & - bems,month,g0) + erod, ilwi, w10m, gwet, airden, & + dz_lowest,bems,month,g0) ! **************************************************************************** ! * Evaluate the source of each dust particles size classes (kg/m3) @@ -173,87 +136,38 @@ SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, & ! * Input: ! * EROD Fraction of erodible grid cell (-) ! * for 1: Sand, 2: Silt, 3: Clay -! * DUSTDEN Dust density (kg/m3) -! * DXY Surface of each grid cell (m2) ! * AIRVOL Volume occupy by each grid boxes (m3) -! * NDT1 Time step (s) +! * DT1 Time step (s) ! * W10m Velocity at the anemometer level (10meters) (m/s) ! * u_tresh Threshold velocity for particule uplifting (m/s) ! * CH_dust Constant to fudge the total emission of dust (s2/m2) +! * dz_lowest heigth of the lowest layer (m) ! * ! * Output: -! * DSRC Source of each dust type (kg/timestep/cell) +! * DSRC Source of each dust type (kg/timestep/m2) +! * BEMS Source of each dust type (kg/timestep/m2) ! * ! * Working: ! * SRC Potential source (kg/m/timestep/cell) ! * ! **************************************************************************** -! USE module_data_gocart -! USE module_data_gocart_dust - - - - INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx + INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx REAL*8, INTENT(IN) :: erod(imx,jmx,ndcls,ndsrc) - INTEGER, INTENT(IN) :: ilwi(imx,jmx),month + INTEGER, INTENT(IN) :: ilwi(imx,jmx),month REAL*8, INTENT(IN) :: w10m(imx,jmx), gwet(imx,jmx) - REAL*8, INTENT(IN) :: dxy(jmx) - REAL*8, INTENT(IN) :: airden(imx,jmx,lmx), airmas(imx,jmx,lmx) + REAL*8, INTENT(IN) :: airden(imx,jmx,lmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(OUT) :: bems(imx,jmx,nmx) + REAL, INTENT(IN ) :: dz_lowest REAL*8 :: den(nmx), diam(nmx) - REAL*8 :: tsrc, u_ts0, cw, u_ts, dsrc, srce + REAL*8 :: u_ts0, u_ts, dsrc, srce REAL, intent(in) :: g0 REAL :: rhoa, g,dt1 INTEGER :: i, j, n, m, k - - REAL*8 :: tcmw(nmx), ar(nmx), tcvv(nmx) - REAL*8 :: ar_wetdep(nmx), kc(nmx) - CHARACTER(LEN=20) :: tcname(nmx), tcunits(nmx) - LOGICAL :: aerosol(nmx) - - -! REAL*8 :: tc1(imx,jmx,lmx,nmx) -! REAL*8, TARGET :: tcms(imx,jmx,lmx,nmx) ! tracer mass (kg; kgS for sulfur case) -! REAL*8, TARGET :: tcgm(imx,jmx,lmx,nmx) ! g/m3 - - !----------------------------------------------------------------------- - ! sea salt specific - !----------------------------------------------------------------------- -! REAL*8, DIMENSION(nmx) :: ssaltden, ssaltreff, ra, rb -! REAL*8 :: ch_ss(nmx,12) - - !----------------------------------------------------------------------- - ! emissions (input) - !----------------------------------------------------------------------- -! REAL*8 :: e_an(imx,jmx,2,nmx), e_bb(imx,jmx,nmx), & -! e_ac(imx,jmx,lmx,nmx) - - !----------------------------------------------------------------------- - ! diagnostics (budget) - !----------------------------------------------------------------------- -! ! tendencies per time step and process -! REAL, TARGET :: bems(imx,jmx,nmx), bdry(imx,jmx,nmx), bstl(imx,jmx,nmx) -! REAL, TARGET :: bwet(imx,jmx,nmx), bcnv(imx,jmx,nmx) -! -! ! integrated tendencies per process -! REAL, TARGET :: tems(imx,jmx,nmx), tstl(imx,jmx,nmx) -! REAL, TARGET :: tdry(imx,jmx,nmx), twet(imx,jmx,nmx), tcnv(imx,jmx,nmx) - - ! global mass balance per time step - REAL*8 :: tmas0(nmx), tmas1(nmx) - REAL*8 :: dtems(nmx), dttrp(nmx), dtdif(nmx), dtcnv(nmx) - REAL*8 :: dtwet(nmx), dtdry(nmx), dtstl(nmx) - REAL*8 :: dtems2(nmx), dttrp2(nmx), dtdif2(nmx), dtcnv2(nmx) - REAL*8 :: dtwet2(nmx), dtdry2(nmx), dtstl2(nmx) - - - - ! executable statemenst DO n = 1, nmx @@ -264,7 +178,6 @@ SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, & g = g0*1.0E2 ! Pointer to the 3 classes considered in the source data files m = ipoint(n) - tsrc = 0.0 DO k = 1, ndsrc ! No flux if wet soil DO i = 1,imx @@ -274,8 +187,6 @@ SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, & SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ & SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0) ! write(0,*)u_ts0,den(n),diam(n),rhoa,g - ! Fraction of emerged surfaces (subtract lakes, coastal ocean,..) -! cw = 1.0 - water(i,j) ! Case of surface dry enough to erode IF (gwet(i,j) < 0.5) THEN @@ -284,23 +195,18 @@ SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, & ! Case of wet surface, no erosion u_ts = 100.0 END IF - srce = frac_s(n)*erod(i,j,m,k)*dxy(j) ! (m2) + srce = frac_s(n)*erod(i,j,m,k) ! (kg s^2 m^-5) IF (ilwi(i,j) == 1 ) THEN - dsrc = ch_dust(n,month)*srce*w10m(i,j)**2 & - * (w10m(i,j) - u_ts)*dt1 ! (kg) -! write(0,*)ch_dust(n,month),srce,w10m(i,j),u_ts,gwet(i,j) + !(kg s^2 m^-5)*(m^3 s^-3)*s = (kg/m2) per dt1 + dsrc = ch_dust(n,month)*srce*w10m(i,j)**2 * (w10m(i,j) - u_ts)*dt1 ELSE dsrc = 0.0 END IF -! dsrc = cw*ch_dust(k)*srce*w10m(i,j)**2 & -! * (w10m(i,j) - u_ts)*dt1 ! (kg) -! dsrc = cw*ch_dust(n,dt(1)%mn)*srce*w10m(i,j)**2 & -! * (w10m(i,j) - u_ts)*dt1 ! (kg) IF (dsrc < 0.0) dsrc = 0.0 ! Update dust mixing ratio at first model level. - tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1) - bems(i,j,n) = dsrc + tc(i,j,1,n) = tc(i,j,1,n) + dsrc/dz_lowest/airden(i,j,1) ! (kg/kg) + bems(i,j,n) = dsrc ! diagnostic (kg/m2) per dt1 END DO END DO END DO diff --git a/chem/module_gocart_dust_afwa.F b/chem/module_gocart_dust_afwa.F index 92cfcf4127..f19ee6f349 100755 --- a/chem/module_gocart_dust_afwa.F +++ b/chem/module_gocart_dust_afwa.F @@ -3,6 +3,9 @@ MODULE GOCART_DUST_AFWA ! AFWA dust routine ! Created by Sandra Jones (AER and AFWA) and Glenn Creighton (AFWA). ! +! A. Ukhov, 11 March 2021, Now "emis_dust" is accumulated dust +! emission (kg/m2). Before was instantenious flux (g m^-2 s^-1). +! Code cleanup, remove unused variables. USE module_data_gocart_dust @@ -11,9 +14,9 @@ MODULE GOCART_DUST_AFWA INTRINSIC max, min CONTAINS - SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy, & - v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,erod_dri,dustin,snowh,zs, & - ivgtyp,isltyp,vegfra,lai_vegmask,xland,xlat,xlong,gsw,dx,g,emis_dust, & + SUBROUTINE gocart_dust_afwa_driver(dt,config_flags,alt, & + chem,rho_phy,smois,u10,v10,p8w,dz8w,erod,erod_dri,dustin,snowh, & + isltyp,vegfra,lai_vegmask,xland,dx,g,emis_dust, & ust,znt,clay_wrf,sand_wrf,clay_nga,sand_nga,afwa_dustloft, & tot_dust,tot_edust,vis_dust,alpha,gamma,smtune,ustune, & ids,ide, jds,jde, kds,kde, & @@ -25,39 +28,29 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags - INTEGER, INTENT(IN ) :: julday, ktau, & - ids,ide, jds,jde, kds,kde, & + INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER,DIMENSION( ims:ime , jms:jme ) , & - INTENT(IN ) :: & - ivgtyp, & - isltyp - REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & - INTENT(IN ) :: moist + INTENT(IN ) :: isltyp + REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, & - INTENT(INOUT ) :: & - emis_dust + INTENT(INOUT ) :: emis_dust + REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , & - INTENT(IN ) :: smois - REAL, DIMENSION( config_flags%num_soil_layers ) , & - INTENT(IN ) :: zs + INTENT(IN ) :: smois REAL, DIMENSION( ims:ime , jms:jme, ndcls ) , & INTENT(IN ) :: erod,erod_dri REAL, DIMENSION( ims:ime , jms:jme, 5 ) , & INTENT(INOUT) :: dustin REAL, DIMENSION( ims:ime , jms:jme ) , & - INTENT(IN ) :: & - u10, & + INTENT(IN ) :: u10, & v10, & - gsw, & vegfra, & lai_vegmask, & xland, & - xlat, & - xlong, & ust, & znt, & clay_wrf, & @@ -65,13 +58,11 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u clay_nga, & sand_nga, & snowh - REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & - INTENT(IN ) :: & - alt, & - t_phy, & - dz8w,p8w, & - u_phy,v_phy,rho_phy - REAL, DIMENSION( ims:ime , jms:jme ) , & + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & + INTENT(IN ) :: alt, & + p8w,dz8w, & + rho_phy + REAL, DIMENSION( ims:ime , jms:jme ) , & INTENT( OUT) :: afwa_dustloft, & tot_edust REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & @@ -81,7 +72,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u ! Local variables - INTEGER :: nmx,smx,i,j,k,imx,jmx,lmx,lhave + INTEGER :: nmx,smx,i,j,k,imx,jmx,lmx INTEGER,DIMENSION (1,1) :: ilwi REAL*8, DIMENSION (1,1) :: erodtot REAL*8, DIMENSION (1,1) :: vegmask @@ -89,8 +80,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u REAL, DIMENSION( ims:ime , jms:jme ) :: clay,sand REAL*8, DIMENSION (1,1) :: drylimit REAL*8, DIMENSION (5) :: tc,bems - REAL*8, DIMENSION (1,1) :: airden,airmas,ustar - REAL*8, DIMENSION (1) :: dxy + REAL*8, DIMENSION (1,1) :: airden,ustar REAL*8, DIMENSION (3) :: massfrac REAL*8 :: volsm REAL :: conver,converi @@ -98,6 +88,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u REAL*8 :: zwant REAL, INTENT(IN ) :: alpha, gamma, smtune, ustune INTEGER :: smois_opt + real :: dz_lowest conver=1.e-9 converi=1.e9 @@ -143,10 +134,9 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u ! Air mass and density at lowest model level. - airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g airden(1,1)=rho_phy(i,kts,j) ustar(1,1)=ust(i,j) - dxy(1)=dx*dx + dz_lowest = dz8w(i,1,j) ! Friction velocity tuning constant (Note: recommend 0.7 for PXLSM, ! else use 1.0. This was created due to make the scheme compatible @@ -297,7 +287,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u ! Call dust emission routine. call source_dust(imx, jmx, lmx, nmx, smx, dt, tc, ustar, massfrac, & - erodtot, ilwi, dxy, gravsm, volsm, airden, airmas, & + erodtot, ilwi, gravsm, volsm, airden,dz_lowest, & bems, ustart, g, drylimit, alpha, gamma, smois_opt) IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then @@ -319,7 +309,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u chem(i,kts,j,p_ac0) = chem(i,kts,j,p_ac0) + 0.03*sum(tc(1:5))*converi*factnuma*soilfac chem(i,kts,j,p_corn) = chem(i,kts,j,p_corn) + 0.97*1.02*sum(tc(1:5))*converi*factnumc*soilfac ELSE - chem(i,kts,j,p_dust_1)=tc(1)*converi + chem(i,kts,j,p_dust_1)=tc(1)*converi ! tc(1...5) is (kg/kg), p_dust_1...5 (ug/kg) chem(i,kts,j,p_dust_2)=tc(2)*converi chem(i,kts,j,p_dust_3)=tc(3)*converi chem(i,kts,j,p_dust_4)=tc(4)*converi @@ -337,17 +327,19 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u afwa_dustloft(i,j)=ustune*w10-ustart*(LOG(10.0/znt(i,j))+psi)/0.4 ENDIF - ! For output diagnostics (g m^-2 s^-1) + ! A. Ukhov + ! for output diagnostics + ! bems (kg/m2) per dt + ! p_edust1...5 is accumulated dust emission (kg/m2) + emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1) + emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2) + emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3) + emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4) + emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5) - emis_dust(i,1,j,p_edust1)=bems(1) - emis_dust(i,1,j,p_edust2)=bems(2) - emis_dust(i,1,j,p_edust3)=bems(3) - emis_dust(i,1,j,p_edust4)=bems(4) - emis_dust(i,1,j,p_edust5)=bems(5) - ! Diagnostic total emitted dust (g m^-2 s^-1) - - tot_edust(i,j)=(bems(1)+bems(2)+bems(3)+bems(4)+bems(5)) + ! Diagnostic accumulated total emitted dust (kg/m2) + tot_edust(i,j)=tot_edust(i,j)+(bems(1)+bems(2)+bems(3)+bems(4)+bems(5)) ENDIF @@ -379,7 +371,7 @@ SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u END SUBROUTINE gocart_dust_afwa_driver SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& - erod, ilwi, dxy, gravsm, volsm, airden, airmas, & + erod, ilwi, gravsm, volsm, airden, dz_lowest, & bems, ustart, g0, drylimit, alpha, gamma, smois_opt) ! **************************************************************************** @@ -394,8 +386,6 @@ SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& ! * DRYLIMIT Upper GRAVSM (VOLSM) limit for air-dry soil (g/g) ! * ALPHA Constant to fudge the total emission of dust (1/m) ! * GAMMA Exponential tuning constant for erodibility (-) - ! * DXY Surface of each grid cell (m2) - ! * AIRMAS Mass of air for each grid box (kg) ! * AIRDEN Density of air for each grid box (kg/m3) ! * USTAR Friction velocity (m/s) ! * MASSFRAC Fraction of mass in each of 3 soil classes (-) @@ -405,6 +395,7 @@ SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& ! * IMX Number of I points (-) ! * JMX Number of J points (-) ! * LMX Number of L points (-) + ! * dz_lowest heigth of the lowest layer (m) ! * ! * Data (see module_data_gocart_dust): ! * SPOINT Pointer to 3 soil classes (-) @@ -440,11 +431,11 @@ SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& ! * DLNDP Dust bin width (-) ! * EMIT Total vertical mass flux (kg/m2/s) ! * EMIT_VOL Total vertical volume flux (m/s) - ! * DSRC Mass of emitted dust (kg/timestep/cell) + ! * DSRC Mass of emitted dust (kg/timestep/m2) ! * ! * Output: - ! * TC Total concentration of dust (kg/kg/timestep/cell) - ! * BEMS Source of each dust type (kg/timestep/cell) + ! * TC Total concentration of dust (kg/kg) + ! * BEMS Source of each dust type (kg/timestep/m2) ! * USTART Threshold friction vel. (bin 7) (m/s) ! * ! **************************************************************************** @@ -455,22 +446,22 @@ SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& REAL*8, INTENT(IN) :: ustar(imx,jmx) REAL*8, INTENT(IN) :: gravsm(imx,jmx) REAL*8, INTENT(IN) :: drylimit(imx,jmx) - REAL*8, INTENT(IN) :: dxy(jmx) - REAL*8, INTENT(IN) :: airden(imx,jmx,lmx), airmas(imx,jmx,lmx) + REAL*8, INTENT(IN) :: airden(imx,jmx,lmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(OUT) :: bems(imx,jmx,nmx) REAL, INTENT(IN) :: g0,dt1 REAL, INTENT(OUT) :: ustart INTEGER, INTENT(IN) :: smois_opt REAL*8, INTENT(IN) :: volsm + REAL, INTENT(IN) :: dz_lowest REAL*8 :: den(smx), diam(smx) - REAL*8 :: dvol(nmx), dlndp(nmx) +! REAL*8 :: dvol(nmx), dlndp(nmx) ! REAL*8 :: distr_dust(nmx) REAL*8 :: dsurface(smx), ds_rel(smx) REAL*8 :: massfrac(3) - REAL*8 :: u_ts0, u_ts, dsrc, srce, dmass, dvol_tot - REAL*8 :: emit, emit_vol + REAL*8 :: u_ts0, u_ts, dsrc, dmass!, dvol_tot + REAL*8 :: emit!, emit_vol REAL :: rhoa, g REAL*8 :: salt, stotal INTEGER :: i, j, m, n, s @@ -651,13 +642,12 @@ SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,& ! Calculate total mass emitted - dsrc = emit*distr_dust(n)*dxy(j)*dt1 ! (kg) + dsrc = emit*distr_dust(n)*dt1 ! (kg m^-2) per dt1 IF (dsrc < 0.0) dsrc = 0.0 ! Update dust mixing ratio at first model level. - - tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1) ! (kg/kg) - bems(i,j,n) = 1000.*dsrc/(dxy(j)*dt1) ! diagnostic (g/m2/s) + tc(i,j,1,n) = tc(i,j,1,n) + dsrc/dz_lowest/airden(i,j,1) ! (kg/kg) + bems(i,j,n) = dsrc ! diagnostic (kg/m2) per dt1 END DO END DO END DO diff --git a/chem/module_qf03.F b/chem/module_qf03.F index a5d92f1ba3..7fc24c8150 100644 --- a/chem/module_qf03.F +++ b/chem/module_qf03.F @@ -6,6 +6,10 @@ MODULE qf03 ! Modify the code for WRF_chem ! ! M. Klose, 2010-2015 Modifications + +! A. Ukhov, 11 March 2021, Now "bems" is dust emission (kg/m2) per dt. +! Before was instantenious flux (kg m^-2 s^-1). + !----------------------------------------------------------------------------------- ! Calculate sediment flux for multi-particle size soils as a weighted average of Q(d) ! dust emission F(d) for covered and moisture soil @@ -351,8 +355,8 @@ subroutine qf03_driver ( nmx, idst, g, rhop, rho, dt, & do n = 1, nbins ! fbin : [kg/m2/s], dz_lowest : [m], rho : [kg/m3], dt : [s] -> tc : [kg/kg-dryair] - tc(n) = tc(n) + cell_fbin(n)/dz_lowest/rho*dt ![kg/kg-dryair] - bems(n) = cell_fbin(n) ![kg/m2/s] + tc(n) = tc(n) + cell_fbin(n)/dz_lowest/rho*dt ! (kg/kg) + bems(n) = cell_fbin(n)*dt ! diagnostic (kg/m2) per dt enddo diff --git a/chem/module_uoc_dust.F b/chem/module_uoc_dust.F index 5bf113713f..c0ebcd8edc 100644 --- a/chem/module_uoc_dust.F +++ b/chem/module_uoc_dust.F @@ -7,6 +7,9 @@ MODULE uoc_dust ! ! For references and available schemes, see module_qf03.F ! Martina Klose, 29 May 2013 +! +! A. Ukhov, 11 March 2021, Now "emis_dust" is accumulated dust +! emission (kg/m2). Before was instantenious flux (ug m^-2 s^-1). ! !---------------------------------------------------------------------------- USE module_data_gocart_dust @@ -17,7 +20,7 @@ MODULE uoc_dust USE module_sf_ruclsm, ONLY:DRYSMC_ruc => DRYSMC CONTAINS - subroutine uoc_dust_driver(ktau,dt,config_flags, & + subroutine uoc_dust_driver(dt,config_flags, & chem,rho_phy,dz8w,smois,ust, & isltyp,vegfra,g,emis_dust, & ust_t_min, imod, rough_cor, smois_cor, & @@ -31,7 +34,7 @@ subroutine uoc_dust_driver(ktau,dt,config_flags, & IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags - INTEGER, INTENT(IN ) :: ktau, imod, & + INTEGER, INTENT(IN ) :: imod, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte @@ -235,13 +238,17 @@ subroutine uoc_dust_driver(ktau,dt,config_flags, & chem(i,kts,j,p_dust_3)=tc(3)*converi chem(i,kts,j,p_dust_4)=tc(4)*converi chem(i,kts,j,p_dust_5)=tc(5)*converi -! for output diagnostics - emis_dust(i,1,j,p_edust1)=bems(1)*converi - emis_dust(i,1,j,p_edust2)=bems(2)*converi - emis_dust(i,1,j,p_edust3)=bems(3)*converi - emis_dust(i,1,j,p_edust4)=bems(4)*converi - emis_dust(i,1,j,p_edust5)=bems(5)*converi ![kg/m2/s] -> [ug/m2/s] - + + ! A. Ukhov + ! for output diagnostics + ! bems (kg/m2) per dt + ! p_edust1...5 is accumulated dust emission (kg/m2) + emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1) + emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2) + emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3) + emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4) + emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5) + else ! no dust source emis_dust(i,1,j,p_edust1)=0. emis_dust(i,1,j,p_edust2)=0.