te_map.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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) = 0.5*te_sp(k)/float(im) + t2_sp(ifirst,k)*pkz(ifirst,1,k)
            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) = 0.5*te_np(k)/float(im) + t2_np(ifirst,k)*pkz(ifirst,jm,k)
            do i=ifirst,ilast
              te(i,jm,k) = te_np(k)
            enddo
         enddo
      endif

#if defined( SPMD )
      if (itot .ne. im) then
        call bufferpack3d(pe,ifirst,ilast,1,km+1,jfirst,jlast,     &
          ilast,ilast,1,km+1,jfirst,jlast,buff_s)

        dest = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
        src  = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
        call mp_barrier()
        call mp_send(dest,src,(km+1)*(jlast-jfirst+1),          &
                     (km+1)*(jlast-jfirst+1),buff_s,buff_r)
        call mp_barrier()
        call mp_recv(src,(km+1)*(jlast-jfirst+1),buff_r)

        call bufferunpack2d(pewest, 1, km+1, jfirst, jlast,        &
                            1, km+1, jfirst, jlast, buff_r)
      endif
#endif

      it = itot / nxu
      jp = nxu * ( jlast - jfirst + 1 )

!$omp  parallel do           &
!$omp  default(shared)       &
!$omp  private(i,j,k,ic,i1w,pe0,pe1,pe2,pe3,ratio)   &
!$omp  private(dak,bkh,rdt5,phis,krd, ixj,i1,i2) &
!$omp  private(te2, dp2, pe1w, pe2w )

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

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

! Copy data to local 2D arrays.
        i1w = i1-1
        if (i1 .eq. 1) i1w = im
        do k=1,km+1
           do i=i1,i2
              pe1(i,k) = pe(i,k,j)
           enddo
           if( itot == im ) then
               pe1w(k) = pe(i1w,k,j)
           else
               pe1w(k) = pewest(k,j)
           endif
        enddo

        do k=1,ks+1
           do i=i1,i2
              pe0(i,k) = ak(k)
              pe2(i,k) = ak(k)
              pe3(i,k) = ak(k)
            enddo
        enddo

        do k=ks+2,km
           do i=i1,i2
              pe0(i,k) = ak(k) + bk(k)* ps(i,j)
              pe2(i,k) = ak(k) + bk(k)*pe1(i,km+1)
           enddo
        enddo

        do i=i1,i2
           pe0(i,km+1) =  ps(i,j)
           pe2(i,km+1) = pe1(i,km+1)
        enddo

! Ghosting for v mapping
        do k=ks+2,km
           pe2w(k) = ak(k) + bk(k)*pe1w(km+1)
        enddo
        pe2w(km+1) = pe1w(km+1)

! Compute omga (dp/dt)
        rdt5 = 0.5 / float(mdt)
        do k=2,km+1
           do i=i1,i2
              pe0(i,k) = pe1(i,k) - pe0(i,k)
           enddo
        enddo

        do i=i1,i2
! update ps
          ps(i,j)   = pe1(i,km+1)
          omga(i,1,j) = rdt5 * pe0(i,2)
        enddo

        do k=2,km
          do i=i1,i2
             omga(i,k,j) = rdt5 * ( pe0(i,k) + pe0(i,k+1) )
          enddo
        enddo

        if(ks .ne. 0) then
           do k=1,ks
             dak = ak(k+1) - ak(k)
             do i=i1,i2
                delp(i,j,k) = dak
             enddo
           enddo
        endif

        do k=ks+1,km
          do i=i1,i2
             delp(i,j,k) = pe2(i,k+1) - pe2(i,k)
          enddo
        enddo

! Compute correction terms to Total Energy
        do i=i1,i2
           phis(i,km+1) = hs(i,j)      
        enddo

        do k=km,1,-1
          do i=i1,i2
             phis(i,k) = phis(i,k+1) + dz(i,j,k)   
          enddo
        enddo

        do k=1,km+1
          do i=i1,i2
             phis(i,k) = phis(i,k) * pe1(i,k)
          enddo
        enddo

! <<< Compute Total Energy >>>
        do k=1,km
          do i=i1,i2
            dp2(i,k) = pe2(i,k+1) - pe2(i,k)
            te2(i,k) =  te(i,j,k)+(phis(i,k+1)-phis(i,k))/(pe1(i,k+1)-pe1(i,k))
          enddo
        enddo

! Map Total Energy
        call map3_ppm (km,   pe1,   te2,                &
                       km,   pe2,   te2,  dp2,          &
                       itot, i1-ifirst+1, i2-ifirst+1,  &
                       1,    kord )
        do k=1,km
           do i=i1,i2
              te(i,j,k) =  te2(i,k)
           enddo
        enddo

! Map constituents

       if( nq .ne. 0 ) then
          if(kord == 8) then
             krd = 8
          else
             krd = 7
          endif

          do ic=1,nq
!
! Review this code!
!
             do k=1,km
                do i=i1,i2
                   te2(i,k) = q3(i,j,k,ic)
                enddo
             enddo
             call map3_ppm (km,   pe1,   te2,          &
                            km,   pe2,   te2,  dp2,    &
                            itot, i1-ifirst+1,       i2-ifirst+1,    &
                            0,    krd )
             do k=1,km
                do i=i1,i2
                   q3(i,j,k,ic) = te2(i,k)
                enddo
             enddo
          enddo

! Ensure that there is absolutely no negatives; this should only be useful
! for during initialization or from a badly interpolated IC
!         call fillz( itot, i1-ifirst+1, i2-ifirst+1, km, nq,       &
!                     q3(ifirst,j,1,1), dp2 )

       endif

! map u
        if(j .ne. 1) then

! WS 99.07.29 : protect j==jfirst case
          if (j > jfirst) then
            do k=2,km+1
              do i=i1,i2
                pe0(i,k) = 0.5*(pe1(i,k)+pe(i,k,j-1))
              enddo
            enddo

            do k=ks+2,km+1
              bkh = 0.5*bk(k)
              do i=i1,i2
                pe3(i,k) = ak(k) + bkh*(pe1(i,km+1)+pe(i,km+1,j-1))
              enddo
            enddo

#if defined( SPMD )
          else
!  WS 99.10.01 : Read in pe(:,:,jfirst-1) from the pesouth buffer
            do k=2,km+1
              do i=i1,i2
                pe0(i,k) = 0.5*(pe1(i,k)+pesouth(i,k))
              enddo
            enddo

            do k=ks+2,km+1
              bkh = 0.5*bk(k)
              do i=i1,i2
                pe3(i,k) = ak(k) + bkh*(pe1(i,km+1)+pesouth(i,km+1))
              enddo
            enddo
#endif
          endif

          call map1_ppm ( km,   pe0,    u,                     &
                          km,   pe3,    u,                     &
                          itot, i1-ifirst+1, i2-ifirst+1,      &
                          j,    jfirst, jlast,  ngus, ngun,    &
                          -1,    kord)
        endif

! map v
        if(j .ne. 1 .and. j .ne. jm) then
          do k=2,km+1
! pe1(i1-1,1:km+1) must be ghosted
            pe0(i1,k) = 0.5*(pe1(i1,k)+pe1w(k))
            do i=i1+1,i2
               pe0(i ,k) = 0.5*(pe1(i,k)+pe1(i-1,k))
            enddo
          enddo

          do k=ks+2,km+1
! pe2(i1-1,ks+2:km+1) must be ghosted
            pe3(i1,k) = 0.5*(pe2(i1,k)+pe2w(k))
            do i=i1+1,i2
               pe3(i,k) = 0.5*(pe2(i,k)+pe2(i-1,k))
            enddo
          enddo

          call map1_ppm ( km,   pe0,    v,                      &
                          km,   pe3,    v,                      &
                          itot, i1-ifirst+1, i2-ifirst+1,       &
                          j,    jfirst, jlast,  ngvs, ngvn,     &
                          -1,    kord)
        endif

! Save new PE to temp storage peln
        do k=2,km
          do i=i1,i2
             peln(i,k,j) = pe2(i,k)
          enddo
        enddo

! Check deformation.
       if( diag ) then
          rmax(ixj) = 0.
          rmin(ixj) = 1.
          do k=1,km
             do i=i1,i2
              ratio(i) = (pe1(i,k+1)-pe1(i,k)) / (pe2(i,k+1)-pe2(i,k))
             enddo

             do i=i1,i2
              if(ratio(i) > rmax(ixj)) then
                 rmax(ixj) = ratio(i)
              elseif(ratio(i) < rmin(ixj)) then
                 rmin(ixj) = ratio(i)
              endif
            enddo
          enddo
       endif
2000  continue

      if( diag ) then
        qmin = rmin(1)
        do ixj=2, jp
          if(rmin(ixj) < qmin) then
            qmin = rmin(ixj)
          endif
        enddo
        CPP_PRT_PREFIX write(6,*) 'rmin=', qmin

        qmax = rmax(1)
        do ixj=2, jp
          if(rmax(ixj) > qmax) then
            qmax = rmax(ixj)
          endif
        enddo
        CPP_PRT_PREFIX write(6,*) 'rmax=', qmax
      endif

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

      do j=jfirst,jlast
        do k=2,km
          do i=ifirst,ilast
            pe(i,k,j) = peln(i,k,j)
          enddo
        enddo
      enddo