Cam sima dycore update

cime_config/namelist_definition_cam.xml

@@ -189,7 +189,27 @@
       <value>UNSET_PATH</value>
     </values>
   </entry>
-
+  <entry id="scale_dry_air_mass">
+    <type>real</type>
+    <category>initial_conditions</category>
+    <group>cam_initfiles_nl</group>
+    <desc>
+      Specify whether and how to perform
+      dry surface pressure scaling. If less than or equal to 0.0,
+      do not perform scaling. If greater than 0.0, perform scaling to scale_dry_air_mass
+      value (in Pa) as the average dry surface pressure target.
+      Default: set by build-namelist.
+    </desc>
+    <values>
+      <value>0.0D0</value>
+      <value aquaplanet="1"> 101080.0D0 </value>
+      <!-- Scale Dry Air Mass: for cases with topography -->
+      <value phys_suite="cam4"   >  98288.0D0 </value>
+      <value phys_suite="cam5"   >  98288.0D0 </value>
+      <value phys_suite="cam6"   >  98288.0D0 </value>
+      <value phys_suite="cam7"   >  98288.0D0 </value>
+    </values>
+  </entry>
 <!-- Topography -->
 
   <entry id="bnd_topo">

src/control/cam_initfiles.F90

@@ -44,6 +44,8 @@ module cam_initfiles
    ! cam_branch_file: Filepath of primary restart file for a branch run
    character(len=cl) :: cam_branch_file = ' '
 
+   real(r8), public, protected :: scale_dry_air_mass = 0.0_r8 ! Toggle and target avg air mass
+
    ! rest_pfile: The restart pointer file contains name of most recently
    !             written primary restart file.
    !             The contents of this file are updated by cam_write_restart
@@ -89,7 +91,7 @@ subroutine cam_initfiles_readnl(nlfile)
       character(len=*), parameter :: subname = 'cam_initfiles_readnl'
 
       namelist /cam_initfiles_nl/ ncdata, bnd_topo, pertlim, cam_branch_file, &
-         unset_path_str
+         unset_path_str, scale_dry_air_mass
       !------------------------------------------------------------------------
 
       if (masterproc) then
@@ -121,7 +123,11 @@ subroutine cam_initfiles_readnl(nlfile)
            mstrid, mpicom, ierr)
       if (ierr /= 0) then
          call endrun(subname//": ERROR: mpi_bcast: cam_branch_file")
-      end if
+       end if
+       call mpi_bcast(scale_dry_air_mass, 1, mpi_real8, mstrid, mpicom, ierr)
+       if (ierr /= 0) then
+         call endrun(subname//": ERROR: mpi_bcast: scale_dry_air_mass")
+       endif
       call mpi_bcast(unset_path_str, len(unset_path_str), mpi_character,      &
            mstrid, mpicom, ierr)
       if (ierr /= 0) then
@@ -198,7 +204,13 @@ subroutine cam_initfiles_readnl(nlfile)
 
          write(iulog,*) '  Maximum abs value of scale factor used to ',       &
               'perturb initial conditions, pertlim= ', pertlim
-
+         if (scale_dry_air_mass > 0) then
+           write(iulog,*) &
+                '  Initial condition dry mass will be scaled to: ',scale_dry_air_mass,' Pa'
+         else
+           write(iulog,*) &
+              '  Initial condition dry mass will not be scaled.'
+         end if
       end if
 
    end subroutine cam_initfiles_readnl

src/dynamics/se/advect_tend.F90

@@ -29,7 +29,7 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
     use cam_abortutils,         only: check_allocate
 
     ! SE dycore:
-    use dimensions_mod,         only: nc,np,nlev,ntrac
+    use dimensions_mod,         only: nc,np,nlev,ntrac,use_cslam
     use element_mod,            only: element_t
     use fvm_control_volume_mod, only: fvm_struct
 
@@ -45,7 +45,7 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
 
     character(len=*), parameter :: subname = 'compute_adv_tends_xyz'
 
-    if (ntrac>0) then
+    if (use_cslam) then
       nx=nc
     else
       nx=np
@@ -65,7 +65,7 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
       adv_tendxyz(:,:,:,:,:) = 0._r8
     endif
 
-    if (ntrac>0) then
+    if (use_cslam) then
       do ie=nets,nete
         do ic = 1, num_advected
           adv_tendxyz(:,:,:,ic,ie) = fvm(ie)%c(1:nc,1:nc,:,ic) - adv_tendxyz(:,:,:,ic,ie)
@@ -105,5 +105,174 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
     deallocate(ftmp)
 #endif
   end subroutine compute_adv_tends_xyz
+#ifdef SCAM
+  !----------------------------------------------------------------------
+  ! computes camiop specific tendencies
+  ! and writes these to the camiop file
+  ! called twice each time step:
+  !   - first call sets the initial mixing ratios/state
+  !   - second call computes and outputs the tendencies
+  !----------------------------------------------------------------------
+  subroutine compute_write_iop_fields(elem,fvm,nets,nete,qn0,n0)
+    use cam_abortutils,         only: endrun
+    use cam_history,            only: outfld, hist_fld_active
+    use time_manager,           only: get_step_size
+    use constituents,           only: pcnst,cnst_name
+    use dimensions_mod,         only: nc,np,nlev,use_cslam,npsq
+    use element_mod,            only: element_t
+    use fvm_control_volume_mod, only: fvm_struct
+    implicit none
+
+    type (element_t), intent(inout) :: elem(:)
+    type(fvm_struct), intent(inout) :: fvm(:)
+    integer,          intent(in) :: nets,nete,qn0,n0
+    real(r8) :: dt
+    real(r8), allocatable        :: q_new(:,:,:)
+    real(r8), allocatable        :: q_adv(:,:,:)
+    real(r8), allocatable        :: t_adv(:,:)
+    real(r8), allocatable        :: out_q(:,:)
+    real(r8), allocatable        :: out_t(:,:)
+    real(r8), allocatable        :: out_u(:,:)
+    real(r8), allocatable        :: out_v(:,:)
+    real(r8), allocatable        :: out_ps(:)
+
+    integer  :: i,j,ic,nx,ie,nxsq,p
+    integer  :: ierr
+    logical  :: init
+    character(len=*), parameter :: sub = 'compute_write_iop_fields:'
+    !----------------------------------------------------------------------------
+
+    if (use_cslam) then
+      nx=nc
+    else
+      nx=np
+    endif
+    nxsq=nx*nx
 
+    init = .false.
+    dt = get_step_size()
+
+    if ( .not. allocated( iop_qtendxyz ) ) then
+      init = .true.
+
+      allocate( iop_qtendxyz(nx,nx,nlev,pcnst,nets:nete),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate iop_qtendxyz' )
+      iop_qtendxyz = 0._r8
+      allocate( derivedfq(nx,nx,nlev,pcnst,nets:nete),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate derivedfq' )
+      derivedfq = 0._r8
+      allocate( iop_qtendxyz_init(nx,nx,nlev,pcnst,nets:nete),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate iop_qtendxyz' )
+      iop_qtendxyz_init = 0._r8
+      allocate( iop_ttendxyz(nx,nx,nlev,nets:nete),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate iop_ttendxyz' )
+      iop_ttendxyz = 0._r8
+      allocate( iop_ttendxyz_init(nx,nx,nlev,nets:nete),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate iop_ttendxyz_init' )
+      iop_ttendxyz_init = 0._r8
+    endif
+
+    ! save initial/calc tendencies on second call to this routine.
+    if (use_cslam) then
+      do ie=nets,nete
+        do ic=1,pcnst
+          iop_qtendxyz(:,:,:,ic,ie) = fvm(ie)%c(1:nc,1:nc,:,ic) - iop_qtendxyz(:,:,:,ic,ie)
+        end do
+      end do
+    else
+      do ie=nets,nete
+        do ic=1,pcnst
+          iop_qtendxyz(:,:,:,ic,ie) = elem(ie)%state%Qdp(:,:,:,ic,qn0)/elem(ie)%state%dp3d(:,:,:,n0)  - iop_qtendxyz(:,:,:,ic,ie)
+       enddo
+      end do
+    end if
+    do ie=nets,nete
+       iop_ttendxyz(:,:,:,ie) = elem(ie)%state%T(:,:,:,n0)  - iop_ttendxyz(:,:,:,ie)
+    end do
+
+    if (init) then
+       do ie=nets,nete
+          iop_ttendxyz_init(:,:,:,ie)   = iop_ttendxyz(:,:,:,ie)
+          iop_qtendxyz_init(:,:,:,:,ie) = iop_qtendxyz(:,:,:,:,ie)
+          derivedfq(:,:,:,:,ie)=elem(ie)%derived%FQ(:,:,:,:)/dt
+       end do
+    end if
+
+    if ( .not. init ) then
+      allocate( q_adv(nxsq,nlev,pcnst),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate q_adv' )
+      q_adv = 0._r8
+      allocate( t_adv(npsq,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate t_adv' )
+      t_adv = 0._r8
+      allocate( q_new(nx,nx,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate q_new' )
+      q_new = 0._r8
+      allocate( out_q(npsq,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate out_q' )
+      out_q = 0._r8
+      allocate( out_t(npsq,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate out_t' )
+      out_t = 0._r8
+      allocate( out_u(npsq,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate out_u' )
+      out_u = 0._r8
+      allocate( out_v(npsq,nlev),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate out_v' )
+      out_v = 0._r8
+      allocate( out_ps(npsq),stat=ierr )
+      if (ierr/=0) call endrun( sub//': not able to allocate out_ps' )
+      out_ps = 0._r8
+      do ie=nets,nete
+         do j=1,nx
+            do i=1,nx
+               t_adv(i+(j-1)*np,:) = iop_ttendxyz(i,j,:,ie)/dt - elem(ie)%derived%FT(i,j,:)
+               out_u(i+(j-1)*np,:) = elem(ie)%state%v(i,j,1,:,n0)
+               out_v(i+(j-1)*np,:) = elem(ie)%state%v(i,j,2,:,n0)
+               out_ps(i+(j-1)*np) = elem(ie)%state%psdry(i,j)
+
+               ! to retain bfb, replace state q and t with roundoff version calculated using the ordering and tendencies of the
+               ! scam prognostic equation
+               elem(ie)%state%T(i,j,:,n0) =  iop_ttendxyz_init(i,j,:,ie) + dt*(elem(ie)%derived%FT(i,j,:) + t_adv(i+(j-1)*np,:))
+               out_t(i+(j-1)*np,:) = elem(ie)%state%T(i,j,:,n0)
+               do p=1,pcnst
+                  q_adv(i+(j-1)*nx,:,p) = iop_qtendxyz(i,j,:,p,ie)/dt - derivedfq(i,j,:,p,ie)
+                  q_new(i,j,:) = iop_qtendxyz_init(i,j,:,p,ie) + dt*(derivedfq(i,j,:,p,ie) + q_adv(i+(j-1)*nx,:,p))
+                  if (use_cslam) then
+                     fvm(ie)%c(i,j,:,p)=q_new(i,j,:)
+                  else
+                     elem(ie)%state%Qdp(i,j,:,p,qn0)=q_new(i,j,:)*elem(ie)%state%dp3d(i,j,:,n0)
+                  end if
+               enddo
+               out_q(i+(j-1)*nx,:) = elem(ie)%state%Qdp(i,j,:,1,qn0)/elem(ie)%state%dp3d(i,j,:,n0)
+            end do
+         end do
+         call outfld('Ps',out_ps,npsq,ie)
+         call outfld('t',out_t,npsq,ie)
+         call outfld('q',out_q,nxsq,ie)
+         call outfld('u',out_u,npsq,ie)
+         call outfld('v',out_v,npsq,ie)
+         call outfld('divT3d',t_adv,npsq,ie)
+         do p=1,pcnst
+            call outfld(trim(cnst_name(p))//'_dten',q_adv(:,:,p),nxsq,ie)
+         enddo
+      end do
+
+      deallocate(iop_ttendxyz)
+      deallocate(iop_ttendxyz_init)
+      deallocate(iop_qtendxyz)
+      deallocate(iop_qtendxyz_init)
+      deallocate(derivedfq)
+      deallocate(out_t)
+      deallocate(out_q)
+      deallocate(out_u)
+      deallocate(out_v)
+      deallocate(out_ps)
+      deallocate(t_adv)
+      deallocate(q_adv)
+      deallocate(q_new)
+
+    endif
+  end subroutine compute_write_iop_fields
+#endif
 end module advect_tend

src/dynamics/se/dp_coupling.F90

@@ -61,7 +61,7 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out)
 
    !SE dycore:
    use fvm_mapping,            only: dyn2phys_vector, dyn2phys_all_vars
-   use time_mod,               only: timelevel_qdp
+   use se_dyn_time_mod,        only: timelevel_qdp
    use control_mod,            only: qsplit
 
    ! arguments
@@ -329,7 +329,7 @@ subroutine p_d_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_in, tl_f, t
    use bndry_mod,        only: bndry_exchange
    use edge_mod,         only: edgeVpack, edgeVunpack
    use fvm_mapping,      only: phys2dyn_forcings_fvm
-
+   use dimensions_mod,   only: use_cslam
    ! arguments
    type(runtime_options), intent(in)     :: cam_runtime_opts ! Runtime settings object
    type(physics_state),   intent(inout)  :: phys_state
@@ -521,9 +521,15 @@ subroutine p_d_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_in, tl_f, t
       call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
       kptr = kptr + 2*nlev
       call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
-      kptr = kptr + nlev
-      call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
-   end do
+      if (.not. use_cslam) then
+        !
+        ! if using CSLAM qdp is being overwritten with CSLAM values in the dynamics
+        ! so no need to do boundary exchange of tracer tendency on GLL grid here
+        !
+        kptr = kptr + nlev
+        call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
+      end if
+    end do
 
    if (iam < par%nprocs) then
      call bndry_exchange(par, edgebuf, location='p_d_coupling')
@@ -534,8 +540,10 @@ subroutine p_d_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_in, tl_f, t
       call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
       kptr = kptr + 2*nlev
       call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
-      kptr = kptr + nlev
-      call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
+      if (.not. use_cslam) then
+        kptr = kptr + nlev
+        call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
+      end if
       if (fv_nphys > 0) then
          do k = 1, nlev
             dyn_in%elem(ie)%derived%FM(:,:,1,k) =                             &
@@ -743,6 +751,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
    ! Ensure N2 = 1 - (O2 + O + H) mmr is greater than 0
    ! Check for unusually large H2 values and set to lower value.
    !------------------------------------------------------------
+   !xxx this code is NOT in cam_development?
    if (cam_runtime_opts%waccmx_option() == 'ionosphere' .or. &
        cam_runtime_opts%waccmx_option() == 'neutral')  then
 
@@ -853,44 +862,4 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
                                      errflg, errmsg)
 
 end subroutine derived_phys_dry
-
-!=========================================================================================
-
-subroutine thermodynamic_consistency(phys_state, const_data_ptr, phys_tend, ncols, pver)
-  !
-   ! Adjust the physics temperature tendency for thermal energy consistency with the
-   ! dynamics.
-   ! Note: mixing ratios are assumed to be dry.
-   !
-   use physconst,         only: cpair
-   use air_composition,   only: get_cp
-
-   ! SE dycore:
-   use dimensions_mod,    only: lcp_moist
-   use control_mod,       only: phys_dyn_cp
-
-   type(physics_state), intent(in)    :: phys_state
-   real(kind_phys), pointer           :: const_data_ptr(:,:,:)
-   type(physics_tend ), intent(inout) :: phys_tend
-   integer,  intent(in)               :: ncols, pver
-
-   real(kind_phys) :: inv_cp(ncols,pver)
-   !----------------------------------------------------------------------------
-
-   if (lcp_moist.and.phys_dyn_cp==1) then
-     !
-     ! scale temperature tendency so that thermal energy increment from physics
-     ! matches SE (not taking into account dme adjust)
-     !
-     ! note that if lcp_moist=.false. then there is thermal energy increment
-     ! consistency (not taking into account dme adjust)
-     !
-     call get_cp(const_data_ptr(1:ncols,1:pver,1:num_advected),.true.,inv_cp)
-
-     phys_tend%dTdt_total(1:ncols,1:pver) = phys_tend%dTdt_total(1:ncols,1:pver)*cpair*inv_cp
-   end if
-end subroutine thermodynamic_consistency
-
-!=========================================================================================
-
 end module dp_coupling