stepon.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  stepon --- Loop over time, perform physics and dynamics
!
! !INTERFACE:
subroutine stepon

! !USES:
   use precision
   use history, only: wshist, wrapup, inithist
   use pmgrid
   use comsrf, only: surface_state2d
   use rgrid
   use prognostics
   use buffer
   use commap, only: gw  ! actually this should be from dynamics_vars,
                         ! but, if so, non-zero diffs occur
   use dynamics_vars, only: cosp, sinp, cose, sine, coslon, sinlon,    &
                      ak, bk, ks, rfac, ptop
   use fv_prints, only: fv_out
   use restart, only: write_restart
   use dynconst, only: omega, rearth
   use physconst, only: gravit, rair

#if (defined COUP_CSM)
   use ccsm_msg, only: csmstop, ccsmfin
#endif

   use ppgrid,         only: begchunk, endchunk
   use physics_types,  only: physics_state, physics_tend
   use dp_coupling
   use physconst, only: zvir, cappa, rga, cpair
   use time_manager, only: advance_timestep, get_step_size, get_nstep, &
                           get_curr_date, is_first_restart_step, &
                           is_last_step, is_end_curr_day

#if defined( SPMD )
   use spmd_dyn, only: comm_z, inter_ijk, inter_ikj
   use parutilitiesmodule, only: pargatherreal, parcollective3d, sumop
   use redistributemodule, only : redistributestart, redistributefinish
#endif

   implicit none

#include <comctl.h>
#include <comhyb.h>
#include <comlun.h>
#include <comsta.h>

!
! !DESCRIPTION:
!
!   Loop over time, calling driving routines for physics, dynamics, 
!   transport
!
!   Warning: The Lin-Rood dynamical core ghost points are hidden from
!   the user. If plon is not equal to plond there may be unexpected 
!   results (especially since many physics arrays are defined with 
!   plond). (Glenn Grant)

! !REVISION HISTORY: 
!   00.01.15    Grant      Creation
!   00.08.23    Lin        Modifications
!   01.03.26    Sawyer     Added ProTeX documentation
!   01.06.12    Sawyer     Use dynamics_vars for LR-specific static variables
!   01.06.28    Mirin      Many changes for multi-2D decomposition;
!                          see comments near beginning of main time loop.
!   01.07.13    Mirin      Hid multi-2D decomposition coding
!   01.12.20    Mirin      Changes to fv_out interface
!   02.02.13    Eaton      Pass surface_state2d through fv_out interface.
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
   integer i,k,j,m             ! longitude, level, latitude and tracer indices
   integer t                   ! history tape indices
   integer :: ncdate           ! current date in integer format [yyyymmdd]
   integer :: ncsec            ! time of day relative to current date [seconds]
   integer :: yr, mon, day     ! year, month, day components of a date

   real(r8) te0          ! Total energy before dynamics

   real (r8) pe11k(plev+1), pe11kln(plev+1) ! Pres. & log for Rayleigh fric.

! Velocity tendencies on a-grid
   real(r8), allocatable :: dudt(:,:,:)
   real(r8), allocatable :: dvdt(:,:,:)
! Moisture at beginning of timestep - temporary storage
   real(r8), allocatable :: qtmp(:,:,:)
! Work array
   real(r8), allocatable :: dummy3(:,:,:)

!-----------------------------------------------------------------------
! The following arrays are for secondary 2D x-y decomposition
!-----------------------------------------------------------------------

   real(r8), allocatable :: phisxy(:,:)    ! Surface geopotential
   real(r8), allocatable ::   psxy(:,:)    ! Surface pressure
   real(r8), allocatable ::   t3xy(:,:,:)  ! (virtual) Temperature
   real(r8), allocatable :: omgaxy(:,:,:)  ! vertical pressure velocity
   real(r8), allocatable ::  u3sxy(:,:,:)  ! Staggered grid winds, latitude
   real(r8), allocatable ::  v3sxy(:,:,:)  ! Satggered grid winds, longitude
   real(r8), allocatable :: delpxy(:,:,:)  ! delta pressure
   real(r8), allocatable ::   ptxy(:,:,:)  ! virtual potential temperature
   real(r8), allocatable ::   q3xy(:,:,:,:)! Moisture and constituents
   real(r8), allocatable ::  pkzxy(:,:,:)  ! finite-volume mean pk
   real(r8), allocatable ::   pkxy(:,:,:)  ! pe**cappa
   real(r8), allocatable ::   pexy(:,:,:)  ! edge pressure
   real(r8), allocatable :: pilnxy(:,:,:)  ! ln(pe)
   real(r8), allocatable :: dudtxy(:,:,:)
   real(r8), allocatable :: dvdtxy(:,:,:)
   real(r8), allocatable :: qtmpxy(:,:,:)
   real(r8), allocatable :: dummy3xy(:,:,:)

!-----------------------------------------------------------------------

! Other local variables

   type(physics_state), dimension(begchunk:endchunk) :: phys_state
   type(physics_tend ), dimension(begchunk:endchunk) :: phys_tend

   real(r8) :: wcstart, wcend   ! wallclock timestamp at start, end of timestep
   real(r8) :: usrstart, usrend ! user timestamp at start, end of timestep
   real(r8) :: sysstart, sysend ! sys timestamp at start, end of timestep

   logical rayf         ! Rayleigh friction flag (off by default)
   logical dcaf         ! Dry convection flag (use only in ideal_phys case)
   logical full_phys    ! flag to convert pt output from fvcore: on output
                        ! pt will be virtual temp (deg K) if full_phys is .T.
                        ! virtual potential temperature if false.
   integer pdt          ! Physics time step
   real(r8) :: dtime    ! Physics time step

! for fv_out
   integer freq_diag          ! Output frequency in seconds
   logical fv_monitor         ! Monitor Mean/Max/Min fields every time step
   data freq_diag  / 21600 /  ! time interval (sec) for calling fv_out
   data fv_monitor / .true. / ! This is CPU-time comsuming; set it to false for
                              ! production runs

   real (r8), allocatable :: delpz(:,:,:), pez(:,:,:)

   parameter (rayf = .false.)         ! off
   parameter (dcaf = .false.)         ! off 
!
! Externals
!
   logical, external :: rstwr  ! whether or not to write restart files

!-----------------------------------------------------------------------

   if ( use_eta ) then
!-----------------------------------------------
! Use ak and bk from the internal set_eta routine
!-----------------------------------------------
      do k = 1, plev+1
         if (iam == 0) write(6,*) k, (hyai(k)-ak(k)*1.e-5),(hybi(k)-bk(k))
         hyai(k) = ak(k) * 1.e-5
         hybi(k) = bk(k)
      end do
   else
!-----------------------------------------
! Use ak and bk as specified by CAM IC
!-----------------------------------------
      do k = 1, plev+1
         ak(k) = hyai(k) * 1.e5
         bk(k) = hybi(k)
         if( bk(k) == 0. ) ks = k-1
      end do
      ptop = ak(1)
      if ( iam == 0 ) then
         write(*,*) 'Using hyai & hybi from IC:', 'KS=',ks,' PTOP=',ptop
      endif
   endif

! Compute gw to be used in physpkg (move to other place?)

   do j=2,plat-1
      gw(j) = sine(j+1) - sine(j)
   end do

      gw(   1) =  1. + sine(2)
      gw(plat) =  1. - sine(plat)

!----------------------------------------------------------
! Lin-Rood dynamical core initialization
!----------------------------------------------------------

   pdt = get_step_size()    ! Physics time step
   dtime = pdt

   if (plon .ne. plond) then
      print *, "STEPON: PLOND must be set equal to PLON when using"
      print *, "the Lin-Rood dynamical core. Stopping."
      print *, "plond = ",plond
      print *, "plon  = ",plon
      call endrun
   endif

#if (!defined STAGGERED)
   print *,"STEPON: pre-processor variable STAGGERED must be set"
   print *,"in you misc.h file. Enter: #define STAGGERED"
   print *,"Then recompile CAM. Quitting."
   call endrun
#endif

! full_phys is set to true, the output from fvcore is virtual temperature
! (deg K), as needed by the physics.  If full_phys is false, the output
! is virtual potential temperature, as needed in the idealized case.
! Note: zvir=0 in idealized case. Therefore, pt is simply the scaled
! potential temp.

   full_phys = .true.

   if ( ideal_phys .or. adiabatic )  then
       full_phys = .false.
       zvir = 0.
   endif

   if (myid_z .eq. 0) then

!$omp parallel do private(i,j)

      do j = beglat, endlat
         do i=1, plon
            pe(i,1,j) = ptop
         enddo
      enddo

   endif

   if ( nlres ) then
!
! Do not recalculate delta pressure (delp) if this is a restart run.
! Re. SJ Lin: The variable "delp" (pressure thikness for a Lagrangian
! layer) must be in the restart file. This is because delp will be
! modified "after" the physics update (to account for changes in water
! vapor), and it can not be reproduced by surface pressure and the
! ETA coordinate's a's and b's.

      if (npr_z .eq. 1) then

!$omp parallel do private(i,j,k)
         do j = beglat, endlat
            do k=1, plev
               do i=1, plon
                  pe(i,k+1,j) = pe(i,k,j) + delp(i,j,k)
               enddo
            enddo
         enddo

      else

#if defined( SPMD )
         allocate (delpz(plon, beglat:endlat, plev))
         allocate (pez(plon, plevp, beglat:endlat))
         delpz(:,:,:) = 0.
         pez(:,:,:) = 0.
!$omp parallel do private(i,j,k)
         do j = beglat, endlat
            do k=1, plevp
               do i=1, plon
                  pez(i,k,j) = 0.
               enddo
            enddo
         enddo
!$omp parallel do private(i,j,k)
         do k=1, plev
            do j = beglat, endlat
               do i=1, plon
                  delpz(i,j,k) = 0.
               enddo
            enddo
         enddo
         call pargatherreal(comm_z, 0, delp, strip3zaty, delpz)
         call parcollective3d(comm_z, sumop, plon, endlat-beglat+1, plev, delpz)

!$omp parallel do private(i,j)
   do j = beglat, endlat
      do i=1, plon
         pez(i,1,j) = ptop
      enddo
   enddo

!$omp parallel do private(i,j,k)
         do j = beglat, endlat
            do k=1, plev
               do i=1, plon
                  pez(i,k+1,j) = pez(i,k,j) + delpz(i,j,k)
               enddo
            enddo
         enddo

!$omp parallel do private(i,j,k)
         do j = beglat, endlat
            do k=beglev, endlevp1
               do i=1, plon
                  pe(i,k,j) = pez(i,k,j)
               enddo
            enddo
         enddo

#endif
      endif

   else
 
! Initial run --> generate pe and delp from the surface pressure
 
!$omp parallel do private(i,j,k)
         do j = beglat, endlat
            do k=max(2,beglev),endlevp1
               do i=1, plon
                  pe(i,k,j) = ak(k) + bk(k) * ps(i,j)
               enddo
            enddo
         enddo

!$omp parallel do private(i,j,k)
         do k = beglev,endlev
            do j = beglat, endlat
               do i=1, plon
                  delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
               enddo
            enddo
         enddo
   endif

!----------------------------------------------------------
! Check total dry air mass; set to 982.22 mb if initial run
! Print out diagnostic message if restart run
!----------------------------------------------------------

  if ( full_phys ) then
     call dryairm( plon,   plat,  plev, beglat, endlat,       &
                   ng_d,   beglev,endlev,                     &
                   .true., ptop,  ps,   q3,     pcnst+pnats,  &
                   pcnst,  delp,  pe,   nlres )
  endif

! Initialize pk, edge pressure to the cappa power.