Merge GCM v12 codebase into main

.circleci/config.yml

@@ -1,11 +1,11 @@
 version: 2.1
 
 # Anchors in case we need to override the defaults from the orb
-#baselibs_version: &baselibs_version v7.27.0
-#bcs_version: &bcs_version v11.6.0
+#baselibs_version: &baselibs_version v8.24.0
+#bcs_version: &bcs_version v12.0.0
 
 orbs:
-  ci: geos-esm/circleci-tools@4
+  ci: geos-esm/circleci-tools@5
 
 workflows:
   build-test:
@@ -21,7 +21,10 @@ workflows:
           #baselibs_version: *baselibs_version
           repo: GEOSgcm
           checkout_fixture: true
-          mepodevelop: true
+          # V12 code uses a special branch for now.
+          fixture_branch: feature/sdrabenh/gcm_v12
+          # We comment out this as it will "undo" the fixture_branch
+          #mepodevelop: true
           persist_workspace: true # Needs to be true to run fv3/gcm experiment, costs extra
 
       # Run AMIP GCM (1 hour, no ExtData)

.github/workflows/workflow.yml

@@ -18,6 +18,7 @@ concurrency:
 jobs:
   build_gcm:
     strategy:
+      fail-fast: false
       matrix:
         compiler: [ifort, gfortran-14, gfortran-15]
         build-type: [Debug]
@@ -26,6 +27,7 @@ jobs:
       compiler: ${{ matrix.compiler }}
       cmake-build-type: ${{ matrix.build-type }}
       fixture-repo: GEOS-ESM/GEOSgcm
+      fixture-ref: feature/sdrabenh/gcm_v12
 
   spack_build:
     uses: GEOS-ESM/CI-workflows/.github/workflows/spack_gcc_build.yml@project/geosgcm
@@ -34,4 +36,6 @@ jobs:
       BUILDCACHE_TOKEN: ${{ secrets.BUILDCACHE_TOKEN }}
     with:
       fixture-repo: GEOS-ESM/GEOSgcm
+      fixture-ref: feature/sdrabenh/gcm_v12
+      load-fms: true
 

GEOS_Shared/CMakeLists.txt

@@ -4,7 +4,7 @@ set (srcs
   GEOS_Utilities.F90 DragCoefficients.F90 GEOS.F90
   tropovars.F90 surfacelayer.F90 GEOS_TopoGet.F90
   GEOS_GmapMod.F90 GEOS_RemapMod.F90 G3_MPI_Util_Mod.F90 windfix.F90
-  PWSSSP.F sphere.F calcdbz.F OVP.F90
+  PWSSSP.F sphere.F calcdbz.F compute_radar_reflectivity.F90 OVP.F90
   getco2.F90 atmOceanIntLayer.F90
   )
 

GEOS_Shared/GEOS_Utilities.F90

@@ -546,7 +546,7 @@ function QSAT0(TL,PL,RAMP,PASCALS,DQSAT) result(QSAT)
        end if
 
        TI = (TI - TMINTBL)*DEGSUBS+1
-       IT = int(TI)
+       IT = min(max(1,int(TI)),TABLESIZE-1)
 
        if(URAMP==TMIX) then
           DQ    = ESTBLX(IT+1) - ESTBLX(IT)
@@ -764,7 +764,7 @@ function DQSAT0(TL,PL,RAMP,PASCALS,QSAT) result(DQSAT)
          end if
 
          TT = (TI - TMINTBL)*DEGSUBS+1
-         IT = int(TT)
+         IT = min(max(1,int(TT)),TABLESIZE-1)
 
          if(URAMP==TMIX) then
             DQQ =  ESTBLX(IT+1) - ESTBLX(IT)

GEOS_Shared/calcdbz.F

@@ -25,7 +25,10 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
 c    calculated in Thompson, Rasmussen, and Manning (2004, Monthly
 c    Weather Review, Vol. 132, No. 2, pp. 519-542.)
 c
-c     4. If iliqskin=1, frozen particles that are at a temperature above
+c     4. If iliqskin=1, snow particles that are at a temperature above
+c     freezing are assumed to scatter as a liquid particle.
+c
+c     4. If iliqskin=2, snow & graupel particles that are at a temperature above
 c     freezing are assumed to scatter as a liquid particle.
 c
 c     More information on the derivation of simulated reflectivity in
@@ -51,8 +54,10 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
       REAL TEMP_C,VIRTUAL_T
       REAL GONV,RONV,SONV
       REAL FACTOR_G,FACTOR_R,FACTOR_S
-      REAL FACTORB_G,FACTORB_R,FACTORB_S
+      REAL FACTORB_G,FACTORB_S
       REAL RHOAIR,Z_E
+      REAL ZE_R,ZE_S,ZE_G
+
 
 c   Constants used to calculate variable intercepts
       REAL R1,RON,RON2,SON,GON
@@ -65,13 +70,11 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
       REAL GAMMA_SEVEN,ALPHA
       REAL RHOWAT,CELKEL,PI,RD
 
-
 c   Constants used to calculate variable intercepts
       R1 = 1.D-15
       RON = 8.D6
       RON2 = 1.D10
       SON = 2.D7
-      GON = 5.D7
       RON_MIN = 8.D6
       RON_QR0 = 0.00010D0
       RON_DELQR0 = 0.25D0*RON_QR0
@@ -81,19 +84,28 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
 c   Constant intercepts
       RN0_R = 8.D6
       RN0_S = 2.D7
-      RN0_G = 4.D6
 
 c   Other constants
       GAMMA_SEVEN = 720.D0
       RHOWAT = 1000.D0
       RHO_R = RHOWAT
       RHO_S = 100.D0
-      RHO_G = 400.D0
       ALPHA = 0.224D0
       CELKEL = 273.15D0
       PI = 3.141592653589793D0
       RD = 287.04D0
 
+c   --- NEW: HAIL VS GRAUPEL CONSTANTS ---
+      IF (IVARINT .GE. 3) THEN
+          RHO_G = 917.D0   ! Hail density (kg/m^3)
+          GON   = 4.D4     ! Hail max intercept (fewer, much larger stones)
+          RN0_G = 4.D4     ! Constant hail intercept
+      ELSE
+          RHO_G = 400.D0   ! Graupel density (kg/m^3)
+          GON   = 5.D7     ! Graupel max intercept (many, smaller stones)
+          RN0_G = 4.D6     ! Constant graupel intercept
+      END IF
+
 c   Force all Q arrays to be 0.0 or greater.
       DO K = 1,BTDIM
          DO J = 1,SNDIM
@@ -150,6 +162,8 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
 
                   IF (ILIQSKIN.EQ.1 .AND. TMK(I,J,K).GT.CELKEL) THEN
                       FACTORB_S = FACTOR_S/ALPHA
+                  ELSEIF (ILIQSKIN.EQ.2 .AND. TMK(I,J,K).GT.CELKEL) THEN
+                      FACTORB_S = FACTOR_S/ALPHA
                       FACTORB_G = FACTOR_G/ALPHA
                   ELSE
                       FACTORB_S = FACTOR_S
@@ -158,11 +172,19 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
 
 c      Calculate variable intercept parameters
 
-                  IF (IVARINT.EQ.1) THEN
+                  IF (IVARINT.GE.1) THEN
 
                       TEMP_C = DMIN1(-0.001D0,TMK(I,J,K)-CELKEL)
                       SONV = DMIN1(2.0D8,2.0D6*EXP(-0.12D0*TEMP_C))
 
+c                     simulate large snow aggregates by halving the intercept parameter.
+c                     Apply this to both option 2 and 3!
+                      IF (IVARINT.GE.2) THEN
+                          IF (QSN(I,J,K).GT.R1) THEN
+                              SONV = SONV * 0.5D0
+                          END IF
+                      END IF
+
                       GONV = GON
                       IF (QGR(I,J,K).GT.R1) THEN
                           GONV = 2.38D0* (PI*RHO_G/
@@ -184,14 +206,31 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
 
                   END IF
 
+
 c      Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is
 c      the sum of z_e for each hydrometeor species:
 
-                  Z_E = FACTOR_R* (RHOAIR*QRA(I,J,K))**1.75D0/
-     +                  RONV**.75D0 + FACTORB_S*
-     +                  (RHOAIR*QSN(I,J,K))**1.75D0/SONV**.75D0 +
-     +                  FACTORB_G* (RHOAIR*QGR(I,J,K))**1.75D0/
-     +                  GONV**.75D0
+c      Calculate individual equivalent reflectivity factors
+                  ZE_R = (RHOAIR*QRA(I,J,K))**1.75D0
+                  ZE_R = FACTOR_R * ZE_R / RONV**.75D0
+
+                  ZE_S = (RHOAIR*QSN(I,J,K))**1.75D0
+                  ZE_S = FACTORB_S * ZE_S / SONV**.75D0
+
+                  ZE_G = (RHOAIR*QGR(I,J,K))**1.75D0
+                  ZE_G = FACTORB_G * ZE_G / GONV**.75D0
+
+c                 --- PSEUDO-MIE PLATEAU FOR WET HAIL ---
+c                 If using Hail (IVARINT>=3) with a liquid skin, Rayleigh math overestimates.
+c                 Cap the hail's specific Ze contribution at 10^7 (approx 70 dBZ).
+                  IF (IVARINT.GE.3 .AND. ILIQSKIN.EQ.2 .AND. 
+     +                TMK(I,J,K).GT.CELKEL) THEN
+                      ZE_G = MIN(ZE_G, 1.0D7)
+                  END IF
+
+c      Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is
+c      the sum of z_e for each hydrometeor species:
+                  Z_E = ZE_R + ZE_S + ZE_G
 
 c      Adjust small values of Z_e so that dBZ is no lower than -30
                   Z_E = MAX(Z_E,.001D0)