#include <misc.h>#include <params.h>subroutine stepon!-----------------------------------------------------------------------!! Purpose:! Loop over time, calling driving routines for physics, dynamics, transport!! Original version:  CCM1!!-----------------------------------------------------------------------   use precision   use history, only: wshist, wrapup   use pmgrid   use pspect   use comslt   use rgrid   use prognostics   use buffer   use commap   use restart, only: write_restart   use physconst, only: cappa, gravit#if ( defined SPMD )   use mpishorthand   use spmd_dyn, only: npes, compute_gsfactors#endif#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,    only: d_p_coupling, p_d_coupling   use time_manager, only: advance_timestep, get_step_size, get_nstep, &                           is_first_step, is_first_restart_step, &                           is_last_step, is_end_curr_day   implicit none#include <comctl.h>#include <comhyb.h>#include <comlun.h>#include <comqfl.h>!------------------------------Parameters-------------------------------!   integer, parameter :: pmap = 20000 ! max dimension of evenly spaced vert. grid used !                                     ! by SLT code to map the departure pts into true !                                     ! model levels.!!---------------------------Local workspace-----------------------------!   integer kdpmpf  (pmap)             ! artificial full vert grid indices   integer kdpmph  (pmap)             ! artificial half vert grid indices   type(physics_state), dimension(begchunk:endchunk) :: phys_state   type(physics_tend ), dimension(begchunk:endchunk) :: phys_tend   real(r8) hyad   (plev)             ! del (A)   real(r8) lam    (plond,platd)      ! longitude coords of extended grid   real(r8) phi    (platd)            ! latitude  coords of extended grid   real(r8) dphi   (platd)            ! latitude intervals (radians)   real(r8) gw     (plat)             ! Gaussian weights   real(r8) sinlam (plond,platd)      ! sin(lam) model domain only              real(r8) coslam (plond,platd)      ! cos(lam) model domain only              real(r8) lbasdy (4,2,platd)        ! latitude derivative weights             real(r8) lbasdz (4,2,plev)         ! vert (full levels) deriv wghts    real(r8) lbassd (4,2,plevp)        ! vert (half levels) deriv wghts    real(r8) lbasiy (4,2,platd)        ! Lagrange cubic interp wghts (lat.)    real(r8) lbasiz (4,2,plev)         ! Lagrange cubic interp wghts (vert)   real(r8) lbassi (4,2,plevp)        ! Lagrange cubic interp wghts (vert)   real(r8) detam  (plev)             ! intervals between vert full levs.   real(r8) detai  (plevp)            ! intervals between vert half levs.   real(r8) dlam   (platd)            ! longitudinal grid interval (radians)   real(r8) cwava  (plat)             ! weight applied to global integrals   real(r8) etamid (plev)             ! vertical coords at midpoints    real(r8) etaint (plevp)            ! vertical coords at interfaces   real(r8) t2     (plond,plev,beglat:endlat)  ! temp tendency   real(r8) flx_net(plond,beglat:endlat)       ! net flux from physics   real(r8) fu     (plond,plev,beglat:endlat)  ! u wind tendency   real(r8) fv     (plond,plev,beglat:endlat)  ! v wind tendency   real(r8) coslat(plond)   real(r8) rcoslat(plond)   real(r8) rpmid(plond,plev)   real(r8) pdel(plond,plev)   real(r8) pint(plond,plevp)   real(r8) pmid(plond,plev)   real(r8) ztodt                     ! time step   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   integer kk                         ! |   integer kk1                        ! |   integer kk2                        ! | - indices   integer kstep                      ! |   integer l                          ! |   integer lvsum                      ! counter   integer lsum                       ! counter   real(r8) vsum                      ! accumulator for SLD binning statistics   logical vflag                      ! logical flag to indicate that binning has been done!                                     ! correctly   integer i,k,lat,j,begj             ! longitude,level,latitude indices   real(r8) tmp1                      ! temp space   integer year,month,day,tod         ! date info   integer ntspdy                     ! Number of timesteps per day   #ifdef SPMD   integer :: numsend          ! number of items to be sent   integer :: numrecv(0:npes-1)! number of items to be received   integer :: displs(0:npes-1) ! displacement array   integer :: numperlat        ! number of items per latitude band#endif!! Externals!   logical, external :: rstwr  ! whether or not to write restart files!!-----------------------------------------------------------------------!! Define eta coordinates: Used for calculation etadot vertical velocity ! for slt.!   call t_startf ('stepon_startup')   do k=1,plev      etamid(k) = hyam(k) + hybm(k)   end do   do k=1,plevp      etaint(k) = hyai(k) + hybi(k)   end do!! Initialize matrix gamma (used in sld T computation)!   gamma(:,:) = 0.   do k = 1,plev      tmp1 = cappa*t0(k)*hypi(plevp)/hypm(k)      gamma(k,k) = 0.5*tmp1      do l=1,k-1         gamma(l,k) = tmp1      end do   end do!! Set slt common block variables!   call grdini(pmap    ,etamid  ,etaint  ,gravit  ,dlam    , &               lam     ,phi     ,dphi    ,gw      ,sinlam  , &               coslam  ,lbasdy  ,lbasdz  ,lbassd  ,lbasiy  , &               lbasiz  ,lbassi  ,detam   ,detai   ,kdpmpf  , &               kdpmph  ,cwava   ,phigs   )   if (is_first_step()) then      do lat=beglat,endlat         j = j1 - 1 + lat         do i=1,nlon(lat)            coslat(i) = cos(clat(lat))            rcoslat(i) = 1./coslat(i)         end do!     ! Set current time pressure arrays for model levels etc.!         call plevs0(nlon(lat), plond, plev, ps(1,lat,n3), pint, pmid, pdel)!         do k=1,plev            do i=1,nlon(lat)               rpmid(i,k) = 1./pmid(i,k)            end do         end do!! Compute appropriate (1/ps)etadot(dp/deta)!         call etadt0 (lat, nlon(lat), &                      rcoslat ,div(1,1,lat,n3), u3(i1,1,j,n3), v3(i1,1,j,n3), dpsl(1,lat), &                      dpsm(1,lat), pdel, ps(1,lat,n3), ed1(1,1,lat))!! Calculate vertical motion field!         call omcalc (rcoslat, div(1,1,lat,n3), u3(i1,1,j,n3), v3(i1,1,j,n3), dpsl(1,lat), &                      dpsm(1,lat), pmid, pdel, rpmid, pint(1,plevp), &                      omga(1,1,lat), nlon(lat))      end do   end if!! Compute pdel from "A" portion of hybrid vertical grid!   do k=1,plev      hyad(k) = hyai(k+1) - hyai(k)   end do   do k=1,plev      do i=1,plon         pdela(i,k) = hyad(k)*ps0      end do   end do!!----------------------------------------------------------! Initialize departure point bin arrays to 0.!----------------------------------------------------------!   kstep = 0   do lat = beglat,endlat      do k = 1,plev         do kk = 1,plev            levkntl(kk,k,lat) = 0.         end do      end do   end do!! Beginning of basic time step loop!   call t_stopf ('stepon_startup')   do      if (masterproc .and. print_step_cost) then         call t_stampf (wcstart, usrstart, sysstart)      end if      ztodt = get_step_size()!!----------------------------------------------------------! PHYSPKG  Call the Physics package!----------------------------------------------------------!      begj = beglatex + numbnd      call t_startf('d_p_coupling')      call d_p_coupling (ps(1,beglat,n3), t3(i1,1,begj,n3), u3(i1,1,begj,n3), &                         v3(i1,1,begj,n3), q3(i1,1,1,begj,n3), &                         omga, phis, phys_state)      call t_stopf('d_p_coupling')      call t_startf('phys_driver')      if (ideal_phys) then         call phys_idealized(phys_state, phys_tend, ztodt, etamid)      else if (adiabatic) then         call phys_adiabatic(phys_state, phys_tend)      else         call physpkg (                     &            phys_state, gw, ztodt, &            phys_tend, cld(1,1,begchunk,n3m1), cld(1,1,begchunk,n3), &            tcwat(1,1,begchunk,n3m1), tcwat(1,1,begchunk,n3), &            qcwat(1,1,begchunk,n3m1), qcwat(1,1,begchunk,n3), &            lcwat(1,1,begchunk,n3m1), lcwat(1,1,begchunk,n3) )      end if      call t_stopf('phys_driver')      call t_startf('p_d_coupling')      call p_d_coupling (phys_state, phys_tend, t2, fu, fv, flx_net,q3(i1,1,1,begj,n3))      call t_stopf('p_d_coupling')      if (is_first_restart_step()) then         call print_memusage      end if!----------------------------------------------------------! DYNPKG Call the Dynamics Package!----------------------------------------------------------      call t_startf ('dynpkg')      call dynpkg (pmap    ,t2      ,fu      ,fv      ,etamid  , &                   etaint  ,cwava   ,detam   ,dlam    ,lam     , &                   phi     ,dphi    ,sinlam  ,coslam  ,lbasdy  , &                   lbasdz  ,lbasiy  ,lbassi  ,lbasiz  ,detai   , &                   kdpmpf  ,kdpmph  ,flx_net , ztodt   )      call t_stopf ('dynpkg')      if (is_first_restart_step()) then         call print_memusage      end if      call t_startf('stepon_single')! 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')!! 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 nstep before returning to top of loop!      call advance_timestep()      kstep = kstep + 1!! Compute some statistical info!      if (nlend) then#if ( defined SPMD )         numperlat = plev*plev         call compute_gsfactors (numperlat, numsend, numrecv, displs)         call mpigatherv (levkntl(1,1,beglat), numsend, mpir8, &                          levkntl            , numrecv, displs, mpir8, 0, mpicom)#endif         if (masterproc) then            do k = 1,plev               do kk = 1,plev                  levknt(kk,k) = 0.                  do lat = 1,plat                     levknt(kk,k) = levknt(kk,k) + levkntl(kk,k,lat)                  end do               end do            end do            lsum = 0            do lat = 1,plat               lsum = lsum + nlon(lat)            end do            vflag = .false.            do k = 1,plev               vsum = 0.               do kk = 1,plev                  vsum = vsum + levknt(kk,k)               end do               lvsum = vsum + 0.01                if( lvsum .ne. lsum*kstep ) vflag = .true.            end do                     do k = 1,plev               do kk = 1,plev                  levknt(kk,k) = levknt(kk,k)/float(lsum*kstep)                  if (levknt(kk,k) .eq. 0.) levknt(kk,k) = 1.e+36               end do            end do                     if(vflag) then               write(6,*) '********************************************'               write(6,1000)               write(6,*) '********************************************'            else               write(6,2000)               k = 1               write(6,3001) hypm(k)/100.,(levknt(kk,k),kk = 1, 6)               k = 2               write(6,3002) hypm(k)/100.,(levknt(kk,k),kk = 1, 7)               k = 3               write(6,3003) hypm(k)/100.,(levknt(kk,k),kk = 1, 8)               k = 4               write(6,3004) hypm(k)/100.,(levknt(kk,k),kk = 1, 9)               k = 5               write(6,3005) hypm(k)/100.,(levknt(kk,k),kk = 1,10)               do k = 6,plev                  kk1 = k-5                  kk2 = k+5                  if(kk2 .gt. plev) kk2 = plev                  write(6,3000) hypm(k)/100.,(levknt(kk,k),kk = kk1,kk2)               end do            end if         end if      end if      call t_stopf('stepon_single')!! Check for end of run!      if (nlend) then#ifdef COUP_CSM         call ccsmfin#endif         return      end if   end do  ! End of timestep loop1000 format(' ERROR in binning departure points; sums were done'// &            ' incorrectly. not printing results.')2000 format(/40x,' BINNING STATISTICS FOR VERTICAL DEPARTURE POINTS'// &                 '    level          -5        -4        -3        -2     ', &                 '   -1    0 (arr pt)     1         2         3         4        ', &                 ' 5'/)3000 format(' ',f9.4, 4x,11f10.5)3001 format(' ',f9.4,54x,11f10.5)3002 format(' ',f9.4,44x,11f10.5)3003 format(' ',f9.4,34x,11f10.5)3004 format(' ',f9.4,24x,11f10.5)3005 format(' ',f9.4,14x,11f10.5)end subroutine stepon