Merge branch 'develop' into release/MAPL-v3

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt

@@ -29,7 +29,7 @@ endif ()
 #      and 10 minutes at O2. But only 7 seconds with O1. So we compile at O1
 if (CMAKE_Fortran_COMPILER_ID MATCHES Intel AND CMAKE_BUILD_TYPE MATCHES Release)
   set_source_files_properties(GEOS_BACM_1M_InterfaceMod.F90 PROPERTIES COMPILE_OPTIONS ${FOPT1})
-  set_source_files_properties(GEOS_MGB2_2M_InterfaceMod.F90 PROPERTIES COMPILE_OPTIONS ${FOPT1})
+ # set_source_files_properties(GEOS_MGB2_2M_InterfaceMod.F90 PROPERTIES COMPILE_OPTIONS ${FOPT1})
 endif ()
 
 esma_add_library (${this}

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90

@@ -4778,7 +4778,7 @@ subroutine SetServices ( GC, RC )
     VERIFY_(STATUS)
 
     call MAPL_AddExportSpec(GC,                             &
-         SHORT_NAME='QCVAR_EXP',                                   &
+         SHORT_NAME='QCVAR',                                   &
          LONG_NAME ='inverse relative variance of cloud water',        &
          UNITS      = '1',                                    &
          DIMS      = MAPL_DimsHorzOnly,                            &
@@ -5222,7 +5222,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
     real, allocatable, dimension(:,:)   :: TMP2D
     ! Internals
     real, pointer, dimension(:,:,:) :: Q, QLLS, QLCN, CLLS, CLCN, QILS, QICN
-    real, pointer, dimension(:,:,:) :: NACTL, NACTI
+    real, pointer, dimension(:,:,:) :: NACTL, NACTI, NCPL, NCPI
     ! Imports
     real, pointer, dimension(:,:,:) :: ZLE, PLE, T, U, V, W
     real, pointer, dimension(:,:)   :: FRLAND, FRLANDICE, FRACI, SNOMAS
@@ -5295,6 +5295,8 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
        call MAPL_GetPointer(INTERNAL, QICN,     'QICN'    , RC=STATUS); VERIFY_(STATUS)
        call MAPL_GetPointer(INTERNAL, NACTL,   'NACTL'    , RC=STATUS); VERIFY_(STATUS)
        call MAPL_GetPointer(INTERNAL, NACTI,   'NACTI'    , RC=STATUS); VERIFY_(STATUS)
+       call MAPL_GetPointer(INTERNAL, NCPL,   'NACTL'    , RC=STATUS); VERIFY_(STATUS)
+       call MAPL_GetPointer(INTERNAL, NCPI,   'NACTI'    , RC=STATUS); VERIFY_(STATUS)
 
        ! Import State
        call MAPL_GetPointer(IMPORT, PLE,     'PLE'     , RC=STATUS); VERIFY_(STATUS)
@@ -5559,7 +5561,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
            end where
          endif
 
-         if (.FALSE.) then
+         if (adjustl(CLDMICR_OPTION)=="MGB2_2M") then
           QST3  = GEOS_QsatLQU (T, PLmb*100.0, DQ=DQST3) !clean up only with respect to liquid water
          else
           DQST3 = GEOS_DQSAT   (T, PLmb, QSAT=QST3)      ! this qsat function expects hPa...

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90

@@ -1115,18 +1115,21 @@ subroutine fix_up_clouds_2M( &
          QG, &
          NR, &
          NS, &
-         NG)
+         NG, &
+         MASS, & 
+         TMP2D)
 
       real, intent(inout), dimension(:,:,:) :: TE,QV,QLC,CF,QLA,AF,QIC,QIA, QR, QS, QG
       real, intent(inout), dimension(:,:,:) :: NI, NL, NS, NR, NG
+      real, dimension(:,:,:),   intent(in)     :: MASS
+      real, dimension(:,:),     intent(  out)  :: TMP2D
+      integer :: IM, JM, LM
 
       real, parameter  :: qmin  = 1.0e-12
       real, parameter :: cfmin  = 1.0e-4
       real, parameter :: nmin  = 100.0
 
 
-
-
       ! Fix if Anvil cloud fraction too small
       where (AF < cfmin)
          QV  = QV + QLA + QIA
@@ -1187,9 +1190,22 @@ subroutine fix_up_clouds_2M( &
          QLC = 0.
          QIC = 0.
       end where
-
-
-
+
+        IM = SIZE( QV, 1 )
+    	JM = SIZE( QV, 2 )
+    	LM = SIZE( QV, 3 )
+
+
+      !make sure QI , NI stay within T limits 
+         call meltfrz_inst2M  ( IM, JM, LM,    &
+              TE              , &
+              QLC          , &
+              QLA         , &
+              QIC           , &
+              QIA          , &               
+              NL         , &
+              NI          )
+
       !make sure no negative number concentrations are passed
       !and that N goes to minimum defaults in the microphysics when mass is too small
 
@@ -1208,6 +1224,18 @@ subroutine fix_up_clouds_2M( &
       where (QS .le. qmin) NS = 0.
 
       where (QG .le. qmin) NG = 0.
+
+      ! need to clean up small negative values. MG does can't handle them
+          call FILLQ2ZERO( QV, MASS, TMP2D) 
+          call FILLQ2ZERO( QG, MASS, TMP2D) 
+          call FILLQ2ZERO( QR, MASS, TMP2D) 
+          call FILLQ2ZERO( QS, MASS, TMP2D) 
+          call FILLQ2ZERO( QLC, MASS, TMP2D)
+          call FILLQ2ZERO( QLA, MASS, TMP2D)  
+          call FILLQ2ZERO( QIC, MASS, TMP2D)
+          call FILLQ2ZERO( QIA, MASS, TMP2D)
+          call FILLQ2ZERO( CF, MASS, TMP2D)
+          call FILLQ2ZERO( AF, MASS, TMP2D)
 
    end subroutine fix_up_clouds_2M
 
@@ -2279,17 +2307,17 @@ end subroutine hystpdf
 
 !==========Estimate RHcrit========================
 !==============================
- subroutine pdf_alpha(PP,P_LM, ALPHA, FRLAND, MINRHCRIT, TURNRHCRIT, EIS, RHC_OPTION)
+ subroutine pdf_alpha(PP,P_LM, ALPHA, FRLAND, MINRHCRIT, TURNRHCRIT, TURNRHCRIT_UPPER, EIS, RHC_OPTION)
 
       real,    intent(in)  :: PP, P_LM !mbar
       real,    intent(out) :: ALPHA
       real,    intent(in)  :: FRLAND
-      real,    intent(in)  :: MINRHCRIT, TURNRHCRIT, EIS
+      real,    intent(in)  :: MINRHCRIT, TURNRHCRIT, EIS, TURNRHCRIT_UPPER
       integer, intent(in)  :: RHC_OPTION !0-Slingo(1985), 1-QUAAS (2012)
       real :: dw_land = 0.20 !< base value for subgrid deviation / variability over land
       real :: dw_ocean = 0.10 !< base value for ocean
       real :: sloperhcrit =20.
-      real :: TURNRHCRIT_UPPER = 300.
+      !real :: TURNRHCRIT_UPPER = 300.
       real ::  aux1, aux2, maxalpha
 
       IF (RHC_OPTION .lt. 1) then
@@ -2307,9 +2335,13 @@ subroutine pdf_alpha(PP,P_LM, ALPHA, FRLAND, MINRHCRIT, TURNRHCRIT, EIS, RHC_OPT
              aux1 = 1.0/(1.0+exp(aux1)) !this function reproduces the old Sligo function.
           end if
 
-          !aux2= 1.0/(1.0+exp(aux2)) !this function would reverse the profile P< TURNRHCRIT_UPPER
-           aux2=1.0
-           ALPHA  = min(maxalpha*aux1*aux2, 0.3)
+           if (TURNRHCRIT_UPPER .gt. 0.0) then 
+          	 aux2= 1.0/(1.0+exp(aux2)) !this function reverses the profile P< TURNRHCRIT_UPPER
+           else
+            aux2=1.0
+           end if 
+
+           ALPHA  = min(maxalpha*aux1*aux2, 0.4)
 
        ELSE
            ! based on Quass 2012 https://doi.org/10.1029/2012JD017495
@@ -3316,9 +3348,12 @@ subroutine update_cld( &
       QSLIQ  = GEOS_QsatLQU( TE, PL*100.0 , DQ=DQx )
       QSICE  = GEOS_QsatICE( TE, PL*100.0 , DQ=DQX )
 
-      if ((QC+QA) .gt. 1.0e-13) then
-         QSx=((QCl+QAl)*QSLIQ + QSICE*(QCi+QAi))/(QC+QA)
-      else
+
+      IF (QCl + QAl .gt. 0.) then 
+        QSx =  QSLIQ
+      ELSEIF (QCi + QAi.gt. 0.) then 
+        QSx =  QSICE        
+      ELSE
 		 DQSx  = GEOS_DQSAT( TE, PL, QSAT=QSx )
       end if
 
@@ -3355,8 +3390,10 @@ subroutine update_cld( &
       if  (QSx .gt. tiny(1.0)) then
          RHCmicro = SCICE - 0.5*DELQ/Qsx
       else
-         RHCmicro = 0.0
+         RHCmicro = 1.0-ALPHA
       end if
+
+      RHCmicro =  max(min(RHCmicro, 0.99), 0.6)
 
       CFALL   = max(CFo, 0.0)
       CFALL   = min(CFo, 1.0)
@@ -3370,8 +3407,7 @@ end subroutine update_cld
 
 
 
-   subroutine meltfrz_inst2M  (     &
-         IM,JM,LM , &
+   subroutine meltfrz_inst2M  ( IM, JM, LM,    &
          TE       , &
          QCL       , &
          QAL       , &
@@ -3380,8 +3416,8 @@ subroutine meltfrz_inst2M  (     &
          NL       , &
          NI            )
 
-      integer, intent(in)                             :: IM,JM,LM
       real ,   intent(inout), dimension(:,:,:)   :: TE,QCL,QCI, QAL, QAI, NI, NL
+      integer, intent(in) :: IM, JM, LM
 
       real ,   dimension(im,jm,lm)              :: dQil, DQmax, QLTOT, QITOT, dNil, FQA
       real :: T_ICE_ALL =  240.
@@ -3391,8 +3427,7 @@ subroutine meltfrz_inst2M  (     &
       QLTOT=QCL + QAL
       FQA = 0.0
 
-
-      where (QITOT+QLTOT .gt. 0.0)
+      where (QITOT+QLTOT .gt. tiny(0.0))
          FQA= (QAI+QAL)/(QITOT+QLTOT)
       end where
 
@@ -3411,7 +3446,7 @@ subroutine meltfrz_inst2M  (     &
          dNil = NL
       end where
 
-      where ((dQil .gt. DQmax) .and. (dQil .gt. 0.0))
+      where ((dQil .gt. DQmax) .and. (dQil .gt. tiny(0.0)))
          dNil  =  NL*DQmax/dQil
       end where
 
@@ -3435,7 +3470,7 @@ subroutine meltfrz_inst2M  (     &
       where ((dQil .le. DQmax) .and. (dQil .gt. 0.0))
          dNil = NI
       end where
-      where ((dQil .gt. DQmax) .and. (dQil .gt. 0.0))
+      where ((dQil .gt. DQmax) .and. (dQil .gt. tiny(0.0)))
          dNil  =  NI*DQmax/dQil
       end where
       dQil = max(  0., dQil )

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/WSUB_ExtData.yaml

@@ -1,7 +1,7 @@
 Collections:
   WSUB_SWclim_2005%m2.nc4:
     template: /discover/nobackup/dbarahon/DEV/SWclim/L72/SWclim_2005%m2.nc4
-    valid_range: "2005-01-01/2005-12-31"
+    valid_range: "2005-01-01/2005-12-01"
 
 Samplings:
   WSUB_sample_0: