Merge GCM v12 codebase into develop

.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.29.0
-#bcs_version: &bcs_version v11.6.0
+#baselibs_version: &baselibs_version v8.5.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)
@@ -45,7 +48,7 @@ workflows:
             - docker-hub-creds
           matrix:
             parameters:
-              compiler: [gfortran, ifort]
+              compiler: [ifort]
           requires:
             - build-GEOSgcm-on-<< matrix.compiler >>
           repo: GEOSgcm

GEOS_RadiationGridComp.F90

@@ -241,6 +241,15 @@ subroutine SetServices ( GC, RC )
                                                               RC=STATUS  )
     VERIFY_(STATUS)
 
+    call MAPL_AddExportSpec ( GC,                                   &
+         SHORT_NAME = 'RADLWNA',                                         &
+         LONG_NAME  = 'air_temperature_tendency_due_to_longwave_no_aerosol', &
+         UNITS      = 'K s-1',                                           &
+         DIMS       = MAPL_DimsHorzVert,                                 &
+         VLOCATION  = MAPL_VLocationCenter,                              &
+                                                              RC=STATUS  )
+    VERIFY_(STATUS)
+
     call MAPL_AddExportSpec ( GC,                                   &
          SHORT_NAME = 'RADSWNA',                                         &
          LONG_NAME  = 'air_temperature_tendency_due_to_shortwave_no_aerosol', &
@@ -639,6 +648,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
   real, pointer, dimension(:,:  )     :: DSFDTS
   real, pointer, dimension(:,:  )     :: TRD
   real, pointer, dimension(:,:,:)     :: FLW
+  real, pointer, dimension(:,:,:)     :: FLWNA
   real, pointer, dimension(:,:,:)     :: FLWCLR
   real, pointer, dimension(:,:,:)     :: FLA
 
@@ -660,6 +670,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
   real, pointer, dimension(:,:,:)     :: RADSW
   real, pointer, dimension(:,:,:)     :: RADLWC
   real, pointer, dimension(:,:,:)     :: RADSWC
+  real, pointer, dimension(:,:,:)     :: RADLWNA
   real, pointer, dimension(:,:,:)     :: RADSWNA
   real, pointer, dimension(:,:,:)     :: RADLWCNA
   real, pointer, dimension(:,:,:)     :: RADSWCNA
@@ -730,6 +741,8 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
     VERIFY_(STATUS)
     call MAPL_GetPointer ( EXPORT, RADSWC  , 'RADSWC'  ,  RC=STATUS )
     VERIFY_(STATUS)
+    call MAPL_GetPointer ( EXPORT, RADLWNA , 'RADLWNA' ,  RC=STATUS )
+    VERIFY_(STATUS)
     call MAPL_GetPointer ( EXPORT, RADSWNA , 'RADSWNA' ,  RC=STATUS )
     VERIFY_(STATUS)
     call MAPL_GetPointer ( EXPORT, RADLWCNA, 'RADLWCNA',  RC=STATUS )
@@ -770,6 +783,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
        VERIFY_(STATUS)
     end if
 
+    if (associated(RADLWNA)                           ) then
+       call MAPL_GetPointer ( GEX(IRR), FLWNA,  'FLXA'   ,  alloc=.TRUE.,RC=STATUS )
+       VERIFY_(STATUS)
+    end if
+
     if (associated(RADLWCNA)                          ) then
        call MAPL_GetPointer ( GEX(IRR), FLA   , 'FLA'    ,  alloc=.TRUE.,RC=STATUS )
        VERIFY_(STATUS)
@@ -804,7 +822,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
     if( associated (RADLW ) .or. associated (RADSW )   .or. &
         associated (RADLWC) .or. associated (RADSWC)   .or. &
         associated (RADSWNA).or. associated (RADSWCNA) .or. &
-        associated (RADLWCNA)                             ) then
+        associated (RADLWNA).or. associated (RADLWCNA) ) then
 
        allocate(DMI(IM,JM,LM),stat=STATUS)
        VERIFY_(STATUS)
@@ -814,6 +832,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC )
        if( associated (RADSW   ) ) RADSW    = (FSW   (:,:,0:LM-1) - FSW   (:,:,1:LM))*DMI
        if( associated (RADLWC  ) ) RADLWC   = (FLWCLR(:,:,0:LM-1) - FLWCLR(:,:,1:LM))*DMI
        if( associated (RADSWC  ) ) RADSWC   = (FSWCLR(:,:,0:LM-1) - FSWCLR(:,:,1:LM))*DMI
+       if( associated (RADLWNA ) ) RADLWNA  = (FLWNA (:,:,0:LM-1) - FLWNA (:,:,1:LM))*DMI
        if( associated (RADSWNA ) ) RADSWNA  = (FSWNA (:,:,0:LM-1) - FSWNA (:,:,1:LM))*DMI
        if( associated (RADLWCNA) ) RADLWCNA = (FLA   (:,:,0:LM-1) - FLA   (:,:,1:LM))*DMI
        if( associated (RADSWCNA) ) RADSWCNA = (FSCNA (:,:,0:LM-1) - FSCNA (:,:,1:LM))*DMI

GEOS_RadiationShared/CMakeLists.txt

@@ -3,6 +3,8 @@ esma_set_this ()
 set (srcs
   radconstants.F90
   gettau.F90
+  rad_types.F90
+  rad_utils.F90
   cloud_condensate_inhomogeneity.F90
   cloud_subcol_gen.F90
   )

GEOS_RadiationShared/RRTMGP_cmake/CMakeLists.txt

@@ -57,9 +57,12 @@ if (MKL_FOUND)
     rng/mo_rng_mklvsl.F90
     rng/mo_rng_mklvsl_plus.F90
     )
-
+  message(STATUS "MKL found. The mklvsl rng will be compiled.")
 else ()
-  ecbuild_warn("MKL was not found. The mklvsl rng will not be compiled. RRTMGP should not be used.")
+  list (APPEND SRCS
+    rng/mo_rng_mt19937.F90
+    )
+  message(STATUS "MKL not found. The mt19937 rng will be compiled.")
 endif()
 
 esma_add_library(${this}
@@ -69,9 +72,9 @@ esma_add_library(${this}
 # Use of the GEOS Vectorized flags in RRTMGP caused a segfault leading to
 # line 726 of mo_gas_optics_kernels.F90. Moving to the "old" no-vect flags
 # allows RRTMGP to run again.
-if (CMAKE_Fortran_COMPILER_ID MATCHES Intel AND CMAKE_BUILD_TYPE MATCHES Release)
-  set (CMAKE_Fortran_FLAGS_RELEASE "${GEOS_Fortran_FLAGS_NOVECT}")
-endif ()
+#if (CMAKE_Fortran_COMPILER_ID MATCHES Intel AND CMAKE_BUILD_TYPE MATCHES Release)
+#  set (CMAKE_Fortran_FLAGS_RELEASE "${GEOS_Fortran_FLAGS_NOVECT}")
+#endif ()
 
 if (MKL_FOUND)
   target_include_directories (${this} PRIVATE ${MKL_INCLUDE_DIRS})

GEOS_RadiationShared/rad_types.F90

@@ -0,0 +1,9 @@
+module rad_types
+
+  ! Handy wrappers to pointer arrays.
+  ! Used to provide ragged arrays, etc.
+  type :: rptr1d_wrap
+    real, pointer, dimension(:) :: p
+  end type rptr1d_wrap
+
+end module rad_types

GEOS_RadiationShared/rad_utils.F90

@@ -0,0 +1,234 @@
+#include "MAPL_Generic.h"
+
+module rad_utils
+
+  use ESMF
+  use MAPL
+
+  implicit none
+  private
+
+  public :: Tbr_from_band_flux
+  public :: choose_solar_scheme, choose_irrad_scheme
+
+contains
+
+  ! -----------------------------------------------------
+  ! Routines to estimate brightness temperature from flux
+  ! -----------------------------------------------------
+
+  ! estimate brightness temperature from a band flux
+  subroutine Tbr_from_band_flux(IM, JM, Fband_, wn1, wn2, Tbr_, RC)
+
+    ! input arguments
+    integer, intent(in ) :: IM, JM
+    real,    intent(in ) :: Fband_(IM,JM) ! band flux [W/m2]
+    real,    intent(in ) :: wn1, wn2      ! bounds of band [m-1]
+
+    ! output arguments
+    real,    intent(out) :: Tbr_(IM,JM)   ! brightness temp [K]
+
+    ! error code
+    integer, optional, intent(out) :: RC
+
+    ! fundamental constants
+    double precision, parameter :: h  = 6.626070040d-34  ! Plancks constant         [J.s]
+    double precision, parameter :: c  = 2.99792458d8     ! Speed of light in vacuum [m/s]
+    double precision, parameter :: kB = 1.38064852d-23   ! Boltzmann constant       [J/K]
+    double precision, parameter :: pi = MAPL_PI_R8
+
+    ! other constants
+    double precision, parameter :: alT = h * c / kB
+    double precision, parameter :: bigS = 2.0d0 * kB**4 * pi / (h**3 * c**2)
+    double precision, parameter :: bigC = 2.0d0 * h * c**2
+
+    ! locals
+    integer :: STATUS
+    double precision, dimension(IM,JM) :: Fband, Tbr, Bmean
+    real :: wnMid
+
+    _ASSERT(wn2 > wn1,'band wavenumber bounds mis-ordered!')
+
+    ! calculations done in double precision
+    Fband = dble(Fband_)
+
+    ! first guess Tbr from narrow band approximation ...
+    ! (1) estimate mean Planck function for a narrow band
+    Bmean = Fband / (pi * (wn2 - wn1))
+    ! (2) invert Planck function for temp at mid-point wavenumber
+    wnMid = (wn1 + wn2) / 2.0d0
+    call invert_Planck_for_T(IM, JM, Bmean, wnMid, bigC, alT, Tbr, __RC__)
+
+    ! now refine with a wide band esimate
+    ! PMN: Iterative routine not ready for prime time
+    !      Produces erroneously large Tbr in cloudy regions
+    !call Tbr_wide_band(IM, JM, Fband, wn1, wn2, bigS, alT, Tbr, __RC__)
+
+    ! put output back in real
+    where (Tbr > 0.0d0)
+      Tbr_ = real(Tbr)
+    elsewhere
+      Tbr_ = MAPL_UNDEF
+    endwhere
+
+    _RETURN(_SUCCESS)
+
+  end subroutine Tbr_from_band_flux
+
+  ! invert Planck function for temperature
+  subroutine invert_Planck_for_T(IM, JM, Bwn, wn, bigC, alT, T, RC)
+
+    ! input arguments
+    integer,          intent(in ) :: IM, JM
+    double precision, intent(in ) :: Bwn(IM,JM)  ! PlanckFn(wavenumber)
+    real,             intent(in ) :: wn          ! wavenumber [m-1]
+    double precision, intent(in ) :: bigC, alT   ! necessary constants
+
+    ! output arguments
+    double precision, intent(out) :: T(IM,JM)    ! temperature [K]
+
+    ! error code
+    integer, optional, intent(out) :: RC
+
+    ! error checking
+    _ASSERT(wn > 0.,'non-positive wavenumber!')
+
+    ! invert Planck function for temp
+    where (Bwn > 0.0d0)
+      T = (alT * wn) / log((bigC * wn**3) / Bwn + 1.0d0)
+    elsewhere
+      T = 0.0d0
+    endwhere
+
+    _RETURN(_SUCCESS)
+
+  end subroutine invert_Planck_for_T
+
+  ! Tbr from wide band approximation
+  subroutine Tbr_wide_band(IM, JM, Fband, wn1, wn2, bigS, alT, Tbr, RC)
+
+    ! input arguments
+    integer,          intent(in   ) :: IM, JM
+    double precision, intent(in   ) :: Fband(IM,JM)  ! band flux [W/m2]
+    real,             intent(in   ) :: wn1, wn2      ! bounds of band [m-1]
+    double precision, intent(in   ) :: bigS, alT     ! necessary constant
+
+    ! Tbr inputs first guess and outputs better estimate
+    double precision, intent(inout) :: Tbr(IM,JM)    ! brightness temp [K]
+
+    ! error code
+    integer, optional, intent(out) :: RC
+
+    ! number of iterations for wide band estimate (converges slowly)
+    integer, parameter :: Nits = 16
+
+    ! locals
+    integer :: n
+    real    :: alTwn1, alTwn2
+
+    ! error checking
+    _ASSERT(Nits >= 1,'must have at least one iteration!')
+
+    ! iterate from first guess Tbr to better estimate
+    alTwn1 = alT * wn1
+    alTwn2 = alT * wn2
+    do n = 1, Nits
+      where (Tbr > 0.0d0) &
+        Tbr = ( Fband / (bigS * (Tfunc(alTwn1/Tbr) - Tfunc(alTwn2/Tbr))) ) ** 0.25d0
+    end do
+
+    _RETURN(_SUCCESS)
+
+  end subroutine Tbr_wide_band
+
+  elemental double precision function Tfunc(x)
+
+    double precision, intent(in) :: x
+
+    ! maximum number of terms in series (converges quickly)
+    integer, parameter :: nmax = 4
+
+    ! locals
+    integer :: n, n2, n3
+    double precision :: emx, cx0, cx1, cx2, cx3, zn
+
+    ! setup
+    emx = exp(-x)
+    cx0 = 6.0d0
+    cx1 = 6.0d0 * x
+    cx2 = 3.0d0 * x**2
+    cx3 =         x**3
+
+    ! do at least 1st order
+    Tfunc = (cx3 + cx2 + cx1 + cx0) * emx
+    if (nmax <= 1) return
+
+    ! higher orders
+    zn = emx
+    do n = 2, nmax
+      n2 = n * n
+      n3 = n * n2
+      zn = zn * emx
+      Tfunc = Tfunc + (cx3 + cx2/n + cx1/n2 + cx0/n3) * zn / n
+    end do
+
+  end function Tfunc
+
+
+  ! ----------------------------------------------------------------------
+  ! Decide which radiation to use for thermodynamics state evolution.
+  ! RRTMGP dominates RRTMG dominates Chou-Suarez.
+  ! Chou-Suarez is the default if nothing else asked for in Resource file.
+  ! ----------------------------------------------------------------------
+
+  subroutine choose_solar_scheme (MAPL, &
+    USE_RRTMGP, USE_RRTMG, USE_CHOU, &
+    RC)
+
+    type (MAPL_MetaComp), pointer, intent(in) :: MAPL
+    logical, intent(out) :: USE_RRTMGP, USE_RRTMG, USE_CHOU
+    integer, optional, intent(out) :: RC  ! return code
+
+    real :: RFLAG
+    integer :: STATUS
+
+    USE_RRTMGP = .false.
+    USE_RRTMG  = .false.
+    USE_CHOU   = .false.
+    call MAPL_GetResource (MAPL, RFLAG, LABEL='USE_RRTMGP_SORAD:', DEFAULT=0., __RC__)
+    USE_RRTMGP = RFLAG /= 0.
+    if (.not. USE_RRTMGP) then
+      call MAPL_GetResource (MAPL, RFLAG, LABEL='USE_RRTMG_SORAD:', DEFAULT=0., __RC__)
+      USE_RRTMG = RFLAG /= 0.
+      USE_CHOU  = .not.USE_RRTMG
+    end if
+
+    _RETURN(_SUCCESS)
+  end subroutine choose_solar_scheme
+
+  subroutine choose_irrad_scheme (MAPL, &
+    USE_RRTMGP, USE_RRTMG, USE_CHOU, &
+    RC)
+
+    type (MAPL_MetaComp), pointer, intent(in) :: MAPL
+    logical, intent(out) :: USE_RRTMGP, USE_RRTMG, USE_CHOU
+    integer, optional, intent(out) :: RC  ! return code
+
+    real :: RFLAG
+    integer :: STATUS
+
+    USE_RRTMGP = .false.
+    USE_RRTMG  = .false.
+    USE_CHOU   = .false.
+    call MAPL_GetResource (MAPL, RFLAG, LABEL='USE_RRTMGP_IRRAD:', DEFAULT=0., __RC__)
+    USE_RRTMGP = RFLAG /= 0.
+    if (.not. USE_RRTMGP) then
+      call MAPL_GetResource (MAPL, RFLAG, LABEL='USE_RRTMG_IRRAD:', DEFAULT=0., __RC__)
+      USE_RRTMG = RFLAG /= 0.
+      USE_CHOU  = .not.USE_RRTMG
+    end if
+
+    _RETURN(_SUCCESS)
+  end subroutine choose_irrad_scheme
+
+end module rad_utils