diff --git a/aqm_files.cmake b/aqm_files.cmake index b23a84f..394512d 100644 --- a/aqm_files.cmake +++ b/aqm_files.cmake @@ -237,7 +237,7 @@ list(APPEND aqm_CCTM_files ${VDIFF}/VDIFF_DATA.F ${VDIFF}/VDIFF_DIAG.F ${VDIFF}/VDIFF_MAP.F - ${VDIFF}/vdiffacmx.F + #${VDIFF}/vdiffacmx.F #${VDIFF}/vdiffproc.F ${VDIFF}/../../biog/megan3/BDSNP_MOD.F ${localCCTM}/o3totcol.f @@ -251,7 +251,8 @@ list(APPEND aqm_CCTM_files ${localCCTM}/ASX_DATA_MOD.F ${localCCTM}/DUST_EMIS.F ${localCCTM}/AERO_PHOTDATA.F - ${localCCTM}/phot.F + ${localCCTM}/vdiffacmx.F + #${localCCTM}/phot.F ${localCCTM}/centralized_io_module.F ${localCCTM}/centralized_io_util_module.F ) diff --git a/src/model/src/AERO_PHOTDATA.F b/src/model/src/AERO_PHOTDATA.F index e7acd25..41c708d 100644 --- a/src/model/src/AERO_PHOTDATA.F +++ b/src/model/src/AERO_PHOTDATA.F @@ -1396,8 +1396,8 @@ SUBROUTINE FASTER_OPTICS ( NR, NI, ALPHV, XLNSIG, BETA_EXT, BETA_SCAT, G ) REAL T1F1, T2F1, T1F2, T2F2, T1F3, T2F3 REAL T1G1, T2G1, T1G2, T2G2, T1G3, T2G3, T1G4, T2G4 - REAL T1G5, T2G5 - REAL(8) T1P1, T2P1 !(Wei Li) + REAL T1G5, T2G5 + REAL(8) T1P1, T2P1 C***the following are for calculating the Penndorff Coefficients @@ -1486,7 +1486,7 @@ SUBROUTINE FASTER_OPTICS ( NR, NI, ALPHV, XLNSIG, BETA_EXT, BETA_SCAT, G ) REAL QQSUM, QQF1,QQF2, QQF3, QQCORR REAL, PARAMETER :: DEGTORAD = PI180 - REAL(8), PARAMETER :: THREE_PI_TWO = 3.0 * PI / 2.0 !(Wei Li) + REAL(8), PARAMETER :: THREE_PI_TWO = 3.0 * PI / 2.0 C***FSB start calculation SIGMA_G = EXP( XLNSIG ) @@ -1534,11 +1534,11 @@ SUBROUTINE FASTER_OPTICS ( NR, NI, ALPHV, XLNSIG, BETA_EXT, BETA_SCAT, G ) ALPHA_I = F2 BEXT = B BSCAT = B - PENN1 = DBLE(0.0) !(Wei Li) - PENN2 = DBLE(0.0) !(Wei Li) + PENN1 = DBLE(0.0) + PENN2 = DBLE(0.0) - ALPHV2 = DBLE(ALPHV * ALPHV) !(Wei Li) - ALPHV3 = DBLE(ALPHV2 * ALPHV) !(Wei Li) + ALPHV2 = DBLE(ALPHV * ALPHV) + ALPHV3 = DBLE(ALPHV2 * ALPHV) IF ( NI .GT. 0.0 ) THEN @@ -1605,14 +1605,14 @@ SUBROUTINE FASTER_OPTICS ( NR, NI, ALPHV, XLNSIG, BETA_EXT, BETA_SCAT, G ) EXPFAC2 = EXP( 2.0 * XLNSIG2 ) EXPFAC3 = EXP( 4.5 * XLNSIG2 ) - T1P1 = DBLE(A1 + A2 * ALPHV2 * EXPFAC2) !(Wei Li) - T2P1 = DBLE(A3 * ALPHV3 * EXPFAC3) !(Wei Li) + T1P1 = DBLE(A1 + A2 * ALPHV2 * EXPFAC2) + T2P1 = DBLE(A3 * ALPHV3 * EXPFAC3) C***PENN1 is the analytic integral of the Pendorff formulae over C*** a log normal particle size distribution. - PENN1 = DBLE(THREE_PI_TWO * ( T1P1 + T2P1 ) ) !(Wei Li) - PENN2 = DBLE(THREE_PI_TWO * T2P1) !(Wei Li) + PENN1 = DBLE(THREE_PI_TWO * ( T1P1 + T2P1 )) + PENN2 = DBLE(THREE_PI_TWO * T2P1) END IF ! test for ni > 0.0 @@ -1852,7 +1852,7 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) REAL C, CEXT, CSCAT REAL B, BEXT, BSCAT REAL BBFAC - REAL(8) ALPHV !(Wei Li) + REAL(8) ALPHV REAL ALPHA_I REAL A, LOGX2, XLNSIG, XLNSIG2, MM1 @@ -1866,7 +1866,7 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) REAL LARGEEXT ! large sphere limit for extinction REAL SMALL_G, LARGE_G - REAL(8) ALPHV2, ALPHV3 !(Wei Li) + REAL(8) ALPHV2, ALPHV3 REAL X_ALPHA, X_ALPHA2, X_ALPHA3 REAL FCORR REAL EXPFAC2, EXPFAC3 @@ -1879,12 +1879,12 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) REAL T1F1, T2F1, T1F2, T2F2, T1F3, T2F3 REAL T1G1, T2G1, T1G2, T2G2, T1G3, T2G3, T1G4, T2G4 REAL T1G5, T2G5 - REAL(8) T1P1, T2P1 !(Wei Li) + REAL(8) T1P1, T2P1 C***the following are for calculating the Penndorff Coefficients REAL A1, A2, A3 - REAL(8) PENN1, PENN2 !(Wei Li) + REAL(8) PENN1, PENN2 REAL XNR, XNI, XNR2, XNI2, XNRI, XNRI2, XNRMI REAL XRI, XRI2, XRI36, XNX, XNX2 REAL Z1, Z12, Z2, XC1 @@ -1969,7 +1969,7 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) REAL QQSUM, QQF1,QQF2, QQF3, QQCORR REAL, PARAMETER :: DEGTORAD = PI180 - REAL(8), PARAMETER :: THREE_PI_TWO = 3.0 * PI / 2.0 !(Wei Li) + REAL(8), PARAMETER :: THREE_PI_TWO = 3.0 * PI / 2.0 REAL, PARAMETER :: SCALE = 1.00E+9 @@ -1982,9 +1982,9 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) C***FSB start calculation XLNSIG = LOG( SIGMA_G ) - ALPHV = DBLE(SCALE * PI * DGN * EXP( 3.0 * XLNSIG * XLNSIG ) / LAMBDA) !(Wei Li) - ALPHV2 = DBLE(ALPHV * ALPHV) !(Wei Li) - ALPHV3 = DBLE(ALPHV * ALPHV * ALPHV) !(Wei Li) + ALPHV = DBLE(SCALE * PI * DGN * EXP( 3.0 * XLNSIG * XLNSIG ) / LAMBDA) + ALPHV2 = DBLE(ALPHV * ALPHV) + ALPHV3 = DBLE(ALPHV * ALPHV * ALPHV) XLNSIG2 = XLNSIG * XLNSIG A = 0.5 / XLNSIG2 @@ -2024,8 +2024,8 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) ALPHA_I = F2 BEXT = B BSCAT = B - PENN1 = DBLE(0.0) !(Wei Li) - PENN2 = DBLE(0.0) !(Wei Li) + PENN1 = DBLE(0.0) + PENN2 = DBLE(0.0) IF ( NI .GT. 0.0 ) THEN @@ -2079,25 +2079,25 @@ SUBROUTINE FAST_OPTICS( NR, NI, LAMBDA, DGN, SIGMA_G, BETA_EXT, BETA_SCAT, G ) Z12 = Z1 * Z1 Z2 = 4.0 * XNRI2 + 12.0 * XNRMI + 9.0 XC1 = 8.0 / ( 3.0 * Z12 ) - A1 = DBLE(24.0 * XRI / Z1) !(Wei Li) + A1 = DBLE(24.0 * XRI / Z1) - A2 = DBLE(44.0 * XRI / 15.0 + 20.0 * XRI / ( 3.0 * Z2 ) + + A2 = DBLE(4.0 * XRI / 15.0 + 20.0 * XRI / ( 3.0 * Z2 ) + & 4.8 * XRI * ( 7.0 * XNRI2 + - & 4.0 * ( XNRMI - 5.0 ) ) / Z12 ) !(Wei Li) + & 4.0 * ( XNRMI - 5.0 ) ) / Z12) - A3 = DBLE(XC1 * ( XNX2 - XRI36 )) !(Wei Li) + A3 = DBLE(XC1 * ( XNX2 - XRI36 )) EXPFAC2 = EXP( 2.0 * XLNSIG2 ) EXPFAC3 = EXP( 4.5 * XLNSIG2 ) - T1P1 = DBLE(A1 + A2 * ALPHV2 * EXPFAC2) !(Wei Li) - T2P1 = DBLE(A3 * ALPHV3 * EXPFAC3) !(Wei Li) + T1P1 = DBLE(A1 + A2 * ALPHV2 * EXPFAC2) + T2P1 = DBLE(A3 * ALPHV3 * EXPFAC3) C***PENN1 is the analytic integral of the Pendorff formulae over C*** a log normal particle size distribution. - PENN1 = DBLE(THREE_PI_TWO * ( T1P1 + T2P1 )) !(Wei Li) - PENN2 = DBLE(THREE_PI_TWO * T2P1 ) !(Wei Li) + PENN1 = DBLE(THREE_PI_TWO * ( T1P1 + T2P1 )) + PENN2 = DBLE(THREE_PI_TWO * T2P1) END IF ! test of ni > 0.0 diff --git a/src/model/src/ASX_DATA_MOD.F b/src/model/src/ASX_DATA_MOD.F index 8270f49..0a03e63 100644 --- a/src/model/src/ASX_DATA_MOD.F +++ b/src/model/src/ASX_DATA_MOD.F @@ -130,10 +130,8 @@ Module ASX_DATA_MOD Logical, Allocatable :: CONVCT ( :,: ) ! convection flag Real, Allocatable :: PBL ( :,: ) ! pbl height (m) ! Real, Allocatable :: NACL_EMIS( :,: ) ! NACL mass emission rate of particles with d <10 um (g/m2/s) - Real, Allocatable :: COSZEN ( :,: ) ! Cosine of the zenith angle - Real, Allocatable :: CFRAC ( :,: ) ! cloud fraction -!> Inline Canopy Processes (Wei Li) +!> Inline Canopy Processes Real, Allocatable :: FCH ( :,: ) ! Forest Canopy Height (m) Real, Allocatable :: FRT ( :,: ) ! Forest Fraction Real, Allocatable :: CLU ( :,: ) ! Clumping Index @@ -142,21 +140,23 @@ Module ASX_DATA_MOD Real, Allocatable :: C1R ( :,: ) ! cumulative LAI fraction hc to 0.75 * hc Real, Allocatable :: C2R ( :,: ) ! cumulative LAI fraction hc to 0.50 * hc Real, Allocatable :: C3R ( :,: ) ! cumulative LAI fraction hc to 0.35 * hc - Real, Allocatable :: C4R ( :,: ) ! cumulative LAI fraction hc to 0.20 * hc - -!> FENGSHA option (Wei Li) - Real, Allocatable :: CLAYF ( :,: ) ! Fractional Clay Content - Real, Allocatable :: SANDF ( :,: ) ! Fractional Sand Content + Real, Allocatable :: C4R ( :,: ) ! cumulative LAI fraction hc to 0.20 * hc +!> FENGSHA option + Real, Allocatable :: CLAYF ( :,: ) ! Fractional Clay Content + Real, Allocatable :: SANDF ( :,: ) ! Fractional Sand Content Real, Allocatable :: DRAG ( :,: ) ! Drag Partion Real, Allocatable :: UTHR ( :,: ) ! Dry Threshold Friction Velocity + + Real, Allocatable :: COSZEN ( :,: ) ! Cosine of the zenith angle + Real, Allocatable :: CFRAC ( :,: ) ! cloud fraction !> U and V wind components on the cross grid points Real, Allocatable :: UWIND ( :,:,: ) ! [m/s] Real, Allocatable :: VWIND ( :,:,: ) ! [m/s] !> 3-D meteorological fields: Real, Allocatable :: KZMIN ( :,:,: ) ! minimum Kz [m**2/s] Real, Allocatable :: PRES ( :,:,: ) ! pressure [Pa] - Real, Allocatable :: PRESF ( :,:,: ) ! full layer pressure [Pa] (Wei Li) + Real, Allocatable :: PRESF ( :,:,: ) ! full layer pressure [Pa] Real, Allocatable :: QV ( :,:,: ) ! water vapor mixing ratio Real, Allocatable :: QC ( :,:,: ) ! cloud water mixing ratio Real, Allocatable :: THETAV ( :,:,: ) ! potential temp @@ -169,8 +169,8 @@ Module ASX_DATA_MOD Real, Allocatable :: RJACM ( :,:,: ) ! reciprocal mid-layer Jacobian Real, Allocatable :: RJACF ( :,:,: ) ! reciprocal full-layer Jacobian Real, Allocatable :: RRHOJ ( :,:,: ) ! reciprocal density X Jacobian - Real, Allocatable :: UWINDA ( :,:,: ) ! [m/s] (Wei Li) - Real, Allocatable :: VWINDA ( :,:,: ) ! [m/s] (Wei Li) + Real, Allocatable :: UWINDA ( :,:,: ) ! [m/s] + Real, Allocatable :: VWINDA ( :,:,: ) ! [m/s] End Type MET_Type Type :: GRID_Type @@ -261,16 +261,22 @@ Module ASX_DATA_MOD Real, allocatable, private :: BUFF2D( :,: ) ! 2D temp var Real, allocatable, private :: BUFF3D( :,:,: ) ! 3D temp var -! Canopy option control (Wei Li) - CHARACTER( 20 ), SAVE :: CTM_CANOPY_SHADE = 'CTM_CANOPY_SHADE '! env var for in-line +! Canopy option control + CHARACTER( 20 ), SAVE :: CTM_CANOPY_SHADE = 'CTM_CANOPY_SHADE'! env var for in-line LOGICAL, PUBLIC, SAVE :: CANOPY_SHADE ! flag in-lining canopy shading -! FENGSHA option control (Wei Li) - CHARACTER( 18 ), SAVE :: CTM_WBDUST_FENGSHA = 'CTM_WBDUST_FENGSHA' ! env var for in-line - LOGICAL, PUBLIC, SAVE :: FENGSHA ! flag for fengsha option - +! FENGSHA option control + CHARACTER( 18 ), SAVE :: CTM_WBDUST_FENGSHA ='CTM_WBDUST_FENGSHA' ! env var for in-line + LOGICAL, PUBLIC, SAVE :: FENGSHA ! flag for fengsha option INTEGER IOSX ! i/o and allocate memory status + + + + + + + DATA subname( 1), dif0( 1), ar( 1), meso( 1), lebas( 1) / 'SO2 ', 0.1089, 10.0, 0.0, 35.0 / DATA subname( 2), dif0( 2), ar( 2), meso( 2), lebas( 2) / 'H2SO4 ', 0.1091, 8000.0, 0.0, 49.0 / DATA subname( 3), dif0( 3), ar( 3), meso( 3), lebas( 3) / 'NO2 ', 0.1361, 2.0, 0.1, 21.0 / @@ -623,15 +629,15 @@ Subroutine INIT_MET ( JDATE, JTIME ) & Met_Data%CONVCT ( NCOLS,NROWS ), & Met_Data%PBL ( NCOLS,NROWS ), ! & Met_Data%NACL_EMIS( NCOLS,NROWS ), + & Met_Data%UWINDA ( NCOLS,NROWS,NLAYS ), + & Met_Data%VWINDA ( NCOLS,NROWS,NLAYS ), & Met_Data%COSZEN ( NCOLS,NROWS ), & Met_Data%CFRAC ( NCOLS,NROWS ), - & Met_Data%UWINDA ( NCOLS,NROWS,NLAYS ), !(Wei Li) - & Met_Data%VWINDA ( NCOLS,NROWS,NLAYS ), !(Wei Li) & Met_Data%UWIND ( NCOLS+1,NROWS+1,NLAYS ), & Met_Data%VWIND ( NCOLS+1,NROWS+1,NLAYS ), & Met_Data%KZMIN ( NCOLS,NROWS,NLAYS ), - & Met_Data%PRESF ( NCOLS,NROWS,1:NLAYS+1 ), !(Wei Li) & Met_Data%PRES ( NCOLS,NROWS,NLAYS ), + & Met_Data%PRESF ( NCOLS,NROWS,1:NLAYS+1 ), & Met_Data%QV ( NCOLS,NROWS,NLAYS ), & Met_Data%QC ( NCOLS,NROWS,NLAYS ), & Met_Data%THETAV ( NCOLS,NROWS,NLAYS ), @@ -716,7 +722,7 @@ Subroutine INIT_MET ( JDATE, JTIME ) Grid_Data%BSLP = 0.0 End If -!> ccccccccccccccccccccc canopy shade option!ccccccccccccccccccccc (Wei Li) +!> ccccccccccccccccccccc canopy shade option!ccccccccccccccccccccc CANOPY_SHADE = ENVYN( 'CTM_CANOPY_SHADE', & 'Flag for in-line canopy shading', & .FALSE., IOSX ) @@ -725,7 +731,7 @@ Subroutine INIT_MET ( JDATE, JTIME ) ! XMSG = 'Using in-line canopy shading option' ! CALL M3MSG2( XMSG ) ! END IF - If ( CANOPY_SHADE ) Then + If ( CANOPY_SHADE ) Then ALLOCATE( Met_Data%FCH ( NCOLS,NROWS ), & Met_Data%FRT ( NCOLS,NROWS ), & Met_Data%CLU ( NCOLS,NROWS ), @@ -740,9 +746,9 @@ Subroutine INIT_MET ( JDATE, JTIME ) XMSG = 'Failure allocating Canopy Shade variables' Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) End If - End If + End If -!> ccccccccccccccccccccc Fengsha option!ccccccccccccccccccccc (Wei Li) +!> ccccccccccccccccccccc Fengsha option!ccccccccccccccccccccc FENGSHA = ENVYN( CTM_WBDUST_FENGSHA, & 'Flag for in-line fengsha ', & .FALSE., IOSX ) @@ -755,7 +761,7 @@ Subroutine INIT_MET ( JDATE, JTIME ) If ( FENGSHA ) Then ALLOCATE( Met_Data%CLAYF ( NCOLS,NROWS ), & Met_Data%SANDF ( NCOLS,NROWS ), - & Met_Data%DRAG ( NCOLS,NROWS ), + & Met_Data%DRAG ( NCOLS,NROWS ), & Met_Data%UTHR ( NCOLS,NROWS ), & STAT = ALLOCSTAT ) If ( ALLOCSTAT .Ne. 0 ) Then @@ -764,7 +770,6 @@ Subroutine INIT_MET ( JDATE, JTIME ) End If End If - !> ccccccccccccccccccccc enable backward compatiblity ccccccccccccccccccccc IF (RCA_AVAIL) THEN @@ -911,11 +916,7 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP ) call interpolate_var ('PRES', jdate, jtime, Met_Data%PRES) - call interpolate_var ('PRESF', jdate, jtime, Met_Data%PRESF) !(Wei Li) - - call interpolate_var ('UWINDA', jdate, jtime, Met_Data%UWINDA) !(Wei Li) - - call interpolate_var ('VWINDA', jdate, jtime, Met_Data%VWINDA) !(Wei Li) + call interpolate_var ('PRESF', jdate, jtime, Met_Data%PRESF) call interpolate_var ('ZF', jdate, jtime, Met_Data%ZF) @@ -941,6 +942,10 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP ) call interpolate_var ('QC', jdate, jtime, Met_Data%QC) + call interpolate_var ('UWINDA', jdate, jtime, Met_Data%UWINDA) + + call interpolate_var ('VWINDA', jdate, jtime, Met_Data%VWINDA) + C-------------------------------- MET_CRO_2D -------------------------------- C Vegetation and surface vars call interpolate_var ('LAI', jdate, jtime, Met_Data%LAI) @@ -949,28 +954,37 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP ) call interpolate_var ('ZRUF', jdate, jtime, Met_Data%Z0) -C Canopy vars (Wei Li) - If ( CANOPY_SHADE ) Then - call interpolate_var ('FCH', jdate, jtime, Met_Data%FCH) - call interpolate_var ('FRT', jdate, jtime, Met_Data%FRT) - call interpolate_var ('CLU', jdate, jtime, Met_Data%CLU) - call interpolate_var ('POPU', jdate, jtime, Met_Data%POPU) - call interpolate_var ('LAIE', jdate, jtime, Met_Data%LAIE) - call interpolate_var ('C1R', jdate, jtime, Met_Data%C1R) - call interpolate_var ('C2R', jdate, jtime, Met_Data%C2R) - call interpolate_var ('C3R', jdate, jtime, Met_Data%C3R) - call interpolate_var ('C4R', jdate, jtime, Met_Data%C4R) - End If +C Canopy vars + If ( CANOPY_SHADE ) Then + call interpolate_var ('FCH', jdate, jtime, Met_Data%FCH) -C FENGSHA vars (Wei Li) - If ( CANOPY_SHADE ) Then - call interpolate_var ('CLAYF', jdate, jtime, Met_Data%CLAYF) - call interpolate_var ('SANDF', jdate, jtime, Met_Data%SANDF) - call interpolate_var ('DRAG', jdate, jtime, Met_Data%DRAG) - call interpolate_var ('UTHR', jdate, jtime, Met_Data%UTHR) - End If + call interpolate_var ('FRT', jdate, jtime, Met_Data%FRT) + + call interpolate_var ('CLU', jdate, jtime, Met_Data%CLU) + + call interpolate_var ('POPU', jdate, jtime, Met_Data%POPU) + call interpolate_var ('LAIE', jdate, jtime, Met_Data%LAIE) + call interpolate_var ('C1R', jdate, jtime, Met_Data%C1R) + + call interpolate_var ('C2R', jdate, jtime, Met_Data%C2R) + + call interpolate_var ('C3R', jdate, jtime, Met_Data%C3R) + + call interpolate_var ('C4R', jdate, jtime, Met_Data%C4R) + + End if !(canopy option) +C FENGSHA vars + If (FENGSHA ) Then + call interpolate_var ('CLAYF', jdate, jtime, Met_Data%CLAYF) + + call interpolate_var ('SANDF', jdate, jtime, Met_Data%SANDF) + + call interpolate_var ('DRAG', jdate, jtime, Met_Data%DRAG) + + call interpolate_var ('UTHR', jdate, jtime, Met_Data%UTHR) + End if C Soil vars call interpolate_var ('SOIM1', jdate, jtime, Met_Data%SOIM1) diff --git a/src/model/src/DUST_EMIS.F b/src/model/src/DUST_EMIS.F index fe3352e..0a6d6fe 100644 --- a/src/model/src/DUST_EMIS.F +++ b/src/model/src/DUST_EMIS.F @@ -17,7 +17,7 @@ ! subject to their copyright restrictions. ! !------------------------------------------------------------------------! -!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: module dust_emis C----------------------------------------------------------------------- @@ -112,7 +112,7 @@ module dust_emis C Number of soil types: For WRF there are 16 types; integer, parameter :: nsltyp = 16 -C Variables for FENGSHA dust scheme (beiming tang) +C Variables for FENGSHA dust scheme real, save :: dust_alpha ! tuning parameter for FENGSHA dust emission flux C Variables for the windblown dust diagnostic file: @@ -273,13 +273,14 @@ function dust_emis_init( jdate, jtime, tstep ) result( success ) success = .false.; return end if - if ( fengsha ) then !(beiming tang) + if ( fengsha ) then C Disable diagnostic output if FENGSHA is used dustem_diag = .false. C Allocate private arrays - allocate( tfb( ncols,nrows ), stat = status ) + !allocate( tfb( ncols,nrows ), stat = status ) !tfa and tfb + !were removed in most recent note. if ( status .ne. 0 ) then xmsg = '*** Failure allocating WMAX or TFB' @@ -471,7 +472,7 @@ function dust_emis_init( jdate, jtime, tstep ) result( success ) end if ! dustem_diag - end if !dust scheme (beiming tang) + end if !dust scheme l2sgj = log( sigj ) * log( sigj ) l2sgk = log( sigk ) * log( sigk ) @@ -635,7 +636,7 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, real, parameter :: sigb_mb = sigb * mb ! = 0.5 real, parameter :: betab_mb = betab * mb ! = 45.0 - character( 24 ) :: ctm_wbdust_fengsha_alpha = 'CTM_WBDUST_FENGSHA_ALPHA' ! env var to (beiming tang) + character( 24 ) :: ctm_wbdust_fengsha_alpha = 'CTM_WBDUST_FENGSHA_ALPHA' ! env var to ! retrieve FENGSHA scaling factor character( 16 ) :: pname = 'GET_DUST_EMIS' @@ -675,7 +676,7 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, real :: lambdav ! vegetation roughness density - Shao et. al [Aus. J. Soil Res., 1996] real :: flxfac1, flxfac2 ! combined soli type mapping factors real :: hflux, vflux ! horizontal and vertical dust flux - real :: v2h ! vertical/horizontal dust flux ratio (beiming tang) + real :: v2h ! vertical/horizontal dust flux ratio real :: wm ! max adsorb water [%] real :: jday integer :: emap( n_dlcat+1 ) @@ -762,9 +763,9 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, if ( firstime ) then firstime = .false. if ( fengsha ) then - dust_alpha = 0.05 ! default (beiming tang) + dust_alpha = 0.05 ! default dust_alpha = envreal( ctm_wbdust_fengsha_alpha, - & 'Emission global scaling factor for FENGSHA dust scheme', + & 'Emission global scaling factor for FENGSHA dust scheme', & dust_alpha, status ) if ( status .ne. 0 ) then xmsg = '*** Failure retrieving FENGSHA scaling factor' @@ -772,20 +773,20 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, end if write(xmsg,'("Using FENGSHA alpha = ",g12.5)') dust_alpha call m3msg2 ( xmsg ) - else + else !this is default case allocate ( ustr( ncols,nrows,n_dlcat+1 ), - & qam( ncols,nrows,n_dlcat+1 ), - & fruf( ncols,nrows,n_dlcat+1 ), - & kvh( ncols,nrows,n_dlcat+1 ), - & elus( ncols,nrows,n_dlcat+1 ), stat = status ) + & qam( ncols,nrows,n_dlcat+1 ), + & fruf( ncols,nrows,n_dlcat+1 ), + & kvh( ncols,nrows,n_dlcat+1 ), + & elus( ncols,nrows,n_dlcat+1 ), stat = status ) if ( status .ne. 0 ) then xmsg = '*** Failure allocating USTR, QAM, FRUF, KVH, or ELUS' call m3exit( pname, jdate, jtime, xmsg, xstat1 ) end if - end if + end if !end fengsha end if -C---Select dust scheme (beiming tang) +C---Select dust scheme if ( fengsha ) then do r = 1, nrows @@ -837,7 +838,7 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, rlay1hgt = rjacm ( c,r ) / cellhgt - dust_em( c,r ) = dust_alpha * vflux * rlay1hgt *tfa(c,r) * tfb(c,r) + dust_em( c,r ) = dust_alpha * vflux * rlay1hgt! *tfa(c,r) * tfb(c,r) end if ! if rain & land & snow & drag @@ -1128,7 +1129,7 @@ subroutine get_dust_emis( jdate, jtime, tstep, rjacm, cellhgt, & out of total cells:', & dryhit, (c-1)*(r-1) #endif - end if !dust scheme (beiming tang) + end if !dust scheme do r = 1, nrows do c = 1, ncols @@ -1374,7 +1375,7 @@ function dust_hflux( ndp, dp, soiltxt, fmoit, fruf, ustr, sd_ep, dens ) end function dust_hflux - function dust_hflux_fengsha( ustar, fmoit, drag, uthr,ssm, dens ) !beiming tang + function dust_hflux_fengsha( ustar, fmoit, drag, uthr,ssm, dens ) & result( hflux ) implicit none diff --git a/src/model/src/o3totcol.f b/src/model/src/o3totcol.f index a500cee..0e069e9 100644 --- a/src/model/src/o3totcol.f +++ b/src/model/src/o3totcol.f @@ -96,12 +96,11 @@ subroutine o3totcol ( latitude, longitude, jdate, jtime, ozone ) real, allocatable, save :: lat( : ) real, allocatable, save :: lon( : ) real, allocatable, save :: oz( :, :, : ) ! two timesteps for interpolation - + + character( 8 ) :: label logical, save :: firsttime = .true. - character( 8 ) :: label !(Wei li) - real, external :: yr2day !(Wei Li: from io/ioapi/yr2day.F) - character*24, external :: dt2str !(Wei Li) - + real, external :: yr2day + character*24, external :: dt2str !---------------------------------------------------------------------- if ( firsttime ) then @@ -115,10 +114,10 @@ subroutine o3totcol ( latitude, longitude, jdate, jtime, ozone ) call m3exit ( pname, jdate, 0, xmsg, xstat1 ) end if - ! read nlat, nlon (Wei Li) + ! read nlat, nlon rewind( tmunit ) - read( tmunit, *) label, nlat - read( tmunit, *) label, nlon + read( tmunit, *) label,nlat + read( tmunit, *) label,nlon write(logdev,'(a,i7,a,i7)')'OMI Ozone column data has Lat by Lon Resolution: ', & nlat,'X',nlon @@ -135,17 +134,16 @@ subroutine o3totcol ( latitude, longitude, jdate, jtime, ozone ) call m3exit ( pname, jdate, 0, xmsg, xstat1 ) end if -! Assign values to array of longitudes: lon +!! Assign values to array of longitudes: lon ! x2 = 360.0 / real( nlon - 1 ) ! do ilon = 1, nlon ! lon( ilon ) = -180.0 + x2 * real( ilon - 1 ) ! end do - !read in longitudes instead (Wei Li) - read( tmunit, * ) label, label, lon + read( tmunit, * ) label, label, lon ! read in longitudes nrecs = 0 - ! read( tmunit, * ) !skip header record. Wei Li:no need here +! read( tmunit, * ) ! skip header record do read( tmunit, *, iostat=ios ) if ( ios .ne. 0 ) exit @@ -166,7 +164,6 @@ subroutine o3totcol ( latitude, longitude, jdate, jtime, ozone ) end if rewind( tmunit ) - ! skip header records read( tmunit, * ) read( tmunit, * ) read( tmunit, * ) @@ -310,7 +307,6 @@ subroutine o3totcol ( latitude, longitude, jdate, jtime, ozone ) end do x1loop ! Determine the corresponding bounding ozone values for all lats and lons rewind( tmunit ) - ! skip header records read( tmunit,* ) read( tmunit,* ) read( tmunit,* ) diff --git a/src/model/src/vdiffacmx.F b/src/model/src/vdiffacmx.F index 06954c4..3ec8114 100644 --- a/src/model/src/vdiffacmx.F +++ b/src/model/src/vdiffacmx.F @@ -16,13 +16,8 @@ ! subject to their copyright restrictions. ! !------------------------------------------------------------------------! - -C RCS file, release, date & time of last delta, author, state, [and locker] -C $Header: /project/yoj/arc/CCTM/src/vdiff/acm2/vdiffacm2.F,v 1.13 2012/01/19 14:37:47 yoj Exp $ - C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -! SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, CNGRD, DDEPJ, DDEPJ_FST ) - SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) + SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, CNGRD ) C----------------------------------------------------------------------- C Asymmetric Convective Model v2 (ACM2/ACM1) -- Pleim(2006/2014) @@ -30,7 +25,7 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) C calculates vertical diffusion C Subroutines and Functions Called: -C INIT3, SEC2TIME, TIME2SEC, WRITE3, NEXTIME, +C SEC2TIME, TIME2SEC, WRITE3, NEXTIME, C M3EXIT, EDDYX, TRI, MATRIX, PA_UPDATE_EMIS, PA_UPDATE_DDEP C Revision History: @@ -46,17 +41,39 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) C matrix solver C 30 May 14 J.Young: split vdiff calculation out of vdiff proc. C 07 Nov 14 J.Bash: Updated for the ASX_DATA_MOD shared data module. +C 02 Nov 2018: L.Zhou, S.Napelenok: isam implementation +C May 2019 J.Pleim Changed from sigma coords to Z coords for compatability w/ MPAS and WRF +C 12 Dec 19 S.L.Napelenok: ddm-3d implementation for version 5.3.1 +C 15 Jun 21 J. Pleim: implemented HONO fix for dry depsotion flux C----------------------------------------------------------------------- USE CGRID_SPCS ! CGRID mechanism species USE GRID_CONF - USE EMIS_DEFN + USE DESID_VARS, ONLY : VDEMIS_DIFF,DESID_LAYS + USE DESID_PARAM_MODULE, ONLY : DESID_N_SRM USE DEPV_DEFN USE ASX_DATA_MOD USE VDIFF_MAP USE UTILIO_DEFN - USE BIDI_MOD +! USE BIDI_MOD +! USE LSM_MOD, ONLY: N_LUFRAC +! USE VDIFF_DIAG, NLPCR => NLPCR_MEAN C USE PT3D_EMIS_DEFN + USE HGRD_DEFN,only : COLSX_PE, ROWSX_PE + USE BDSNP_MOD, only: GET_N_DEP +#ifdef isam + USE SA_DEFN, ONLY: N_SPCTAG, ISAM, VNAM_SPCTAG, TRANSPORT_SPC, + & SA_VDEMIS_DIFF, ITAG, NTAG_SA, NSPC_SA, + & S_SPCTAG, T_SPCTAG, SA_DDEP, OTHRTAG, ISAM_SPEC, + & L_NO3, SA_BIDI, BIDITAG, L_NH4 +#endif + +#ifdef sens + USE DDM3D_DEFN, ONLY : NP, NPMAX, SNGRD, S_DDEP, S_PLDV, SVDEMIS_DIFF, + & SENS, S_EMIS, S_DD, S_UU, S_DDBF, S_POL, + & S_DELC, S_PLDV_HONO +#endif + IMPLICIT NONE @@ -70,8 +87,6 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) REAL, INTENT( INOUT ) :: SEDDY ( :,:,: ) ! flipped EDDYV REAL, INTENT( INOUT ) :: DDEP ( :,:,: ) ! ddep accumulator REAL, INTENT( INOUT ) :: ICMP ( :,:,: ) ! component flux accumlator - REAL, INTENT( INOUT ), OPTIONAL :: DDEPJ ( :,:,:,: ) ! ddep for mosaic - REAL, INTENT( INOUT ), OPTIONAL :: DDEPJ_FST( :,:,:,: ) ! ddep for stomtal/cuticular pathway REAL, INTENT( INOUT ) :: CNGRD ( :,:,:,: ) ! cgrid replacement C Parameters: @@ -80,6 +95,8 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) REAL, PARAMETER :: THETA = 0.5, & THBAR = 1.0 - THETA +! REAL, PARAMETER :: EPS = 1.0E-06 + C External Functions: None C Local Variables: @@ -87,39 +104,128 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) CHARACTER( 16 ), SAVE :: PNAME = 'VDIFFACMX' LOGICAL, SAVE :: FIRSTIME = .TRUE. + LOGICAL, SAVE :: SPECLOG = .TRUE. ! For BDSNP REAL, ALLOCATABLE, SAVE :: DD_FAC ( : ) ! combined subexpression REAL, ALLOCATABLE, SAVE :: DDBF ( : ) ! secondary DDEP - REAL, ALLOCATABLE, SAVE :: CMPF ( : ) ! intermediate CMP + REAl, ALLOCATABLE, SAVE :: CMPF ( : ) ! intermediate CMP REAL, ALLOCATABLE, SAVE :: CONC ( :,: ) ! secondary CGRID expression REAL, ALLOCATABLE, SAVE :: EMIS ( :,: ) ! emissions subexpression REAL DTDENS1 ! DT * layer 1 air density C ACM Local Variables + REAL :: EDDY ( NLAYS ) ! local converted eddyv + REAL MEDDY ! ACM2 intermediate var + REAL MBAR ! ACM2 mixing rate (S-1) + REAL :: MBARKS( NLAYS ) ! by layer + REAL :: MDWN ( NLAYS ) ! ACM down mix rate + REAL :: MFAC ( NLAYS ) ! intermediate loop factor + REAL :: AA ( NLAYS ) ! matrix column one + REAL :: BB1 ( NLAYS ) ! diagonal for MATRIX1 + REAL :: BB2 ( NLAYS ) ! diagonal for TRI + REAL :: CC ( NLAYS ) ! subdiagonal + REAL :: EE1 ( NLAYS ) ! superdiagonal for MATRIX1 + REAL :: EE2 ( NLAYS ) ! superdiagonal for TRI + REAL, ALLOCATABLE, SAVE :: DD ( :,: ) ! R.H.S + REAL, ALLOCATABLE, SAVE :: UU ( :,: ) ! returned solution REAL DFACP, DFACQ +! REAL :: DFSP( NLAYS ), DFSQ( NLAYS ) ! intermediate loop factors +! REAL DELC, DELP, RP, RQ REAL RP, RQ +! REAL :: LFAC1( NLAYS ) ! intermediate factor for CONVT +! REAL :: LFAC2( NLAYS ) ! intermediate factor for CONVT +! REAL :: LFAC3( NLAYS ) ! intermediate factor for eddy +! REAL :: LFAC4( NLAYS ) ! intermediate factor for eddy REAL, ALLOCATABLE, SAVE :: DEPVCR ( : ) ! dep vel in one cell + ! one cell for each landuse category REAL, ALLOCATABLE, SAVE :: EFAC1 ( : ) REAL, ALLOCATABLE, SAVE :: EFAC2 ( : ) REAL, ALLOCATABLE, SAVE :: POL ( : ) ! prodn/lossrate = PLDV/DEPV REAL PLDV_HONO ! PLDV for HONO REAL DEPV_NO2 ! dep vel of NO2 REAL DEPV_HNO3 ! dep vel of HNO3 - INTEGER, SAVE :: NO2_HIT, HONO_HIT, HNO3_HIT, NO2_MAP, HNO3_MAP - INTEGER, SAVE :: NH3_HIT +! REAL FNL ! ACM2 Variable +! INTEGER NLP, NL, LCBL + INTEGER, SAVE :: NO2_HIT = 0, HONO_HIT = 0, HNO3_HIT = 0, NO2_MAP= 0, HONO_MAP = 0, HNO3_MAP = 0 + INTEGER, SAVE :: O3_HIT = 0, O3_MAP = 0 + INTEGER, SAVE :: NH3_HIT = 0 +! REAL DTLIM, DTS, DTACM, RZ REAL DTS - - INTEGER, SAVE :: LOGDEV INTEGER ASTAT INTEGER C, R, L, S, V, I, J ! loop induction variables INTEGER MDATE, MTIME ! internal simulation date&time +!--Local Arrays for Z-coord implimentation + REAL :: DZH ( NLAYS ) ! ZF(L) - ZF(L-1) + REAL :: DZHI ( NLAYS ) ! 1/DZH + REAL :: DZFI ( NLAYS ) ! ZH(L+1) - ZH(L) + integer gl_c, gl_r + +#ifdef isam + REAL :: TOTAL_SA_NO2 + REAL, ALLOCATABLE, SAVE :: SA_DDBF( : ) + INTEGER IBGN, JSPCTAG + + REAL, ALLOCATABLE,SAVE :: SAEMIS( :,: ) + REAL, ALLOCATABLE,SAVE :: SACONC( :,: ) + REAL, ALLOCATABLE,SAVE :: SA_DD( :,: ) + REAL, ALLOCATABLE,SAVE :: SA_UU( :,: ) + + REAL, ALLOCATABLE,SAVE :: SAFRAC( : ) + REAL, ALLOCATABLE,SAVE :: SA_NO2( : ) + REAL, ALLOCATABLE,SAVE :: SA_SUM( : ) + + INTEGER, SAVE :: ISAM_INDEX_NO2 = 0 ! ...Index locating NO2 in ISAM + INTEGER, SAVE :: ISAM_INDEX_NH3 = 0 ! ...Index locating NH3 in ISAM + INTEGER, SAVE :: PLDV_INDEX_NH3 = 0 ! ...Index locating NH3 in PLDV + + INTEGER, ALLOCATABLE, SAVE :: ISAM_DEPV( : ) + INTEGER, ALLOCATABLE, SAVE :: INDEX_SA_HONO( : ) + INTEGER, ALLOCATABLE, SAVE :: INDEX_SA_NH3( : ) + + + CHARACTER( 16 ) :: ISAM_SPECIES + + INTEGER TOP, BOT + + REAL NH3_SUM +#endif + + +! INTERFACE +! SUBROUTINE MATRIX1 ( KL, A, B, E, D, X ) +! INTEGER, INTENT( IN ) :: KL +! REAL, INTENT( IN ) :: A( : ), B( : ), E( : ) +! REAL, INTENT( IN ) :: D( :,: ) +! REAL, INTENT( OUT ) :: X( :,: ) +! END SUBROUTINE MATRIX1 +! SUBROUTINE TRI ( L, D, U, B, X ) +! REAL, INTENT( IN ) :: L( : ), D( : ), U( : ) +! REAL, INTENT( IN ) :: B( :,: ) +! REAL, INTENT( OUT ) :: X( :,: ) +! END SUBROUTINE TRI + +#ifdef isam + SUBROUTINE SA_MATRIX1 ( KL, A, B, E, D, X ) + INTEGER, INTENT( IN ) :: KL + REAL, INTENT( IN ) :: A( : ), B( : ), E(: ) + REAL, INTENT( IN ) :: D( :,: ) + REAL, INTENT( OUT ) :: X( :,: ) + END SUBROUTINE SA_MATRIX1 + + SUBROUTINE SA_TRI ( L, D, U, B, X ) + REAL, INTENT( IN ) :: L( : ), D( : ), U( : ) + REAL, INTENT( IN ) :: B( :,: ) + REAL, INTENT( OUT ) :: X( :,: ) + END SUBROUTINE SA_TRI +#endif +! END INTERFACE + C----------------------------------------------------------------------- IF ( FIRSTIME ) THEN FIRSTIME = .FALSE. - LOGDEV = INIT3() MDATE = 0; MTIME = 0 @@ -150,48 +256,371 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) END IF CONC = 0.0; EMIS = 0.0 ! array assignment +! ALLOCATE ( DD( N_SPC_DIFF,NLAYS ), +! & UU( N_SPC_DIFF,NLAYS ), STAT = ASTAT ) +! IF ( ASTAT .NE. 0 ) THEN +! XMSG = 'Failure allocating DD or UU' +! CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) +! END IF +! DD = 0.0; UU = 0.0 ! array assignment + HONO_HIT = 0; HNO3_HIT = 0; NO2_HIT = 0; NH3_HIT = 0 - HNO3_MAP = 0; NO2_MAP = 0 + HONO_MAP = 0; HNO3_MAP = 0; NO2_MAP = 0 DO V = 1, N_SPC_DEPV IF ( DV2DF_SPC( V ) .EQ. 'NO2' ) THEN NO2_HIT = V NO2_MAP = DV2DF( V ) ELSE IF ( DV2DF_SPC( V ) .EQ. 'HONO' ) THEN HONO_HIT = V + HONO_MAP = DV2DF( V ) ELSE IF ( DV2DF_SPC( V ) .EQ. 'HNO3' ) THEN HNO3_HIT = V HNO3_MAP = DV2DF( V ) ELSE IF ( DV2DF_SPC( V ) .EQ. 'NH3' ) THEN NH3_HIT = V + ELSE IF ( DV2DF_SPC( V ) .EQ. 'O3' ) THEN + O3_HIT = V + O3_MAP = DV2DF( V ) + END IF + END DO + +#ifdef isam + ALLOCATE ( SA_DDBF( N_SPCTAG ), + & SACONC( N_SPCTAG, NLAYS ), + & SAEMIS( N_SPCTAG, NLAYS ), + & SA_DD ( N_SPCTAG, NLAYS ), + & SA_UU ( N_SPCTAG, NLAYS ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure ISAM diffusion variables' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + ALLOCATE ( SAFRAC ( N_SPCTAG ), + & SA_SUM ( NSPC_SA ), + & ISAM_DEPV( N_SPCTAG ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure ISAM depv variables' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE ( SA_NO2( NTAG_SA ), + & INDEX_SA_HONO( NTAG_SA ), + & INDEX_SA_NH3 ( NTAG_SA ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure ISAM HONO/NH3 variables' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + + SACONC = 0.0 + SAEMIS = 0.0 + SA_DD = 0.0 + SA_UU = 0.0 + + SAFRAC = 0.0 + ISAM_DEPV = 0 + SA_NO2 = 1.0 / REAL( NTAG_SA ) + INDEX_SA_HONO = 0 + INDEX_SA_NH3 = 0 + +! set default partitioning of surface fluxes + DO JSPCTAG = 1, N_SPCTAG + IF ( T_SPCTAG( JSPCTAG ) .EQ. OTHRTAG ) THEN + SAFRAC( JSPCTAG ) = 1.0 + ELSE + SAFRAC( JSPCTAG ) = 0.0 END IF END DO +! find NO2 in tracked species + DO S = 1, NSPC_SA + IF( ISAM_SPEC( S,1 ) .EQ. 'NO2' )THEN + ISAM_INDEX_NO2 = S + EXIT + END IF + END DO + +! find NH3 in tracked species + IF ( SA_BIDI ) THEN + DO S = 1, NSPC_SA + IF( ISAM_SPEC( S,1 ) .EQ. 'NH3' )THEN + ISAM_INDEX_NH3 = S + EXIT + END IF + END DO + S = -1 + S = INDEX1( 'NH3', N_SPC_DEPV, DEPV_SPC ) + IF ( S .GT. 0 ) THEN + PLDV_INDEX_NH3 = S + ELSE + XMSG = 'NH3 not found in DEPV_SPC array' + CALL M3EXIT( 'PLDV_INDEX_NH3', 0, 0, XMSG, XSTAT1 ) + END IF + END IF + +! set indices determining depv treatment, equals zero if none + ITAG = 0 + WRITE(LOGDEV,'(/,A7,1X,2(A16,1X))')'JSPCTAG','ISAM_SPECIES','DEPV Value' + DO JSPCTAG = 1, N_SPCTAG + ISAM_SPECIES = ISAM_SPEC( S_SPCTAG( JSPCTAG ),1 ) + IF( TRIM( ISAM_SPECIES ) .EQ. 'HONO' )THEN + ITAG = ITAG + 1 + INDEX_SA_HONO( ITAG ) = JSPCTAG + END IF + DO V = 1, N_SPC_DEPV + IF ( TRIM( ISAM_SPECIES ) .EQ. DV2DF_SPC( V ) ) THEN + ISAM_DEPV( JSPCTAG ) = V + END IF + END DO + END DO + WRITE(LOGDEV,'(/,A4,1X,A13,1X,A16)')'ITAG','INDEX_SA_HONO','ISAM_SPECIES' + DO ITAG = 1, NTAG_SA + JSPCTAG = INDEX_SA_HONO( ITAG ) + IF ( JSPCTAG .GT. 0 ) THEN + WRITE(LOGDEV,'(I2,3X,I4,8X,A16)')ITAG,JSPCTAG,VNAM_SPCTAG( JSPCTAG ) + ELSE + WRITE(LOGDEV,'(I2,3X,I4,8X,A16)')ITAG,JSPCTAG,'MISSING' + END IF + END DO +c WRITE(LOGDEV,* )'TAG_species, Default Partitioning Coeff.' +c DO JSPCTAG = 1, N_SPCTAG +c WRITE(LOGDEV,*)VNAM_SPCTAG( JSPCTAG ),' ,',SAFRAC( JSPCTAG ) +c END DO + IF ( ABFLUX .AND. L_NH4 ) THEN + ITAG = 0 + DO JSPCTAG = 1, N_SPCTAG + ISAM_SPECIES = ISAM_SPEC( S_SPCTAG( JSPCTAG ),1 ) + IF( TRIM( ISAM_SPECIES ) .EQ. 'NH3' )THEN + ITAG = ITAG + 1 + INDEX_SA_NH3( ITAG ) = JSPCTAG + END IF + END DO + +c DO ITAG = 1, NTAG_SA +c JSPCTAG = INDEX_SA_NH3( ITAG ) +c SAFRAC( JSPCTAG ) = 0.0 ! to not double count the bi-di emmissions +c END DO + END IF +#endif + +#ifdef sens + ALLOCATE ( S_POL ( NPMAX, N_SPC_DEPV ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure allocating S_POL' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + S_POL = 0.0 + + ALLOCATE ( S_DDBF( N_SPC_DEPV, NPMAX ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure allocating S_DBF' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + S_DDBF = 0.0 + + ALLOCATE ( SENS( N_SPC_DIFF,NLAYS,NPMAX ), + & S_EMIS( N_SPC_DIFF,NLAYS,NPMAX ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure allocating SENS or S_EMIS' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + SENS = 0.0; S_EMIS = 0.0 ! array assignment + + ALLOCATE ( S_DD( N_SPC_DIFF,NLAYS,NPMAX ), + & S_UU( N_SPC_DIFF,NLAYS,NPMAX ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure allocating S_DD or S_UU' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + S_DD = 0.0; S_UU = 0.0 ! array assignment + + ALLOCATE ( S_PLDV_HONO( NPMAX ), STAT = ASTAT ) + IF ( ASTAT .NE. 0 ) THEN + XMSG = 'Failure allocating S_PLDV_HONO' + CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + END IF + S_PLDV_HONO = 0.0 ! array assignment + + +#endif END IF ! if Firstime -C CALL GET_PT3D_EMIS () C ------------------------------------------- Row, Col LOOPS ----------- DTS = DTSEC + DO 345 R = 1, NROWS + DO 344 C = 1, NCOLS + + DZH(1) = Met_Data%ZF( C,R,1 ) + DZHI(1) = 1./DZH(1) + DO L = 2, NLAYS + DZH(L) = Met_Data%ZF( C,R,L ) - Met_Data%ZF( C,R,L-1 ) + DZHI(L) = 1./DZH(L) + ENDDO + DO L = 1, NLAYS - 1 + DZFI(L) = 1. / ( Met_Data%ZH( C,R,L+1 ) - Met_Data%ZH( C,R,L ) ) + ENDDO + DZFI(NLAYS) = DZFI(NLAYS-1) + +!C for ACM time step +! DTLIM = DTSEC +! +!C dt = .75 dzf*dzh / Kz +! DO L = 1, NLAYS - 1 +! DTLIM = MIN( DTLIM, 0.75 / ( SEDDY( L,C,R ) * DZHI(L)*DZFI(L) ) ) +! END DO +! MBARKS = 0.0 ! array assignment +! MDWN = 0.0 ! array assignment +! +!C conjoin ACM & EDDY --------------------------------------------------- +! +! MBAR = 0.0 +! FNL = 0.0 +! +! IF ( Met_Data%CONVCT( C,R ) ) THEN ! Do ACM for this column +! LCBL = Met_Data%LPBL( C,R ) +! MEDDY = SEDDY( 1,C,R ) * DZFI(1) / (Met_Data%PBL( C,R ) - Met_Data%ZF(C,R,1)) +! FNL = 1.0 / ( 1.0 + ( ( KARMAN / ( -Met_Data%HOL( C,R ) ) ) ** 0.3333 ) +! & / ( 0.72 * KARMAN ) ) +! MBAR = MEDDY * FNL +! IF ( MEDDY .LT. EPS ) THEN +! gl_c = c + COLSX_PE(1,mype+1) -1 +! gl_r = r + ROWSX_PE(1,mype+1) -1 +! WRITE( LOGDEV,* ) ' Warning --- MEDDY < 1e-6 s-1' +! WRITE( LOGDEV,* ) ' SEDDY, MEDDY, FNL, HOL = ', +! & SEDDY( 1,C,R ), MEDDY, FNL, Met_Data%HOL( C,R ) +! XMSG = '*** ACM fails ***' +! WRITE( LOGDEV,*)' c,r=', gl_c,gl_r,' pbl,ust=',Met_Data%PBL( C,R ),Met_Data%USTAR( C,R ) +! CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT2 ) +! END IF - DO 345 R = 1, MY_NROWS - DO 344 C = 1, MY_NCOLS +! IF ( ( FNL .LE. 0.0 ) .OR. ! never gonna happen for CONVCT +! & ( LCBL .GE. NLAYS-1 ) .OR. ! .GT. never gonna happen +! & ( Met_Data%HOL( C,R ) .GT. -0.00001 ) ) ! never gonna happen +! & WRITE( LOGDEV,1015 ) LCBL, MBAR, FNL, SEDDY( 1,C,R ), Met_Data%HOL( C,R ) +!1015 FORMAT( ' LCBL, MBAR, FNL, SEDDY1, HOL:', I3, 1X, 4(1PE13.5) ) +! +! DO L = 2, LCBL +! SEDDY( L,C,R ) = ( 1.0 - FNL ) * SEDDY( L,C,R ) +! MBARKS( L ) = MBAR +! MDWN( L ) = MBAR * (Met_Data%PBL( C,R ) - Met_Data%ZF(C,R,L-1)) * DZHI(L) +! END DO +! SEDDY( 1,C,R ) = ( 1.0 - FNL ) * SEDDY( 1,C,R ) +! MBARKS(1) = MBAR +! MBARKS(LCBL) = MDWN(LCBL) +! MDWN(LCBL+1) = 0.0 +!C Modify Timestep for ACM2 +! RZ = (Met_Data%ZF(C,R,LCBL) - Met_Data%ZF(C,R,1)) * DZHI(1) +! DTACM = 1.0 / ( MBAR * RZ ) +! DTLIM = MIN( 0.75 * DTACM, DTLIM ) +! ELSE +! LCBL = 1 +! END IF + +!C----------------------------------------------------------------------- + +! NLP = INT( DTSEC / DTLIM + 0.99 ) +! IF ( VDIFFDIAG ) NLPCR( C,R ) = REAL( NLP ) +! DTS = DTSEC / REAL( NLP ) DTDENS1 = DTS * Met_Data%DENS1( C,R ) DFACP = THETA * DTS DFACQ = THBAR * DTS +!#ifdef Verbose_Vdiff +! IF ( R .EQ. NROWS / 2 .AND. C .EQ. NCOLS / 2 ) +! & WRITE( LOGDEV,1021 ) Met_Data%CONVCT( C,R ), DTS, EDDYV( C,R,1 ), MBAR, FNL +!1021 FORMAT( ' CONVCT, DTS, EDDYV, MBAR, FNL: ', L3, 1X, 4(1PE13.5) ) +!#endif + DO L = 1, NLAYS DO V = 1, N_SPC_DIFF CONC( V,L ) = CNGRD( DIFF_MAP( V ),L,C,R ) END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SACONC( JSPCTAG,L ) = ISAM( C,R,L,S_SPCTAG( JSPCTAG ),T_SPCTAG( JSPCTAG ) ) + END DO +#endif END DO +#ifdef sens + DO NP = 1, NPMAX + DO L = 1, NLAYS + DO V = 1, N_SPC_DIFF + SENS( V,L,NP ) = SNGRD( DIFF_MAP( V ),L,C,R,NP ) + END DO + END DO + END DO +#endif + +#ifdef isam + SA_SUM = 0.0 + DO V = 1, NSPC_SA + DO ITAG = 1, NTAG_SA + SA_SUM( V ) = SA_SUM( V ) + ISAM( C,R,1,V,ITAG ) + END DO + SA_SUM( V ) = MAX ( 1.0E-30, SA_SUM( V ) ) + END DO + + SAFRAC = 0.0 + DO JSPCTAG = 1, N_SPCTAG + SAFRAC( JSPCTAG ) = SACONC( JSPCTAG,1 ) / SA_SUM ( S_SPCTAG( JSPCTAG ) ) + END DO + + IF ( ABFLUX .AND. L_NH4 ) THEN + DO ITAG = 1, NTAG_SA + JSPCTAG = INDEX_SA_NH3( ITAG ) + SAFRAC( JSPCTAG ) = 0.0 ! to not double count the bi-di emmissions + END DO + END IF + + IF( L_NO3 .AND. SFC_HONO ) THEN +! compute the flux partitioning for HONO from NO2 surface reaction + DO ITAG = 1, NTAG_SA + SA_NO2( ITAG ) = MAX( ISAM( C,R,1,ISAM_INDEX_NO2,ITAG ), 1.0E-30 ) + END DO + TOTAL_SA_NO2 = 1.0 / SUM( SA_NO2 ) + DO ITAG = 1, NTAG_SA + JSPCTAG = INDEX_SA_HONO( ITAG ) + SAFRAC( JSPCTAG ) = SA_NO2( ITAG ) * TOTAL_SA_NO2 + END DO + END IF +#endif + EMIS = 0.0 ! array assignment - DO L = 1, EMLAYS - DO V = 1, N_SPC_DIFF - EMIS( V,L ) = DTS * VDEMIS( DF2EM( V ),L,C,R ) + IF ( DESID_N_SRM .GE. 1 ) + & EMIS( :,1:DESID_LAYS ) = DTS * VDEMIS_DIFF( :,:,C,R ) + +#ifdef isam + SAEMIS = 0.0 + +! modify ground emissions for bidirectional species (for bidi, PLVD > 0.0) + IF ( SA_BIDI ) THEN + SA_VDEMIS_DIFF( ISAM_INDEX_NH3,1,C,R,BIDITAG ) + & = PLDV(PLDV_INDEX_NH3,C,R) * Met_Data%RDEPVHT( C,R ) + END IF + +! collapse ISAM emissions array + DO L = 1, DESID_LAYS + DO ITAG = 1, NTAG_SA + BOT = (ITAG-1)*NSPC_SA+1 + TOP = NSPC_SA*ITAG + SAEMIS( BOT:TOP,L ) = DTS * SA_VDEMIS_DIFF( :,L,C,R,ITAG ) END DO END DO +#endif + +#ifdef sens + S_EMIS = 0.0 + DO NP = 1, NPMAX + S_EMIS( :,1:DESID_LAYS,NP ) = DTS * SVDEMIS_DIFF( :,:,C,R,NP ) + END DO +#endif + +! DO L = 1, NLAYS +! DFSP( L ) = DFACP * DZHI( L ) +! DFSQ( L ) = DFACQ * DZHI( L ) +! EDDY( L ) = SEDDY( L,C,R ) * DZFI(L) +! END DO RP = DFACP * Met_Data%RDEPVHT( C,R ) RQ = DFACQ * Met_Data%RDEPVHT( C,R ) @@ -202,17 +631,64 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) EFAC1 ( V ) = EXP( -DEPVCR( V ) * RP ) EFAC2 ( V ) = EXP( -DEPVCR( V ) * RQ ) POL ( V ) = PLDV( V,C,R ) / DEPVCR( V ) - IF ( ABFLUX .AND. V .EQ. NH3_HIT ) THEN - DO I = 1, LCMP - CMPF( I ) = ICMP( I,C,R ) - END DO - END IF +#ifdef sens + DO NP = 1, NPMAX + S_POL ( NP,V ) = S_PLDV( NP,V,C,R ) / DEPVCR( V ) + S_DDBF( V, NP ) = S_DDEP( V, C, R, NP ) + END DO +#endif END DO PLDV_HONO = PLDV( HONO_HIT,C,R ) +#ifdef sens + DO NP = 1, NPMAX + S_PLDV_HONO( NP ) = S_PLDV( NP,HONO_HIT,C,R ) + END DO +#endif -C----------------------------------------------------------------------- +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DDBF( JSPCTAG ) = SA_DDEP( C,R,JSPCTAG ) + END DO +#endif - DO V = 1, N_SPC_DEPV +!C These don`t change in the NLP sub-time step loop:--------------------- +! DO L = 1, NLAYS +! AA ( L ) = 0.0 +! BB1( L ) = 0.0 +! EE1( L ) = 0.0 +! CC ( L ) = 0.0 +! EE2( L ) = 0.0 +! BB2( L ) = 0.0 +! END DO +! IF ( Met_Data%CONVCT( C,R ) ) THEN +! L = 1 +! DELP = Met_Data%PBL( C,R ) - Met_Data%ZF( C,R,L ) +! BB1( L ) = 1.0 + DELP * DFSP( L ) * MBARKS( L ) +! LFAC1( L ) = DFSQ( L ) * DELP * MBARKS( L ) +! LFAC2( L ) = DFSQ( L ) * MDWN( L+1 ) * DZH( L+1 ) +! DO L = 2, LCBL +! AA ( L ) = -DFACP * MBARKS( L ) +! BB1( L ) = 1.0 + DFACP * MDWN( L ) +! EE1( L ) = -DFSP( L-1 ) * DZH( L ) * MDWN( L ) +! MFAC( L ) = DZH( L+1 ) * DZHI( L ) * MDWN( L+1 ) +! END DO +! END IF +! +! DO L = 1, NLAYS +! EE2( L ) = - DFSP( L ) * EDDY( L ) +! LFAC3( L ) = DFSQ( L ) * EDDY( L ) +! END DO +! +! BB2( 1 ) = 1.0 - EE2( 1 ) +! DO L = 2, NLAYS +! CC ( L ) = - DFSP( L ) * EDDY( L-1 ) +! BB2( L ) = 1.0 - CC( L ) - EE2( L ) +! LFAC4( L ) = DFSQ( L ) * EDDY( L-1 ) +! END DO + +! DO 301 NL = 1, NLP ! loop over sub time + + DO V = 1, N_SPC_DEPV C --------- HET HONO RX ----------------- @@ -225,6 +701,12 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) DEPV_HNO3 = DEPVCR( V ) + PLDV_HONO / CONC( NO2_MAP,1 ) DD_FAC( V ) = DTDENS1 * DD_CONV( V ) * DEPV_HNO3 DDBF( V ) = DDBF( V ) + THETA * DD_FAC( V ) * CONC( S,1 ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC1( V ) + S_DDBF( V, NP ) = S_DDBF( V, NP ) + THETA * DD_FAC( V ) * SENS ( S,1,NP ) + END DO +#endif C Use special treatment for NO2 C Loss of NO2 via the heterogeneous reaction is accounted for as an additional @@ -239,32 +721,347 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) POL ( V ) = PLDV( V,C,R ) / DEPV_NO2 CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC1( V ) DDBF( V ) = DDBF( V ) + THETA * DD_FAC( V ) * CONC( S,1 ) +#ifdef sens + DO NP = 1, NPMAX + S_POL( NP, V ) = S_PLDV( NP,V,C,R ) / DEPV_NO2 + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC1( V ) + S_DDBF( V, NP ) = S_DDBF( V, NP ) + THETA * DD_FAC( V ) * SENS ( S,1,NP ) + END DO +#endif + + ELSE IF ( V .EQ. HONO_HIT ) THEN + S = HONO_MAP + CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC1( V ) + DDBF( V ) = DDBF( V ) + THETA * DD_FAC( V ) * CONC( S,1 ) +C Don't add HONO emissions as negative dep flux +! & - DTDENS1 * DD_CONV( V ) * PLDV( V,C,R ) ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC1( V ) + S_DDBF( V, NP ) = S_DDBF( V, NP ) + THETA * DD_FAC( V ) * SENS ( S,1,NP ) + END DO +#endif C --------- END of HET HONO RX ---------- ELSE + +C Pass selected N species to the BDSNP Soil NO emissions scheme + + IF ( MGN_ONLN_DEP ) THEN + + IF(SPECLOG) then + IF( V .eq. N_SPC_DEPV) THEN + SPECLOG = .false. ! no need to do any species more than once + WRITE( LOGDEV,*) 'BDSNP Species list complete', speclog + END IF + END IF + + IF ( (INDEX(TRIM( DV2DF_SPC( V ) ), 'NH3') .NE. 0) .OR. + & (INDEX(TRIM( DV2DF_SPC( V ) ), 'NH4') .NE. 0) .OR. + & (INDEX(TRIM( DV2DF_SPC( V ) ), 'HNO3').NE. 0) .OR. + & (INDEX(TRIM( DV2DF_SPC( V ) ), 'NO3') .NE. 0) .OR. + & (INDEX(TRIM( DV2DF_SPC( V ) ), 'NO2') .NE. 0) .OR. + & (INDEX(TRIM( DV2DF_SPC( V ) ), 'PAN') .NE. 0)) THEN + + + IF( SPECLOG ) THEN !write species each time it is used + WRITE( LOGDEV,*) 'BDSNP Dry Species Used:', TRIM(DV2DF_SPC( V ) ), V, N_SPC_DEPV + END IF + + IF ( ( DDBF(V)- DDEP( V,C,R) ) .LT. 0.0 ) THEN !negative error checking + + XMSG = 'Negative Deposition' +! WRITE( LOGDEV,*) 'BDSNP Negative Deposition vdiff, variable:', +! & TRIM( DV2DF_SPC( V )), ( DDBF(V)- DDEP( V,C,R) ), C, R +! CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) + CALL GET_N_DEP (DV2DF_SPC( V ), 0/ + & DTSEC, C, R ) + else + CALL GET_N_DEP (DV2DF_SPC( V ), ( DDBF(V)- DDEP( V,C,R) )/ + & DTSEC, C, R ) + END IF !end negative error checking + + + END IF !end species check + + END IF !end BDSNP check + S = DV2DF( V ) CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC1( V ) - DDBF( V ) = DDBF( V ) + THETA * DD_FAC( V ) * CONC( S,1 ) + DDBF( V ) = DDBF( V ) + THETA * ( DD_FAC( V ) * CONC( S,1 ) +C Add evasion as negative dep flux + & - DTDENS1 * DD_CONV( V ) * PLDV( V,C,R ) ) + +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC1( V ) + S_DDBF( V, NP ) = S_DDBF( V, NP ) + THETA * ( DD_FAC( V ) * SENS( S, 1,NP ) + & - DTDENS1 * DD_CONV( V ) * S_PLDV( NP,V,C,R ) ) + END DO +#endif - IF ( ABFLUX .AND. V .EQ. NH3_HIT ) THEN - DO I = 1, LCMP - CMPF( I ) = CMPF( I ) + THETA * CMP( I,C,R ) * DD_CONV( V ) * DTDENS1 - END DO - END IF END IF END DO -C --------- ADD EMISSIONS --------------- - DO L = 1, NLAYS DO V = 1, N_SPC_DIFF +! DD( V,L ) = 0.0 +! UU( V,L ) = 0.0 CONC( V,L ) = CONC( V,L ) + EMIS( V,L ) +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP) = 0.0 + S_UU( V,L,NP) = 0.0 + END DO +#endif + END DO + END DO + +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .GT. 0 ) THEN + SACONC( JSPCTAG,1 ) = SACONC( JSPCTAG,1 ) * EFAC1( S ) + & + SAFRAC( JSPCTAG ) * POL( S ) * ( 1.0 - EFAC1( S ) ) + SA_DDBF( JSPCTAG ) = SA_DDBF( JSPCTAG ) + & + THETA * DD_FAC( S ) * SACONC( JSPCTAG,1 ) + END IF + END DO + +c Recalculate bidi NH3 deposition + IF ( ABFLUX .AND. SA_BIDI ) THEN + NH3_SUM = 0.0 + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .EQ. NH3_HIT ) NH3_SUM = NH3_SUM + SACONC( JSPCTAG,1 ) + END DO + + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .EQ. NH3_HIT ) THEN + IF( NH3_SUM .GT. 1.0E-25 ) THEN + SA_DDBF( JSPCTAG ) = DDBF(NH3_HIT) * ( SACONC( JSPCTAG,1 ) / NH3_SUM ) + ELSE + SA_DDBF( JSPCTAG ) = 0.0 + END IF + END IF + END DO + END IF + + DO L = 1, NLAYS + DO V = 1, N_SPCTAG + SA_DD( V,L ) = 0.0 + SA_UU( V,L ) = 0.0 + END DO + END DO +#endif + +!C Compute tendency of CBL concentrations - semi-implicit solution +!C Set MATRIX1 elements A (col 1), B (diag.), E (superdiag.) and D (RHS) +! +! IF ( Met_Data%CONVCT( C,R ) ) THEN +! +! L = 1 +! DO V = 1, N_SPC_DIFF +! DD( V,L ) = CONC( V,L ) +! & - LFAC1( L ) * CONC( V,L ) +! & + LFAC2( L ) * CONC( V,L+1 ) +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP ) = SENS( V,L,NP ) + & - LFAC1( L ) * SENS( V,L,NP ) + & + LFAC2( L ) * SENS( V,L+1,NP ) + END DO +#endif +! END DO + +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DD( JSPCTAG, L) = SACONC( JSPCTAG,L ) + & - LFAC1( L ) * SACONC( JSPCTAG,L ) + & + LFAC2( L ) * SACONC( JSPCTAG,L+1 ) + ENDDO +#endif + +! DO L = 2, LCBL +! DO V = 1, N_SPC_DIFF +! DELC = MBARKS( L ) * CONC( V,1 ) +! & - MDWN( L ) * CONC( V,L ) +! & + MFAC( L ) * CONC( V,L+1 ) +! DD( V,L ) = CONC( V,L ) + DFACQ * DELC +#ifdef sens + DO NP = 1, NPMAX + S_DELC = MBARKS( L ) * SENS( V,1,NP ) + & - MDWN( L ) * SENS( V,L,NP ) + & + MFAC( L ) * SENS( V,L+1,NP ) + S_DD( V,L,NP ) = SENS( V,L,NP ) + DFACQ * S_DELC + END DO +#endif +! END DO + +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + DELC = MBARKS( L ) * SACONC( JSPCTAG,1 ) + & - MDWN( L ) * SACONC( JSPCTAG,L ) + & + MFAC( L ) * SACONC( JSPCTAG,L+1 ) + SA_DD( JSPCTAG,L ) = SACONC( JSPCTAG,L ) + DFACQ * DELC + END DO +#endif + +! END DO + +! CALL MATRIX1 ( LCBL, AA, BB1, EE1, DD, UU ) +#ifdef isam + CALL SA_MATRIX1( LCBL, AA, BB1,EE1, SA_DD, SA_UU) +#endif + +#ifdef sens + DO NP = 1, NPMAX + CALL MATRIX1 ( LCBL, AA, BB1, EE1, S_DD(:,:,NP), S_UU(:,:,NP) ) END DO +#endif + + +!C update conc +! DO L = 1, LCBL +! DO V = 1, N_SPC_DIFF +! CONC( V,L ) = UU( V,L ) +#ifdef sens + DO NP = 1, NPMAX + SENS( V,L,NP ) = S_UU( V,L,NP ) + END DO +#endif +! END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SACONC ( JSPCTAG,L ) = SA_UU( JSPCTAG,L ) + ENDDO +#endif +! END DO + +!C reinitialize for TRI solver +! DO L = 1, NLAYS +! DO V = 1, N_SPC_DIFF +! DD( V,L ) = 0.0 +! UU( V,L ) = 0.0 +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP ) = 0.0 + S_UU( V,L,NP ) = 0.0 + END DO +#endif +! END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DD( JSPCTAG,L ) = 0.0 + SA_UU( JSPCTAG,L ) = 0.0 + ENDDO +#endif +! END DO + +! END IF + +! L = 1 +! DO V = 1, N_SPC_DIFF +! DD( V,L ) = CONC( V,L ) +! & + LFAC3( L ) * ( CONC( V,L+1 ) - CONC( V,L ) ) +! & + EMIS( V,L ) +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP ) = SENS( V,L,NP ) + & + LFAC3( L ) * ( SENS( V,L+1,NP ) - SENS( V,L,NP ) ) + & + S_EMIS( V,L,NP ) + END DO +#endif +! END DO + +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DD( JSPCTAG,L ) = SACONC( JSPCTAG,L ) + & + LFAC3( L ) * ( SACONC( JSPCTAG,L+1 ) - SACONC( JSPCTAG,L ) ) + & + SAEMIS( JSPCTAG,L ) END DO +#endif + +! DO L = 2, NLAYS-1 +! DO V = 1, N_SPC_DIFF +! DD( V,L ) = CONC( V,L ) +! & + LFAC3( L ) * ( CONC( V,L+1 ) - CONC( V,L ) ) +! & - LFAC4( L ) * ( CONC( V,L ) - CONC( V,L-1 ) ) +! & + EMIS( V,L ) +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP ) = SENS( V,L,NP ) + & + LFAC3( L ) * ( SENS( V,L+1,NP ) - SENS( V,L,NP ) ) + & - LFAC4( L ) * ( SENS( V,L,NP ) - SENS( V,L-1,NP ) ) + & + S_EMIS( V,L,NP ) + END DO +#endif +! END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DD( JSPCTAG,L ) = SACONC( JSPCTAG,L ) + & + LFAC3( L ) * ( SACONC( JSPCTAG,L+1 ) - SACONC( JSPCTAG,L ) ) + & - LFAC4( L ) * ( SACONC( JSPCTAG,L ) - SACONC( JSPCTAG,L-1 ) ) + & + SAEMIS( JSPCTAG,L ) + END DO +#endif +! END DO + +! L = NLAYS +! DO V = 1, N_SPC_DIFF +! DD( V,L ) = CONC( V,L ) +! & - LFAC4( L ) * ( CONC( V,L ) - CONC( V,L-1 ) ) +#ifdef sens + DO NP = 1, NPMAX + S_DD( V,L,NP ) = SENS( V,L,NP ) + & - LFAC4( L ) * ( SENS( V,L,NP ) - SENS( V,L-1,NP ) ) + END DO +#endif +! END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DD( JSPCTAG,L ) = SACONC( JSPCTAG,L ) + & - LFAC4( L ) * ( SACONC( JSPCTAG,L ) - SACONC( JSPCTAG,L-1 ) ) + END DO +#endif + +! CALL TRI ( CC, BB2, EE2, DD, UU ) +#ifdef isam + CALL SA_TRI ( CC, BB2, EE2, SA_DD, SA_UU ) +#endif + +#ifdef sens + DO NP = 1, NPMAX + CALL TRI ( CC, BB2, EE2, S_DD(:,:,NP), S_UU(:,:,NP) ) + END DO +#endif + + +!C Load into CGRID +! DO L = 1, NLAYS +! DO V = 1, N_SPC_DIFF +! CONC( V,L ) = UU( V,L ) +#ifdef sens + DO NP = 1, NPMAX + SENS( V,L,NP ) = S_UU( V,L,NP ) + END DO +#endif +! END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SACONC( JSPCTAG,L ) = SA_UU( JSPCTAG,L ) + END DO +#endif +! END DO + + + + + -C --------- END EMISSIONS --------------- DO V = 1, N_SPC_DEPV @@ -275,43 +1072,136 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC2( V ) DDBF( V ) = DDBF( V ) + THBAR * DD_FAC( V ) * CONC( S,1 ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC2( V ) + S_DDBF( V,NP ) = S_DDBF( V,NP ) + THBAR * DD_FAC( V ) * SENS( S,1,NP ) + END DO +#endif + ELSE IF ( V .EQ. NO2_HIT ) THEN S = NO2_MAP CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC2( V ) DDBF( V ) = DDBF( V ) + THBAR * DD_FAC( V ) * CONC( S,1 ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC2( V ) + S_DDBF( V,NP ) = S_DDBF( V,NP ) + THBAR * DD_FAC( V ) * SENS( S,1,NP ) + END DO +#endif + + ELSE IF ( V .EQ. HONO_HIT ) THEN + S = HONO_MAP + CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC2( V ) + DDBF( V ) = DDBF( V ) + THBAR * DD_FAC( V ) * CONC( S,1 ) +C Don't add HONO emissions as negative dep flux +! & - DTDENS1 * DD_CONV( V ) * PLDV( V,C,R ) ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC2( V ) + S_DDBF( V,NP ) = S_DDBF( V,NP ) + THBAR * DD_FAC( V ) * SENS( S,1,NP ) + END DO +#endif + C --------- END of HET HONO RX ---------- ELSE S = DV2DF( V ) CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC2( V ) - DDBF( V ) = DDBF( V ) + THBAR * DD_FAC( V ) * CONC( S,1 ) - - IF ( ABFLUX .AND. V .EQ. NH3_HIT ) THEN - DO I = 1, LCMP - CMPF( I ) = CMPF( I ) + THBAR * CMP( I,C,R ) * DD_CONV( V ) * DTDENS1 - END DO - END IF + DDBF( V ) = DDBF( V ) + THBAR * ( DD_FAC( V ) * CONC( S,1 ) +C Add evasion as negative dep flux + & - DTDENS1 * DD_CONV( V ) * PLDV( V,C,R ) ) +#ifdef sens + DO NP = 1, NPMAX + SENS( S,1,NP ) = S_POL( NP,V ) + ( SENS( S,1,NP ) - S_POL( NP,V ) ) * EFAC2( V ) + S_DDBF( V,NP ) = S_DDBF( V,NP ) + THBAR * ( DD_FAC( V ) * SENS( S,1,NP ) + & - DTDENS1 * DD_CONV( V ) * S_PLDV( NP,V,C,R ) ) + END DO +#endif + END IF END DO +#ifdef isam +C Update ISAM Dry Deposition + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .GT. 0 ) THEN + SACONC( JSPCTAG,1 ) = SACONC( JSPCTAG,1 ) * EFAC2( S ) + & + SAFRAC( JSPCTAG ) * POL( S ) * ( 1.0 - EFAC2( S ) ) + SA_DDBF( JSPCTAG ) = SA_DDBF( JSPCTAG ) + & + THBAR * DD_FAC( S ) * SACONC( JSPCTAG,1 ) + END IF + END DO + +c Recalculate bidi NH3 deposition + IF ( ABFLUX .AND. SA_BIDI ) THEN + NH3_SUM = 0.0 + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .EQ. NH3_HIT ) NH3_SUM = NH3_SUM + SACONC( JSPCTAG,1 ) + END DO + + DO JSPCTAG = 1, N_SPCTAG + S = ISAM_DEPV( JSPCTAG ) + IF ( S .EQ. NH3_HIT ) THEN + IF( NH3_SUM .GT. 1.0E-25 ) THEN + SA_DDBF( JSPCTAG ) = DDBF(NH3_HIT) * ( SACONC( JSPCTAG,1 ) / NH3_SUM ) + ELSE + SA_DDBF( JSPCTAG ) = 0.0 + END IF + END IF + END DO + END IF +#endif + +301 CONTINUE ! end sub time loop + + DO L = 1, NLAYS DO V = 1, N_SPC_DIFF CNGRD( DIFF_MAP( V ),L,C,R ) = CONC( V,L ) END DO +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + IF( TRANSPORT_SPC( JSPCTAG ) )THEN + ISAM( C,R,L,S_SPCTAG( JSPCTAG ),T_SPCTAG( JSPCTAG ) ) = SACONC( JSPCTAG,L ) + END IF + END DO + +#endif END DO DO V = 1, N_SPC_DEPV DDEP( V,C,R ) = DDBF( V ) END DO - - IF ( ABFLUX ) THEN - DO I = 1, LCMP - ICMP( I,C,R ) = CMPF( I ) + +#ifdef isam + DO JSPCTAG = 1, N_SPCTAG + SA_DDEP( C,R,JSPCTAG ) = SA_DDBF( JSPCTAG ) + END DO +#endif + +#ifdef sens + DO NP = 1, NPMAX + + DO L = 1, NLAYS + DO V = 1, N_SPC_DIFF + SNGRD( DIFF_MAP( V ),L,C,R,NP ) = SENS( V,L,NP ) + END DO END DO - END IF + + DO V = 1, N_SPC_DEPV + S_DDEP( V,C,R,NP ) = S_DDBF( V,NP ) + END DO + + END DO +#endif + + 344 CONTINUE ! end loop on col C 345 CONTINUE ! end loop on row R