dynpkg.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

      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')

   endif

2000  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

!EOC
end subroutine dynpkg
!-----------------------------------------------------------------------