#include <misc.h>!-----------------------------------------------------------------------!BOP! !ROUTINE:  dynpkg --- Driver for the NASA finite-volume dynamical core!! !INTERFACE:subroutine dynpkg( im,      jm,      km,      u,       v,      &                   ps,      delp,    pe,      pk,      pkz,    &                   peln,    ptop,    jfirst,  jlast,           &                   kfirst,  klast,   klastp,  phis,            &                   omga,    cp,      nq,      q3,      ndt,    &                   om,      rair,    ae,      ns0,     nc,     &                   pt,      tvm,     convt,   iord,    jord,   &                   kord,    te0,     uxy,     vxy,     psxy,   &                   delpxy,  pexy,    pkxy,    pkzxy,   pelnxy, &                   phisxy,  omgaxy,  q3xy,    ptxy,    tvmxy,  &                   ifirstxy,ilastxy, jfirstxy,jlastxy, dcaf,   &                   rayf,    ideal )! !USES:   use precision   use dynamics_vars, only : coslon, sinlon, cosl5, sinl5, acosp,  &                       acap, rcap, sine, cosp, sinp, cose,   &                       ns, icd, jcd, ng_c, ng_d, ng_s, q3t, yzt, xyt   use pmgrid, only: twod_decomp, myid_z, npr_z, npr_y, iam   use physconst, only: cappa, gravit#if defined( SPMD )   use spmd_dyn, only : comm_z, inter_ijk, inter_ikjp, inter_q3   use mod_comm, only : bufferpack3d, bufferunpack3d, buff_s, buff_r,   &                        mp_send, mp_recv, mp_barrier   use redistributemodule, only : redistributetype, redistributestart,  &                  redistributefinish#endif   implicit none! !INPUT PARAMETERS:   integer, intent(in):: im        ! dimension in east-west   integer, intent(in):: jm        ! dimension in North-South   integer, intent(in):: km        ! number of Lagrangian layers   integer, intent(in):: jfirst    ! starting latitude index for MPI   integer, intent(in):: jlast     ! ending latitude index for MPI   integer, intent(in):: kfirst    ! starting vertical index for MPI   integer, intent(in):: klast     ! ending vertical index for MPI   integer, intent(in):: klastp    ! klast, except km+1 when klast=km   integer, intent(in):: nq        ! total # of tracers to be advected   integer, intent(in):: nc        ! declared dimension of q3   integer, intent(in):: ndt       ! the large time step in seconds                                   ! Also the mapping time step in this setup   integer, intent(in):: iord      ! parameter controlling monotonicity in E-W                                   ! recommendation: iord=4   integer, intent(in):: jord      ! parameter controlling monotonicity in N-S                                    ! recommendation: jord=4   integer, intent(in):: kord      ! parameter controlling monotonicity in mapping   integer, intent(in):: ifirstxy, ilastxy, jfirstxy, jlastxy  ! xy decomposition!                                    index ranges                                   ! recommendation: kord=4   real(r8), intent(in):: phis(im,jfirst:jlast)   ! surface geopotential (grav*zs)   real(r8), intent(in):: om       ! angular velocity of earth's rotation     real(r8), intent(in):: cp       ! heat capacity of air at constant pressure   real(r8), intent(in):: ae       ! radius of the earth (m)   real(r8), intent(in):: rair     ! Gas constant of the air   real(r8), intent(in):: ptop     ! Pressure at model top (interface pres)   logical,  intent(in):: convt    ! Flag to control pt/tvm output                                   ! If true: output tvm, the virtual temperature                                   ! If false: pt is updated   logical,  intent(in):: dcaf     ! Dry convection flag (use only in ideal_phys case)   logical,  intent(in):: rayf     ! Rayleigh friction flag (off by default)   logical,  intent(in):: ideal    ! True for ideal physics   real(r8), intent(in):: phisxy(ifirstxy:ilastxy,jfirstxy:jlastxy)  ! xy-decomposed version of phis! !INPUT/OUTPUT PARAMETERS:   integer, intent(inout):: ns0    ! total number of splits for Lagrangian                                   ! dynamics; a suitable value will be automatically                                   ! determined if ns0 = 0 (i.e., if you don't know what value                                   ! to be used try ns0=0   real(r8), intent(inout):: pt(im,jfirst-ng_d:jlast+ng_d,kfirst:klast)                            ! scaled (virtual) potential temperature   real(r8), intent(inout):: ps(im,jfirst:jlast)                ! Surface pressure (pa)    real(r8), intent(inout):: delp(im,jfirst:jlast,kfirst:klast) ! Pressure thickness   real(r8), intent(inout):: u(im,jfirst-ng_d:jlast+ng_s,kfirst:klast)    ! u wind velocities, staggered grid   real(r8), intent(inout):: v(im,jfirst-ng_s:jlast+ng_d,kfirst:klast)    ! v wind velocities, staggered grid   real(r8), intent(inout):: q3(im,jfirst-ng_d:jlast+ng_d,kfirst:klast,nc) ! Tracers!--------------------------------------------------------------------------------------! The following three arrays must be pre-computed as input to benergy(). They are NOT! needed if consv=.F.; updated on output (to be used by physdrv)! Please refer to routine pkez on the algorithm for computing pkz! from pe and pk!--------------------------------------------------------------------------------------   real(r8), intent(inout):: pe(im,kfirst:klast+1,jfirst:jlast)    ! Pres at layer edges    real(r8), intent(inout):: pk(im,jfirst:jlast,kfirst:klast+1)    ! pe**cappa   real(r8), intent(inout):: pkz(im,jfirst:jlast,kfirst:klast)     ! finite-volume mean of pk!--------------------------------------------------------------------------------------! !OUTPUT (input values are not used):!--------------------------------------------------------------------------------------   real(r8), intent(out):: te0                   ! Total energy before dynamics   real(r8), intent(out):: peln(im,kfirst:klast+1,jfirst:jlast) ! log pressure (pe) at layer edges   real(r8), intent(out):: omga(im,kfirst:klast,jfirst:jlast)   ! vertical pressure velocity (pa/sec)   real(r8), intent(out):: tvm(im,kfirst:klast,jfirst:jlast)   ! virtual temperature if convt! The following variables are xy-decomposed instanciations:   real(r8), intent(out):: uxy(ifirstxy:ilastxy,jfirstxy:jlastxy+1,km)  ! ghosted N1   real(r8), intent(out):: vxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km)   real(r8), intent(out):: psxy(ifirstxy:ilastxy,jfirstxy:jlastxy)   real(r8), intent(out):: delpxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km)   real(r8), intent(out):: pexy(ifirstxy:ilastxy,km+1,jfirstxy:jlastxy)   real(r8), intent(out):: pkxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km+1)   real(r8), intent(out):: pkzxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km)   real(r8), intent(out):: pelnxy(ifirstxy:ilastxy,km+1,jfirstxy:jlastxy)   real(r8), intent(out):: omgaxy(ifirstxy:ilastxy,km,jfirstxy:jlastxy)   real(r8), intent(out):: q3xy(ifirstxy:ilastxy,jfirstxy:jlastxy,km,nc)   real(r8), intent(out):: ptxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km)   real(r8), intent(out):: tvmxy(ifirstxy:ilastxy,km,jfirstxy:jlastxy)! !DESCRIPTION:!! Developer: Shian-Jiann Lin, NASA/GSFC; email: lin@dao.gsfc.nasa.gov!! Top view of D-grid prognostatic variables: u, v, and delp (and other scalars)!!               u(i,j+1)!                 |!      v(i,j)---delp(i,j)---v(i+1,j)!                 |!               u(i,j)!! External routine required: the user needs to supply a subroutine to set up!                            "Eulerian vertical coordinate" for remapping purpose.!                             Currently this routine is named as set_eta()!                             In principle any terrian following vertical!                             coordinate can be used. The input to fvcore!                             need not be on the same vertical coordinate!                             as the output.!                             If SPMD is defined the Pilgrim communication!                             library developed by Will Sawyer will be needed.!! Remarks: values at poles for both u and v need not be defined; but values for!          all other scalars needed to be defined at both poles (as polar cap mean!          quantities). Tracer advection is done "off-line" using the!          large time step. Consistency is maintained by using the time accumulated!          Courant numbers and horizontal mass fluxes for the FFSL algorithm.!          The input "pt" can be either dry potential temperature!          defined as T/pkz (adiabatic case) or virtual potential temperature!          defined as T*/pkz (full phys case). IF convt is true, pt is not updated.!          Instead, tvm, virtual temperature is ouput.!          ipt is updated if convt is false.!          The user may set the value of nx to optimize the SMP performance!          The optimal valuse of nx depends on the total number of available!          shared memory CPUs per node (NS). Assuming the maximm MPI !          decomposition is used in the y-direction, set nx=1 if the!          NS <=4; nx=4 if NS=16.!! !REVISION HISTORY:!   SJL 99.04.13:  Initial SMP version delivered to Will Sawyer!   WS  99.10.03:  1D MPI completed and tested; !   WS  99.10.11:  Additional documentation!   WS  99.10.28:  benergy and te_map added; arrays pruned!   SJL 00.01.01:  SMP and MPI enhancements; documentation!   WS  00.07.13:  Changed PILGRIM API!   WS  00.08.28:  SPMD instead of MPI_ON!   AAM 00.08.10:  Add kfirst:klast!   WS  00.12.19:  phis now distr., LLNL2DModule initialized here!   WS  01.02.02:  bug fix: parsplit only called for FIRST time!   WS  01.04.09:  Added initialization of ghost regions!   WS  01.06.10:  Removed if(first) section; use module!   AAM 01.06.27:  Extract te_map call into separate routine!   AAM 01.07.13:  Get rid of dynpkg2; recombine te_map;!                  perform forward transposes for 2D decomposition!   WS  01.12.10:  Ghosted PT (changes benergy, cd_core, te_map, hswf)!   WS  02.05.03:  Allocated 2D arrays consistently to pass strict check!!EOP!-----------------------------------------------------------------------!BOC! Local variables   integer i, j, k, iq          ! Loop indicies   real(r8) umax                ! Maximum winds, m/s   parameter (umax = 300.0)   logical consv                ! Flag to force conservation of toal energy   parameter (consv = .true.)   integer    nx          ! # of split pieces in x-direction; for performance, the   parameter (nx = 4)     ! user may set nx=1 if there is NO shared memory multitasking   integer ipe, it   integer nsplit, n, n2   integer incount, outcount! Geometric arrays! Move the following 3D arrays to an initialization routine?   real(r8), allocatable :: worka(:,:,:),dp0(:,:,:),cx(:,:,:),cy(:,:,:)   real(r8), allocatable :: mfx(:,:,:), mfy(:,:,:)   real(r8), allocatable :: delpf(:,:,:), uc(:,:,:), vc(:,:,:)   real(r8), allocatable :: dwz(:,:,:), pkc(:,:,:), wz(:,:,:)   real(r8), allocatable :: dpt(:,:,:)   real(r8), allocatable :: pkcc(:,:,:), wzc(:,:,:)! The following variables are work arrays for xy=>yz transpose   real(r8), allocatable :: pkkp(:,:,:), wzkp(:,:,:), pekp(:,:,:)! The following variables are xy instanciations   real(r8), allocatable :: mfxxy(:,:,:), dp0xy(:,:,:), wzxy(:,:,:)! ps3 and psxy3 are dummy 3d variants of ps and psxy, resp.   real(r8), allocatable :: ps3(:,:,:), psxy3(:,:,:)   double precision zamda, zam5   logical first, fill   integer imh   real(r8) dt   real(r8) bdt   data first /.true./   data fill  /.true./              ! perform a simple filling algorithm                                        ! in case negatives were found! Allocate temporary work arrays! Change later to use pointers for SMP performance???! (prime candidates: uc, vc, delpf)      allocate( worka(im,jfirst:     jlast,     kfirst:klast) )      allocate(   dp0(im,jfirst:     jlast,     kfirst:klast) )      allocate(   mfx(im,jfirst:     jlast,     kfirst:klast) )      allocate(   mfy(im,jfirst:     jlast+1,   kfirst:klast) )      allocate(    cx(im,jfirst-ng_d:jlast+ng_d,kfirst:klast) )      allocate(    cy(im,jfirst:     jlast+1,   kfirst:klast) )      allocate( delpf(im,jfirst-ng_d:jlast+ng_d,kfirst:klast) )      allocate(    uc(im,jfirst-ng_d:jlast+ng_d,kfirst:klast) )      allocate(    vc(im,jfirst-2:   jlast+2,   kfirst:klast) )      allocate(   dpt(im,jfirst-1:   jlast+1,   kfirst:klast) )      allocate(   dwz(im,jfirst-1:    jlast,    kfirst:klast+1) )      allocate(   pkc(im,jfirst-1:   jlast+1,   kfirst:klast+1) )       allocate(    wz(im,jfirst-1:   jlast+1,   kfirst:klast+1) )      allocate(  pkcc(im,jfirst  :   jlast  ,   kfirst:klast+1) )       allocate(   wzc(im,jfirst  :   jlast  ,   kfirst:klast+1) ) !! WS: 02.05.03:  2D arrays must be allocated consistently in non-2D case!      allocate(pkkp(im,jfirst:jlast,kfirst:klastp))      allocate(wzkp(im,jfirst:jlast,kfirst:klastp))      allocate(pekp(im,kfirst:klastp,jfirst:jlast))      allocate(wzxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km+1))      allocate(mfxxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km))      allocate(dp0xy(ifirstxy:ilastxy,jfirstxy:jlastxy,km))      allocate(ps3(im,jfirst:jlast,kfirst:klast))      allocate(psxy3(ifirstxy:ilastxy,jfirstxy:jlastxy,km))! First touch pkc and wz??? (bufferpack is multitask in vertical but geop! computations are parallel in j-loop)   if ( km > 1 ) then         ! not shallow water equations      if( consv ) then! Compute globally integrated Total Energy (te0)         call t_startf('benergy')         call benergy(im,     jm,     km,     u,      v,         &                      pt,     ng_d,   ng_s,   delp,   pe,        &                      pk,     pkz,    phis,   cp,     te0,       &                      mfx,    dp0,    jfirst, jlast,  kfirst,    &                      klast,  klastp)         call t_stopf('benergy')       endif   endif! Determine splitting      bdt = ndt! Second level splitting      n2 = max ( 1, ns/4 )      nsplit = (ns+n2-1) / n2      dt = bdt / float(nsplit*n2)  do 2000 n=1, n2   if( nq > 0 ) then!$omp parallel do private(i, j, k)      do k=kfirst,klast         do j=jfirst,jlast            do i=1,im! Save initial delp field before the small-time-step! Initialize the CFL number accumulators: (cx, cy)! Initialize total mass fluxes: (mfx, mfy)               dp0(i,j,k) = delp(i,j,k)                cx(i,j,k) = 0.                cy(i,j,k) = 0.               mfx(i,j,k) = 0.               mfy(i,j,k) = 0.            enddo         enddo      enddo   endif   do it=1, nsplit      if(it == nsplit .and. n == n2) then         ipe = 1                     ! end of fvcore; output pe for te_map