te_map.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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
        enddo
1000  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)