update to cam6_4_009

.gitmodules

@@ -36,15 +36,15 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/ESCOMP/atmospheric_physics
-        fxtag = atmos_phys0_02_006
+        fxtag = atmos_phys0_03_000
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
 
 [submodule "fv3"]
 	path = src/dynamics/fv3
 	url = https://github.com/ESCOMP/CAM_FV3_interface.git
         fxrequired = AlwaysRequired
-        fxtag = fv3int_053124
+        fxtag = fv3int_061924
 	fxDONOTUSEurl = https://github.com/ESCOMP/CAM_FV3_interface.git
 
 [submodule "geoschem"]

bld/configure

@@ -2185,6 +2185,7 @@ sub write_filepath
       print $fh "$camsrcdir/src/physics/simple\n";
       print $fh "$camsrcdir/src/atmos_phys/kessler\n";
       print $fh "$camsrcdir/src/atmos_phys/held_suarez\n";
+      print $fh "$camsrcdir/src/atmos_phys/tj2016\n";
     }
 
     # Weak scaling fix. This has to come before physics/cam and before dycores

src/dynamics/fv/dp_coupling.F90

@@ -576,7 +576,7 @@ subroutine d_p_coupling(grid, phys_state, phys_tend,  pbuf2d, dyn_out)
        ! (note: cam_thermo_dry_air_update assumes dry unless optional conversion factor provided)
        !
        call set_state_pdry(phys_state(lchnk))    ! First get dry pressure to use for this timestep
-       call set_wet_to_dry(phys_state(lchnk))    ! Dynamics had moist, physics wants dry
+       call set_wet_to_dry(phys_state(lchnk), convert_cnst_type='dry') ! Dynamics had moist, physics wants dry
        if (dry_air_species_num>0) then
 !------------------------------------------------------------
 ! Apply limiters to mixing ratios of major species

src/dynamics/se/dycore/global_norms_mod.F90

@@ -639,7 +639,11 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,&
     !
     ! if user or namelist is not specifying sponge del4 settings here are best guesses (empirically determined)
     !
-    if (top_090_140km.or.top_140_600km) then ! defaults for waccm(x)
+    if (top_042_090km) then
+       if (sponge_del4_lev       <0) sponge_del4_lev        = 4
+       if (sponge_del4_nu_fac    <0) sponge_del4_nu_fac     = 3.375_r8 !max value without having to increase subcycling of div4
+       if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 3.375_r8 !max value without having to increase subcycling of div4
+    else if (top_090_140km.or.top_140_600km) then ! defaults for waccm(x)
       if (sponge_del4_lev       <0) sponge_del4_lev        = 20
       if (sponge_del4_nu_fac    <0) sponge_del4_nu_fac     = 5.0_r8
       if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 10.0_r8

src/dynamics/se/dycore/prim_advance_mod.F90

@@ -444,7 +444,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2,
     use dimensions_mod, only: nu_scale_top,nu_lev,kmvis_ref,kmcnd_ref,rho_ref,km_sponge_factor
     use dimensions_mod, only: nu_t_lev
     use control_mod,    only: nu, nu_t, hypervis_subcycle,hypervis_subcycle_sponge, nu_p, nu_top
-    use control_mod,    only: molecular_diff
+    use control_mod,    only: molecular_diff,sponge_del4_lev
     use hybrid_mod,     only: hybrid_t!, get_loop_ranges
     use element_mod,    only: element_t
     use derivative_mod, only: derivative_t, laplace_sphere_wk, vlaplace_sphere_wk, vlaplace_sphere_wk_mol
@@ -668,7 +668,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2,
       call tot_energy_dyn(elem,fvm,nets,nete,nt,qn0,'dCH')
       do ie=nets,nete
         !$omp parallel do num_threads(vert_num_threads), private(k,i,j,v1,v2,heating)
-        do k=ksponge_end,nlev
+        do k=sponge_del4_lev+2,nlev
           !
           ! only do "frictional heating" away from sponge
           !

src/ionosphere/waccmx/wei05sc.F90

@@ -1126,7 +1126,7 @@ real(r8) function km_n(m,rn)
         return
      end if
      rm = real(m, r8)
-     km_n = sqrt(2._r8*exp(lngamma(rn+rm+1._r8)-lngamma(rn-rm+1._r8))) / &
+     km_n = sqrt(2._r8*exp(log_gamma(rn+rm+1._r8)-log_gamma(rn-rm+1._r8))) / &
           (2._r8**m*factorial(m))
   end function km_n
 !-----------------------------------------------------------------------
@@ -1296,32 +1296,6 @@ integer function value_locate(vec,val)
         end if
      end do
   end function value_locate
-!-----------------------------------------------------------------------
-  real(r8) function lngamma(xx)
-     !
-     ! This is an f90 translation from C code copied from
-     ! gammln routine from "Numerical Recipes in C" Chapter 6.1.
-     ! see: http://numerical.recipes
-     !
-
-     real(r8), intent(in) :: xx
-     real(r8) :: x,y,tmp,ser
-     real(r8) :: cof(6) = (/76.18009172947146_r8, -86.50532032941677_r8,      &
-          24.01409824083091_r8, -1.231739572450155_r8,                        &
-          0.1208650973866179e-2_r8, -0.5395239384953e-5_r8/)
-     integer :: j
-     !
-     y = xx
-     x = xx
-     tmp = x+5.5_r8
-     tmp = tmp-(x + 0.5_r8) * log(tmp)
-     ser = 1.000000000190015_r8
-     do j = 1, 5
-        y = y + 1
-        ser = ser + (cof(j) / y)
-     end do
-     lngamma = -tmp+log(2.5066282746310005_r8*ser/x)
-  end function lngamma
 !-----------------------------------------------------------------------
   real(r8) function factorial(n)
      integer,intent(in) :: n

src/physics/cam/clubb_intr.F90

@@ -2549,8 +2549,9 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
     ! Copy the state to state1 array to use in this routine
     call physics_state_copy(state, state1)
 
-    ! constituents are all treated as dry mmr by clubb
-    call set_wet_to_dry(state1)
+    ! Constituents are all treated as dry mmr by clubb.  Convert the water species to
+    ! a dry basis.
+    call set_wet_to_dry(state1, convert_cnst_type='wet')
 
     if (clubb_do_liqsupersat) then
       call pbuf_get_field(pbuf, npccn_idx, npccn)

src/physics/cam/gw_drag.F90

@@ -1396,7 +1396,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
   call physics_state_copy(state, state1)
 
   ! constituents are all treated as wet mmr
-  call set_dry_to_wet(state1)
+  call set_dry_to_wet(state1, convert_cnst_type='dry')
 
   lchnk = state1%lchnk
   ncol  = state1%ncol

src/physics/cam/hb_diff.F90

@@ -262,7 +262,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     !-----------------------------------------------------------------------
 
-    use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm
+    use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm_free
 
     !------------------------------Arguments--------------------------------
     !
@@ -321,7 +321,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     ! Get free atmosphere exchange coefficients
     !
-    call austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
+    call austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
          ml2, ri, s2, kvf)
 
     kvq(:ncol,:) = kvf(:ncol,:)

src/physics/cam/pbl_utils.F90

@@ -27,7 +27,7 @@ module pbl_utils
 public calc_obklen
 public virtem
 public compute_radf
-public austausch_atm
+public austausch_atm, austausch_atm_free
 
 real(r8), parameter :: ustar_min = 0.01_r8
 
@@ -408,4 +408,62 @@ subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
 
 end subroutine austausch_atm
 
+subroutine austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
+
+  !---------------------------------------------------------------------- !
+  !                                                                       !
+  ! same as austausch_atm but only mixing for Ri<0                        !
+  ! i.e. no background mixing and mixing for Ri>0                         !
+  !                                                                       !
+  !---------------------------------------------------------------------- !
+
+  ! --------------- !
+  ! Input arguments !
+  ! --------------- !
+
+  integer,  intent(in)  :: pcols                ! Atmospheric columns dimension size
+  integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
+  integer,  intent(in)  :: pver                 ! Number of atmospheric layers
+  integer,  intent(in)  :: ntop                 ! Top layer for calculation
+  integer,  intent(in)  :: nbot                 ! Bottom layer for calculation
+
+  real(r8), intent(in)  :: ml2(pver+1)          ! Mixing lengths squared
+  real(r8), intent(in)  :: s2(pcols,pver)       ! Shear squared
+  real(r8), intent(in)  :: ri(pcols,pver)       ! Richardson no
+
+  ! ---------------- !
+  ! Output arguments !
+  ! ---------------- !
+
+  real(r8), intent(out) :: kvf(pcols,pver+1)    ! Eddy diffusivity for heat and tracers
+
+  ! --------------- !
+  ! Local Variables !
+  ! --------------- !
+
+  real(r8)              :: fofri                ! f(ri)
+  real(r8)              :: kvn                  ! Neutral Kv
+
+  integer               :: i                    ! Longitude index
+  integer               :: k                    ! Vertical index
+
+  ! ----------------------- !
+  ! Main Computation Begins !
+  ! ----------------------- !
+
+  kvf(:ncol,:)           = 0.0_r8
+  ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
+  do k = ntop, nbot - 1
+     do i = 1, ncol
+        if( ri(i,k) < 0.0_r8 ) then
+           fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
+        else
+           fofri = 0.0_r8
+        end if
+        kvn = ml2(k) * sqrt(s2(i,k))
+        kvf(i,k+1) = kvn * fofri
+     end do
+  end do
+end subroutine austausch_atm_free
+
 end module pbl_utils

src/physics/cam/physics_types.F90

@@ -1481,40 +1481,72 @@ end subroutine set_state_pdry
 
 !===============================================================================
 
-subroutine set_wet_to_dry (state)
+subroutine set_wet_to_dry(state, convert_cnst_type)
+
+  ! Convert mixing ratios from a wet to dry basis for constituents of type
+  ! convert_cnst_type.  Constituents are given a type when they are added
+  ! to the constituent array by a call to cnst_add during the register
+  ! phase of initialization.  There are two constituent types: 'wet' for
+  ! water species and 'dry' for non-water species.
 
   use constituents,  only: pcnst, cnst_type
 
   type(physics_state), intent(inout) :: state
+  character(len=3),    intent(in)    :: convert_cnst_type
 
+  ! local variables
   integer m, ncol
+  character(len=*), parameter :: sub = 'set_wet_to_dry'
+  !-----------------------------------------------------------------------------
+
+  ! check input
+  if (.not.(convert_cnst_type == 'wet' .or. convert_cnst_type == 'dry')) then
+    write(iulog,*) sub//': FATAL: convert_cnst_type not recognized: '//convert_cnst_type
+    call endrun(sub//': FATAL: convert_cnst_type not recognized: '//convert_cnst_type)
+  end if
 
   ncol = state%ncol
 
-  do m = 1,pcnst
-     if (cnst_type(m).eq.'dry') then
+  do m = 1, pcnst
+     if (cnst_type(m) == convert_cnst_type) then
         state%q(:ncol,:,m) = state%q(:ncol,:,m)*state%pdel(:ncol,:)/state%pdeldry(:ncol,:)
-     endif
+     end if
   end do
 
 end subroutine set_wet_to_dry
 
 !===============================================================================
 
-subroutine set_dry_to_wet (state)
+subroutine set_dry_to_wet(state, convert_cnst_type)
+
+  ! Convert mixing ratios from a dry to wet basis for constituents of type
+  ! convert_cnst_type.  Constituents are given a type when they are added
+  ! to the constituent array by a call to cnst_add during the register
+  ! phase of initialization.  There are two constituent types: 'wet' for
+  ! water species and 'dry' for non-water species.
 
   use constituents,  only: pcnst, cnst_type
 
   type(physics_state), intent(inout) :: state
+  character(len=3),    intent(in)    :: convert_cnst_type
 
+  ! local variables
   integer m, ncol
+  character(len=*), parameter :: sub = 'set_dry_to_wet'
+  !-----------------------------------------------------------------------------
+
+  ! check input
+  if (.not.(convert_cnst_type == 'wet' .or. convert_cnst_type == 'dry')) then
+    write(iulog,*) sub//': FATAL: convert_cnst_type not recognized: '//convert_cnst_type
+    call endrun(sub//': FATAL: convert_cnst_type not recognized: '//convert_cnst_type)
+  end if
 
   ncol = state%ncol
 
-  do m = 1,pcnst
-     if (cnst_type(m).eq.'dry') then
+  do m = 1, pcnst
+     if (cnst_type(m) == convert_cnst_type) then
         state%q(:ncol,:,m) = state%q(:ncol,:,m)*state%pdeldry(:ncol,:)/state%pdel(:ncol,:)
-     endif
+     end if
   end do
 
 end subroutine set_dry_to_wet

src/physics/cam/physpkg.F90

@@ -1948,7 +1948,7 @@ subroutine tphysac (ztodt,   cam_in,  &
       !
       ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call
       !
-      call set_dry_to_wet(state)
+      call set_dry_to_wet(state, convert_cnst_type='dry')
 
       if (trim(cam_take_snapshot_before) == "physics_dme_adjust") then
         call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&

src/physics/cam/vertical_diffusion.F90

@@ -911,7 +911,7 @@ subroutine vertical_diffusion_tend( &
   ! ----------------------- !
 
   ! Assume 'wet' mixing ratios in diffusion code.
-  call set_dry_to_wet(state)
+  call set_dry_to_wet(state, convert_cnst_type='dry')
 
   rztodt = 1._r8 / ztodt
   lchnk  = state%lchnk
@@ -1384,7 +1384,7 @@ subroutine vertical_diffusion_tend( &
      endif
   end do
   ! convert wet mmr back to dry before conservation check
-  call set_wet_to_dry(state)
+  call set_wet_to_dry(state, convert_cnst_type='dry')
 
   if (.not. do_pbl_diags) then
      slten(:ncol,:)         = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt

src/physics/cam7/physpkg.F90

@@ -2395,7 +2395,7 @@ subroutine tphysac (ztodt,   cam_in,  &
       !
       ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call
       !
-      call set_dry_to_wet(state)
+      call set_dry_to_wet(state, convert_cnst_type='dry')
 
       if (trim(cam_take_snapshot_before) == "physics_dme_adjust") then
          call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&

src/physics/carma/cam/carma_intr.F90

@@ -1036,7 +1036,7 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli
     ! The CARMA interface assumes that mass mixing ratios are relative to a
     ! wet atmosphere, so convert any dry mass mixing ratios to wet.
     call physics_state_copy(state, state_loc)
-    call set_dry_to_wet(state_loc)
+    call set_dry_to_wet(state_loc, convert_cnst_type='dry')
 
     spdiags(:, :, :)       = 0.0_r8
     gpdiags(:, :, :, :)    = 0.0_r8

src/physics/carma/models/cirrus/carma_cloudfraction.F90

@@ -24,8 +24,7 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri
     use carma_mod
 
     use shr_kind_mod,   only: r8 => shr_kind_r8
-    use physics_types,  only: physics_state, physics_ptend, set_wet_to_dry, &
-                              set_dry_to_wet
+    use physics_types,  only: physics_state, physics_ptend
     use constituents,   only: cnst_get_ind
     use cam_abortutils, only: endrun
 

src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90

@@ -24,8 +24,7 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri
     use carma_mod
 
     use shr_kind_mod,   only: r8 => shr_kind_r8
-    use physics_types,  only: physics_state, physics_ptend, set_wet_to_dry, &
-                              set_dry_to_wet
+    use physics_types,  only: physics_state, physics_ptend
     use constituents,   only: cnst_get_ind
     use cam_abortutils, only: endrun
 

src/physics/simple/physpkg.F90

@@ -654,7 +654,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf)
         !
         ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call
         !
-        call set_dry_to_wet(state)
+        call set_dry_to_wet(state, convert_cnst_type='dry')
         call physics_dme_adjust(state, tend, qini, totliqini, toticeini, ztodt)
         call tot_energy_phys(state, 'phAM')
         call tot_energy_phys(state, 'dyAM', vc=vc_dycore)