#include <misc.h>
#include <params.h>
!-----------------------------------------------------------------------
!BOP
! !ROUTINE: stepon --- Loop over time, perform physics and dynamics
!
! !INTERFACE:
subroutine stepon 1,44
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use history, only: wshist, wrapup, write_inithist, inithist
use pmgrid
use constituents, only: pcnst, pnats, cnst_get_type_byind
use comsrf, only: surface_state2d
use rgrid
use prognostics
use commap, only: w
use dynamics_vars, only: cosp, sinp, cose, sine, coslon, sinlon, &
ak, bk, ks, ptop
use fv_prints, only: fv_out
use restart, only: write_restart
use dynconst, only: omega, rearth
use physconst, only: gravit, rair
use abortutils, only: endrun
#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 phys_buffer, only: pbuf
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_step, is_first_restart_step, &
is_last_step, is_end_curr_day
#if defined( SPMD )
use spmd_dyn, only: comm_z, ijk_yz_to_xy
use parutilitiesmodule, only: pargatherreal, parcollective3d, sumop
use mod_comm, only: mp_sendirr, mp_recvirr
#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.
! 03.08.05 Sawyer Removed Rayleigh Friction, related variables
! 03.08.13 Sawyer Ghost zone in u3sxy removed.
!
!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
logical l_write_inithist
real(r8) te0 ! Total energy before dynamics
! Velocity tendencies on a-grid
real(r8), allocatable :: dudt(:,:,:)
real(r8), allocatable :: dvdt(:,:,:)
!delta pressure dry
real(r8), allocatable :: delpdry(:,:,:)
! 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 :: 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 :: dummy3xy(:,:,:)
!-----------------------------------------------------------------------
! Other local variables
type(physics_state), allocatable :: phys_state(:)
type(physics_tend ), allocatable :: 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 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(:,:,:)
!
! 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
!----------------------------------------------------------
! 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.
!$omp parallel do private(i,j,k)
do k = beglev, endlevp1
do j = beglat, endlat
do i = 1, plon
pk(i,j,k) = pe(i,k,j)**cappa
enddo
enddo
enddo
! Generate pkz, the conversion factor betw pt and t3
call pkez(1, plon, plev, beglat, endlat, &
beglev, endlev, 1, plon, pe, pk, cappa, &
ks, piln, pkz, .false. )
if ( .not. nlres ) then
! Compute pt for initial run: scaled virtual potential temperature
! defined as (virtual temp deg K)/pkz. pt will be written to restart (SJL)
!$omp parallel do private(i,j,k)
do k = beglev, endlev
do j = beglat, endlat
do i = 1, plon
pt(i,j,k) = t3(i,k,j)*(1.+zvir*q3(i,j,k,1))/pkz(i,j,k)
enddo
enddo
enddo
endif
!----------------------------------------------------------------
! Convert mixing ratios initialized as dry to moist for dynamics
!----------------------------------------------------------------
if ( .not. nlres ) then
! on initial time step, dry mixing ratio advected constituents have been
! initialized to dry mixing ratios. dynpkg expects moist m.r. so convert here.
! first calculate delpdry. The set_pdel_state subroutine
! is called after the dynamics in d_p_coupling to set more variables.
! This is not in tracers.F90 because it is only used by LR dynamics.
allocate (delpdry(plon,beglat:endlat,beglev:endlev))
do i =1, plon
do j = beglat, endlat
do k = beglev,endlev
delpdry(i,j,k) = delp(i,j,k)*(1.-q3(i,j,k,1))
enddo
enddo
enddo
do m = 1,ppcnst
if (cnst_get_type_byind(m).eq.'dry') then
do i=1, plon
do j = beglat, endlat
do k = beglev,endlev
q3(i,j,k,m) = q3(i,j,k,m)*delpdry(i,j,k)/delp(i,j,k)
end do
end do
end do
end if
end do
end if ! .not. nlres
allocate(phys_state(begchunk:endchunk))
allocate(phys_tend(begchunk:endchunk))
!----------------------------------------------------------
! Allocate auxiliary variables
!----------------------------------------------------------
allocate (dudt(plon,beglev:endlev,beglat:endlat))
allocate (dvdt(plon,beglev:endlev,beglat:endlat))
allocate (dummy3(plon,beglat:endlat,beglev:endlev))
!----------------------------------------------------------
! Allocate variables for secondary 2D xy decomposition
!----------------------------------------------------------
allocate (phisxy(beglonxy:endlonxy , beglatxy:endlatxy))
allocate ( psxy(beglonxy:endlonxy , beglatxy:endlatxy))
allocate (omgaxy(beglonxy:endlonxy, plev, beglatxy:endlatxy))
allocate ( u3sxy(beglonxy:endlonxy, beglatxy:endlatxy, plev)) ! now unghosted, was N1
allocate ( v3sxy(beglonxy:endlonxy, beglatxy:endlatxy, plev))
allocate (delpxy(beglonxy:endlonxy, beglatxy:endlatxy, plev))
allocate ( ptxy(beglonxy:endlonxy, beglatxy:endlatxy, plev))
allocate ( q3xy(beglonxy:endlonxy, beglatxy:endlatxy, plev, pcnst+pnats))
allocate ( pkzxy(beglonxy:endlonxy, beglatxy:endlatxy, plev))
allocate ( pkxy(beglonxy:endlonxy, beglatxy:endlatxy, plevp))
allocate ( pexy(beglonxy:endlonxy, plevp, beglatxy:endlatxy))
allocate (pilnxy(beglonxy:endlonxy, plevp, beglatxy:endlatxy))
allocate (dudtxy(beglonxy:endlonxy, plev, beglatxy:endlatxy))
allocate (dvdtxy(beglonxy:endlonxy, plev, beglatxy:endlatxy))
allocate (dummy3xy(beglonxy:endlonxy,beglatxy:endlatxy,plev))
!-----------------------------------------------------------------------
! Transpose phis if 2D decomposition
! The transposed phis (called phisxy) is needed for the geopotential
! calculation (when using transpose option) and for remapping
! Embed in 3D array since transpose machinery cannot handle 2D arrays
!-----------------------------------------------------------------------
if (twod_decomp .eq. 1) then
#if defined( SPMD )
!$omp parallel do private(i,j,k)
do k = beglev,endlev
do j = beglat,endlat
do i = 1,plon
dummy3(i,j,k) = phis(i,j)
enddo
enddo
enddo
call mp_sendirr(dummy3, ijk_yz_to_xy%SendDesc, &
ijk_yz_to_xy%RecvDesc, dummy3xy)
call mp_recvirr(dummy3xy, ijk_yz_to_xy%RecvDesc )
do j = beglatxy,endlatxy
do i = beglonxy,endlonxy
phisxy(i,j) = dummy3xy(i,j,1)
enddo
enddo
#endif
else
do j = beglat,endlat
do i = 1,plon
phisxy(i,j) = phis(i,j)
enddo
enddo
endif
!----------------------------------------------------------
! Beginning of basic time step loop
!----------------------------------------------------------
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! ATTENTION *** ATTENTION *** ATTENTION *** ATTENTION *** ATTENTION
!
!
! A 2D xy decomposition is used for handling the Lagrangian surface
! remapping, the ideal physics, and (optionally) the geopotential
! calculation.
!
! The transpose from yz to xy decomposition takes place within dynpkg.
! The xy decomposed variables are then transposed directly to the
! physics decomposition within d_p_coupling.
!
! The xy decomposed variables have names corresponding to the
! yz decomposed variables: simply append "xy". Thus, "uxy" is the
! xy decomposed version of "u".
!
! To assure that the latitudinal decomposition operates
! as efficiently as before, a separate parameter "twod_decomp" has
! been defined; a value of 1 refers to the multi-2D decomposition with
! transposes; a value of 0 means that the decomposition is effectively
! one-dimensional, thereby enabling the transpose logic to be skipped;
! there is an option to force computation of transposes even for case
! where decomposition is effectively 1-D.
!
! For questions/comments, contact Art Mirin, mirin@llnl.gov
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
do
if (masterproc .and. print_step_cost) then
call t_stampf (wcstart, usrstart, sysstart)
end if
!--------------------------------------------------------------------------
! Perform finite-volume dynamics -- this dynamical core contains some
! yet to be published algorithms. Its use in the CAM is
! for software development purposes only.
! Please contact S.-J. Lin (slin@dao.gsfc.nasa.gov)
! if you plan to use this mudule for scientific purposes. Contact S.-J. Lin
! or Will Sawyer (sawyer@dao.gsfc.nasa.gov) if you plan to modify the
! software.
!--------------------------------------------------------------------------
call t_startf('dynpkg')
!----------------------------------------------------------
! For 2-D decomposition, phisxy is input to dynpkg, and the other
! xy variables are output. Some are computed through direct
! transposes, and others are derived.
!----------------------------------------------------------
call dynpkg(plon, plat, plev, u3s, v3s, &
ps, delp, pe, pk, pkz, &
piln, ptop, beglat, endlat, &
beglev, endlev, endlevp, phis, &
omga, cpair, pcnst, q3, pdt, &
omega, rair, rearth, nsplit, pcnst+pnats, &
pt, full_phys, iord, jord, &
kord, te0, &
u3sxy, v3sxy, psxy, delpxy, pexy, &
pkxy, pkzxy, pilnxy, phisxy, omgaxy, &
q3xy, ptxy, &
beglonxy, endlonxy, beglatxy, endlatxy )
call t_stopf('dynpkg')
call t_startf('physics')
call t_startf('dp_coupling_total')
call t_startf('d_p_coupling')
!----------------------------------------------------------
! Move data into phys_state structure.
!----------------------------------------------------------
call d_p_coupling (ps, u3s, v3s, pt, coslon, sinlon, &
q3, omga, phis, pe, piln, pk, &
pkz, phys_state, phys_tend, pbuf, full_phys, &
psxy, u3sxy, v3sxy, ptxy, &
q3xy, omgaxy, phisxy,pexy, pilnxy, &
pkxy, pkzxy )
call t_stopf('d_p_coupling')
call t_stopf('dp_coupling_total')
!----------------------------------------------------------
! Call full physics
!----------------------------------------------------------
if ( full_phys ) then
!
call t_startf('physpkg')
call physpkg (phys_state, w, dtime, phys_tend, pbuf)
call t_stopf('physpkg')
!----------------------------------------------------------
! Otherwise, adiabatic physics to update calendar and history
!----------------------------------------------------------
else
call phys_adiabatic (phys_state, phys_tend)
endif
!----------------------------------------------------------
! Update dynamics variables using phys_state & phys_tend.
! 2-D decomposition: Compute dudtxy, dvdtxy, ptxy and q3xy; for ideal
! physics, scale ptxy by (old) pkzxy; then transpose to yz variables
! 1-D decomposition: Compute dudt, dvdt, pt and q3; for ideal physics,
! scale pt by old pkz.
! Call uv3s_update to update u3s and v3s from dudt and dvdt.
! Call p_d_adjust to update pt, q3, pe, delp, ps, piln, pkz and pk.
! For adiabatic case, transpose to yz variables.
!----------------------------------------------------------
call t_startf('dp_coupling_total')
call t_startf('p_d_coupling')
call p_d_coupling(phys_state, phys_tend, full_phys, &
adiabatic, t3, &
q3, pt, dudt, dvdt, pkz, &
q3xy, ptxy, dudtxy, dvdtxy, pkzxy, &
dtime, u3s, v3s, u3sxy, v3sxy, &
zvir, cappa, ptop, pk, &
piln, ps, &
pe, pexy, delp, delpxy )
call t_stopf('p_d_coupling')
call t_stopf('dp_coupling_total')
!----------------------------------------------------------
! Monitor max/min/mean of selected fields
!
! SEE BELOW **** SEE BELOW **** SEE BELOW
! Beware that fv_out uses both dynamics and physics instanciations.
! However, I think that they are used independently, so that the
! answers are correct. Still, this violates the notion that the
! physics state is no longer active after p_d_coupling.
!----------------------------------------------------------
call get_curr_date(yr, mon, day, ncsec)
ncdate = yr*10000 + mon*100 + day
i = pdt + ncsec ! step complete, but nstep not incremented yet
if ( fv_monitor .and. mod(i, freq_diag) == 0 ) then
call fv_out( plon, plat , plev, beglat, endlat, &
ng_d, beglev, endlev, pk, pt , &
ptop, ps, q3, pcnst+pnats, pcnst, &
delp, pe, surface_state2d, phys_state, ncdate, &
i, full_phys )
endif
call t_stopf('physics')
if (is_first_step() .or. is_first_restart_step()) then
call print_memusage ('stepon after dynpkg')
end if
! Set end of run flag.
#if ( ! defined COUP_CSM )
if (is_last_step()) nlend = .true.
#else
if (csmstop) then
if ( masterproc ) write(6,*)'atm: Stopping at the end of this day'
if (is_end_curr_day()) nlend = .true.
end if
#endif
!----------------------------------------------------------
! History and restart logic
!----------------------------------------------------------
! Write and/or dispose history tapes if required
call t_startf ('wshist')
call wshist ()
call t_stopf ('wshist')
!
! Shift time pointers (used in cloud water). For LR, this must be
! done after physpkg and before writing the restart file.
!
call shift_time_indices ()
! Write restart file
if (rstwr() .and. nrefrq /= 0) then
call t_startf ('write_restart')
call write_restart
call t_stopf ('write_restart')
end if
! Dispose necessary files
call t_startf ('wrapup')
call wrapup
call t_stopf ('wrapup')
if (masterproc .and. print_step_cost) then
call t_stampf (wcend, usrend, sysend)
write(6,'(a,3f8.3,a)')'Prv timestep wallclock, usr, sys=', &
wcend-wcstart, usrend-usrstart, sysend-sysstart, ' seconds'
end if
! Increment timestep before returning to top of loop
call advance_timestep()
! Write initial file if requested.
call get_curr_date(yr, mon, day, ncsec)
l_write_inithist = .false.
if (inithist == '6-HOURLY' ) then
if (mod(ncsec, 21600) == 0) l_write_inithist = .true.
end if
if (ncsec == 0 .and. inithist == 'DAILY' ) l_write_inithist = .true.
if (ncsec == 0 .and. day == 1 .and. inithist == 'MONTHLY') l_write_inithist = .true.
if (ncsec == 0 .and. day == 1 .and. mon == 1 .and. inithist == 'YEARLY' ) l_write_inithist = .true.
if (l_write_inithist) call write_inithist
! Check for end of run
if (nlend) then
call print_memusage ('End stepon')
deallocate(phys_state)
deallocate(phys_tend)
deallocate(dudt)
deallocate(dvdt)
deallocate(dummy3)
deallocate(phisxy)
deallocate( psxy)
deallocate(omgaxy)
deallocate( u3sxy)
deallocate( v3sxy)
deallocate(delpxy)
deallocate( ptxy)
deallocate( q3xy)
deallocate( pkzxy)
deallocate( pkxy)
deallocate( pexy)
deallocate(pilnxy)
deallocate(dudtxy)
deallocate(dvdtxy)
deallocate(dummy3xy)
#ifdef COUP_CSM
call ccsmfin
#endif
return
end if
end do ! End of timestep loop
!EOC
end subroutine stepon
!-----------------------------------------------------------------------