benergy.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
!-----------------------------------------------------------------------
!BOP
! !ROUTINE: benergy --- Calculate the total energy
!
! !INTERFACE:

      subroutine benergy(im, jm, km, u, v, pt, ng_d, ng_s, delp, pe, pk, pkz, hs, &
                        cp,  te0, te, dz, jfirst, jlast,              &
                        kfirst, klast, klastp)

! !USES:

      use precision, only : r8
      use dynamics_vars, only : cosp, acap
      use pmgrid, only : npr_z, myid_z, npr_y, myid_y, strip3zaty, strip3zatypt
#if defined( SPMD )
      use spmd_dyn, only : comm_z, comm_y
      use mod_comm, only: mp_send_s, mp_recv_n, mp_barrier
      use parutilitiesmodule, only : sumop, pargatherreal, parcollective
#endif
      implicit none

! !INPUT PARAMETERS:
      integer im, jm, km, jfirst, jlast                  ! Dimensions
      integer kfirst, klast, klastp
      integer ng_d                                       ! Ghosting zone of pt
      integer ng_s                                       ! max(ng_c+1,ng_d)
      real (r8)  u(im, jfirst-ng_d:jlast+ng_s, kfirst:klast)   ! Winds x
      real (r8)  v(im, jfirst-ng_s:jlast+ng_d, kfirst:klast)   ! Winds y
      real (r8) pt(im, jfirst-ng_d:jlast+ng_d, kfirst:klast)   ! Potential temperature
      real (r8) delp(im, jfirst:jlast, kfirst:klast)           ! Delta pressure
      real (r8) pkz(im, jfirst:jlast, kfirst:klast)
      real (r8) pe(im, kfirst:klast+1, jfirst:jlast)           ! Edge pressure
      real (r8) pk(im, jfirst:jlast, kfirst:klast+1)
      real (r8) hs(im,jfirst:jlast)
      real (r8) cp

! !INPUT work arrays:
      real (r8) te(im, jfirst:jlast, kfirst:klast)             ! Work array (cache perf.)
      real (r8) dz(im, jfirst:jlast, kfirst:klast)             ! Work array (cache perf.)

! !OUTPUT PARAMETERS:
      real (r8) te0              ! globally integrated total energy

! !DESCRIPTION:
!    Determines the globally integrated total energy, if jfirst == 1
!    and jlast == jnp, otherwise it calculates the total energy in
!    the latitude slice jfirst:jlast. 
!
! !REVISION HISTORY:
!
! SJL 99.04.13 : Delivered as release 0.9.8
! WS  99.05.18 : Added im, jm, km, te, dz as arguments
! WS  99.05.25 : Replaced IMR by IM, JMR by JM-1; removed fvcore.h
! WS  99.10.11 : Ghosted U, now fully limited to jfirst:jlast
! WS  99.11.23 : Pruned te, additional cleaning
! WS  00.05.14 : Renamed ghost indices as per Kevin's definitions
! WS  00.07.13 : Changed PILGRIM API
! WS  00.08.28 : Cosmetic changes
! AAM 00.08.28 : Added kfirst,klast
! WS  00.12.01 : Replaced MPI_ON with SPMD; hs now distributed
! AAM 01.06.15 : Changes for zero diff
! WS  01.12.10 : Ghosted PT
! WS  01.12.31 : Ghosted U,V
!
!EOP
!---------------------------------------------------------------------
!BOC

! Local
      real (r8) u2(im,jfirst:jlast+1)
      real (r8) v2(im,jfirst:jlast)
      real (r8) tte(jfirst:jlast)
      real (r8) bte(im)
      real (r8) tteglob(jm)
      real (r8) te_sp
      real (r8) te_np
      real (r8) gztop(im,jfirst:jlast)

      integer i, j, k, js1g0, js2g0, jn1g0, jn1g1, jn2g0, ktot, jtot
      real (r8) xsum

! Geometric arrays

#if !defined( SPMD )
      integer comm_y
#endif
      real (r8), allocatable :: dzte_tmp(:,:,:)
      real (r8), allocatable :: pez(:,:,:), pey(:,:,:)


! WS 99.07.27 :  Set loop limits appropriately

      js1g0  = max(1,jfirst)
      js2g0  = max(2,jfirst)
      jn2g0  = min(jm-1,jlast)
      jn1g0  = min(jm,jlast)
      jn1g1  = min(jm,jlast+1)
      ktot = klast-kfirst+1
      jtot = jlast-jfirst+1
 
#if defined(SPMD)
!
! Send one latitude to south; ghost U one latitude at jlast+1
!
      call mp_send_s(im, jm, jfirst, jlast, kfirst, klast, &
                     ng_d, ng_s, u)
      call mp_barrier()
      call mp_recv_n(im, jm, jfirst, jlast, kfirst, klast, &
                     ng_d, ng_s, u)
#endif

!$omp  parallel do         &
!$omp  default(shared)     &
!$omp  private(i,j,k, u2, v2,te_sp,te_np)

      do k=kfirst,klast
        do j=js2g0,jn1g1
           do i=1,im
              u2(i,j) = u(i,j,k)**2
           enddo
        enddo

        do j=js2g0,jn2g0
           do i=1,im
              v2(i,j) = v(i,j,k)**2
           enddo
        enddo

        do j=js2g0,jn2g0
           do i=1,im-1
           te(i,j,k) = 0.25 * ( u2(i,j) + u2(i,j+1) +         &
                                v2(i,j) + v2(i+1,j)  )
           enddo
! i=im
           te(im,j,k) = 0.25 * ( u2(im,j) + u2(im,j+1) +      &
                                 v2(im,j) + v2(1,j)  )
        enddo


        do j=js2g0,jn2g0
           do i=1,im
              te(i,j,k) = delp(i,j,k) * ( te(i,j,k) +         &
                             cp*pt(i,j,k)*pkz(i,j,k)  ) 
           enddo
        enddo


! poles

        if ( jfirst .EQ. 1 ) then
          te_sp = 0.
          do i=1,im
            te_sp = te_sp + u2(i,  2) + v2(i,  2)
          enddo
          te_sp =   delp(1,  1,k) * (0.5*te_sp/float(im) +    &
                    cp*pt(1, 1,k)*pkz(1,1,k)      )
          do i=1,im
            te(i,  1,k) = te_sp
          enddo
        endif
  
        if ( jlast .EQ. jm ) then
          te_np = 0.
          do i=1,im
            te_np = te_np + u2(i,jm) + v2(i,jm-1)
          enddo
          te_np =   delp(1,jm,k) * (0.5*te_np/float(im) +     &
                    cp*pt(1,jm,k)*pkz(1,jm,k)   )
          do i=1,im
            te(i,jm,k) = te_np
          enddo
        endif

! Compute dz
        do j=js1g0,jn1g0
           do i=1,im
              dz(i,j,k) = cp*pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))
           enddo
        enddo

      enddo

! Perform vertical integration

      allocate(dzte_tmp(im,jfirst:jlast,km))
      allocate(pez(im,km+1,jfirst:jlast))
      allocate(pey(im,kfirst:klastp,jfirst:jlast))

! Compute gztop

#if defined( SPMD )
      if (myid_z .eq. 0) then
#endif
!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k)
        do k = kfirst, klast
         do j = jfirst, jlast
          do i = 1, im
           dzte_tmp(i,j,k) = dz(i,j,k)
          enddo
         enddo
        enddo
#if defined( SPMD )
      endif
#endif
#if defined( SPMD )
      if (npr_z .gt. 1) then
        call pargatherreal(comm_z, 0, dz, strip3zaty, dzte_tmp)
      endif
#endif

#if defined( SPMD )
      if (myid_z .eq. 0) then
#endif

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k)
        do j = jfirst, jlast
          do i=1,im
            gztop(i,j) = hs(i,j)
          enddo
          do k=1,km
            do i=1,im
              gztop(i,j) = gztop(i,j) + dzte_tmp(i,j,k)
            enddo
          enddo
        enddo

#if defined( SPMD )
      endif
#endif

#if defined( SPMD )
      if (myid_z .eq. 0) then
#endif
!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k)
        do j = jfirst, jlast
         do k = kfirst, klastp
          do i = 1, im
           pez(i,k,j) = pe(i,k,j)
          enddo
         enddo
        enddo
#if defined( SPMD )
      endif
#endif
#if defined( SPMD )
      if (npr_z .gt. 1) then
!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k)
        do j = jfirst, jlast
         do k = kfirst, klastp
          do i = 1, im
           pey(i,k,j) = pe(i,k,j)
          enddo
         enddo
        enddo
        call pargatherreal(comm_z, 0, pey, strip3zatypt, pez)
      endif
#endif

#if defined( SPMD )
      if (myid_z .eq. 0) then
#endif
!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k)
        do k = kfirst, klast
         do j = jfirst, jlast
          do i = 1, im
           dzte_tmp(i,j,k) = te(i,j,k)
          enddo
         enddo
        enddo
#if defined( SPMD )
      endif
#endif
#if defined( SPMD )
      if (npr_z .gt. 1) then
        call pargatherreal(comm_z, 0, te, strip3zaty, dzte_tmp)
      endif
#endif

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

      do j=js1g0, jn1g0
         tte(j) = 0.
      enddo

#if defined( SPMD )
      if (myid_z .eq. 0) then
#endif

!$omp  parallel do       &
!$omp  default(shared)   &
!$omp  private(i,j,k,bte)
         do j=js1g0, jn1g0
            tte(j) = 0.
            if  (j .eq. 1) then
               tte(1) = tte(1) - pez(1,1,1)*gztop(1,1)
               tte(1) = tte(1) + pez(1,km+1,1)*hs(1,1)
               do k=1,km
                  tte(1) = tte(1) + dzte_tmp(1,1,k)
               enddo
               tte(1) = acap * tte(1)
            elseif (j .eq. jm) then
               tte(jm) = tte(jm) - pez(1,1,jm)*gztop(1,jm)
               tte(jm) = tte(jm) + pez(1,km+1,jm)*hs(1,jm)
               do k=1,km
                  tte(jm) = tte(jm) + dzte_tmp(1,jm,k)
               enddo
               tte(jm) = acap * tte(jm)
            else
               do i = 1,im
                  bte(i) = 0.
               enddo
               do i = 1,im
                  bte(i) = bte(i) - pez(i,1,j)*gztop(i,j)
                  bte(i) = bte(i) + pez(i,km+1,j)*hs(i,j)
                  do k=1,km
                     bte(i) = bte(i) + dzte_tmp(i,j,k)
                  enddo
               enddo
               do i = 1,im
                  tte(j) = tte(j) + bte(i)
               enddo
               tte(j) = tte(j) * cosp(j)
            endif
         enddo
#if defined( SPMD )
      endif
#endif

#if defined( SPMD )
      if (npr_z .gt. 1) then
        call parcollective(comm_z, sumop, jtot, tte)
      endif
#endif

      deallocate(dzte_tmp)
      deallocate(pez)
      deallocate(pey)

! Sum over j

      te0=0.
      do j = jfirst, jlast
        tteglob(j) = tte(j)
      enddo
      call par_vecsum(jm, jfirst, jlast, tteglob, te0, comm_y, npr_y)

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