te_map.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

      call pkez(nxu, im, km, jfirst, jlast, 1, km, ifirst, ilast,  &
                pe, pk, akap, ks, peln, pkz, .true.)

! ((((((((((((((((( compute globally integrated TE >>>>>>>>>>>>>>>>

      if( consv ) then

!$omp  parallel do         &
!$omp  default(shared)     &
!$omp  private(i,j,k)

        do k=1,km
          do j=jfirst,jlast
             do i=ifirst,ilast
                dz(i,j,k) = te(i,j,k) * delp(i,j,k)
             enddo
          enddo
        enddo

!$omp  parallel do        &
!$omp  default(shared)    &
!$omp  private(i,j,k,bte)

! Perform vertical integration

        do 4000 j=jfirst,jlast

          if ( j == 1 ) then
! SP
            tte(1) = 0.
  
            do k=1,km
              tte(1) = tte(1) + dz(ifirst,1,k)
            enddo

          elseif ( j .eq. jm) then
! NP
            tte(jm) = 0.

            do k=1,km
              tte(jm) = tte(jm) + dz(ifirst,jm,k)
            enddo

          else
! Interior
            do i=ifirst,ilast
              bte(i) = 0.
            enddo

            do k=1,km
              do i=ifirst,ilast
                bte(i) = bte(i) + dz(i,j,k)
              enddo
            enddo

            xysum(j,1) = 0.
            do i=ifirst,ilast
              xysum(j,1) = xysum(j,1) + bte(i)

#if !defined (OLDWAY)
              tmpij(i,j,1) = bte(i)
#endif

            enddo

          endif
4000    continue

#if defined (SPMD)
        if (nprxy_x .gt. 1) then
# if defined (OLDWAY)
          call parcollective(commxy_x, sumop, jlast-jfirst+1, 1, xysum)
# else
          call par_xsum(tmpij, ifirst, ilast, im, jlast-jfirst+1, xysum)
# endif
        endif
#endif

!$omp  parallel do        &
!$omp  default(shared)    &
!$omp  private(j)

        do j = max(jfirst,2), min(jlast,jm-1)
           tte(j) = xysum(j,1)*cosp(j)
        enddo

        if ( jfirst == 1 ) tte(1)  = acap * tte(1)
        if ( jlast == jm ) tte(jm) = acap * tte(jm)

        te1 = 0.
        call par_vecsum(jm, jfirst, jlast, tte, te1, comm_use, npry_use)

      endif   ! consv


#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)
      call mp_barrier()
      call mp_recv(iam+nprxy_x, outcount, buff_r)

      if ( jlast < jm ) then
        call bufferunpack3d( u,ifirst,ilast,jfirst-ngus,jlast+ngun,1,km,    &
                                ifirst,ilast,jlast+1,jlast+1,1,km,buff_r )
      endif
#endif

      if( consv ) then

!$omp  parallel do       &
!$omp& default(shared)   &
!$omp& private(i,j)
 
       do j=js2g0, jn2g0
          xysum(j,1) = 0.
          xysum(j,2) = 0.
        do i=ifirst,ilast
          xysum(j,1) = xysum(j,1) + ps(i,j)
          xysum(j,2) = xysum(j,2) + peln(i,km+1,j)
#if !defined (OLDWAY)
          tmpij(i,j,1) = ps(i,j)
          tmpij(i,j,2) = peln(i,km+1,j) 
#endif
        enddo
       enddo

#if defined( SPMD )
       if (nprxy_x .gt. 1) then
# if defined (OLDWAY)
          call parcollective(commxy_x, sumop, jlast-jfirst+1, 2, xysum)
# else
          call par_xsum(tmpij, ifirst, ilast, im, 2*(jlast-jfirst+1), xysum)
# endif
       endif
#endif

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(j)
 
       do j=js2g0, jn2g0
        tte(j) = cp*cosp(j)*(xysum(j,1) - ptop*float(im) -           &
                 akap*ptop*(xysum(j,2) - peln(ifirst,1,j)*float(im)) )
! peln(i,1,j) should be independent of i (AAM)
       enddo

       if ( jfirst .eq. 1 ) tte(1) = acap*cp * (ps(ifirst,1) - 2.*ptop -    &
                akap*ptop*(peln(ifirst,km+1,1) - peln(ifirst,1,1) ) )
       if ( jlast .eq. jm ) tte(jm)= acap*cp * (ps(ifirst,jm) -             &
                akap*ptop*(peln(ifirst,km+1,jm) - peln(ifirst,1,jm) ) )
      endif ! consv

      if (consv) then

       sum=0.
       call par_vecsum(jm, jfirst, jlast, tte, sum, comm_use, npry_use)

       dtmp = (te0 - te1) / sum
       if( diag ) then
         CPP_PRT_PREFIX write(6,*) 'te=',te0, ' Energy deficit in T = ', dtmp
       endif

      endif              ! end consv check

! Single x-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)

      do 8000 k=1,km
! Compute KE
        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=js2g0,jn2g0
          do i=ifirst,ilast
            te(i,j,k) = te(i,j,k) - 0.25 * ( u2(i,j) + u2(i,j+1)     &
                                            +v2(i,j) + v2(i+1,j) )
          enddo
        enddo

        if ( jfirst .eq. 1 ) then
! South pole
          do i=ifirst,ilast
            u2_sp(i,k) = u2(i,2)
            v2_sp(i,k) = v2(i,2)
          enddo
        endif

        if ( jlast .eq. jm ) then
! North pole
          do i=ifirst,ilast
            u2_np(i,k) = u2(i,jm)
            v2_np(i,k) = v2(i,jm-1)
          enddo
        endif

8000  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)
#endif
            enddo
         enddo

#if defined( SPMD )
         if (nprxy_x .gt. 1) then
# if defined (OLDWAY)
            call parcollective(commxy_x, sumop, km, te_sp)
# else
            call par_xsum(tmpik, ifirst, ilast, im, km, te_sp)
# endif
         endif
#endif

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i, k)

         do k = 1, km
            te_sp(k) = te(ifirst,1,k) - 0.5*te_sp(k)/float(im)
            do i=ifirst,ilast
              te(i,  1,k) = te_sp(k)
            enddo
         enddo
      endif

      if ( jlast .eq. jm ) then

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i, k)

         do k = 1, km
            te_np(k) = 0.
            do i=ifirst,ilast
#if defined (OLDWAY)
              te_np(k) = te_np(k) + u2_np(i,k) + v2_np(i,k)
#else
              tmpik(i,k) = u2_np(i,k) + v2_np(i,k)
              te_np(k) = te_np(k) + tmpik(i,k)
#endif
            enddo
         enddo

#if defined( SPMD )
         if (nprxy_x .gt. 1) then
# if defined (OLDWAY)
            call parcollective(commxy_x, sumop, km, te_np)
# else
            call par_xsum(tmpik, ifirst, ilast, im, km, te_np)
# endif
         endif
#endif

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i, k)

         do k = 1, km
            te_np(k) = te(ifirst,jm,k) - 0.5*te_np(k)/float(im)
            do i=ifirst,ilast
              te(i,jm,k) = te_np(k)
            enddo
         enddo
      endif

! Recover (virtual) temperature
!$omp  parallel do        &
!$omp  default(shared)    &
!$omp  private(ixj, i1, i2, i, j, k, rg, gz, dlnp)

!     do 9000 j=jfirst,jlast
      do 9000 ixj=1,jp

         j  = jfirst + (ixj-1) / nxu
         i1 = ifirst + it * mod(ixj-1, nxu)
         i2 = i1 + it - 1

         rg = akap * cp
         do i=i1,i2
            gz(i) = hs(i,j)      
         enddo

        do k=km,1,-1
          do i=i1,i2
            dlnp  = rg*(peln(i,k+1,j) - peln(i,k,j))
            tvm(i,k,j)  = delp(i,j,k)*(te(i,j,k) - gz(i)) /     &
                        ( cp*delp(i,j,k) - pe(i,k,j)*dlnp )
! Update phis
            gz(i) = gz(i) + dlnp*tvm(i,k,j)
          enddo

          if( consv ) then
              do i=i1,i2
                 tvm(i,k,j) = tvm(i,k,j) + dtmp
              enddo
          endif

          if( .not. convt ) then
              do i=i1,i2
                 pt(i,j,k) = tvm(i,k,j) / pkz(i,j,k)
              enddo
          endif
        enddo           ! end k-loop
9000  continue

      return
!EOC
      end
!-----------------------------------------------------------------------