Merge branch 'develop' into feature/mmanyin/restore_blend_option

Author: sdrabenh

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/ConvPar_GF2020.F90

@@ -3514,9 +3514,9 @@ SUBROUTINE CUP_gf(its,ite,kts,kte ,itf,ktf, mtp, nmp &
                            zdo(i,k  )*(hcdo(i,k  )-heo_cup(i,k  ) ) )*g/dp*edto(i)
 
             !---meltglac-------------------------------------------------
-            dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))*0.5*(qrco(i,k+1)+qrco(i,k)) &
-                                              - melting(i,k))*g/dp
-
+            !dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))*0.5*(qrco(i,k+1)+qrco(i,k)) &
+            !                                  - melting(i,k))*g/dp
+            dellah(i,k) = dellah(i,k) - xlf*melting(i,k)*g/dp
             !-- for output only
             subten_H(i,k) = -(zuo(i,k+1)*(-heo_cup(i,k+1)) - zuo(i,k)*(-heo_cup(i,k)))*g/dp       &
                             +(zdo(i,k+1)*(-heo_cup(i,k+1)) - zdo(i,k)*(-heo_cup(i,k)))*g/dp*edto(i)
@@ -3678,9 +3678,9 @@ SUBROUTINE CUP_gf(its,ite,kts,kte ,itf,ktf, mtp, nmp &
                          +(zdo(i,k+1)*(hcdo(i,k+1)-heo_cup(i,k+1) ) -                    &
                            zdo(i,k  )*(hcdo(i,k  )-heo_cup(i,k  ) ) )*g/dp*edto(i)
 
-            dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))* &
-                                       0.5*(qrco(i,k+1)+qrco(i,k)) - melting(i,k))*g/dp
-
+            !dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))* &
+            !                           0.5*(qrco(i,k+1)+qrco(i,k)) - melting(i,k))*g/dp
+            dellah(i,k) = dellah(i,k) - xlf*melting(i,k)*g/dp
             !--- for output only
             subten_H(i,k) = -(zuo(i,k+1)*(-heo_cup(i,k+1)) - zuo(i,k)*(-heo_cup(i,k)))*g/dp       &
                             +(zdo(i,k+1)*(-heo_cup(i,k+1)) - zdo(i,k)*(-heo_cup(i,k)))*g/dp*edto(i)
@@ -3707,8 +3707,9 @@ SUBROUTINE CUP_gf(its,ite,kts,kte ,itf,ktf, mtp, nmp &
             dellah(i,k) =-( zuo(i,k+1)*hco (i,k+1) - zuo(i,k)*hco (i,k) )*g/dp           &
                          +( zdo(i,k+1)*hcdo(i,k+1) - zdo(i,k)*hcdo(i,k) )*g/dp*edto(i)
 
-            dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))* &
-                                        0.5*(qrco(i,k+1)+qrco(i,k)) - melting(i,k))*g/dp
+            !dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))* &
+            !                            0.5*(qrco(i,k+1)+qrco(i,k)) - melting(i,k))*g/dp
+            dellah(i,k) = dellah(i,k) - xlf*melting(i,k)*g/dp
             !- update with subsidence term from the FCT scheme
             dellah(i,k) = dellah(i,k) + sub_tend(1,k)
             !--- for output only

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90

@@ -32,6 +32,13 @@ module GEOSmoist_Process_Library
     module procedure ICE_FRACTION_1D
     module procedure ICE_FRACTION_SC
   end interface ICE_FRACTION
+
+  ! SRF_TYPE constants
+  integer, parameter :: SRF_TYPE_LAND  = 1
+  integer, parameter :: SRF_TYPE_SNOW  = 2
+  integer, parameter :: SRF_TYPE_ICE   = 3
+  integer, parameter :: SRF_TYPE_OCEAN = 0
+
   ! ICE_FRACTION constants
    ! In anvil/convective clouds
    real, parameter :: aT_ICE_ALL = 252.16
@@ -403,7 +410,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
       ICEFRCT_C = MAX(ICEFRCT_C,0.00)
       ICEFRCT_C = ICEFRCT_C**aICEFRPWR
      ! Sigmoidal functions like figure 6b/6c of Hu et al 2010, doi:10.1029/2009JD012384
-      if (SRF_TYPE >= 2.0) then
+      select case (nint(SRF_TYPE))
+      case (SRF_TYPE_SNOW, SRF_TYPE_ICE)
         ! Over snow (SRF_TYPE == 2.0) and ice (SRF_TYPE == 3.0)
         if (ICE_RADII_PARAM == 1) then
           ! Jason formula
@@ -424,8 +432,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**iICEFRPWR
-      else if (SRF_TYPE > 1.0) then
-        ! Over Land
+      case (SRF_TYPE_LAND)
+        ! Over Land (SRF_TYPE == 1)
         ICEFRCT_M  = 0.00
         if ( TEMP <= lT_ICE_ALL ) then
            ICEFRCT_M = 1.000
@@ -435,8 +443,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**lICEFRPWR
-      else
-        ! Over Oceans
+      case (SRF_TYPE_OCEAN)
+        ! Over Oceans (SRF_TYPE == 0)
         ICEFRCT_M  = 0.00
         if ( TEMP <= oT_ICE_ALL ) then
            ICEFRCT_M = 1.000
@@ -446,7 +454,11 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**oICEFRPWR
-      endif
+      case default
+        ! You should not be here
+        print *, 'ICE_FRACTION_SC: Unknown SRF_TYPE = ',SRF_TYPE
+        error stop
+      end select
       ! Combine the Convective and MODIS functions
         ICEFRCT  = ICEFRCT_M*(1.0-CNV_FRACTION) + ICEFRCT_C*(CNV_FRACTION)
 #endif