#include <misc.h>!-----------------------------------------------------------------------!BOP! !ROUTINE: te_map --- Map vertical Lagrangian coordinates to normal grid!! !INTERFACE:   subroutine te_map(consv,   convt,  ps,      omga,    pe,         &                     delp,    pkz,    pk,      mdt,     im,         &                     jm,      km,     nx,      jfirst,  jlast,      &                     ng_d,    ngus,   ngun,    ngvs,    ngvn,       &                     ifirst,  ilast,   nq,         &                     u,       v,      pt,      q3,      hs,         &                     cp,      akap,   kord,    peln,    te0,        &                     te,      dz,     tvm,     nc )!! !USES:   use precision, only : r8   use dynamics_vars, only : acap, ks, ptop, cosp, ak, bk!  use fill_module   use pmgrid, only: npr_y, nprxy_y, nprxy_x, myid_y, myidxy_y,      &                     myidxy_x, twod_decomp, iam#if defined( SPMD )   use parutilitiesmodule, only: sumop,  parcollective   use spmd_dyn, only: comm_y, commxy_x, commxy_y   use mod_comm, only: bufferpack2d, bufferunpack2d,                    &                       bufferpack3d, bufferunpack3d,                    &                       mp_send, mp_recv, mp_barrier,                    &                       buff_s, buff_r#endif!-----------------------------------------------------------------------! #define OLDWAY!! Activate the above line to attain zero diff with earlier code version.! The new method is needed for answers to be independent of MPI configuration.!-----------------------------------------------------------------------   implicit none#if defined( SPMD )#define CPP_PRT_PREFIX  if(iam.eq.0)#else#define CPP_PRT_PREFIX#endif! !INPUT PARAMETERS:   logical consv                 ! flag to force TE conservation   logical convt                 ! flag to control pt output (see below)   integer mdt                   ! mapping time step (same as phys)   integer im, jm, km            ! x, y, z dimensions   integer nq                    ! number of tracers (including h2o)   integer nc                    !    integer ng_d                  ! number of ghost latitudes on D grid   integer ngus                  ! number of ghost latitudes for U south   integer ngun                  ! number of ghost latitudes for U north   integer ngvs                  ! number of ghost latitudes for V south   integer ngvn                  ! number of ghost latitudes for V north   integer nx                    ! number of SMP "decomposition" in x   integer jfirst, jlast         ! starting & ending latitude index   integer ifirst, ilast         ! starting & ending longitude index   real(r8) hs(ifirst:ilast,jfirst:jlast) ! surface geopotential   real(r8) cp   real(r8) te0 ! !INPUT/OUTPUT PARAMETERS:   real(r8) pk(ifirst:ilast,jfirst:jlast,km+1) ! pe to the kappa   real(r8) u(ifirst:ilast,jfirst-ngus:jlast+ngun,km)    ! u-wind (m/s)   real(r8) v(ifirst:ilast,jfirst-ngvs:jlast+ngvn,km)    ! v-wind (m/s)   real(r8) q3(ifirst:ilast,jfirst-ng_d:jlast+ng_d,km,nc)! tracers including specific humidity   real(r8) pe(ifirst:ilast,km+1,jfirst:jlast) ! pressure at layer edges   real(r8) ps(ifirst:ilast,jfirst:jlast)      ! surface pressure   real(r8) pt(ifirst:ilast,jfirst-ng_d:jlast+ng_d,km)   ! virtual potential temperature as input                                     ! Output: virtual temperature if convt is true                                     ! false: output is (virtual) potential temperature    real(r8)  te(ifirst:ilast,jfirst:jlast,km)  ! Work array (cache performance)   real(r8)  dz(ifirst:ilast,jfirst:jlast,km)     ! Work array (cache performance) ! !OUTPUT PARAMETERS:   real(r8) delp(ifirst:ilast,jfirst:jlast,km) ! pressure thickness   real(r8)  tvm(ifirst:ilast,km,jfirst:jlast)    !   real(r8) omga(ifirst:ilast,km,jfirst:jlast)    ! vertical press. velocity (pascal/sec)   real(r8) peln(ifirst:ilast,km+1,jfirst:jlast)  ! log(pe)   real(r8) pkz(ifirst:ilast,jfirst:jlast,km)     ! layer-mean pk for converting t to pt! !DESCRIPTION:!! !REVISION HISTORY:!! WS 99.05.19 : Replaced IMR, JMR, JNP and NL with IM, JM-1, JM and KM! WS 99.05.25 : Revised conversions with IMR and JMR; removed fvcore.h! WS 99.07.29 : Reduced computation region to jfirst:jlast! WS 99.07.30 : Tested concept with message concatenation of te_map calls! WS 99.10.01 : Documentation; indentation; cleaning! SJL 99.12.31: SMP "decomposition" in-E-W direction! WS 00.05.14 : Renamed ghost indices as per Kevin's definitions! WS 00.07.13 : Changed PILGRIM API! AM 00.08.29 : Variables in this routine will ultimately be decomposed in (i,j).! AM 01.06.13 : 2-D decomposition; reordering summation causes roundoff difference.! WS 01.06.10 : Removed "if(first)" section in favor of a variable module! AM 01.06.27 : Merged yz decomposition technology into ccm code.! WS 02.01.14 : Upgraded to mod_comm!!EOP!-----------------------------------------------------------------------!BOC! Local arrays:      real(r8) rmin(nx*jm), rmax(nx*jm)      real(r8) tte(jm)! x-y      real(r8)  u2(ifirst:ilast,jfirst:jlast+1)      real(r8)  v2(ifirst:ilast+1,jfirst:jlast)      real(r8)  t2(ifirst:ilast,jfirst:jlast)      real(r8)  veast(jfirst:jlast,km)      real(r8)  pewest(km+1,jfirst:jlast)! y-z      real(r8)  pe0(ifirst:ilast,km+1)      real(r8)  pe1(ifirst:ilast,km+1)      real(r8)  pe2(ifirst:ilast,km+1)      real(r8)  pe3(ifirst:ilast,km+1)      real(r8) phis(ifirst:ilast,km+1)      real(r8) u2_sp(ifirst:ilast,km)      real(r8) v2_sp(ifirst:ilast,km)      real(r8) t2_sp(ifirst:ilast,km)      real(r8) u2_np(ifirst:ilast,km)      real(r8) v2_np(ifirst:ilast,km)      real(r8) t2_np(ifirst:ilast,km)      real(r8) dp2(ifirst:ilast,km)      real(r8) te2(ifirst:ilast,km)! x      real(r8)     gz(ifirst:ilast)      real(r8)  ratio(ifirst:ilast)      real(r8)    bte(ifirst:ilast)! z      real(r8) pe1w(km+1)      real(r8) pe2w(km+1)      integer i1w, nxu      integer i, j, k, ic, js2g0, jn2g0, jn1g1      integer kord      integer krd      real(r8) akap, dak, bkh, qmax, qmin      real(r8) te_sp(km), te_np(km)      real(r8) xysum(jfirst:jlast,2)      real(r8) tmpik(ifirst:ilast,km)      real(r8) tmpij(ifirst:ilast,jfirst:jlast,2)      real(r8) dtmp      real(r8) sum      real(r8) rdt5      real(r8) rg      real(r8) te1      real(r8) dlnp      integer ixj, jp, it, i1, i2#if defined( SPMD )      integer :: incount, outcount, incount2, outcount2, dest, src      real(r8), allocatable, save :: pesouth(:,:)#endif      integer comm_use, npry_use, itot      logical diag      logical first      data diag    /.false./      data first   /.true./! WS 99.07.27 :  Set loop limits appropriately      js2g0  = max(2,jfirst)      jn1g1  = min(jm,jlast+1)      jn2g0  = min(jm-1,jlast)      do j=jfirst,jlast         xysum(j,1) = 0.         xysum(j,2) = 0.      enddo      do j=jfirst,jlast         do i=ifirst,ilast            tmpij(i,j,1) = 0.            tmpij(i,j,2) = 0.         enddo      enddo      do k=1,km         do i=ifirst,ilast            tmpik(i,k) = 0.         enddo      enddo      itot = ilast-ifirst+1      nxu = 1      if (itot .eq. im) nxu = nx      if (first) then#if defined( SPMD )        if (.not. allocated(pesouth)) then          allocate( pesouth(ifirst:ilast,km+1) )        endif#endif        first = .false.      endif#if defined( SPMD )      if (twod_decomp .eq. 1) then         comm_use = commxy_y         npry_use = nprxy_y      else         comm_use = comm_y         npry_use = npr_y      endif#endif      call pkez(nxu, im, km, jfirst, jlast, 1, km, ifirst, ilast,      &                pe, pk, akap, ks, peln, pkz, .false.) #if defined( SPMD )      incount  = 0      outcount = 0! Send u southward      if ( jfirst > 1 ) then        call bufferpack3d( u,ifirst,ilast,jfirst-ngus,jlast+ngun,1,km,    &                              ifirst,ilast,jfirst,jfirst,1,km,buff_s )        incount  = itot*km      endif      if ( jlast < jm ) then        outcount = itot*km      endif      call mp_barrier()      call mp_send(iam-nprxy_x, iam+nprxy_x, incount, outcount, buff_s, buff_r)! Should split send/receive of pe for performance! Send pe northward      incount2 = 0      outcount2= 0      if ( jfirst > 1 ) then        outcount2 = itot*(km+1)      endif      if ( jlast .lt. jm ) then        call bufferpack3d( pe, ifirst,ilast,1,km+1,jfirst,jlast,ifirst,ilast,1,km+1,    &                              jlast,jlast,buff_s(incount+1) )        incount2 = itot*(km+1)      endif      call mp_send(iam+nprxy_x, iam-nprxy_x, incount2, outcount2, &                   buff_s(incount+1), buff_r(outcount+1) )      call mp_barrier()      call mp_recv(iam+nprxy_x, outcount, buff_r)      call mp_recv(iam-nprxy_x, outcount2, buff_r(outcount+1))      if ( jlast .lt. jm ) then        call bufferunpack3d( u,ifirst,ilast,jfirst-ngus,jlast+ngun,1,km,    &                                ifirst,ilast,jlast+1,jlast+1,1,km,buff_r )      endif      if ( jfirst > 1 ) then        call bufferunpack2d( pesouth, ifirst,ilast,1,km+1,     &                             ifirst,ilast,1,km+1,buff_r(outcount+1) )      endif#endif! Single subdomain case (periodic)      do k=1,km         do j=jfirst,jlast            veast(j,k) = v(ifirst,j,k)         enddo      enddo! Nontrivial x decomposition#if defined (SPMD)      if (itot .ne. im) then        call bufferpack3d(v,ifirst,ilast,jfirst-ngvs,jlast+ngvn,1,km,    &          ifirst,ifirst,jfirst,jlast,1,km,buff_s)        dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)        src  = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)        call mp_barrier()        call mp_send(dest,src,km*(jlast-jfirst+1),km*(jlast-jfirst+1),&                     buff_s,buff_r)        call mp_barrier()        call mp_recv(src,km*(jlast-jfirst+1),buff_r)        call bufferunpack2d( veast, jfirst, jlast, 1, km,      &                             jfirst, jlast, 1, km, buff_r )      endif#endif!$omp  parallel do        &!$omp  default(shared)    &!$omp  private(i,j, k, u2, v2, t2)! Compute cp*T + KE      do 1000 k=1,km        do j=js2g0,jn1g1           do i=ifirst,ilast              u2(i,j) = u(i,j,k)**2           enddo        enddo        do j=js2g0,jn2g0          do i=ifirst,ilast             v2(i,j) = v(i,j,k)**2          enddo          v2(ilast+1,j) = veast(j,k)**2        enddo        do j=jfirst,jlast           do i=ifirst,ilast              t2(i,j) = cp*pt(i,j,k)           enddo        enddo        do j=js2g0,jn2g0          do i=ifirst,ilast            te(i,j,k) = 0.25 * ( u2(i,j) + u2(i,j+1) +        &                                 v2(i,j) + v2(i+1,j) ) +      &                        t2(i,j)*pkz(i,j,k)          enddo        enddo! WS 99.07.29 : Restructuring creates small round-off.  Not clear why...! Do collective Mpisum (in i) for te_sp and te_np below (AAM)!        if ( jfirst == 1 ) then! South pole          do i=ifirst,ilast             u2_sp(i,k) = u2(i,2)             v2_sp(i,k) = v2(i,2)             t2_sp(i,k) = t2(i,1)          enddo        endif        if ( jlast == jm ) then! North pole          do i=ifirst,ilast             u2_np(i,k) = u2(i,jm)             v2_np(i,k) = v2(i,jm-1)             t2_np(i,k) = t2(i,jm)          enddo        endif! Compute dz; geo-potential increments        do j=jfirst,jlast           do i=ifirst,ilast              dz(i,j,k) = t2(i,j)*(pk(i,j,k+1)-pk(i,j,k))           enddo        enddo1000  continue      if ( jfirst .eq. 1 ) then!$omp  parallel do        &!$omp  default(shared)    &!$omp  private(i, k)         do k = 1, km            te_sp(k) = 0.            do i=ifirst,ilast#if defined (OLDWAY)              te_sp(k) = te_sp(k) + u2_sp(i,k) + v2_sp(i,k)#else              tmpik(i,k) = u2_sp(i,k) + v2_sp(i,k)              te_sp(k) = te_sp(k) + tmpik(i,k)