GEOS-ESM · mmanyin · Jan 21, 2025 · Jun 12, 2023 · Jun 12, 2023 · Jun 12, 2023
@@ -16,7 +16,7 @@ setupcoag.o smallconc.o step.o supersat.o vaporp.o vertadv.o \
 vertdif.o vertical.o versol.o rhopart.o psolve.o zeromicro.o \
 nsubsteps.o setupgrow.o setupgkern.o setupnuc.o growevapl.o microfast.o \
 gsolve.o actdropl.o freezglaerl_murray2010.o growp.o downgxfer.o \
-gasexchange.o melticel.o upgxfer.o freezdropl.o \
+gasexchange.o melticel.o upgxfer.o freezdropl.o pfastdmdt.o \
 downgevapply.o evapp.o evap_poly.o evap_mono.o \
 evap_ingrp.o tsolve.o miess.o vaporp_h2o_buck1981.o \
 vaporp_h2o_murphy2005.o maxconc.o setupbdif.o setupvf_std.o \
@@ -159,6 +159,9 @@ growp.o : growp.F90 carma_globaer.h carma_mod.mod carmastate_mod.mod carma_types
 gsolve.o : gsolve.F90 carma_globaer.h carma_mod.mod carmastate_mod.mod carma_types_mod.mod carma_constants_mod.mod carma_enums_mod.mod carma_precision_mod.mod 
 	$(FORTRAN) $(FFLAGS) -c $<
 
+pfastdmdt.o : pfastdmdt.F90 carma_globaer.h carma_mod.mod carmastate_mod.mod carma_types_mod.mod carma_constants_mod.mod carma_enums_mod.mod carma_precision_mod.mod 
+	$(FORTRAN) $(FFLAGS) -c $<
+
 fractal_meanfield_mod.o : fractal_meanfield_mod.F90 carma_globaer.h carma_mod.mod carmastate_mod.mod carma_types_mod.mod carma_constants_mod.mod carma_enums_mod.mod carma_precision_mod.mod adgaquad_mod.mod lusolvec_mod.mod 
 	$(FORTRAN) $(FFLAGS) -c $<
 

@@ -85,11 +85,11 @@ module carma_enums_mod
   integer, public, parameter :: I_SWG_RURAL     = 14  !! Rural
 
   ! Routines to calculate gas vapor pressures
+  integer, public, parameter :: I_VAPRTN_NULL              = -1  !! For non-condensing gases
   integer, public, parameter :: I_VAPRTN_H2O_BUCK1981      = 1   !! H2O, Buck[1981]
   integer, public, parameter :: I_VAPRTN_H2O_MURPHY2005    = 2   !! H2O, Murphy & Koop [2005]
   integer, public, parameter :: I_VAPRTN_H2O_GOFF1946      = 3   !! H2O, Goff & Gratch [1946], used in CAM
   integer, public, parameter :: I_VAPRTN_H2SO4_AYERS1980   = 4   !! H2SO4, Ayers [1980] & Kumala [1990]
-  integer, public, parameter :: I_VAPRTN_NULL              = 5   !! For non-condensing gases
 
   ! Routines to calculate fall velocities
   integer, public, parameter :: I_FALLRTN_STD              = 1   !! Standard CARMA 2.3 routine (spherical only)
@@ -118,6 +118,7 @@ module carma_enums_mod
   integer, public, parameter :: RC_ERROR          = -1  !! Failure
   integer, public, parameter :: RC_WARNING        = 1   !! Warning
   integer, public, parameter :: RC_WARNING_RETRY  = 2   !! Warning, Retry Suggested
+  integer, public, parameter :: RC_WARNING_PFAST  = 3   !! Warning, Parameterization Suggested
 
 
   !  Define values of symbols used to specify horizontal & vertical grid type.

@@ -145,6 +145,7 @@
 #define do_print_init carma%f_do_print_init
 #define do_step       carma%f_do_step
 #define do_substep    carma%f_do_substep
+#define do_pfast      carma%f_do_pfast
 #define do_thermo     carma%f_do_thermo
 #define do_vdiff      carma%f_do_vdiff
 #define do_vtran      carma%f_do_vtran

@@ -333,7 +333,7 @@ end subroutine CARMA_Create
   !!  @version Feb-2009 
   !!  @author  Chuck Bardeen 
   subroutine CARMA_Initialize(carma, rc, do_cnst_rlh, do_coag, do_detrain, do_fixedinit, &
-      do_grow, do_incloud, do_explised, do_print_init, do_substep, do_thermo, do_vdiff, &
+          do_grow, do_incloud, do_explised, do_print_init, do_substep, do_pfast, do_thermo, do_vdiff, &
       do_vtran, do_drydep, vf_const, minsubsteps, maxsubsteps, maxretries, conmax, &
       do_pheat, do_pheatatm, dt_threshold, cstick, gsticki, gstickl, tstick, do_clearsky, &
       do_partialinit)
@@ -348,6 +348,7 @@ subroutine CARMA_Initialize(carma, rc, do_cnst_rlh, do_coag, do_detrain, do_fixe
     logical, intent(in), optional       :: do_incloud    !! do incloud growth and coagulation?
     logical, intent(in), optional       :: do_explised   !! do sedimentation with substepping
     logical, intent(in), optional       :: do_substep    !! do substepping
+    logical, intent(in), optional       :: do_pfast      !! do parameterized microfast
     logical, intent(in), optional       :: do_print_init !! do prinit initializtion information
     logical, intent(in), optional       :: do_thermo     !! do thermodynamics
     logical, intent(in), optional       :: do_vdiff      !! do Brownian diffusion
@@ -384,6 +385,7 @@ subroutine CARMA_Initialize(carma, rc, do_cnst_rlh, do_coag, do_detrain, do_fixe
     carma%f_do_pheatatm   = .FALSE.
     carma%f_do_print_init = .FALSE.
     carma%f_do_substep    = .FALSE.
+    carma%f_do_pfast      = .FALSE.
     carma%f_do_thermo     = .FALSE.
     carma%f_do_vdiff      = .FALSE.
     carma%f_do_vtran      = .FALSE.
@@ -408,6 +410,7 @@ subroutine CARMA_Initialize(carma, rc, do_cnst_rlh, do_coag, do_detrain, do_fixe
     if (present(do_pheatatm))   carma%f_do_pheatatm   = do_pheatatm
     if (present(do_print_init)) carma%f_do_print_init = (do_print_init .and. carma%f_do_print)
     if (present(do_substep))    carma%f_do_substep    = do_substep
+    if (present(do_pfast))      carma%f_do_pfast      = do_pfast
     if (present(do_thermo))     carma%f_do_thermo     = do_thermo
     if (present(do_vdiff))      carma%f_do_vdiff      = do_vdiff
     if (present(do_vtran))      carma%f_do_vtran      = do_vtran 

@@ -353,6 +353,7 @@ module carma_types_mod
     !   do_print_init If .true. then do print initializtion info          {init}
     !   do_step     if .true. then varstepping succeeded                  {init}
     !   do_substep  if .true. then use substepping                        {init}
+    !   do_pfast    if .true. then use microfast parameterization         {init}
     !   do_thermo   if .true. then do solve thermodynamic equation        {init}
     !   do_vdiff    If .true. then do Brownian diffusion                  {init}
     !   do_vtran    If .true. then do vertical transport                  {init}
@@ -394,6 +395,7 @@ module carma_types_mod
     logical                                       :: f_do_print_init
     logical                                       :: f_do_step
     logical                                       :: f_do_substep
+    logical                                       :: f_do_pfast
     logical                                       :: f_do_thermo
     logical                                       :: f_do_cnst_rlh
     logical, allocatable, dimension(:,:)          :: f_if_nuc       !(NELEM,NELEM)

@@ -75,8 +75,9 @@ subroutine freezaerl_koop2000(carma, cstate, iz, rc)
   real(kind=f)                         :: fkelv
   real(kind=f)                         :: fkelvi
 
-
-  rc = RC_OK
+  ! PAC: Assuming success inside of a subroutine is a problem, any errors
+  !       are erased by this.
+  !rc = RC_OK
 
   !  Aerosol freezing limited to T < 240K
   if (t(iz) <= 240._f) then

@@ -76,7 +76,9 @@ subroutine freezaerl_mohler2010(carma, cstate, iz, rc)
   real(kind=f)                         :: fkelvi
 
 
-  rc = RC_OK
+  ! PAC: Assuming success inside of a subroutine is a problem, any errors
+  !       are erased by this.
+  !rc = RC_OK
 
   !  Aerosol freezing limited to T < 240K
   if (t(iz) <= 240._f) then

@@ -60,7 +60,9 @@ subroutine freezglaerl_murray2010(carma, cstate, iz, rc)
   real(kind=f)                         :: ssi, ssiold
 
   ! Assume success.
-  rc = RC_OK
+  ! PAC: Assuming success inside of a subroutine is a problem, any errors
+  !       are erased by this.
+  !rc = RC_OK
 
   ! Loop over particle groups.
   do igroup = 1,NGROUP

@@ -64,7 +64,7 @@ subroutine gsolve(carma, cstate, iz, previous_ice, previous_liquid, scale_thresh
 
     if (gc(iz,igas) < 0.0_f) then
       if (do_substep) then
-        if (nretries == maxretries) then 
+        if (nretries == maxretries .and. .not. do_pfast) then 
           if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, lat, lon, gc(iz,igas), gasprod(igas), &
             supsati(iz,igas), supsatl(iz,igas), t(iz)
           if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz,igas), supsatlold(iz,igas), told(iz), d_gc(iz,igas), d_t(iz)

@@ -151,14 +151,17 @@ subroutine microfast(carma, cstate, iz, scale_threshold, rc)
     call downgevapply(carma, cstate, iz, rc)
     if (rc < RC_OK) return
 
+    ! If this is a parameterized timestep, scale bin values to prevent gas overshoot
+    if (rc == RC_WARNING_PFAST) call pfastdmdt(carma, cstate, iz, rc)
+
     call gsolve(carma, cstate, iz, previous_ice, previous_liquid, scale_threshold, rc)
-    if (rc /=RC_OK) return
+    if (rc /= RC_OK .and. rc /= RC_WARNING_PFAST) return
   endif
 
   ! Update temperature if thermal processes requested
   if (do_thermo) then
     call tsolve(carma, cstate, iz, scale_threshold, rc)
-    if (rc /= RC_OK) return
+    if (rc /= RC_OK .and. rc /= RC_WARNING_PFAST) return
   endif
 
   !  Update saturation ratios

@@ -188,16 +188,29 @@ subroutine newstate_calc(carma, cstate, scale_threshold, rc)
               nretries = nretries + 1
 
               if (nretries > maxretries) then
-                if (do_print) write(LUNOPRT,1) iz, isubstep, ntsubsteps, nretries - 1._f
-                rc = RC_ERROR
-                exit
+                if (do_pfast) then
+                    if (nretries > maxretries+1) then
+                        if (do_print) write(LUNOPRT,1) iz, isubstep, ntsubsteps, -1._f*nretries
+                        rc = RC_ERROR
+                        exit
+                    endif
+                    rc = RC_WARNING_PFAST
+                else
+                    if (do_print) write(LUNOPRT,1) iz, isubstep, ntsubsteps, nretries - 1._f
+                    rc = RC_ERROR
+                    exit
+                end if
               end if
 
               ! Try twice the substeps
               !
               ! NOTE: We are going to rely upon retries, so don't clutter the log
               ! with retry print statements. They slow down the run.
-              ntsubsteps = ntsubsteps * 2
+              if (rc == RC_WARNING_PFAST) then
+                  ntsubsteps = minsubsteps
+              else
+                  ntsubsteps = ntsubsteps * 2
+              endif
 
 !              if (do_print) write(LUNOPRT,*) "newstate::WARNING - Substep failed, retrying with ", ntsubsteps, " substeps."
 
@@ -213,7 +226,7 @@ subroutine newstate_calc(carma, cstate, scale_threshold, rc)
                 if (rc < RC_OK) return
               end do
 
-              rc = RC_OK
+              if (rc == RC_WARNING_RETRY) rc = RC_OK
               takeSteps = .true.
               exit
 
@@ -238,6 +251,8 @@ subroutine newstate_calc(carma, cstate, scale_threshold, rc)
       nsubstep = nsubstep + ntsubsteps
       nretry   = nretry   + nretries
 
+      if (rc == RC_WARNING_PFAST) rc = RC_OK
+
       if (do_substep) zsubsteps(iz) = ntsubsteps
     end do
 

@@ -0,0 +1,68 @@
+! Include shortname defintions, so that the F77 code does not have to be modified to
+! reference the CARMA structure.
+#include "carma_globaer.h"
+
+!!  This routine normalizes dmdt as to not make gas concentration go negative.
+!!
+!! @author Parker Case
+!! @version Mar-2024
+subroutine pfastdmdt(carma, cstate, iz, rc)
+
+  ! types
+  use carma_precision_mod
+  use carma_enums_mod
+  use carma_constants_mod
+  use carma_types_mod
+  use carmastate_mod
+  use carma_mod
+
+  implicit none
+
+  type(carma_type), intent(in)         :: carma       !! the carma object
+  type(carmastate_type), intent(inout) :: cstate      !! the carma state object
+  integer, intent(in)                  :: iz          !! z index
+  integer, intent(inout)               :: rc          !! return code, negative indicates failure
+
+  ! Local declarations
+  real(kind=f)  :: rvap
+  real(kind=f)  :: gc_cgs
+  real(kind=f)  :: gc_target
+  real(kind=f)  :: scalefactor
+  integer       :: igroup
+  integer       :: ielem
+  integer       :: igas
+  integer       :: ibin
+
+  do igroup = 1,NGROUP
+
+    ielem = ienconc(igroup)     ! element of particle number concentration
+
+    igas = igrowgas(ielem)      ! condensing gas
+
+    if ((itype(ielem) == I_VOLATILE) .and. (igas /= 0)) then
+
+      ! Calculate vapor pressures.
+      call vaporp(carma, cstate, iz, igas, rc)
+
+      ! Define gas constant for this gas
+      rvap = RGAS / gwtmol(igas)
+
+      ! Current gas concentration
+      gc_cgs = gc(iz,igas) / (zmet(iz)*xmet(iz)*ymet(iz))
+
+      ! Gas concentration at equilibrium
+      gc_target = pvapl(iz,igas) / (rvap * t(iz))
+
+      ! Determine total mass added to bins in implicit timestep
+      scalefactor = sum((pc(iz,:,ielem) - pcl(iz,:,ielem))*rmass(:,igroup))/(gc_cgs - gc_target)
+
+      ! Loop through bins and apply correction
+      do ibin = 1,NBIN
+        pc(iz,ibin,ielem) = pcl(iz,ibin,ielem) + (pc(iz,ibin,ielem) - pcl(iz,ibin,ielem)) / scalefactor
+      end do
+    end if
+
+  end do
+
+  return
+end
@@ -1088,7 +1088,7 @@ SUBROUTINE CARMA_GridCompRun ( gcCARMA, qa, impChem, expChem, nymd, nhms, &
 ! This is set in the internal_spec but possibly is bootstrapped.
 ! If bootstrapped set to current temperature.
   n = nCARMAbegin + reg%NBIN*reg%NELEM + reg%NGAS
-  if(qa(n)%data3d(i1,j2,km) < 0.) qa(n)%data3d = tmpu
+  if(qa(n)%data3d(i1,j2,km) <= 0.) qa(n)%data3d = tmpu
 ! And same for gases -- first, initialize water vapor to current
   do igas = 1, reg%NGAS
    n  = nCARMAbegin + reg%NBIN*reg%NELEM - 1 + igas

@@ -310,12 +310,8 @@ SUBROUTINE CARMA_Emissions ( gcCARMA, qa, impChem, expChem, nymd, nhms, cdt, &
 !    In this implementation we compute f(u,v,...) in a separate routine
 !    and do the rest here.
 
-!    Read dust source function from file if not already allocated
-     if( .not. associated(gcCARMA%dust_source)) then
-      allocate(gcCARMA%dust_source(i1:i2,j1:j2),__STAT__)
-      gcCARMA%dust_source = 0.
-      call MAPL_GetPointer( impChem, gcCARMA%dust_source, 'CARMA_DU_SRC', __RC__)
-     endif
+!    Read dust source function from file
+     call MAPL_GetPointer( impChem, gcCARMA%dust_source, 'CARMA_DU_SRC', __RC__)
 
 !    Do the emission calculation
 !    The DEAD emission calculation occurs outside the size bins, returning

@@ -987,7 +987,6 @@ SUBROUTINE Finalize_ ( gc, import, export, clock, rc )
 
 !  Destroy emissions
 !  -----------------
-   if(associated(gcCARMA%dust_source))       deallocate( gcCARMA%dust_source, __STAT__)
    if(associated(gcCARMA%vLat))              deallocate( gcCARMA%vLat, __STAT__)
    if(associated(gcCARMA%vLon))              deallocate( gcCARMA%vLon, __STAT__)
    if(associated(gcCARMA%vSO2))              deallocate( gcCARMA%vSO2, __STAT__)

@@ -64,8 +64,9 @@ ICOMPOSITION: 1       2      3
 # GASNAME  = name
 # IGCOMP   = gas composition: 1 (I_GCOMP_H2O), 2 (I_GCOMP_H2SO4), 3 (I_GCOMP_SO2),
 #                             4 (I_GCOMP_HNO3)
-# IGVAPREQ = vapor pressure equation: 1 (I_VAPRTN_H2O_BUCK1981), 2 (I_VAPRTN_H2O_MURPHY2005)
-#                                     3 (I_VAPRTN_H2O_GOFF1946), 4 (I_VAPRTN_H2SO4_AYERS1980)
+# IGVAPREQ = vapor pressure equation: -1 (I_VAPRTN_NULL), 1 (I_VAPRTN_H2O_BUCK1981), 
+#                                     2 (I_VAPRTN_H2O_MURPHY2005), 3 (I_VAPRTN_H2O_GOFF1946),
+#                                     4 (I_VAPRTN_H2SO4_AYERS1980)
 GASNAME:      H2O     H2SO4   HNO3
 IGCOMP:       1       2       4
 IGVAPREQ:     2       4       -1