#include <misc.h>
#include <params.h>
#if ( defined SCAM )
#include <max.h>
#endif

subroutine stepon 1,43
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Loop over time, calling driving routines for physics, dynamics, 
! transport
! 
! Method: 
! 
! Author: 
! Original version:  CCM1
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, D. Williamson, August 1992
! Reviewed:          B. Boville, D. Williamson, April 1996
! Restructured:      J. Truesdale, May 1999
!
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use history, only: wshist, wrapup, write_inithist, inithist
   use pmgrid
   use pspect
   use comslt
   use rgrid
   use prognostics
   use restart, only: write_restart
#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,    only: d_p_coupling, p_d_coupling
   use commap
   use physconst, only: gravit
   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, get_curr_date
#if ( defined BFB_CAM_SCAM_IOP )
   use iop, only: init_iop_fields,fixmassav,betasav,alphasav,divt3dsav,divq3dsav,dqfx3savm1
#endif
#if ( defined SCAM )
   use scamMod, only :use_iop,use_pert_frc,doiopupdate
   use iop, only: init_iop_fields,fixmassav,betasav,alphasav,divt3dsav,divq3dsav,dqfx3savm1
#endif
   implicit none

   integer pmap   ! max dimension of evenly spaced vert. grid used 
!                    ! by SLT code to map the departure pts into true 
!                    ! model levels.
   parameter ( pmap = 20000 )
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comhyb.h>
!-----------------------------------------------------------------------
#include <comlun.h>
!-----------------------------------------------------------------------
#include <comqfl.h>
!-----------------------------------------------------------------------
#if ( defined SCAM )
#include <comfrc.h>
!-----------------------------------------------------------------------
#endif
!
   integer kdpmpf(pmap)             ! artificial full vert grid indices
   integer kdpmph(pmap)             ! artificial half vert grid indices

   type(physics_state), allocatable :: phys_state(:)
   type(physics_tend ), allocatable :: 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) 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), allocatable :: t2(:,:,:) ! temp tendency
   real(r8), allocatable :: fu(:,:,:) ! u wind tendency
   real(r8), allocatable :: fv(:,:,:) ! v wind tendency
   real(r8) flx_net(plond,beglat:endlat)       ! net flux from physics
   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) dtime               ! timestep size
   real(r8) ztodt               ! twice time step unless nstep=0
   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 i,k,lat,j,begj       ! longitude,level,latitude indices
   integer iter
   integer :: yr, mon, day    ! year, month, and day components of a date
   integer :: ncsec           ! current time of day [seconds]
   logical l_write_inithist

#if ( defined SCAM )
   integer dummy
   real(r8) omegapert
   real(r8) pertval              ! perturbation valie
   real(r8) drand48
   external drand48
#endif
!
! Externals
!
   logical, external :: rstwr  ! whether or not to write restart files
!
!-----------------------------------------------------------------------
   call t_startf ('stepon_startup')
   dtime = get_step_size()

!
! Define eta coordinates: Used for calculation etadot vertical velocity 
! for slt.
!
   do k=1,plev
      etamid(k) = hyam(k) + hybm(k)
   end do
   do k=1,plevp
      etaint(k) = hyai(k) + hybi(k)
   end do

#if ( defined SCAM )
   call plevs0(nlon(beglat), plond, plev, ps(1,beglat,n3), pint, pmid, pdel)
   etamid(:)=pmid(1,:)	
   etaint(:)=pint(1,:)	
#endif
!
! Set slt common block variables
!
   call grdini(pmap    ,etamid  ,etaint  ,gravit  ,dlam    , &
               lam     ,phi     ,dphi    ,gw      ,sinlam  , &
               coslam  ,lbasdy  ,lbasdz  ,lbassd  ,lbasiy  , &
               detam   ,detai   ,kdpmpf  ,kdpmph  ,cwava   )
!
! Initial guess for trajectory midpoints in spherical coords.
! nstep = 0:  use arrival points as initial guess for trajectory midpoints.
! nstep > 0:  use calculated trajectory midpoints from previous time 
! step as first guess.
! NOTE:  reduce number of iters necessary for convergence after nstep = 1.
!
   if (is_first_step()) then
      do lat=beglat,endlat
         j = j1 - 1 + lat
         do k=1,plev
            do i=1,nlon(lat)
#if ( !defined SCAM )
               lammp(i,k,lat) = float(i-1)*dlam(j1-1+lat)
               phimp(i,k,lat) = clat(lat)
               sigmp(i,k,lat) = etamid(k)
#else
               sigmp(i,k,lat) = pmid(i,k)
#endif
            end do
         end do
#if ( !defined SCAM )
         do i=1,nlon(lat)
            coslat(i) = cos(clat(lat))
            rcoslat(i) = 1./coslat(i)
         end do
#endif
!     
! 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
#if ( ! defined SCAM )
!
! 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))
#else
         
         omga(1,:,lat)=wfld(:)
#endif
      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

   allocate(phys_state(begchunk:endchunk))
   allocate(phys_tend(begchunk:endchunk))
   allocate(t2(plond,plev,beglat:endlat))
   allocate(fu(plond,plev,beglat:endlat))
   allocate(fv(plond,plev,beglat:endlat))
!
! Beginning of basic time step loop
!
   call t_stopf ('stepon_startup')


#if ( defined BFB_CAM_SCAM_IOP )
   if (is_first_step()) then
      begj = beglatex + numbnd
      call init_iop_fields(ps(1,beglat,n3m2), t3(i1,1,begj,n3m2), &
        u3(i1,1,begj,n3m2), v3(i1,1,begj,n3m2), q3(i1,1,1,begj,n3m2))
      fixmassav(:)=0.
      betasav(:)=0.
      alphasav(:,:)=0.
      divt3dsav(:,:,:)=0.
      divq3dsav(:,:,:,:)=0.
      dqfx3savm1(:,:,:,:)=0.
   endif
#endif

! Begin time loop.

   do
      call t_startf('stepon_st')
      if (masterproc .and. print_step_cost) then
         call t_stampf (wcstart, usrstart, sysstart)
      end if

      ztodt = 2.0*dtime
!
! If initial time step adjust dt
!
      if (is_first_step()) ztodt = dtime
!
! adjust hydrostatic matrices if the time step has changed.  This only
! happens on transition from time 0 to time 1. 
      if (get_nstep() == 1) then
         call settau(dtime)
      end if
!
!----------------------------------------------------------
! PHYSPKG  Call the Physics package
!----------------------------------------------------------
!
      begj = beglatex + numbnd
      call t_stopf('stepon_st')
      call t_startf('d_p_coupling')
      call d_p_coupling (ps(1,beglat,n3m2), t3(i1,1,begj,n3m2), u3(i1,1,begj,n3m2), &
                         v3(i1,1,begj,n3m2), q3(i1,1,1,begj,n3m2), &
                         omga, phis, phys_state, phys_tend, pbuf)
      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, pbuf)
      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, &
                         qminus(i1,1,1,begj),  q3(i1,1,1,begj,n3))
      call t_stopf('p_d_coupling')
#if (defined SCAM )
!
! In order to provide forcasted values for u,v, and ps from a ccm derived
! iop dataset we must advance timestep and perform iop read here.  Dynamics
! residules are also read here to allow forcasting of q and t.
 
     call advance_timestep()
     call get_curr_date(yr, mon, day, ncsec)
!
!=====================================================================
!     Determine whether it is time for an IOP update;
!     doiopupdate set to true if model time step > next available IOP
!     dataset time step
!
     if (use_iop) then
        call setiopupdate
     end if
!
!     Update IOP properties e.g. omega, divT, divQ
!
     if (doiopupdate) then
        call readiopdata(dummy) ! not checking error code
!
!  Randomly perturb omega if requested
!         
        omegapert = .30             ! +/- 30 %
!         
        if (use_pert_frc .eqv. .true. ) then
           do k=1,plev
              pertval = omegapert * (2.*(0.5 - drand48()))
              wfld(k) = wfld(k) * (1. + pertval )
           end do
        end if
!
     end if                    ! (do_iopupdate)

!************************************************************************
!     IF surface pressure changes with time we need to remap the vertical
!     coordinate for the slt advection process.  It has been empirically
!     determined that we can get away with 500 for pmap (instead of 20000)
!     This is necessary to make the procedure computationally feasible
!

     call grdini(pmap    ,etamid    ,etaint    ,gravit  ,dlam    , &
          lam     ,phi     ,dphi    ,gw      ,sinlam  , &
          coslam  ,lbasdy  ,lbasdz  ,lbassd  ,lbasiy  , &
          detam   ,detai   ,kdpmpf  ,kdpmph  ,cwava   )
      
#endif
!----------------------------------------------------------
! DYNPKG Call the Dynamics Package
!----------------------------------------------------------

      call t_startf('dynpkg')
      call dynpkg (t2      ,fu      ,fv      ,etamid  ,etaint  , &
                   cwava   ,detam   ,dlam    ,lam     ,phi     , &
                   dphi    ,sinlam  ,coslam  ,lbasdy  ,lbasdz  , &
                   lbassd  ,lbasiy  ,detai   ,kdpmpf  ,kdpmph  , &
                   flx_net ,ztodt   )
      call t_stopf('dynpkg')

      call t_startf('stepon_st')
      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

#if ( ! defined SCAM )
!
!----------------------------------------------------------
! 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
!
! Advance timestep before returning to top of loop
!
      call advance_timestep()
      call get_curr_date(yr, mon, day, ncsec)
      
! Write initial file if requested.

      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

#else                  ! else if running SCAM then end timestep and dealloc
      nlend = .true.
#endif
      call t_stopf('stepon_st')
!
! Check for end of run
!
      if (nlend) then
         call print_memusage ('End stepon')
         deallocate(phys_state)
         deallocate(phys_tend)
         deallocate(t2)
         deallocate(fu)
         deallocate(fv)
#ifdef COUP_CSM
         call ccsmfin
#endif
         return
      end if

   end do  ! End of timestep loop

end subroutine stepon