elseif(it == 1 .and. n == 1) then         ipe = -1                    ! start of cd_core      else         ipe = 0      endif! Call the Lagrangian dynamical core using small tme step      call t_startf('cd_core')      call cd_core(im,     jm,    km,     nq,     nx,              &                   jfirst, jlast, kfirst, klast,                   &                   klastp,  u,    v,      pt,     delp,            &                   pe,      pk,   dt,     ptop,   umax,            &                   ae,      rcap, cp,     cappa,  icd,             &                   jcd,     iord, jord,   ng_c,   ng_d,            &                   ng_s,    ipe,  om,     phis,                    &                   cx,      cy,   mfx,    mfy,    delpf,           &                   uc,      vc,   pkz,    dpt,    worka,           &                   dwz, pkc, wz, phisxy, ptxy, pkxy,               &                   pexy, pkcc, wzc, wzxy, delpxy, pkkp, wzkp,      &                   pekp, ifirstxy, ilastxy, jfirstxy, jlastxy)      call t_stopf('cd_core')   enddo   if( nq .ne. 0 ) then! Perform large-tme-step scalar transport using the accumulated CFL and! mass fluxes       call t_startf('trac2d')      call trac2d( dp0,    q3,     nc,     nq,     cx,             &                   cy,     mfx,    mfy,    iord,   jord,           &                   ng_d,   fill,   im,     jm,     km,             &                   jfirst, jlast,  kfirst, klast,  pkz,            &                   worka  )      call t_stopf('trac2d')   endif2000  continue#if defined (SPMD)   if (twod_decomp .eq. 1) then!! Transpose ps, u, v, and q3 from yz to xy decomposition!! Note: pt, pe and pk will have already been transposed through! call to geopk in cd_core. geopk does not actually require! secondary xy decomposition; direct 16-byte technique works just! as well, perhaps better. However, transpose method is used on last! call to avoid having to compute these three transposes now.!      call t_startf('transpose_fwd')! Embed ps in 3D array, per requirement of Pilgrim!$omp parallel do private(i,j,k)      do k = kfirst,klast         do j = jfirst,jlast            do i = 1,im               mfx(i,j,k) = ps(i,j)            enddo         enddo      enddo      call redistributestart (inter_ijk, .true., mfx)!! TEMPORARY!!!$omp parallel do private(i,j,k,iq)      do k = kfirst,klast         do j = jfirst,jlast            do i = 1,im               yzt(i,j,k) = u(i,j,k)            enddo         enddo      enddo      call redistributefinish(inter_ijk, .true., mfxxy)!$omp parallel do private(i,j)      do j = jfirstxy,jlastxy         do i = ifirstxy,ilastxy            psxy(i,j) = mfxxy(i,j,1)         enddo      enddo      call redistributestart (inter_ijk, .true., yzt) ! send U!! TEMPORARY!!!$omp parallel do private(i,j,k,iq)      do iq = 1,nc         do k = kfirst,klast            do j = jfirst,jlast               do i = 1,im                  q3t(i,j,k,iq) = q3(i,j,k,iq)               enddo            enddo         enddo      enddo      call redistributefinish(inter_ijk, .true., xyt) ! recv UXY      call redistributestart (inter_q3, .true., q3t)!! TEMPORARY!!!$omp parallel do private(i,j,k)      do k = 1,km         do j = jfirstxy,jlastxy            do i = ifirstxy,ilastxy               uxy(i,j,k) = xyt(i,j,k)            enddo         enddo      enddo!$omp parallel do private(i,j,k)      do k = kfirst,klast         do j = jfirst,jlast            do i = 1,im               yzt(i,j,k) = v(i,j,k)            enddo         enddo      enddo      call redistributefinish(inter_q3, .true., q3xy)      call redistributestart (inter_ijk, .true., yzt)  ! send V      call redistributefinish(inter_ijk, .true., vxy)  ! recv VXY      call t_stopf('transpose_fwd')    endif#endif    if ( km > 1 ) then           ! not shallow water equations! Perform vertical remapping from Lagrangian control-volume to! the Eulerian coordinate as specified by the routine set_eta.! Note that this finite-volume dycore is otherwise independent of the vertical! Eulerian coordinate.      call t_startf('te_map')      if (twod_decomp .eq. 1) then! ! te_map requires uxy, vxy, psxy, pexy, pkxy, phisxy, q3xy, and ptxy!         call te_map(consv,  convt,  psxy,  omgaxy, pexy,            &                     delpxy, pkzxy,  pkxy,  ndt,    im,              &                     jm,     km,     nx,    jfirstxy, jlastxy,       &                     0,      0,      1,     0,        0,             &                     ifirstxy, ilastxy,              &                     nq,     uxy,    vxy,   ptxy,   q3xy,            &                     phisxy, cp,     cappa, kord,   pelnxy,          &                     te0,    mfxxy,  dp0xy, tvmxy,  nc )!! te_map computes uxy, vxy, tvmxy, psxy, delpxy, pexy, pkxy, pkzxy,! pelnxy, omgaxy, q3xy and ptxy.!      else         call te_map(consv,  convt,  ps,    omga,   pe,              &                     delp,   pkz,    pk,    ndt,    im,              &                     jm,     km,     nx,    jfirst, jlast,           &                     ng_d,   ng_d,   ng_s,  ng_s,   ng_d,            &                     1,      im,                      &                     nq,     u,      v,     pt,     q3,              &                     phis,   cp,     cappa, kord,   peln,            &                     te0,    mfx,    dp0,   tvm,    nc )      endif      call t_stopf('te_map')    endif#if defined( SPMD )    if (twod_decomp .eq. 1) then       call t_startf('transpose_bck1')       if ( .not. convt ) then!! Transpose delpxy to delp for simplified physics (for full_phys,! delp is recomputed after physics advance)!          call redistributestart (inter_ijk, .false., delpxy)          call redistributefinish(inter_ijk, .false., delp)       endif!! Transpose pexy into pekp, then embed in pe and perform boundary update! (pexy is not needed for physics update)!       call redistributestart (inter_ikjp, .false., pexy)       call redistributefinish(inter_ikjp, .false., pekp)!$omp parallel do private(i,j,k)       do j = jfirst,jlast          do k = kfirst,klastp             do i = 1,im                pe(i,k,j) = pekp(i,k,j)             enddo          enddo       enddo       if (npr_z > 1) then          incount = 0          outcount = 0          if (kfirst > 1) then             call bufferpack3d(pe, 1, im, kfirst, klast+1, jfirst, jlast,  &                               1, im, kfirst, kfirst, jfirst, jlast, buff_s )             incount = im * (jlast-jfirst+1)          endif          if (klast < km) then             outcount = im * (jlast-jfirst+1)          endif          call mp_barrier()          call mp_send(iam-npr_y,iam+npr_y,incount,outcount,buff_s,buff_r)          call mp_barrier()          call mp_recv(iam+npr_y,outcount,buff_r)          if (klast < km) then             call bufferunpack3d(pe,1,im,kfirst,klast+1,jfirst,jlast,     &                                  1,im,klast+1,klast+1,jfirst,jlast,buff_r)          endif       endif!! Transpose psxy into ps, using 3D temporary arrays! (psxy is not needed for physics update)!       do k=1,km          do j=jfirstxy,jlastxy             do i=ifirstxy,ilastxy                psxy3(i,j,k) = psxy(i,j)             enddo          enddo       enddo       call redistributestart (inter_ijk, .false., psxy3)       call redistributefinish(inter_ijk, .false., ps3)       do j=jfirst,jlast          do i=1,im             ps(i,j) = ps3(i,j,kfirst)          enddo       enddo       call t_stopf('transpose_bck1')    endif#endif    deallocate( mfy )    deallocate( mfx )    deallocate(  cy )    deallocate(  cx )    deallocate( dp0 )    deallocate( delpf )    deallocate( uc    )    deallocate( vc    )    deallocate( dpt   )    deallocate( pkc   )    deallocate( dwz   )    deallocate(  wz   )    deallocate( worka )    deallocate( pkcc )    deallocate( wzc )    deallocate( pkkp )    deallocate( wzkp )    deallocate( pekp )    deallocate( wzxy )    deallocate( mfxxy )    deallocate( dp0xy )    deallocate( ps3 )    deallocate( psxy3 )!----------------------------------------------------------! Idealized physics: do Held-Suarez-Williamson-Lin forcing.! Since actual variable names depend on whether we are using! 2D decomposition, branching is required.!----------------------------------------------------------    if (ideal) then       call t_startf('ideal_phys')       if (twod_decomp .eq. 1) then!--------------------------------------------------------------------------! For 2D decomposition, hswf requires u3sxy, v3sxy, ptxy, pexy and ! pkzxy, and computes u3sxy, v3sxy and ptxy.!--------------------------------------------------------------------------          call hswf( im,   jm,   km,   jfirstxy, jlastxy,         &                     ifirstxy,         ilastxy,                   &                     uxy,  vxy,  ptxy, 0, 0, 1, 0, 0,             &                     pexy,     pkzxy,           &                     ndt,  cappa,      gravit,   rair,    dcaf,   &                    .true.,      rayf, sinp,     cosp,    sine,   &                     cose, coslon,     sinlon )       else          call hswf( im,   jm,   km,   jfirst,   jlast,           &                     1,    im,   u,    v,        pt,              &                     ng_d, ng_d, ng_s, ng_s,     ng_d,            &                     pe,   pkz,                                   &                     ndt,  cappa,      gravit,   rair,    dcaf,   &                    .true.,      rayf, sinp,     cosp,    sine,   &                     cose, coslon,     sinlon )       endif       call t_stopf('ideal_phys')    endif!EOCend subroutine dynpkg!-----------------------------------------------------------------------