#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       endif2000  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