Cam sima dycore update

cime_config/namelist_definition_cam.xml

@@ -65,6 +65,8 @@
       <value dyn="se" hgrid="ne5np4" nlev="66">${DIN_LOC_ROOT}/atm/waccm/ic/wa3_ne5np4_1950_spinup.cam2.i.1960-01-01-00000_c150810.nc</value>
       <value dyn="se" hgrid="ne30np4" nlev="70" sim_year="1850">${DIN_LOC_ROOT}/atm/waccm/ic/waccm5_1850_ne30np4_L70_0001-01-11-00000_c151217.nc</value>
       <value dyn="se" hgrid="ne30np4" nlev="70">${DIN_LOC_ROOT}/atm/waccm/ic/fw2000_ne30np4_L70_c181221.nc</value>
+      <value dyn="se" hgrid="ne16np4" nlev="58">${DIN_LOC_ROOT}/atm/cam/inic/se/FHISTC_LTso_ne16pg3_ne16pg3_mg17.i.c260220.nc</value>
+      <value dyn="se" hgrid="ne16np4" nlev="93">${DIN_LOC_ROOT}/atm/cam/inic/se/FHISTC_MTso_ne16pg3_ne16pg3_mg17.i.c260220.nc</value>
 
       <!-- MPAS dycore ICs -->
       <value dyn="mpas" hgrid="mpasa120"  nlev="26" analytic_ic="1">${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L26_notopo_coords_c251105.nc</value>
@@ -116,7 +118,28 @@
       <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: 0.0
+    </desc>
+    <values>
+      <value>0.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>
+      <!-- Scale Dry Air Mass: for cases with no topography (aquaplanet) -->
+      <value aquaplanet="1"> 101080.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
@@ -122,6 +124,10 @@ subroutine cam_initfiles_readnl(nlfile)
       if (ierr /= 0) then
          call endrun(subname//": ERROR: mpi_bcast: cam_branch_file")
       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/control/runtime_opts.F90

@@ -25,7 +25,6 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon)
       use cam_abortutils,               only: endrun
       use cam_logfile,                  only: cam_logfile_readnl, iulog
       use cam_initfiles,                only: cam_initfiles_readnl
-      use cam_constituents,             only: cam_constituents_readnl
       use cam_ccpp_scheme_namelists,    only: cam_read_ccpp_scheme_namelists
       use runtime_obj,                  only: cam_set_runtime_opts, unset_str
       use cam_ccpp_cap,                 only: ccpp_physics_suite_schemes
@@ -81,7 +80,6 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon)
       call cam_logfile_readnl(nlfilename)   !The log settings must always be read first
       call physconst_readnl(nlfilename)
       call cam_initfiles_readnl(nlfilename)
-      call cam_constituents_readnl(nlfilename)
       call history_readnl(nlfilename)
       call phys_readnl(nlfilename) ! Should set phys_suite_name
       call vert_coord_readnl(nlfilename)

src/dynamics/mpas/dyn_comp_impl.F90

@@ -130,7 +130,7 @@ module subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
         use cam_abortutils, only: check_allocate
         use cam_constituents, only: const_name, const_is_water_species
         use cam_constituents, only: num_constituents, num_advected
-        use cam_constituents, only: const_is_advected, readtrace
+        use cam_constituents, only: const_is_advected
         use cam_control_mod, only: initial_run
         use cam_initfiles, only: initial_file_get_id, topo_file_get_id
         use cam_logfile, only: debugout_debug, debugout_info
@@ -216,14 +216,8 @@ module subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
                 ! the actual availability (checked and handled by MPAS).
                 call dyn_exchange_constituent_states(direction='e', exchange=.true., conversion=.false.)
 
-                ! Namelist option that controls if constituents are to be read from a file.
-                if (readtrace) then
-                    ! Read variables that belong to the "input" stream in MPAS.
-                    call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input')
-                else
-                    ! Read variables that belong to the "input" stream in MPAS, excluding constituents.
-                    call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input-scalars')
-                end if
+                ! Read variables that belong to the "input" stream in MPAS.
+                call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input')
             end if
         else
             ! Run type is branch or restart run.

src/dynamics/se/advect_tend.F90

@@ -6,7 +6,6 @@ module advect_tend
 
   use shr_kind_mod, only : r8 => shr_kind_r8
 
-  save
   private
 
   public :: compute_adv_tends_xyz
@@ -22,14 +21,16 @@ module advect_tend
   !   - second call computes and outputs the tendencies
   !----------------------------------------------------------------------
   subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
-!    use cam_history,           only: outfld, hist_fld_active
-    use time_manager,           only: get_step_size
+
+    use shr_kind_mod,          only: cl=>shr_kind_cl
+    use cam_history,           only: history_out_field
+    use time_manager,          only: get_step_size
 !    use cam_constituents,      only: tottnam,num_advected
-    use cam_constituents,       only: num_advected
-    use cam_abortutils,         only: check_allocate
+    use cam_constituents,      only: num_advected
+    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
 
@@ -43,29 +44,31 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
     logical  :: init
     real(r8), allocatable, dimension(:,:) :: ftmp
 
+    character(len=cl) :: errmsg
+
     character(len=*), parameter :: subname = 'compute_adv_tends_xyz'
 
-    if (ntrac>0) then
+    if (use_cslam) then
       nx=nc
     else
       nx=np
-    endif
+    end if
 
-    allocate( ftmp(nx*nx,nlev), stat=iret )
+    allocate( ftmp(nx*nx,nlev), stat=iret, errmsg=errmsg )
     call check_allocate(iret, subname, 'ftmp(nx*nx,nlev)', &
-                        file=__FILE__, line=__LINE__)
+                        file=__FILE__, line=__LINE__, errmsg=errmsg)
 
     init = .false.
     if ( .not. allocated( adv_tendxyz ) ) then
       init = .true.
-      allocate( adv_tendxyz(nx,nx,nlev,num_advected,nets:nete), stat=iret )
+      allocate( adv_tendxyz(nx,nx,nlev,num_advected,nets:nete), stat=iret, errmsg=errmsg )
       call check_allocate(iret, subname, 'adv_tendxyz(nx,nx,nlev,num_advected,nets:nete)', &
-                          file=__FILE__, line=__LINE__)
+                          file=__FILE__, line=__LINE__, errmsg=errmsg)
 
       adv_tendxyz(:,:,:,:,:) = 0._r8
-    endif
+    end if
 
-    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)
@@ -75,12 +78,11 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
       do ie=nets,nete
         do ic = 1, num_advected
           adv_tendxyz(:,:,:,ic,ie) = elem(ie)%state%Qdp(:,:,:,ic,qn0)/elem(ie)%state%dp3d(:,:,:,n0)  - adv_tendxyz(:,:,:,ic,ie)
-        enddo
+        end do
       end do
     end if
 
-!Remove once history outputs are enabled:
-#if 0
+#ifdef cam_thermo_history
     if ( .not. init ) then
       dt = get_step_size()
       idt = 1._r8/dt
@@ -92,11 +94,11 @@ subroutine compute_adv_tends_xyz(elem,fvm,nets,nete,qn0,n0)
               ftmp(i+(j-1)*nx,:) = adv_tendxyz(i,j,:,ic,ie)
             end do
           end do
-          call outfld(tottnam(ic), ftmp,nx*nx, ie)
+          call history_out_field(tottnam(ic), ftmp)
         end do
       end do
       deallocate(adv_tendxyz)
-    endif
+    end if
     deallocate(ftmp)
 #else
     if (.not. init) then
@@ -105,5 +107,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
+    end if
+    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
+    end if
+
+    ! 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)
+       end do
+      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
+               end do
+               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)
+         end do
+      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)
+
+    end if
+  end subroutine compute_write_iop_fields
+#endif
 end module advect_tend