stepon.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!$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


!----------------------------------------------------------
! Allocate auxiliary variables
!----------------------------------------------------------

   allocate (dudt(plon,beglev:endlev,beglat:endlat))
   allocate (dvdt(plon,beglev:endlev,beglat:endlat))
   allocate (qtmp(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 (  t3xy(beglonxy:endlonxy, plev, beglatxy:endlatxy))
      allocate (omgaxy(beglonxy:endlonxy, plev, beglatxy:endlatxy))    
      allocate ( u3sxy(beglonxy:endlonxy, beglatxy:endlatxy+1, plev))  ! ghosted 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 (qtmpxy(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
!-----------------------------------------------------------------------

#if defined( SPMD )
   if (twod_decomp .eq. 1) then
!$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 redistributestart (inter_ijk, .true., dummy3)
      call redistributefinish(inter_ijk, .true., dummy3xy)
!$omp parallel do private(i,j)
      do j = beglatxy,endlatxy
         do i = beglonxy,endlonxy
            phisxy(i,j) = dummy3xy(i,j,1)
         enddo
      enddo
   endif
#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".
!
!     The Rayleigh friction has been put into its own routine
!     called rayl_fric.
!
!     Ideal_physics (hswf) is now called directly from dynpkg.
!
!     To assure that the latitudinal decomposition operates
!     as efficiently as before, a separate parameter "twod_decomp" has
!     been defined; a value of 0 means the original 1D decomposition
!     (with the 2D decomposition formalism skipped), and a value of 1
!     utilizes the 2D decomposition formalism.
!
!     For questions/comments, contact Art Mirin, [e-mail]
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

   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 ([e-mail])
! if you plan to use this mudule for scientific purposes. Contact S.-J. Lin
! or Will Sawyer ([e-mail]) 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.
!
! New edge pressure pexy is back-transposed to pe, and new surface
! pressure ps is computed. (For full physics, pe is then used and
! recomputed, along with ps, after physics advance.)
!
! For simplified physics, new delpxy is back-transposed to delp.
! (For full physics, delp is recomputed after physics advance, so
! the current delpxy is temporary and there is no need to back-transpose.)
!
! For ideal physics, hswf is called directly from dynpkg.
!----------------------------------------------------------
      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,      t3,      full_phys, iord,  jord,        &
                  kord,    te0,                                    &
                  u3sxy,   v3sxy,   psxy,    delpxy,  pexy,        &
                  pkxy,    pkzxy,   pilnxy,  phisxy,  omgaxy,      &
                  q3xy,    ptxy,    t3xy,                          &
                  beglonxy, endlonxy, beglatxy, endlatxy,          &
                  dcaf,    rayf,    ideal_phys )

      call t_stopf('dynpkg')

      call t_startf('physics')

      call t_startf('d_p_coupling')
!----------------------------------------------------------
! Move data into phys_state structure.
! For 2-D decompositon, qtmpxy is back-transposed to qtmp.
!----------------------------------------------------------
      call d_p_coupling (ps,   u3s,  v3s,   pt,    coslon, sinlon,    &
                         t3,   q3,  omga,  phis,  pe,     piln,   pk, &
                         pkz,  phys_state,  phys_tend,    full_phys,  &
                         qtmp, psxy, u3sxy, v3sxy, ptxy,  t3xy,       &
                         q3xy, omgaxy,      phisxy,pexy,  pilnxy,     &
                         pkxy, pkzxy,  qtmpxy,     pe11k, pe11kln )
      call t_stopf('d_p_coupling')


!----------------------------------------------------------
! Call full physics
!----------------------------------------------------------
      if ( full_phys ) then
!
         call t_startf('physpkg')
         call physpkg (phys_state,  gw,  dtime,  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)  )
         call t_stopf('physpkg')
!----------------------------------------------------------
! Otherwise, adiabatic physics to update calendar and history
!----------------------------------------------------------
      else
         call phys_adiabatic (phys_state, phys_tend)
      endif
!----------------------------------------------------------
! Note: ideal physics (hswf) is called directly from dynpkg
!----------------------------------------------------------

!----------------------------------------------------------
! Update Rayleigh friction.
!----------------------------------------------------------
      if ( .not. adiabatic ) then
         call t_startf('physpkg')
         call rayl_fric(phys_state, phys_tend, dtime, pe11k, pe11kln,        &
                        cpair, cappa, rfac, rayf )
         call t_stopf('physpkg')
      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('p_d_coupling')
      call p_d_coupling(phys_state,   phys_tend,   full_phys,      &
                        adiabatic,                                 &
                        q3,    pt,    dudt,    dvdt,     pkz,      &
                        q3xy,  ptxy,  dudtxy,  dvdtxy,   pkzxy,    &
                        dtime, u3s,   v3s,     u3sxy,    v3sxy,    &
                        zvir,  cappa, ptop,    pk,                 &
                        piln,  ps,    qtmp,                        &
                        pe,    pexy,  delp,    delpxy     )
      call t_stopf('p_d_coupling')

!----------------------------------------------------------
! 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_restart_step()) call print_memusage

! 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')
#if defined( SPMD )
      if (ncsec==0 .and. day==1) then
         if (inithist=='MONTHLY' .or. inithist=='YEARLY' .and. mon==1) then
            if (twod_decomp .eq. 1) then
                call redistributestart (inter_ikj, .false., t3xy)
                call redistributefinish(inter_ikj, .false., t3)
            endif
         endif
      endif
#endif
      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()

! Check for end of run

      if (nlend) then
#ifdef COUP_CSM
         call ccsmfin
#endif
         return
      end if

   end do  ! End of timestep loop
!EOC
end subroutine stepon
!-----------------------------------------------------------------------