diff --git a/CHANGELOG.md b/CHANGELOG.md index 951e30d15..85655727c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -69,6 +69,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Prevent repeated printing of KPP integrate errors to the stdout stream. - Fixed selection of troposphere-stratosphere boundary in `global_ch4_mod.F90` - Removed operator splitting in CH4 simulation that was biasing diagnostics +- Fixed parallelization in Luo wetdep simulations caused by uninitialized variable ## [14.1.1] - 2023-03-03 ### Added diff --git a/GeosCore/apm_driv_mod.F90 b/GeosCore/apm_driv_mod.F90 index 419c1eb8a..8141715b6 100644 --- a/GeosCore/apm_driv_mod.F90 +++ b/GeosCore/apm_driv_mod.F90 @@ -1036,49 +1036,51 @@ SUBROUTINE APM_DRIV( Input_Opt, State_Chm, State_Diag, & !$OMP END PARALLEL DO ENDIF - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L, N ) & - !$OMP PRIVATE( SIZENUM, PRESS, TK, RHIN ) & - !$OMP PRIVATE( CACID,PACID ) & - !$OMP PRIVATE( MSO4,MSO4BULK,MNIT,MNH4,SOAT) & - !$OMP PRIVATE( MBCS, MOCS, MSULFT,MDSTS,MSALTS) & - !$OMP PRIVATE( MBC, MOC, MMSA) & - !$OMP PRIVATE( XMDST) & - !$OMP PRIVATE( MASS1, MASS2) & - !$OMP PRIVATE( CSO2,CNH3,XN0,CAMINE,CAMINEEMIT,YAMINEEMIT) & - !$OMP PRIVATE( CCO,CNO,CNO2,CNO3,CHNO3,CISOP,CMTPA) & - !$OMP PRIVATE( NH3EMIT) & - !$OMP PRIVATE( CSOG) & - !$OMP PRIVATE( CSOA) & - !$OMP PRIVATE( VOL) & - !$OMP PRIVATE( CLVSOG,MSULFLV,MBCLV,MOCLV,MDSTLV,MSALTLV) & - !$OMP PRIVATE( XM1D,XN1D,XNOLD,TEMPOUT1,ATOM4N,AEROCOMOUT1D) & - !$OMP PRIVATE( XQ,PLVSOG01,PLVSOG1,GFTOT1,GFTOT2,DENWET1,DENWET2) & - !$OMP PRIVATE( IACT10,IACT20,IACT30,FCLOUD1,AERAREA1,AERDRYR1,GAMMAPM1) & - !$OMP PRIVATE( RACT1,RACT2,RACT3) & - !$OMP PRIVATE( NCOAG1,NCOAG2) & - !$OMP PRIVATE( YSPGF,XBCLIFE,XOCLIFE,XCSNH3) &!Yu+ - !$OMP PRIVATE( XOH, XSINK,XAREA,XX0,XX1,DXX,ACS,XLAT, XLON,XAMINE) & - !$OMP PRIVATE( KYEAR,KMON,KDAY,KHOUR,KMIN,ISITE,JSITE,NSITE) & - !$OMP PRIVATE( XU,XV,TOP, TOPP) & - !$OMP PRIVATE( KKOUT) & - !$OMP PRIVATE( ZBEXT,ZW,ZG) &!OPT+ - !$OMP PRIVATE( ZBABS) &!OPT+ - !$OMP PRIVATE( YBEXT,XBEXT1k,YW,YG) &!OPT+ - !$OMP PRIVATE( IWL) &!OPT+ - !$OMP PRIVATE( ITYP) &!mxy+ - !$OMP PRIVATE( YCCN,YCDN,YCDNSP,YCLDF,YCLDLIQ,YCLDICE,YRCLDL,VZ) &!Yu+ 7/2012 - !$OMP PRIVATE( XCDN,XCDNSP ) & - !$OMP PRIVATE( YF,YC,SCOS,LOCALTIME ) & - !$OMP PRIVATE( PRESS0, YSIGMA ) & - !$OMP PRIVATE( wbar,relhum,yqc,yna,YB,YREI,YK) & - !$OMP PRIVATE( dumc, dumnc, pgam, lamc) & - !$OMP PRIVATE( CCLD,CLDLIQ,CLDICE ) & - !$OMP PRIVATE( REL,REI) & - !$OMP PRIVATE( taucloud, taucloudl, taucloudi, ssacloudl, ssacloudi ) & - !!$OMP PRIVATE( nuci, onihf, oniimm, onidep, onimey) & - !$OMP SCHEDULE( DYNAMIC ) +! Disable this parallel loop, which causes differences in output +! for the time being. -- Bob Yantosca (24 May 2023) +! !$OMP PARALLEL DO & +! !$OMP DEFAULT( SHARED ) & +! !$OMP PRIVATE( I, J, L, N ) & +! !$OMP PRIVATE( SIZENUM, PRESS, TK, RHIN ) & +! !$OMP PRIVATE( CACID,PACID ) & +! !$OMP PRIVATE( MSO4,MSO4BULK,MNIT,MNH4,SOAT) & +! !$OMP PRIVATE( MBCS, MOCS, MSULFT,MDSTS,MSALTS) & +! !$OMP PRIVATE( MBC, MOC, MMSA) & +! !$OMP PRIVATE( XMDST) & +! !$OMP PRIVATE( MASS1, MASS2) & +! !$OMP PRIVATE( CSO2,CNH3,XN0,CAMINE,CAMINEEMIT,YAMINEEMIT) & +! !$OMP PRIVATE( CCO,CNO,CNO2,CNO3,CHNO3,CISOP,CMTPA) & +! !$OMP PRIVATE( NH3EMIT) & +! !$OMP PRIVATE( CSOG) & +! !$OMP PRIVATE( CSOA) & +! !$OMP PRIVATE( VOL) & +! !$OMP PRIVATE( CLVSOG,MSULFLV,MBCLV,MOCLV,MDSTLV,MSALTLV) & +! !$OMP PRIVATE( XM1D,XN1D,XNOLD,TEMPOUT1,ATOM4N,AEROCOMOUT1D) & +! !$OMP PRIVATE( XQ,PLVSOG01,PLVSOG1,GFTOT1,GFTOT2,DENWET1,DENWET2) & +! !$OMP PRIVATE( IACT10,IACT20,IACT30,FCLOUD1,AERAREA1,AERDRYR1,GAMMAPM1) & +! !$OMP PRIVATE( RACT1,RACT2,RACT3) & +! !$OMP PRIVATE( NCOAG1,NCOAG2) & +! !$OMP PRIVATE( YSPGF,XBCLIFE,XOCLIFE,XCSNH3) &!Yu+ +! !$OMP PRIVATE( XOH, XSINK,XAREA,XX0,XX1,DXX,ACS,XLAT, XLON,XAMINE) & +! !$OMP PRIVATE( KYEAR,KMON,KDAY,KHOUR,KMIN,ISITE,JSITE,NSITE) & +! !$OMP PRIVATE( XU,XV,TOP, TOPP) & +! !$OMP PRIVATE( KKOUT) & +! !$OMP PRIVATE( ZBEXT,ZW,ZG) &!OPT+ +! !$OMP PRIVATE( ZBABS) &!OPT+ +! !$OMP PRIVATE( YBEXT,XBEXT1k,YW,YG) &!OPT+ +! !$OMP PRIVATE( IWL) &!OPT+ +! !$OMP PRIVATE( ITYP) &!mxy+ +! !$OMP PRIVATE( YCCN,YCDN,YCDNSP,YCLDF,YCLDLIQ,YCLDICE,YRCLDL,VZ) &!Yu+ 7/2012 +! !$OMP PRIVATE( XCDN,XCDNSP ) & +! !$OMP PRIVATE( YF,YC,SCOS,LOCALTIME ) & +! !$OMP PRIVATE( PRESS0, YSIGMA ) & +! !$OMP PRIVATE( wbar,relhum,yqc,yna,YB,YREI,YK) & +! !$OMP PRIVATE( dumc, dumnc, pgam, lamc) & +! !$OMP PRIVATE( CCLD,CLDLIQ,CLDICE ) & +! !$OMP PRIVATE( REL,REI) & +! !$OMP PRIVATE( taucloud, taucloudl, taucloudi, ssacloudl, ssacloudi ) & +! !!$OMP PRIVATE( nuci, onihf, oniimm, onidep, onimey) & +! !$OMP SCHEDULE( DYNAMIC ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -2139,7 +2141,9 @@ SUBROUTINE APM_DRIV( Input_Opt, State_Chm, State_Diag, & ENDDO ENDDO - !$OMP END PARALLEL DO +! Disable this parallel loop, which causes differences in output. +! -- Bob Yantosca (24 May 2023) +! !$OMP END PARALLEL DO write(*,*)'LuoSSA',sum(TCOD3D(:,:,:,1))/size(TCOD3D(:,:,:,1)), & sum(TCOD3D(:,:,:,5))/size(TCOD3D(:,:,:,5)) diff --git a/GeosCore/carbon_mod.F90 b/GeosCore/carbon_mod.F90 index 6706be823..b173ec1d2 100644 --- a/GeosCore/carbon_mod.F90 +++ b/GeosCore/carbon_mod.F90 @@ -1800,20 +1800,35 @@ SUBROUTINE SOA_CHEMISTRY( Input_Opt, State_Chm, State_Diag, & ! add KRO2XXX to private (hotp 5/7/10) ! add JSV to private (hotp 5/13/10) ! remove NOX (hotp 5/22/10) - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L, JHC, IPR, GM0, AM0 ) & + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& #ifdef APM - !$OMP PRIVATE( N, NTEMP) & + !$OMP PRIVATE( IFINORG, OCBIN_SUM, N, NTEMP )& #endif - !$OMP PRIVATE( VOL, FAC, RTEMP, KO3, KOH, KNO3, CAIR ) & - !$OMP PRIVATE( MPRODUCT, MSOA_OLD, VALUE, UPPER, LOWER, MNEW, TOL ) & - !$OMP PRIVATE( ORG_AER, ORG_GAS, KOM, MPOC ) & - !$OMP PRIVATE( KRO2NO, KRO2HO2, JSV ) + !$OMP PRIVATE( I, J, L, JHC, IPR, GM0, AM0 )& + !$OMP PRIVATE( VOL, FAC, RTEMP, KO3, KOH, KNO3, CAIR )& + !$OMP PRIVATE( VALUE, UPPER, LOWER, MNEW, TOL )& + !$OMP PRIVATE( ORG_AER, ORG_GAS, KOM, MPOC )& + !$OMP PRIVATE( KRO2NO, KRO2HO2, JSV ) DO L = 1, State_Grid%MaxChemLev DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero important variables at top of loop + AM0 = 0.0_fp + GM0 = 0.0_fp + KO3 = 0.0_fp + KOH = 0.0_fp + KNO3 = 0.0_fp + KRO2NO = 0.0_fp + KRO2HO2 = 0.0_fp + LOWER = 0.0_fp + MNEW = 0.0_fp + MPOC = 0.0_fp + TOL = 0.0_fp + UPPER = 0.0_fp + VALUE = 0.0_fp + ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE @@ -1821,7 +1836,7 @@ SUBROUTINE SOA_CHEMISTRY( Input_Opt, State_Chm, State_Diag, & VOL = State_Met%AIRVOL(I,J,L) ! conversion factor from kg to ug/m3 - FAC = 1.e+9_fp / VOL + FAC = 1.0e+9_fp / VOL ! air conc. in kg/m3 CAIR = State_Met%AD(I,J,L) / VOL @@ -1832,7 +1847,7 @@ SUBROUTINE SOA_CHEMISTRY( Input_Opt, State_Chm, State_Diag, & ! Get SOA yield parameters ! ALPHA is a module variable now. (ccc, 2/2/10) ! add arguments for RO2+NO, RO2+HO2 rates (hotp 5/7/10) - CALL SOA_PARA( RTEMP, KO3, KOH, KNO3, KOM, & + CALL SOA_PARA( RTEMP, KO3, KOH, KNO3, KOM, & I, J, L, KRO2NO, KRO2HO2, State_Met ) ! Partition mass of gas & aerosol species @@ -1858,8 +1873,8 @@ SUBROUTINE SOA_CHEMISTRY( Input_Opt, State_Chm, State_Diag, & ! Initialize other arrays to be safe (dkh, 11/10/06) ! update dims (hotp 5/22/10) - ORG_AER(:,:) = 0e+0_fp - ORG_GAS(:,:) = 0e+0_fp + ORG_AER = 0.0_fp + ORG_GAS = 0.0_fp ! Individual SOA's: convert from [kg] to [ug/m3] or [kgC] to [ugC/m3] DO JSV = 1, MAXSIMSV @@ -1913,85 +1928,95 @@ SUBROUTINE SOA_CHEMISTRY( Input_Opt, State_Chm, State_Diag, & ! (Colette Heald, 12/3/09) !----------------------------------------------------------- #ifdef APM - OCBIN_SUM = 0.e+0_fp + OCBIN_SUM = 0.0_fp DO N = 1, NBCOC - OCBIN_SUM = OCBIN_SUM + SUM(Spc(APMIDS%id_OCBIN1+N-1)%Conc(:,:,:)) + OCBIN_SUM = OCBIN_SUM + SUM( Spc(APMIDS%id_OCBIN1+N-1)%Conc ) ENDDO - MPOC = FAC * OCBIN_SUM + MPOC = FAC * OCBIN_SUM MPOC = MPOC * 2.1d0 IFINORG = 2 - IF(IFINORG.EQ.1) THEN !Yu+ consider inorg in the SOA partition - - IF ( APMIDS%id_SO4 > 0 .and. & - APMIDS%id_NH4 > 0 .and. & - APMIDS%id_NIT > 0 ) THEN - ! Then compute SOG condensation onto SO4, NH4, NIT aerosols - MPOC = MPOC + ( Spc(APMIDS%id_NH4)%Conc(I,J,L) + & - Spc(APMIDS%id_NIT)%Conc(I,J,L) ) * FAC - - IF(NSO4>=1)THEN - NTEMP=APMIDS%id_SO4BIN1-1 - DO N=(NTEMP+1),(NTEMP+NSO4) - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - IF(NCTSO4>=1)THEN - NTEMP=APMIDS%id_CTSO4-1 - DO N=(NTEMP+1),(NTEMP+NCTSO4) - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - IF(NCTBC>=1)THEN - NTEMP=APMIDS%id_CTBC-1 - DO N=(NTEMP+1),(NTEMP+NCTBC) - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - IF(NCTOC>=1)THEN - NTEMP=APMIDS%id_CTOC-1 - DO N=(NTEMP+1),(NTEMP+NCTOC) - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - IF(NCTSEA>=1)THEN - NTEMP=APMIDS%id_CTSEA-1 - DO N=(NTEMP+1),(NTEMP+NCTSEA) - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - IF(NSEA>=1)THEN - NTEMP=APMIDS%id_SEABIN1-1 - DO N=(NTEMP+1),(NTEMP+16) ! SALTbin16 = 1 um - MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC - ENDDO - ENDIF - - !Add MSA - MPOC = MPOC + Spc(APMIDS%id_MSA)%Conc(I,J,L) * FAC - - ENDIF - - ELSEIF(IFINORG.EQ.2) THEN !Consider SV SOA partition on LV SOA - - MPOC = MPOC + FAC * (Spc(APMIDS%id_CTSO4)%Conc(I,J,L) + & !MSULFLV - Spc(APMIDS%id_CTBC+1)%Conc(I,J,L) + & !MBCLV - Spc(APMIDS%id_CTOC+1)%Conc(I,J,L) + & !MOCLV - Spc(APMIDS%id_CTDST+1)%Conc(I,J,L) + & !MDSTLV - Spc(APMIDS%id_CTSEA+1)%Conc(I,J,L)) !MSALTLV - - ELSE - - ! Compute SOG condensation onto OC aerosol - MPOC = ( Spc(id_OCPI)%Conc(I,J,L) + Spc(id_OCPO)%Conc(I,J,L) ) * FAC - MPOC = MPOC * 2.1d0 - - ENDIF +!############################################################################# +! NOTE: IFINORG is always 2 so the other IF branches never get done. +! Comment these out for now for better numerical efficiency +! -- Bob Yantosca (23 May 2023) +! IF(IFINORG.EQ.1) THEN !Yu+ consider inorg in the SOA partition +! +! IF ( APMIDS%id_SO4 > 0 .and. & +! APMIDS%id_NH4 > 0 .and. & +! APMIDS%id_NIT > 0 ) THEN +! ! Then compute SOG condensation onto SO4, NH4, NIT aerosols +! MPOC = MPOC + ( Spc(APMIDS%id_NH4)%Conc(I,J,L) + & +! Spc(APMIDS%id_NIT)%Conc(I,J,L) ) * FAC +! +! IF(NSO4>=1)THEN +! NTEMP=APMIDS%id_SO4BIN1-1 +! DO N=(NTEMP+1),(NTEMP+NSO4) +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! IF(NCTSO4>=1)THEN +! NTEMP=APMIDS%id_CTSO4-1 +! DO N=(NTEMP+1),(NTEMP+NCTSO4) +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! IF(NCTBC>=1)THEN +! NTEMP=APMIDS%id_CTBC-1 +! DO N=(NTEMP+1),(NTEMP+NCTBC) +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! IF(NCTOC>=1)THEN +! NTEMP=APMIDS%id_CTOC-1 +! DO N=(NTEMP+1),(NTEMP+NCTOC) +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! IF(NCTSEA>=1)THEN +! NTEMP=APMIDS%id_CTSEA-1 +! DO N=(NTEMP+1),(NTEMP+NCTSEA) +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! IF(NSEA>=1)THEN +! NTEMP=APMIDS%id_SEABIN1-1 +! DO N=(NTEMP+1),(NTEMP+16) ! SALTbin16 = 1 um +! MPOC = MPOC + Spc(N)%Conc(I,J,L) * FAC +! ENDDO +! ENDIF +! +! !Add MSA +! MPOC = MPOC + Spc(APMIDS%id_MSA)%Conc(I,J,L) * FAC +! +! ENDIF +! +! ELSEIF(IFINORG.EQ.2) THEN !Consider SV SOA partition on LV SOA +!############################################################################# + + MPOC = MPOC + FAC * (Spc(APMIDS%id_CTSO4 )%Conc(I,J,L) + & !MSULFLV + Spc(APMIDS%id_CTBC +1)%Conc(I,J,L) + & !MBCLV + Spc(APMIDS%id_CTOC +1)%Conc(I,J,L) + & !MOCLV + Spc(APMIDS%id_CTDST+1)%Conc(I,J,L) + & !MDSTLV + Spc(APMIDS%id_CTSEA+1)%Conc(I,J,L)) !MSALTLV + +!############################################################################# +! NOTE: IFINORG is always 2 so the other IF branches never get done. +! Comment these out for now for better numerical efficiency +! -- Bob Yantosca (23 May 2023) +! ELSE +! +! ! Compute SOG condensation onto OC aerosol +! MPOC = ( Spc(id_OCPI)%Conc(I,J,L) + Spc(id_OCPO)%Conc(I,J,L) ) * FAC +! MPOC = MPOC * 2.1d0 +! +! ENDIF +!############################################################################# #else ! Now treat either traditional POA or semivolatile POA (hotp 7/25/10) IF ( id_OCPI > 0 .and. id_OCPO > 0 ) THEN @@ -8357,7 +8382,7 @@ END SUBROUTINE CHEM_OCPINEW !------------------------------------------------------------------------------ !BOP ! -! !IROUTINE: dry_settlingbin +! !IROUTINE: bcdry_settlingbin ! ! !DESCRIPTION: Subroutine DRY\_SETTLINGBIN computes the dry settling of ! aerosol tracers. Modified for APM simulation. (G. Luo) @@ -8466,6 +8491,23 @@ SUBROUTINE BCDRY_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero private loop variables + CONST = 0.0_fp + DELZ = 0.0_fp + DELZ1 = 0.0_fp + DEN = 0.0_fp + DP = 0.0_fp + MASS = 0.0_fp + OLD = 0.0_fp + P = 0.0_fp + PDP = 0.0_fp + REFF = 0.0_fp + SLIP = 0.0_fp + TEMP = 0.0_fp + TC0 = 0.0_fp + VISC = 0.0_fp + VTS = 0.0_fp + DO L = 1, State_Grid%NZ DO N = APMIDS%id_BCBIN1, IDTEMP MASS(L) = Spc(N)%Conc(I,J,L) @@ -8608,7 +8650,7 @@ END SUBROUTINE BCDRY_SETTLINGBIN !------------------------------------------------------------------------------ !BOP ! -! !IROUTINE: dry_settlingbin +! !IROUTINE: ocdry_settlingbin ! ! !DESCRIPTION: Subroutine DRY\_SETTLINGBIN computes the dry settling of ! aerosol tracers. Modified for APM simulation. (G. Luo) @@ -8717,6 +8759,23 @@ SUBROUTINE OCDRY_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero private loop variables + CONST = 0.0_fp + DELZ = 0.0_fp + DELZ1 = 0.0_fp + DEN = 0.0_fp + DP = 0.0_fp + MASS = 0.0_fp + OLD = 0.0_fp + P = 0.0_fp + PDP = 0.0_fp + REFF = 0.0_fp + SLIP = 0.0_fp + TEMP = 0.0_fp + TC0 = 0.0_fp + VISC = 0.0_fp + VTS = 0.0_fp + DO L = 1, State_Grid%NZ DO N = APMIDS%id_OCBIN1, IDTEMP MASS(L) = MASS(L) + Spc(N)%Conc(I,J,L) @@ -8855,4 +8914,4 @@ SUBROUTINE OCDRY_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & END SUBROUTINE OCDRY_SETTLINGBIN #endif !EOC - END MODULE CARBON_MOD +END MODULE CARBON_MOD diff --git a/GeosCore/drydep_mod.F90 b/GeosCore/drydep_mod.F90 index aba461db8..f18bc0fc5 100644 --- a/GeosCore/drydep_mod.F90 +++ b/GeosCore/drydep_mod.F90 @@ -903,14 +903,15 @@ SUBROUTINE METERO( State_Grid, State_Met, CZ1, TC0, OBK, CFRAC, & REAL(f8) :: SP REAL(f8) :: SFCWINDSQR - !================================================================= + !======================================================================== ! METERO begins here! - !================================================================= + !======================================================================== ! Loop over surface grid boxes - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, THIK, SP, SFCWINDSQR ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, THIK, SP, SFCWINDSQR )& + !$OMP COLLAPSE( 2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -928,7 +929,7 @@ SUBROUTINE METERO( State_Grid, State_Met, CZ1, TC0, OBK, CFRAC, & PRESSU(I,J) = SP * 1.e+2_f8 !============================================================== - ! Return meterological quantities as 1-D arrays for DEPVEL + ! Return meterological quantities for DEPVEL !============================================================== ! Roughness height [m] @@ -957,7 +958,7 @@ SUBROUTINE METERO( State_Grid, State_Met, CZ1, TC0, OBK, CFRAC, & ! changed to 98% due to vapor pressure lowering above sea water ! (Lewis & Schwartz, 2004) ! jaegle (5/11/11) - RHB(I,J) = MIN( 0.98e+0_f8, State_Met%RH(I,J,1) * 1.e-2_f8 ) + RHB(I,J) = MIN( 0.98_f8, State_Met%RH(I,J,1) * 1.0e-2_f8 ) ! 10m windspeed [m/s] (jaegle 5/11/11) SFCWINDSQR = State_Met%U10M(I,J)**2 + State_Met%V10M(I,J)**2 @@ -1163,7 +1164,7 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(f8) :: RDC,RLUXX,RGSX,RCL,DTMP1,DTMP2,DTMP3,DTMP4 REAL(f8) :: CZH,CKUSTR,REYNO,CORR1,CORR2,Z0OBK REAL(f8) :: RA,RB,DUMMY1,DUMMY2,DUMMY3,DUMMY4 - REAL(f8) :: XMWH2O,DAIR,TEMPK,TEMPC + REAL(f8) :: XMWH2O,DAIR,TEMPK,TEMPC,F0_K INTEGER :: IOLSON,II,IW INTEGER :: K,LDT REAL(f8) :: RCLX,RIXX !,BIOFIT @@ -1325,48 +1326,90 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & ENDIF #endif - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( CZ, TEMPK, TEMPC, K, VD ) & - !$OMP PRIVATE( LDT, RSURFC, C1, XNU, RT, IOLSON ) & - !$OMP PRIVATE( II, RI, RLU, RAC, RGSS, RGSO ) & - !$OMP PRIVATE( RCLS, RCLO, RAD0, RIX, GFACT, GFACI ) & - !$OMP PRIVATE( RDC, XMWH2O, RIXX, RLUXX, RGSX, RCLX ) & - !$OMP PRIVATE( DTMP1, DTMP2, DTMP3, DTMP4, VDS, CZH ) & - !$OMP PRIVATE( CKUSTR, REYNO, CORR1, CORR2, Z0OBK, RA ) & - !$OMP PRIVATE( Z0OBK_Alt, RA_Alt,DUMMY2_Alt, DUMMY4_Alt ) & - !$OMP PRIVATE( DUMMY1, DUMMY2, DUMMY3, DUMMY4, DAIR, RB ) & - !$OMP PRIVATE( C1X, VK, I, J, IW ) & - !$OMP PRIVATE( DIAM, DEN, XLAI_FP, SUNCOS_FP, CFRAC_FP ) & - !$OMP PRIVATE( N_SPC, alpha, DEPVw ) & + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( CZ, TEMPK, TEMPC, K, VD )& + !$OMP PRIVATE( LDT, RSURFC, C1, XNU, RT, IOLSON )& + !$OMP PRIVATE( II, RI, RLU, RAC, RGSS, RGSO )& + !$OMP PRIVATE( RCLS, RCLO, RAD0, RIX, GFACT, GFACI )& + !$OMP PRIVATE( RDC, XMWH2O, RIXX, RLUXX, RGSX, RCLX )& + !$OMP PRIVATE( DTMP1, DTMP2, DTMP3, DTMP4, VDS, CZH )& + !$OMP PRIVATE( CKUSTR, REYNO, CORR1, CORR2, Z0OBK, RA )& + !$OMP PRIVATE( Z0OBK_Alt, RA_Alt, DUMMY2_Alt, DUMMY4_Alt )& + !$OMP PRIVATE( DUMMY1, DUMMY2, DUMMY3, DUMMY4, DAIR, RB )& + !$OMP PRIVATE( C1X, VK, I, J, IW )& + !$OMP PRIVATE( DIAM, DEN, XLAI_FP, SUNCOS_FP, CFRAC_FP )& + !$OMP PRIVATE( N_SPC, alpha, DEPVw )& #ifdef TOMAS - !$OMP PRIVATE( BIN ) & + !$OMP PRIVATE( BIN )& #endif - !$OMP PRIVATE( SpcId, SpcInfo ) + !$OMP PRIVATE( SpcId, SpcInfo, F0_K )& + !$OMP COLLAPSE( 2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX - ! Initialize + ! Zero variables that aren't zeroed below + VD = 0.0_f8 + RSURFC = 0.0_f8 + II = 0.0 + RI = 0.0_f8 + RLU = 0.0_f8 + RAC = 0.0_f8 + RGSS = 0.0_f8 + RGSO = 0.0_f8 + RCLS = 0.0_f8 + RCLO = 0.0_f8 + RAD0 = 0.0_f8 + RIX = 0.0_f8 + GFACT = 0.0_f8 + GFACI = 0.0_f8 + RDC = 0.0_f8 + XMWH2O = 0.0_f8 + RIXX = 0.0_f8 + RLUXX = 0.0_f8 + RGSX = 0.0_f8 + RCLX = 0.0_f8 + DTMP1 = 0.0_f8 + DTMP2 = 0.0_f8 + DTMP3 = 0.0_f8 + DTMP4 = 0.0_f8 + VDS = 0.0_f8 + CZH = 0.0_f8 + CKUSTR = 0.0_f8 + REYNO = 0.0_f8 + CORR1 = 0.0_f8 + CORR2 = 0.0_f8 + Z0OBK = 0.0_f8 + RA = 0.0_f8 Z0OBK_Alt = 0.0_f8 + RA_Alt = 0.0_f8 DUMMY2_Alt = 0.0_f8 DUMMY4_Alt = 0.0_f8 - RA_Alt = 0.0_f8 + DUMMY1 = 0.0_f8 + DUMMY2 = 0.0_f8 + DUMMY3 = 0.0_f8 + DUMMY4 = 0.0_f8 + DAIR = 0.0_f8 + RB = 0.0_f8 + C1X = 0.0_f8 + VK = 0.0_f8 + DIAM = 0.0_f8 + DEN = 0.0_f8 + XLAI_FP = 0.0_f8 + SUNCOS_FP = 0.0_f8 + CFRAC_FP = 0.0_f8 + N_SPC = 0 + alpha = 0.0_f8 + DEPVw = 0.0_f8 + F0_K = 0.0_f8 !** CZ is Altitude (m) at which deposition velocity is computed - CZ = CZ1(I,J) + CZ = CZ1(I,J) !** TEMPK and TEMPC are surface air temperatures in K and in C TEMPK = TEMP(I,J) TEMPC = TEMP(I,J)-273.15e+0_f8 - !* Initialize variables - DO K = 1,NUMDEP - VD(K) = 0.0e+0_f8 - DO LDT = 1,NTYPE - RSURFC(K,LDT) = 0.e+0_f8 - END DO - END DO - !** Calculate the kinematic viscosity XNU (m2 s-1) of air !** as a function of temperature. !** The kinematic viscosity is used to calculate the roughness heights @@ -1375,7 +1418,7 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & !** The expression for the temperature dependence of XNU !** is from the FORTRAN code in Appendix II of Wesely [1988]; !** I wasn't able to find an original reference but it seems benign enough. - C1 = TEMPK/273.15e+0_f8 + C1 = TEMPK/273.15e+0_f8 XNU = 0.151e+0_f8*(C1**1.77e+0_f8)*1.0e-04_f8 !* Compute bulk surface resistance for gases. @@ -1534,6 +1577,10 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & !* DO 160 K = 1,NUMDEP + ! Save F0(K) in a private variable to avoid diffs due + ! to parallelization -- Bob Yantosca (17 May 2023) + F0_K = F0(K) + !** exit for non-depositing species or aerosols. IF (.NOT. LDEP(K) .OR. AIROSOL(K)) GOTO 155 @@ -1568,18 +1615,27 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & ELSE IF ((N_SPC .EQ. ID_O3) .AND. (State_Met%isSnow(I,J))) THEN RSURFC(K,LDT) = 10000.0_f8 ELSE - ! Check latitude and longitude, alter F0 only for Amazon rainforest for Hg0 - ! (see reference: Feinberg et al., ESPI, 2022: Evaluating atmospheric mercury (Hg) uptake by vegetation in a - ! chemistry-transport model) - IF (N_SPC .EQ. ID_Hg0) THEN ! Check for Hg0 - IF ( II .EQ. 6 .AND. & ! if rainforest land type - State_Grid%XMid(I,J) > -82 .AND. & ! bounding box of Amazon - State_Grid%XMid(I,J) < -33 .AND. & - State_Grid%YMid(I,J) > -34 .AND. & - State_Grid%YMid(I,J) < 14 ) THEN - F0(K) = 2.0e-01_f8 ! increase reactivity, as inferred from observations - ELSE - F0(K) = 3.0e-05_f8 ! elsewhere, lower reactivity + ! Check latitude and longitude, alter F0 only for Amazon + ! rainforest for Hg0 (see reference: Feinberg et al., ESPI, + ! 2022: Evaluating atmospheric mercury (Hg) uptake by + ! vegetation in a chemistry-transport model) + ! + ! Remove IF/ELSE block using never-nesting technique. + ! - Bob Yantosca (17 May 2023) + IF ( N_SPC == ID_Hg0 ) THEN + + ! Assume lower reactivity + F0_K = 3.0e-05_f8 + + ! But if this is the rainforest land type and we fall + ! within the bounding box of the Amazon rainforest, + ! then increase reactivity as inferred from observations. + IF ( II == 6 .AND. & + State_Grid%XMid(I,J) > -82.0_f8 .AND. & + State_Grid%XMid(I,J) < -33.0_f8 .AND. & + State_Grid%YMid(I,J) > -34.0_f8 .AND. & + State_Grid%YMid(I,J) < 14.0_f8 ) THEN + F0_K = 2.0e-01_f8 ENDIF ENDIF @@ -1588,25 +1644,25 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & #ifdef LUO_WETDEP RIXX = RIX*DIFFG(TEMPK,PRESSU(I,J),XMWH2O)/ & DIFFG(TEMPK,PRESSU(I,J),XMW(K)) & - + 1.e+0_f8/(HSTAR3D(I,J,K)/3000.e+0_f8+100.e+0_f8*F0(K)) + + 1.e+0_f8/(HSTAR3D(I,J,K)/3000.e+0_f8+100.e+0_f8*F0_K) #else RIXX = RIX*DIFFG(TEMPK,PRESSU(I,J),XMWH2O)/ & DIFFG(TEMPK,PRESSU(I,J),XMW(K)) & - + 1.e+0_f8/(HSTAR(K)/3000.e+0_f8+100.e+0_f8*F0(K)) + + 1.e+0_f8/(HSTAR(K)/3000.e+0_f8+100.e+0_f8*F0_K) #endif RLUXX = 1.e+12_f8 IF (RLU(LDT).LT.9999.e+0_f8) & #ifdef LUO_WETDEP - RLUXX = RLU(LDT)/(HSTAR3D(I,J,K)/1.0e+05_f8 + F0(K)) + RLUXX = RLU(LDT)/(HSTAR3D(I,J,K)/1.0e+05_f8 + F0_K) #else - RLUXX = RLU(LDT)/(HSTAR(K)/1.0e+05_f8 + F0(K)) + RLUXX = RLU(LDT)/(HSTAR(K)/1.0e+05_f8 + F0_K) #endif ! If POPs simulation, scale cuticular resistances with octanol- ! air partition coefficient (Koa) instead of HSTAR ! (clf, 1/3/2011) IF (IS_POPS) & - RLUXX = RLU(LDT)/(KOA(K)/1.0e+05_f8 + F0(K)) + RLUXX = RLU(LDT)/(KOA(K)/1.0e+05_f8 + F0_K) !* !* To prevent virtually zero resistance to species with huge @@ -1620,14 +1676,14 @@ SUBROUTINE DEPVEL( Input_Opt, State_Chm, State_Diag, State_Grid, & !* #ifdef LUO_WETDEP RGSX = 1.e+0_f8/(HSTAR3D(I,J,K)/1.0e+05_f8/RGSS(LDT) + & - F0(K)/RGSO(LDT)) + F0_K/RGSO(LDT)) RCLX = 1.e+0_f8/(HSTAR3D(I,J,K)/1.0e+05_f8/RCLS(LDT) + & - F0(K)/RCLO(LDT)) + F0_K/RCLO(LDT)) #else RGSX = 1.e+0_f8/(HSTAR(K)/1.0e+05_f8/RGSS(LDT) + & - F0(K)/RGSO(LDT)) + F0_K/RGSO(LDT)) RCLX = 1.e+0_f8/(HSTAR(K)/1.0e+05_f8/RCLS(LDT) + & - F0(K)/RCLO(LDT)) + F0_K/RCLO(LDT)) #endif !* !** Get the bulk surface resistance of the canopy, RSURFC, from diff --git a/GeosCore/dust_mod.F90 b/GeosCore/dust_mod.F90 index a720c88e3..a82001a11 100644 --- a/GeosCore/dust_mod.F90 +++ b/GeosCore/dust_mod.F90 @@ -1046,6 +1046,23 @@ SUBROUTINE DRY_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero private loop variables + CONST = 0.0_fp + DELZ = 0.0_fp + DELZ1 = 0.0_fp + DEN = 0.0_fp + DP = 0.0_fp + MASS = 0.0_fp + OLD = 0.0_fp + P = 0.0_fp + PDP = 0.0_fp + REFF = 0.0_fp + SLIP = 0.0_fp + TEMP = 0.0_fp + TC0 = 0.0_fp + VISC = 0.0_fp + VTS = 0.0_fp + DO L = 1, State_Grid%NZ DO N = APMIDS%id_DSTBIN1, IDTEMP MASS(L) = MASS(L) + Spc(N)%Conc(I,J,L) diff --git a/GeosCore/mercury_mod.F90 b/GeosCore/mercury_mod.F90 index 95cfc222c..63d97477d 100644 --- a/GeosCore/mercury_mod.F90 +++ b/GeosCore/mercury_mod.F90 @@ -679,8 +679,8 @@ SUBROUTINE ChemMercury( Input_Opt, State_Chm, State_Diag, & INTEGER :: TotSteps, TotFuncs, TotJacob, TotAccep INTEGER :: TotRejec, TotNumLU, HCRC, IERR INTEGER :: Day, S, errorCount - REAL(fp) :: REL_HUM, Start, Finish, rtim - REAL(fp) :: itim, TOUT, T, TIN + REAL(fp) :: REL_HUM, rtim, itim, TOUT + REAL(fp) :: T, TIN ! Strings CHARACTER(LEN=16) :: thisName @@ -961,10 +961,10 @@ SUBROUTINE ChemMercury( Input_Opt, State_Chm, State_Diag, & !$OMP PARALLEL DO & !$OMP DEFAULT( SHARED )& - !$OMP PRIVATE( I, J, L, N )& - !$OMP PRIVATE( IERR, RCNTRL, START, FINISH, ISTATUS )& - !$OMP PRIVATE( RSTATE, SpcID, KppID, F, P )& - !$OMP PRIVATE( Vloc, Aout, NN, C_before_integrate )& + !$OMP PRIVATE( I, J, L, N )& + !$OMP PRIVATE( IERR, RCNTRL, ISTATUS, RSTATE )& + !$OMP PRIVATE( SpcID, KppID, F, P )& + !$OMP PRIVATE( Vloc, Aout, NN, C_before_integrate )& !$OMP COLLAPSE( 3 )& !$OMP SCHEDULE ( DYNAMIC, 24 )& !$OMP REDUCTION( +:errorCount ) diff --git a/GeosCore/ocean_mercury_mod.F90 b/GeosCore/ocean_mercury_mod.F90 index 60e479c08..a1013bb95 100644 --- a/GeosCore/ocean_mercury_mod.F90 +++ b/GeosCore/ocean_mercury_mod.F90 @@ -467,9 +467,10 @@ SUBROUTINE READ_HG2_PARTITIONING( Input_Opt, State_Grid, State_Met, & ! convert REAL*4 to REAL(fpp) DST_CONC(:,:,1:State_Grid%MaxChemLev) = ARRAYdst1(:,:,1:State_Grid%MaxChemLev) - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L )& + !$OMP COLLAPSE( 3 ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -1077,62 +1078,72 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ENDIF ! Loop over surface boxes - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, vi, A_M2, Hg2_RED ) & - !$OMP PRIVATE( J, k_ox, OC_tot, Hg2_CONV ) & - !$OMP PRIVATE( N, TK, CHg0, k_red_bio ) & - !$OMP PRIVATE( C, TC, RADz, Hg0_OX, k_red_rad ) & - !$OMP PRIVATE( D, EC, k_red, OLDMLD, TOTDEPall ) & - !$OMP PRIVATE( Y, Ze, ScCO2, FRAC_O, Frac_Hg2, Hg2aq_tot ) & - !$OMP PRIVATE( H, Kw, MLDCM, TOTDEP, OC_tot_kg ) & - !$OMP PRIVATE( X, SPM, CHg0aq, Hg2_GONE ) & - !$OMP PRIVATE( Sc, Usq, C_tot, JorgC_kg ) & - !$OMP PRIVATE( IS_OCEAN_BOX ) & - !$OMP PRIVATE( FRAC_OPEN_OCEAN, FRAC_OCEAN_OR_ICE, Snow_Hg2aq ) & - !$OMP PRIVATE( FRAC_REDUCIBLE, UVI_RATIO ) & - !$OMP PRIVATE( Kd_part ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, A_M2 )& + !$OMP PRIVATE( C, CHg0, CHg0aq )& + !$OMP PRIVATE( C_tot, D, EC )& + !$OMP PRIVATE( Frac_Hg2, Frac_O, FRAC_OCEAN_OR_ICE )& + !$OMP PRIVATE( FRAC_OPEN_OCEAN, FRAC_REDUCIBLE, H )& + !$OMP PRIVATE( Hg0_OX, Hg2aq_tot, Hg2_Conv )& + !$OMP PRIVATE( Hg2_GONE, Hg2_RED, IS_OCEAN_BOX )& + !$OMP PRIVATE( JorgC_kg, k_ox, k_red )& + !$OMP PRIVATE( k_red_bio, k_red_rad, kd_part )& + !$OMP PRIVATE( Kw, MLDcm, N )& + !$OMP PRIVATE( OLDMLD, OC_tot, OC_tot_kg )& + !$OMP PRIVATE( RADz, Sc, ScCO2 )& + !$OMP PRIVATE( SNOW_Hg2aq, SPM, TC )& + !$OMP PRIVATE( TK, TOTDEP, TOTDEPall )& + !$OMP PRIVATE( Usq, UVI_RATIO, vi )& + !$OMP PRIVATE( Ze )& + !$OMP COLLAPSE( 2 )& + !$OMP SCHEDULE( DYNAMIC, 24 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX - ! Grid box surface area [m2] - A_M2 = State_Grid%Area_M2(I,J) - - ! Grid-box latitude [degrees] - Y = State_Grid%YMid(I,J) - - ! Initialize values - Kw = 0e+0_fpp - Hg2_CONV = 0e+0_fpp - TK = 0e+0_fpp - TC = 0e+0_fpp - JorgC_kg = 0e+0_fpp - EC = 0e+0_fpp - RADz = 0e+0_fpp - k_red = 0e+0_fpp - k_red_rad = 0e+0_fpp - k_red_bio = 0e+0_fpp - SPM = 0e+0_fpp - Frac_Hg2 = 0e+0_fpp - Hg2aq_tot = 0e+0_fpp - Hg2_RED = 0e+0_fpp - C_tot = 0e+0_fpp - Ze = 0e+0_fpp - OC_tot = 0e+0_fpp - OC_tot_kg = 0e+0_fpp - Hg0_OX = 0e+0_fpp - D = 0e+0_fpp - TOTDEP = 0e+0_fpp - TOTDEPall = 0e+0_fpp - k_ox = 0e+0_fpp - SNOW_Hg2aq = 0e+0_fpp - UVI_RATIO = 0e+0_fpp - FRAC_REDUCIBLE = 0e+0_fpp - - OLDMLD = MLDav(I,J) - MLDav(I,J) = MLDav(I,J) + dMLD(I,J) * DTSRCE - MLDcm = MLDav(I,J) + ! Initialize variables (now alphabetized) + A_M2 = State_Grid%Area_M2(I,J) ! Grid box area [m2] + CHg0 = 0.0_fpp + CHg0aq = 0.0_fpp + C_tot = 0.0_fpp + D = 0.0_fpp + EC = 0.0_fpp + Frac_Hg2 = 0.0_fpp + Frac_O = 0.0_fpp + FRAC_REDUCIBLE = 0.0_fpp + H = 0.0_fpp + Hg0_OX = 0.0_fpp + Hg2aq_tot = 0.0_fpp + Hg2_CONV = 0.0_fpp + Hg2_GONE = 0.0_fpp + Hg2_RED = 0.0_fpp + JorgC_kg = 0.0_fpp + k_ox = 0.0_fpp + k_red = 0.0_fpp + k_red_bio = 0.0_fpp + k_red_rad = 0.0_fpp + kd_part = 0.0_fpp + Kw = 0.0_fpp + OLDMLD = MLDav(I,J) + MLDav(I,J) = MLDav(I,J) + dMLD(I,J) * DTSRCE + MLDcm = MLDav(I,J) + OC_tot = 0.0_fpp + OC_tot_kg = 0.0_fpp + RADz = 0.0_fpp + Sc = 0.0_fpp + ScCO2 = 0.0_fpp + SNOW_Hg2aq = 0.0_fpp + SPM = 0.0_fpp + TC = 0.0_fpp + TK = 0.0_fpp + TOTDEP = 0.0_fpp + TOTDEPall = 0.0_fpp + Usq = 0.0_fpp + UVI_RATIO = 0.0_fpp + vi = 0.0_fpp + X = State_Grid%XMid(I,J) ! Longitude [deg] + Y = State_Grid%YMid(I,J) ! Laittude [deg] + Ze = 0.0_fpp ! Add error trap to prevent new MLD from being negative ! (jaf, 7/6/11) @@ -1386,9 +1397,6 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !%%% and also removed NN as the 3rd array index !%%% -- Bob Yantosca (23 Jun 2022) - ! Grid box longitude [degrees] - X = State_Grid%XMid(I,J) - !-------------------------------------------------------- ! Atlantic !-------------------------------------------------------- @@ -2124,9 +2132,11 @@ SUBROUTINE OCEAN_MERCURY_READ( Input_Opt, State_Grid, & ! Overwrite fields with Arctic specific parameters !----------------------------------------------------------------- - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, X, Y ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, X, Y )& + !$OMP COLLAPSE( 2 )& + !$OMP SCHEDULE( DYNAMIC, 24 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX diff --git a/GeosCore/seasalt_mod.F90 b/GeosCore/seasalt_mod.F90 index bf1d9c58f..eead8fa74 100644 --- a/GeosCore/seasalt_mod.F90 +++ b/GeosCore/seasalt_mod.F90 @@ -1804,7 +1804,23 @@ SUBROUTINE WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX - MASS = 0.d0 + ! Zero private loop variables + CONST = 0.0_fp + DELZ = 0.0_fp + DELZ1 = 0.0_fp + DEN = 0.0_fp + DP = 0.0_fp + MASS = 0.0_fp + OLD = 0.0_fp + P = 0.0_fp + PDP = 0.0_fp + REFF = 0.0_fp + SLIP = 0.0_fp + TEMP = 0.0_fp + TC0 = 0.0_fp + VISC = 0.0_fp + VTS = 0.0_fp + DO L = 1, State_Grid%NZ DO N = IDTEMP1, IDTEMP2 MASS(L) = MASS(L) + Spc(N)%Conc(I,J,L) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index 7d0729369..ecb76ceb8 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -365,10 +365,9 @@ SUBROUTINE CHEMSULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, & #ifdef APM IDNH3 = HCO_GetHcoID( 'NH3', HcoState ) - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( L, J, I, A_M2 ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( L, J, I, A_M2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -385,10 +384,9 @@ SUBROUTINE CHEMSULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP END PARALLEL DO IDSO2 = HCO_GetHcoID( 'SO2', HcoState ) - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( L, J, I, A_M2 ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( L, J, I, A_M2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -2654,16 +2652,18 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PRIVATE( HCHO0, HMSc, HMS0, OH0, KaqHCHO )& !$OMP PRIVATE( KaqHMS, KaqHMS2, L7, L7S, L7_b )& !$OMP PRIVATE( L7S_b, L8, L8S, LSTOT_HMS )& - !$OMP SCHEDULE( DYNAMIC, 1 ) + !$OMP COLLAPSE( 3 )& + !$OMP SCHEDULE( DYNAMIC, 24 ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX ! Initialize for safety's sake - Ld = 0.0_fp - LSTOT0 = 0.0_fp - LSTOT = 0.0_fp - LSTOT_HMS = 0.0_fp + Ld = 0.0_fp + LSTOT0 = 0.0_fp + LSTOT = 0.0_fp + LSTOT_HMS = 0.0_fp + one_m_KRATE = 0.0_fp ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE @@ -7520,8 +7520,8 @@ END SUBROUTINE HET_DROP_CHEM !\\ ! !INTERFACE: ! - SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & - State_Met, RC ) + SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, & + State_Grid, State_Met, RC ) ! ! !USES: ! @@ -7586,13 +7586,9 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Pointers TYPE(SpcConc), POINTER :: Spc(:) -#ifdef APM - REAL(fp), POINTER :: PSO4_SO2APM2(:,:,:) -#endif #ifdef APM - REAL*8 :: MASS0, MASS, PMASS - REAL*8 :: RKTs, E_RKTs, DTCHEM + REAL(fp), PARAMETER :: AIRMW_96 = AIRMW / 96.0d0 #endif !================================================================= @@ -7612,38 +7608,32 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & !------------------------------------------ ! Call APM size-resolved drydep algorithm !------------------------------------------ - CALL WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, State_Grid, & - State_Met, RC ) - - ! Point to PSO4_SO2APM2 now moved to State_Met - PSO4_SO2APM2 => State_Met%PSO4_SO2APM2 + CALL WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & + State_Grid, State_Met, RC ) #endif ! Point to chemical species array [kg] Spc => State_Chm%Species ! Loop over chemistry grid boxes - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L, SO4, SO4s, SO40, SO40s ) & - !$OMP PRIVATE( SO4d, SO40d, SO40_dust ) & - !$OMP PRIVATE( IBIN, PSO4d, IDTRC ) & -#ifdef APM - !$OMP PRIVATE( N, MASS0, MASS, PMASS, RKTs, E_RKTs ) & -#endif - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L, N, SO4, SO4s, SO40 )& + !$OMP PRIVATE( SO40s, SO4d, SO40d, SO40_dust, IBIN, PSO4d, IDTRC )& + !$OMP SCHEDULE( DYNAMIC ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero loop variables + SO4 = 0.0_fp + SO4s = 0.0_fp + SO40 = 0.0_fp + SO40s = 0.0_fp + ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE - ! Initialize for safety's sake - SO4 = 0e+0_fp - SO4s = 0e+0_fp - SO4d = 0e+0_fp ! tdf 04/07/08 - !============================================================== ! Initial concentrations before chemistry !============================================================== @@ -7676,13 +7666,10 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & ! SO4 production from SO2 [v/v/timestep] SO4s = SO40s + PSO4_ss(I,J,L) - !tdf + !============================================================== + ! SO4d (SO4 w/in dust aerosol) chemistry: tdf 04/07/08 + !============================================================== IF ( LDSTUP ) THEN - - !============================================================== - ! SO4d (SO4 w/in dust aerosol) chemistry: tdf 04/07/08 - !============================================================== - IDTRC(1) = id_SO4d1 IDTRC(2) = id_SO4d2 IDTRC(3) = id_SO4d3 @@ -7717,55 +7704,57 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & !APM_GanLuo+ #ifdef APM - IF(NSO4>=1)THEN - DO N=1,NSO4 + IF ( NSO4 >= 1 ) THEN + DO N = 1, NSO4 ! Updated SO4 (gas phase) [v/v] - Spc(APMIDS%id_SO4BIN1+N-1)%Conc(I,J,L) = & - Spc(APMIDS%id_SO4BIN1+N-1)%Conc(I,J,L) + & - (PSO4_SO2APM(I,J,L)+PSO4_SO2APM2(I,J,L)*(AIRMW/96.D0)/ & - (g0_100*State_Met%DELP_DRY(I,J,L)))* & + Spc(APMIDS%id_SO4BIN1+N-1)%Conc(I,J,L) = & + Spc(APMIDS%id_SO4BIN1+N-1)%Conc(I,J,L) + & + ( PSO4_SO2APM(I,J,L) + & + State_Chm%PSO4_SO2APM2(I,J,L) * AIRMW_96 / & + ( g0_100 * State_Met%DELP_DRY(I,J,L) ) ) * & FCLOUD(I,J,L,N) ENDDO ENDIF - IF(NCTBC>=1)THEN - DO N=1,1 - Spc(APMIDS%id_CTBC+N-1)%Conc(I,J,L) = & - Spc(APMIDS%id_CTBC+N-1)%Conc(I,J,L) + & - (PSO4_SO2APM(I,J,L)+PSO4_SO2APM2(I,J,L)*(AIRMW/96.D0)/ & - (g0_100*State_Met%DELP_DRY(I,J,L)))* & - FCLOUD(I,J,L,(NSO4+N)) - ENDDO + IF ( NCTBC >=1 ) THEN + N = 1 + Spc(APMIDS%id_CTBC)%Conc(I,J,L) = & + Spc(APMIDS%id_CTBC)%Conc(I,J,L) + & + ( PSO4_SO2APM(I,J,L) + & + State_Chm%PSO4_SO2APM2(I,J,L) * AIRMW_96 / & + ( g0_100 * State_Met%DELP_DRY(I,J,L) ) ) * & + FCLOUD(I,J,L,(NSO4+N)) ENDIF - IF(NCTOC>=1)THEN - DO N=1,1 - Spc(APMIDS%id_CTOC+N-1)%Conc(I,J,L) = & - Spc(APMIDS%id_CTOC+N-1)%Conc(I,J,L) + & - (PSO4_SO2APM(I,J,L)+PSO4_SO2APM2(I,J,L)*(AIRMW/96.D0)/ & - (g0_100*State_Met%DELP_DRY(I,J,L)))* & - FCLOUD(I,J,L,(NSO4+N)) - ENDDO + IF ( NCTOC >= 1 ) THEN + N = 1 + Spc(APMIDS%id_CTOC)%Conc(I,J,L) = & + Spc(APMIDS%id_CTOC)%Conc(I,J,L) + & + ( PSO4_SO2APM(I,J,L) + & + State_Chm%PSO4_SO2APM2(I,J,L) * AIRMW_96 / & + ( g0_100 *State_Met%DELP_DRY(I,J,L) ) ) * & + FCLOUD(I,J,L,(NSO4+N)) ENDIF - IF(NCTDST>=1)THEN - DO N=1,1 - Spc(APMIDS%id_CTDST+N-1)%Conc(I,J,L) = & - Spc(APMIDS%id_CTDST+N-1)%Conc(I,J,L) + & - (PSO4_SO2APM(I,J,L)+PSO4_SO2APM2(I,J,L)*(AIRMW/96.D0)/ & - (g0_100*State_Met%DELP_DRY(I,J,L)))* & - FCLOUD(I,J,L,(NSO4+3)) - ENDDO + IF ( NCTDST >= 1 ) THEN + N = 1 + Spc(APMIDS%id_CTDST)%Conc(I,J,L) = & + Spc(APMIDS%id_CTDST)%Conc(I,J,L) + & + ( PSO4_SO2APM(I,J,L) + & + State_Chm%PSO4_SO2APM2(I,J,L) * AIRMW_96 / & + ( g0_100 * State_Met%DELP_DRY(I,J,L) ) ) * & + FCLOUD(I,J,L,(NSO4+3)) ENDIF - IF(NCTSEA>=1)THEN - DO N=1,1 - Spc(APMIDS%id_CTSEA+N-1)%Conc(I,J,L) = & - Spc(APMIDS%id_CTSEA+N-1)%Conc(I,J,L) + & - (PSO4_SO2APM(I,J,L)+PSO4_SO2APM2(I,J,L)*(AIRMW/96.D0)/ & - (g0_100*State_Met%DELP_DRY(I,J,L)))* & - FCLOUD(I,J,L,(NSO4+4)) + PSO4_SO2SEA(I,J,L) - ENDDO + IF ( NCTSEA >= 1 ) THEN + N = 1 + Spc(APMIDS%id_CTSEA)%Conc(I,J,L) = & + Spc(APMIDS%id_CTSEA)%Conc(I,J,L) + & + ( PSO4_SO2APM(I,J,L) + & + State_Chm%PSO4_SO2APM2(I,J,L) * AIRMW_96 / & + ( g0_100 * State_Met%DELP_DRY(I,J,L) ) ) * & + FCLOUD(I,J,L,(NSO4+4)) + & + PSO4_SO2SEA(I,J,L) ENDIF #endif @@ -7784,7 +7773,8 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP END PARALLEL DO #ifdef APM - PSO4_SO2APM2 = 0.D0 + ! Reset the PSO4_SO2APM2 which tracks SO4 from wetdep + State_Chm%PSO4_SO2APM2 = 0.0d0 #endif ! Free pointers @@ -7868,15 +7858,19 @@ SUBROUTINE CHEM_SO4_AQ( Input_Opt, State_Chm, State_Grid, State_Met, RC ) RETURN ENDIF - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L ) & - !$OMP PRIVATE( KMIN, SO4OXID, BINACT1, BINACT2 ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L, KMIN, SO4OXID, BINACT1, BINACT2 ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero private loop variables + BINACT1 = 0.0_fp + BINACT2 = 0.0_fp + KMIN = 0.0_fp + SO4OXID = 0.0_fp + ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE @@ -7887,13 +7881,13 @@ SUBROUTINE CHEM_SO4_AQ( Input_Opt, State_Chm, State_Grid, State_Met, RC ) ! JKodros (6/2/15 - Set activating bin based on which TOMAS bin !length being used) #if defined( TOMAS12 ) - CALL GETACTBIN( I, J, L, id_NK5, .TRUE. , BINACT1, State_Chm, RC ) + CALL GETACTBIN( I, J, L, id_NK5, .TRUE. , BINACT1, State_Chm, RC ) - CALL GETACTBIN( I, J, L, id_NK5, .FALSE., BINACT2, State_Chm, RC ) + CALL GETACTBIN( I, J, L, id_NK5, .FALSE., BINACT2, State_Chm, RC ) #elif defined( TOMAS15 ) - CALL GETACTBIN( I, J, L, id_NK8, .TRUE. , BINACT1, State_Chm, RC ) + CALL GETACTBIN( I, J, L, id_NK8, .TRUE. , BINACT1, State_Chm, RC ) - CALL GETACTBIN( I, J, L, id_NK8, .FALSE., BINACT2, State_Chm, RC ) + CALL GETACTBIN( I, J, L, id_NK8, .FALSE., BINACT2, State_Chm, RC ) #elif defined( TOMAS30 ) CALL GETACTBIN( I, J, L, id_NK10, .TRUE. , BINACT1, State_Chm, RC ) @@ -9242,7 +9236,7 @@ END SUBROUTINE CLEANUP_SULFATE !------------------------------------------------------------------------------ !BOP ! -! !IROUTINE: dry_settlingbin +! !IROUTINE: wet_settlingbin ! ! !DESCRIPTION: Subroutine WET\_SETTLINGBIN computes the dry settling of ! aerosol tracers. Modified for APM simulation. (G. Luo) @@ -9250,8 +9244,8 @@ END SUBROUTINE CLEANUP_SULFATE !\\ ! !INTERFACE: ! - SUBROUTINE WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, State_Grid, & - State_Met, RC ) + SUBROUTINE WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, & + State_Grid, State_Met, RC ) ! ! !USES: ! @@ -9344,17 +9338,34 @@ SUBROUTINE WET_SETTLINGBIN( Input_Opt, State_Chm, State_Diag, State_Grid, & IDTEMP1 = APMIDS%id_SO4BIN1 IDTEMP2 = APMIDS%id_SO4BIN1+NSO4-1 - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L, N, K, DEN, REFF, DP ) & - !$OMP PRIVATE( CONST, VTS, TEMP, P, PDP, SLIP ) & - !$OMP PRIVATE( MASS, OLD, VISC, TC0, DELZ, DELZ1 ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L, N, K, DEN, REFF )& + !$OMP PRIVATE( DP, CONST, VTS, TEMP, P, PDP, SLIP )& + !$OMP PRIVATE( MASS, OLD, VISC, TC0, DELZ, DELZ1 )& + !$OMP SCHEDULE( DYNAMIC, 1 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX + ! Zero private loop variables + CONST = 0.0_fp + DELZ = 0.0_fp + DELZ1 = 0.0_fp + DEN = 0.0_fp + DP = 0.0_fp + MASS = 0.0_fp + OLD = 0.0_fp + P = 0.0_fp + PDP = 0.0_fp + REFF = 0.0_fp + SLIP = 0.0_fp + TEMP = 0.0_fp + TC0 = 0.0_fp + VISC = 0.0_fp + VTS = 0.0_fp + DO L = 1, State_Grid%NZ - MASS(L) = 0.d8 + !MASS(L) = 0.0d0 DO N = IDTEMP1, IDTEMP2 MASS(L) = MASS(L) + Spc(N)%Conc(I,J,L) ENDDO diff --git a/GeosCore/wetscav_mod.F90 b/GeosCore/wetscav_mod.F90 index dee2b03b4..caf725006 100644 --- a/GeosCore/wetscav_mod.F90 +++ b/GeosCore/wetscav_mod.F90 @@ -417,9 +417,10 @@ SUBROUTINE MAKE_QQ( State_Chm, State_Grid, State_Met, LS, RC ) ! Initialize RC = GC_SUCCESS - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED ) & + !$OMP PRIVATE( I, J, L ) & + !$OMP COLLAPSE( 3 ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -1353,6 +1354,9 @@ SUBROUTINE RAINOUT( I, J, L, N, & ENDIF ENDIF + + ! Free pointer + p_pHCloud => NULL() #else !================================================================= ! %%% SPECIAL CASE %%% @@ -1738,6 +1742,15 @@ SUBROUTINE WASHOUT( I, J, L, & ! to be a gas-phase species elsewhere (e.g. dry deposition) !================================================================= #ifdef LUO_WETDEP + + ! Initialize + Hplus = 0.0_f8 + HCSO2 = 0.0_f8 + HCNH3 = 0.0_f8 + Ks1 = 0.0_f8 + Ks2 = 0.0_f8 + T_Term = 0.0_f8 + IF ( SpcInfo%WD_Is_HNO3 ) THEN ! Washout is a kinetic process @@ -1869,7 +1882,8 @@ SUBROUTINE WASHOUT( I, J, L, & ELSE IF ( SpcInfo%MP_SizeResAer .or. SpcInfo%MP_SizeResNum ) THEN #ifdef APM ! Washout is a kinetic process - KIN = .TRUE. + KIN = .TRUE. + RIN = 0.0_fp IF(SpcInfo%Name(1:8)=='APMSPBIN')THEN RIN = RDRY(N-APMIDS%id_SO4BIN1+1) * GFTOT3D(I,J,L,1) @@ -3268,7 +3282,8 @@ SUBROUTINE WETDEP( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PRIVATE( QDOWN, IS_RAINOUT, IS_WASHOUT, N ) & !$OMP PRIVATE( DEP_HG, SpcInfo, Hg_Cat, EC ) & !$OMP PRIVATE( COND_WATER_CONTENT ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP COLLAPSE( 2 ) & + !$OMP SCHEDULE( DYNAMIC, 24 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -4552,8 +4567,8 @@ SUBROUTINE DO_WASHOUT_ONLY( LS, I, J, L, & + GAINED * 96e+0_fp / 64e+0_fp #ifdef APM - State_Met%PSO4_SO2APM2(I,J,L) = State_Met%PSO4_SO2APM2(I,J,L) & - + GAINED * 96e+0_fp / 64e+0_fp + State_Chm%PSO4_SO2APM2(I,J,L) = & + State_Chm%PSO4_SO2APM2(I,J,L) + GAINED * 96e+0_fp / 64e+0_fp #endif Spc(N)%Conc(I,J,L) = Spc(N)%Conc(I,J,L) * ( 1e+0_fp - WASHFRAC ) @@ -4891,8 +4906,8 @@ SUBROUTINE DO_COMPLETE_REEVAP( LS, I, J, & ( WETLOSS * 96e+0_fp / 64e+0_fp ) #ifdef APM - State_Met%PSO4_SO2APM2(I,J,L) = State_Met%PSO4_SO2APM2(I,J,L) - & - ( WETLOSS * 96e+0_fp / 64e+0_fp ) + State_Chm%PSO4_SO2APM2(I,J,L) = & + State_Chm%PSO4_SO2APM2(I,J,L) - ( WETLOSS * 96e+0_fp / 64e+0_fp ) #endif #ifdef TOMAS @@ -5826,10 +5841,11 @@ SUBROUTINE SETUP_WETSCAV( Input_Opt, State_Chm, State_Grid, State_Met, RC ) ! Only do computation if wetdep or convection is turned on IF ( Input_Opt%LWETD .or. Input_Opt%LCONV ) THEN - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L, TK, PL ) & - !$OMP SCHEDULE( DYNAMIC ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L, TK, PL )& + !$OMP COLLAPSE( 3 )& + !$OMP SCHEDULE( DYNAMIC, 24 ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX diff --git a/Headers/state_chm_mod.F90 b/Headers/state_chm_mod.F90 index 7bae81c31..4cb7c8280 100644 --- a/Headers/state_chm_mod.F90 +++ b/Headers/state_chm_mod.F90 @@ -330,6 +330,13 @@ MODULE State_Chm_Mod REAL(fp), POINTER :: CH4_EMIS (:,:,: ) ! CH4 emissions [kg/m2/s]. ! third dim is cat, total 15 +#ifdef APM + !----------------------------------------------------------------------- + ! Fields for APM aerosol microphysics + !----------------------------------------------------------------------- + REAL(fp), POINTER :: PSO4_SO2APM2(:,:,: ) +#endif + !----------------------------------------------------------------------- ! Registry of variables contained within State_Chm !----------------------------------------------------------------------- @@ -612,6 +619,10 @@ SUBROUTINE Zero_State_Chm( State_Chm, RC ) ! Add quantities for coupling to CESM State_Chm%H2SO4_PRDR => NULL() #endif +#ifdef APM + ! Add fields for APM microphysics + State_Chm%PSO4_SO2APM2 => NULL() +#endif END SUBROUTINE Zero_State_Chm !EOC @@ -2268,6 +2279,27 @@ SUBROUTINE Init_State_Chm( Input_Opt, State_Chm, State_Grid, RC ) ENDIF ENDIF +#ifdef APM + !======================================================================= + ! Initialize State_Chm quantities for APM microphysics simulations + !======================================================================= + chmId = 'PSO4SO2APM2' + CALL Init_and_Register( & + Input_Opt = Input_Opt, & + State_Chm = State_Chm, & + State_Grid = State_Grid, & + chmId = chmId, & + Ptr2Data = State_Chm%PSO4_SO2APM2, & + noRegister = .TRUE., & + RC = RC ) + + IF ( RC /= GC_SUCCESS ) THEN + errMsg = TRIM( errMsg_ir ) // TRIM( chmId ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF +#endif + !======================================================================== ! Once we are done registering all fields, we need to define the ! registry lookup table. This algorithm will avoid hash collisions. @@ -3717,6 +3749,15 @@ SUBROUTINE Cleanup_State_Chm( State_Chm, RC ) State_Chm%NOXLAT => NULL() ENDIF +#ifdef APM + IF ( ASSOCIATED( State_Chm%PSO4_SO2APM2 ) ) THEN + DEALLOCATE( State_Chm%PSO4_SO2APM2, STAT=RC ) + CALL GC_CheckVar( 'State_Chm%PSO4_SO2APM2', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Chm%PSO4_SO2APM2 => NULL() + ENDIF +#endif + !----------------------------------------------------------------------- ! Template for deallocating more arrays, replace xxx with field name !----------------------------------------------------------------------- diff --git a/Headers/state_met_mod.F90 b/Headers/state_met_mod.F90 index 8cba2ee70..6da6d4427 100644 --- a/Headers/state_met_mod.F90 +++ b/Headers/state_met_mod.F90 @@ -268,7 +268,6 @@ MODULE State_Met_Mod ! [cm3 H2O/cm2 area/s] REAL(fp), POINTER :: QQ (:,:,:) ! Rate of new precip formation [cm3 H2O/cm3 air/s] REAL(fp), POINTER :: REEVAP (:,:,:) ! Rate of precip reevaporation - REAL(fp), POINTER :: PSO4_SO2APM2 (:,:,:) !---------------------------------------------------------------------- ! Fields for boundary layer mixing @@ -498,7 +497,6 @@ SUBROUTINE Zero_State_Met( State_Met, RC ) State_Met%IMIX => NULL() State_Met%FPBL => NULL() State_Met%REEVAP => NULL() - State_Met%PSO4_SO2APM2 => NULL() State_Met%PBL_MAX_L = 0 END SUBROUTINE Zero_State_Met @@ -3035,26 +3033,7 @@ SUBROUTINE Init_State_Met( Input_Opt, State_Grid, State_Met, RC ) RETURN ENDIF -#ifdef APM - !--------------------------------------------------------------------- - ! PSO4_SO2APM2 - !--------------------------------------------------------------------- - metId = 'PSO4SO2APM2' - CALL Init_and_Register( & - Input_Opt = Input_Opt, & - State_Met = State_Met, & - State_Grid = State_Grid, & - metId = metId, & - Ptr2Data = State_Met%PSO4_SO2APM2, & - noRegister = .TRUE., & - RC = RC ) - IF ( RC /= GC_SUCCESS ) THEN - errMsg = TRIM( errMsg_ir ) // TRIM( metId ) - CALL GC_Error( errMsg, RC, thisLoc ) - RETURN - ENDIF -#endif ENDIF @@ -4448,18 +4427,9 @@ SUBROUTINE Cleanup_State_Met( State_Met, RC ) State_Met%REEVAP => NULL() ENDIF -#ifdef APM - IF ( ASSOCIATED( State_Met%PSO4_SO2APM2 ) ) THEN - DEALLOCATE( State_Met%PSO4_SO2APM2, STAT=RC ) - CALL GC_CheckVar( 'State_Met%PSO4_SO2APM2', 2, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Met%PSO4_SO2APM2 => NULL() - ENDIF -#endif - - !------------------------------------------------------------------------- + !------------------------------------------------------------------------- ! Template for deallocating more arrays, replace xxx with field name - !------------------------------------------------------------------------- + !------------------------------------------------------------------------- !IF ( ASSOCIATED( State_Met%xxx ) ) THEN ! DEALLOCATE( State_Met%xxx, STAT=RC ) ! CALL GC_CheckVar( 'State_Met%xxx', 2, RC ) diff --git a/KPP/fullchem/fullchem_SulfurChemFuncs.F90 b/KPP/fullchem/fullchem_SulfurChemFuncs.F90 index 51d6e73e8..a3cfadd17 100644 --- a/KPP/fullchem/fullchem_SulfurChemFuncs.F90 +++ b/KPP/fullchem/fullchem_SulfurChemFuncs.F90 @@ -1027,7 +1027,11 @@ SUBROUTINE SET_SO2( I, J, L, Input_Opt, & ! Get liquid water content [m3 H2O/m3 air] within cloud from met flds ! Units: [kg H2O/kg air] * [kg air/m3 air] * [m3 H2O/1e3 kg H2O] #ifdef LUO_WETDEP - ! Luo et al wetdep scheme + ! QQ3D and similar arrays are only allocated for wetdep or convection + Is_QQ3D = ( Input_Opt%LWETD .or. Input_Opt%LCONV ) + + ! If QQ3D is allocated, compute LWC according to Luo et al wetdep scheme. + ! Otherwise, compute LWC according to the default wetdep scheme. IF ( Is_QQ3D ) THEN LWC = State_Met%QL(I,J,L) * State_Met%AIRDEN(I,J,L) * 1e-3_fp + & MAX( 0.0_fp, State_Chm%QQ3D(I,J,L) * DTCHEM ) @@ -1035,7 +1039,7 @@ SUBROUTINE SET_SO2( I, J, L, Input_Opt, & LWC = State_Met%QL(I,J,L) * State_Met%AIRDEN(I,J,L) * 1e-3_fp ENDIF #else - ! Default scheme + ! Compute LWC according to the default wetdep scheme LWC = State_Met%QL(I,J,L) * State_Met%AIRDEN(I,J,L) * 1e-3_fp #endif