Add constituent hooks to remaining analytic IC options.

nusbaume · nusbaume · commit 0d7919c524ed · 2024-09-10T15:08:50.000-06:00
diff --git a/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 b/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90
@@ -50,10 +50,13 @@ module ic_baro_dry_jw06
 
   subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
                                 Q, m_cnst, mask, verbose)
-    use dyn_tests_utils,  only: vc_moist_pressure, vc_dry_pressure, vc_height
-    use cam_constituents, only: const_get_index
-    !use constituents,    only: cnst_name
-    !use const_init,      only: cnst_init_default
+    use shr_kind_mod,              only: cx => shr_kind_cx
+    use dyn_tests_utils,           only: vc_moist_pressure, vc_dry_pressure, vc_height
+    use runtime_obj,               only: wv_stdname
+    use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+    use cam_ccpp_cap,              only: cam_model_const_properties
+    use ccpp_kinds,                only: kind_phys
+    use cam_constituents,          only: const_get_index, const_qmin
 
     !-----------------------------------------------------------------------
     !
@@ -79,23 +82,31 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
     logical, allocatable              :: mask_use(:)
     logical                           :: verbose_use
     logical                           :: lu,lv,lt,lq,l3d_vars
+    logical                           :: const_has_default
     integer                           :: i, k, m
     integer                           :: ncol
     integer                           :: nlev
     integer                           :: ncnst
     integer                           :: iret
-    integer                           :: ix_rain, ix_cld_liq
+    integer                           :: ix_q, m_cnst_ix_q
     character(len=*), parameter       :: subname = 'BC_DRY_JW06_SET_IC'
+    character(len=cx)                 :: errmsg !CCPP error message
     real(r8)                          :: tmp
     real(r8)                          :: r(size(latvals))
     real(r8)                          :: eta
     real(r8)                          :: factor
     real(r8)                          :: perturb_lon, perturb_lat
     real(r8)                          :: phi_vertical
     real(r8)                          :: u_wind(size(latvals))
+    real(kind_phys)                   :: const_default_value !Constituent default value
+    real(kind_phys)                   :: const_qmin_value    !Constituent minimum value
 
-       a_omega                = rearth*omega
-       exponent               = rair*gamma/gravit
+    !Private array of constituent properties (for property interface functions)
+    type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:)
+
+    !Set local constants:
+    a_omega                = rearth*omega
+    exponent               = rair*gamma/gravit
 
     allocate(mask_use(size(latvals)), stat=iret)
     call check_allocate(iret, subname, 'mask_use(size(latvals))', &
@@ -116,10 +127,6 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
       verbose_use = .true.
     end if
 
-    !set constituent indices
-    call const_get_index('cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', ix_cld_liq)
-    call const_get_index('rain_mixing_ratio_wrt_moist_air_and_condensed_water', ix_rain)
-
     ncol = size(latvals, 1)
     nlev = -1
 
@@ -236,40 +243,63 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
         end if
       end if
       if (lq) then
+         !Get water vapor constituent index:
+         call const_get_index(wv_stdname, ix_q)
+
+         !Determine which "Q" variable index matches water vapor:
+         m_cnst_ix_q = findloc(m_cnst, ix_q, dim=1)
+
         do k = 1, nlev
           where(mask_use)
-            Q(:,k,1) = 0.0_r8
+            Q(:,k,m_cnst_ix_q) = 0.0_r8
           end where
         end do
-!Un-comment once constituents are working in CAMDEN -JN:
-#if 0
         if(masterproc.and. verbose_use) then
-          write(iulog,*) '         ', trim(cnst_name(m_cnst(1))), ' initialized by "',subname,'"'
+          write(iulog,*) '         ', wv_stdname, ' initialized by "',subname,'"'
         end if
-#endif
       end if
     end if
 
-!Un-comment once constituents are working in CAMDEN -JN:
-#if 0
-    if (lq) then
-      ncnst = size(m_cnst, 1)
-      if ((vcoord == vc_moist_pressure) .or. (vcoord == vc_dry_pressure)) then
-        do m = 2, ncnst
-          call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),&
-               mask=mask_use, verbose=verbose_use, notfound=.false.)
-        end do
-      end if
-    end if
-#else
     if (lq) then
+      ncnst = size(m_cnst)
       if ((vcoord == vc_moist_pressure) .or. (vcoord == vc_dry_pressure)) then
-        !Initialize cloud liquid and rain until constituent routines are enabled:
-        Q(:,:,ix_cld_liq) = 0.0_r8
-        Q(:,:,ix_rain)    = 0.0_r8
-      end if
-    end if
-#endif
+        do m = 1, ncnst
+
+          !Skip water vapor, as it was aleady set above:
+          if (m == m_cnst_ix_q) cycle
+
+          !Extract constituent minimum value:
+          const_qmin_value = const_qmin(m_cnst(m))
+
+          !Initialize constituent to its minimum value:
+          Q(:,:,m) = real(const_qmin_value, r8)
+
+          if (iret /= 0) then
+            call endrun(errmsg, file=__FILE__, line=__LINE__)
+          end if
+
+          if (const_has_default) then
+
+            !If default value exists, then extract default value
+            !from constituent properties object:
+            call const_props(m_cnst(m))%default_value(const_default_value, &
+                                                      iret,        &
+                                                      errmsg)
+            if (iret /= 0) then
+              call endrun(errmsg, file=__FILE__, line=__LINE__)
+            end if
+
+            !Set constituent to default value in masked region:
+            do k=1,nlev
+              where(mask_use)
+                Q(:,k,m) = real(const_default_value, r8)
+              end where
+            end do
+
+          end if !has_default
+        end do   !m_cnst
+      end if     !lq
+    end if       !l3d_vars
 
     deallocate(mask_use)
 
diff --git a/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 b/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90
@@ -7,19 +7,6 @@ module ic_us_standard_atmosphere
 !
 !-------------------------------------------------------------------------------
 
-use shr_kind_mod,        only: r8 => shr_kind_r8
-use spmd_utils,          only: masterproc
-
-use hycoef,              only: ps0, hyam, hybm
-use physconst,           only: gravit
-!use constituents,        only: cnst_name
-!use const_init,          only: cnst_init_default
-
-use std_atm_profile,     only: std_atm_pres, std_atm_height, std_atm_temp
-
-use cam_logfile,         only: iulog
-use cam_abortutils,      only: endrun, check_allocate
-
 implicit none
 private
 save
@@ -33,6 +20,19 @@ module ic_us_standard_atmosphere
 subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, &
        PHIS_OUT, Q, m_cnst, mask, verbose)
 
+   use shr_kind_mod,              only: r8 => shr_kind_r8, cx => shr_kind_cx
+   use ccpp_kinds,                only: kind_phys
+   use spmd_utils,                only: masterproc
+   use hycoef,                    only: ps0, hyam, hybm
+   use physconst,                 only: gravit
+   use std_atm_profile,           only: std_atm_pres, std_atm_height, std_atm_temp
+   use cam_logfile,               only: iulog
+   use cam_abortutils,            only: endrun, check_allocate
+   use runtime_obj,               only: wv_stdname
+   use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+   use cam_ccpp_cap,              only: cam_model_const_properties
+   use cam_constituents,          only: const_get_index, const_qmin
+
    !----------------------------------------------------------------------------
    !
    ! Set initial values for static atmosphere with vertical profile from US
@@ -58,14 +58,22 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, &
    ! Local variables
    logical, allocatable              :: mask_use(:)
    logical                           :: verbose_use
+   logical                           :: const_has_default
    integer                           :: i, k, m
    integer                           :: ncol
    integer                           :: nlev, nlevp
    integer                           :: ncnst
    integer                           :: iret
+   integer                           :: ix_q, m_cnst_ix_q
    character(len=*), parameter       :: subname = 'us_std_atm_set_ic'
+   character(len=cx)                 :: errmsg !CCPP error message
    real(r8)                          :: psurf(1)
    real(r8), allocatable             :: pmid(:), zmid(:), zmid2d(:,:)
+   real(kind_phys)                   :: const_default_value !Constituent default value
+   real(kind_phys)                   :: const_qmin_value    !Constituent minimum value
+
+   !Private array of constituent properties (for property interface functions)
+   type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:)
    !----------------------------------------------------------------------------
 
    ! check input consistency
@@ -208,35 +216,58 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, &
          zmid2d = 0.5_r8*(zint(:,1:nlev) + zint(:,2:nlev+1))
       end if
 
-      ncnst = size(m_cnst, 1)
+      !Get water vapor constituent index:
+      call const_get_index(wv_stdname, ix_q)
+
+      !Determine which "Q" variable index matches water vapor:
+      m_cnst_ix_q = findloc(m_cnst, ix_q, dim=1)
+
+      ncnst = size(m_cnst)
       do m = 1, ncnst
-         if (m_cnst(m) == 1) then
+         if (m_cnst(m) == m_cnst_ix_q) then
             ! No water vapor in profile
             do k = 1, nlev
                where(mask_use)
                   Q(:,k,m_cnst(m)) = 0.0_r8
                end where
             end do
-!Un-comment once constituents are working in CAMDEN -JN:
-#if 0
             if(masterproc .and. verbose_use) then
-               write(iulog,*) '          ', trim(cnst_name(m_cnst(m))), ' initialized by '//subname
-            end if
-         else
-            if (present(zint)) then
-               call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),&
-                  mask=mask_use, verbose=verbose_use, notfound=.false., z=zmid2d)
-            else
-               call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),&
-                  mask=mask_use, verbose=verbose_use, notfound=.false.)
+               write(iulog,*) '          ', wv_stdname, ' initialized by '//subname
             end if
-#else
          else
-           !Initialize cloud liquid and rain until constituent routines are enabled:
-           Q(:,:,m_cnst(m)) = 0.0_r8
-#endif
-         end if
-      end do
+
+           !Extract constituent minimum value:
+           const_qmin_value = const_qmin(m_cnst(m))
+
+           !Initialize constituent to its minimum value:
+           Q(:,:,m) = real(const_qmin_value, r8)
+
+           if (iret /= 0) then
+              call endrun(errmsg, file=__FILE__, line=__LINE__)
+           end if
+
+           if (const_has_default) then
+
+              !If default value exists, then extract default value
+              !from constituent properties object:
+              call const_props(m_cnst(m))%default_value(const_default_value, &
+                                                        iret,        &
+                                                        errmsg)
+              if (iret /= 0) then
+                call endrun(errmsg, file=__FILE__, line=__LINE__)
+              end if
+
+              !Set constituent to default value in masked region:
+              do k=1,nlev
+                 where(mask_use)
+                    Q(:,k,m) = real(const_default_value, r8)
+                 end where
+              end do
+
+           end if !has_default
+
+         end if !water vapor
+      end do    !ncnst
 
       if (allocated(zmid2d)) deallocate(zmid2d)
 
diff --git a/src/dynamics/tests/namelist_definition_analy_ic.xml b/src/dynamics/tests/namelist_definition_analy_ic.xml
@@ -8,7 +8,9 @@
     <type>char*80</type>
     <category>dyn_test</category>
     <group>analytic_ic_nl</group>
-    <valid_values>none,held_suarez_1994,moist_baroclinic_wave_dcmip2016,dry_baroclinic_wave_dcmip2016,dry_baroclinic_wave_jw2006</valid_values>
+    <valid_values>
+      none,held_suarez_1994,moist_baroclinic_wave_dcmip2016,dry_baroclinic_wave_dcmip2016,dry_baroclinic_wave_jw2006,us_standard_atmosphere
+    </valid_values>
     <desc>
       Specify the type of analytic initial conditions for an initial run.
       held_suarez_1994: Initial conditions specified in Held and Suarez (1994)
