cd_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
!-----------------------------------------------------------------------
!BOP
! !ROUTINE: cd_core --- Dynamical core for both C- and D-grid Lagrangian
!                       dynamics
!
! !INTERFACE:
 subroutine cd_core(im,   jm,   km,    nq,     nx,                 &
                    jfirst, jlast , kfirst, klast, klastp, u,      &
                    v,   pt,   delp,   pe,     pk,                 &
                    dt,   ptopin, umax,   ae,     rcap,            &
                    cp,   akap, iord_c, jord_c, iord_d, jord_d,    &
                    ng_c, ng_d, ng_s, ipe,    om,     hs,          &
                    cx3  ,  cy3, mfx, mfy,                         &
                    delpf, uc, vc, ptc, dpt, ptk,                  &
                    wz3, pkc, wz,  hsxy, ptxy, pkxy,               &
                    pexy, pkcc, wzc, wzxy, delpxy, pkkp, wzkp,     &
                    pekp, ifirstxy, ilastxy, jfirstxy, jlastxy )

! !USES:
   use precision, only : r8
   use sw_core, only : c_sw, d_sw
   use pft_module, only : fftfax, pft2d, pft_cf
   use dynamics_vars, only : sinp, cosp, cose, acosp,              &
                       sinlon, coslon, cosl5, sinl5,               &
                       ak, bk, ptop, ns, yzt
   use pmgrid, only: twod_decomp, myid_y, myid_z, npr_y, npr_z, iam

#if defined( SPMD )
   use redistributemodule, only: redistributestart, redistributefinish
   use spmd_dyn, only: comm_y, comm_z, inter_ijk, inter_ijkp

   use mod_comm, only: mp_send3d_ns, mp_recv3d_ns,           &
                       mp_send2_ns, mp_recv2_ns,             &
                       mp_send2_n, mp_recv2_s,               &
                       mp_send_s, mp_recv_n,                 &
                       mp_send, mp_recv, mp_barrier,         &
                       bufferpack3d, bufferunpack3d,         &
                       buff_s, buff_r

#define CPP_PRT_PREFIX  if(iam.eq.0)
#else
#define CPP_PRT_PREFIX
#endif

   implicit none

! !INPUT PARAMETERS:

! Input paraterers:
  integer im, jm, km
  integer nx                 ! # of split pieces in longitude direction
  integer jfirst
  integer jlast
  integer kfirst
  integer klast
  integer klastp             ! klast, except km+1 when klast=km
  integer ifirstxy,ilastxy   ! xy-decomp. longitude ranges
  integer jfirstxy,jlastxy   ! xy-decomp. latitude ranges
  integer ipe                ! ipe=1:  end of cd_core()
                             ! ipe=-1: start of cd_core()
                             ! ipe=0 :
  integer nq                 ! # of tracers to be advected by trac2d
  integer iord_c, jord_c     ! scheme order on C grid in X and Y dir.
  integer iord_d, jord_d     ! scheme order on D grid in X and Y dir.
  integer ng_c               ! ghost latitudes on C grid
  integer ng_d               ! ghost lats on D (Max NS dependencies, ng_d >= ng_c)
  integer ng_s               ! max(ng_c+1,ng_d) significant if ng_c = ng_d

  real(r8) ae                ! Radius of the Earth (m)
  real(r8) om                ! rotation rate
  real(r8) ptopin
  real(r8) umax
  real(r8) dt            !small time step in seconds
  real(r8) rcap
  real(r8) cp
  real(r8) akap

! Input time independent arrays:
  real(r8) hs(im,jfirst:jlast)         !surface geopotential
  real(r8) hsxy(ifirstxy:ilastxy,jfirstxy:jlastxy) !surface geopotential XY-decomp.

! !INPUT/OUTPUT PARAMETERS:
! Prognostic variables:
  real(r8)   u(im,jfirst-ng_d:jlast+ng_s,kfirst:klast)          ! Wind in X
  real(r8)   v(im,jfirst-ng_s:jlast+ng_d,kfirst:klast)          ! Wind in Y
  real(r8) delp(im,jfirst:jlast,kfirst:klast)          ! Delta pressure
  real(r8)   pt(im,jfirst-ng_d:jlast+ng_d,kfirst:klast)! Potential temperature

! Input/output: accumulated winds & mass fluxes on c-grid for large-
!               time-step transport
  real(r8) cx3(im,jfirst-ng_d:jlast+ng_d,kfirst:klast)! Accumulated Courant number in X
  real(r8) cy3(im,jfirst:jlast+1,kfirst:klast)        ! Accumulated Courant number in Y
  real(r8) mfx(im,jfirst:jlast,kfirst:klast)          ! Mass flux in X  (unghosted)
  real(r8) mfy(im,jfirst:jlast+1,kfirst:klast)        ! Mass flux in Y

! !OUTPUT PARAMETERS:
  real(r8) pe(im,kfirst:klast+1,jfirst:jlast)         ! Edge pressure
  real(r8) pk(im,jfirst:jlast,kfirst:klast+1)         ! Pressure to the kappa
  real(r8) ptxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km) ! Potential temperature XY decomp
  real(r8) pkxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km+1) ! P-to-the-kappa XY decomp
  real(r8) pexy(ifirstxy:ilastxy,km+1,jfirstxy:jlastxy) ! Edge pressure XY decomp

! Input work arrays:
  real(r8) delpf(im,jfirst-ng_d:jlast+ng_d,kfirst:klast)
  real(r8)   uc(im,jfirst-ng_d:jlast+ng_d,kfirst:klast)
  real(r8)   vc(im,jfirst-2:   jlast+2,   kfirst:klast)

  real(r8) ptc(im,jfirst:jlast,kfirst:klast)
  real(r8) ptk(im,jfirst:jlast,kfirst:klast)
  real(r8) dpt(im,jfirst-1:jlast+1,kfirst:klast)
  real(r8) wz3(im,jfirst-1:jlast  ,kfirst:klast+1)
  real(r8) pkc(im,jfirst-1:jlast+1,kfirst:klast+1) 
  real(r8)  wz(im,jfirst-1:jlast+1,kfirst:klast+1)
  real(r8) pkcc(im,jfirst:jlast,kfirst:klast+1) 
  real(r8) wzc(im,jfirst:jlast,kfirst:klast+1) 
  real(r8) wzxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km+1)
  real(r8) delpxy(ifirstxy:ilastxy,jfirstxy:jlastxy,km)
  real(r8) pkkp(im,jfirst:jlast,kfirst:klastp)
  real(r8) wzkp(im,jfirst:jlast,kfirst:klastp)
  real(r8) pekp(im,kfirst:klastp,jfirst:jlast)

#if defined ( HIGH_P )
  real(r8) wzz(im,jfirst:jlast+1,kfirst:klast+1)
#endif

! ! !DESCRIPTION:
!    Perform a dynamical update for one small time step; the small
!    time step is limitted by the fastest wave within the Lagrangian control-
!    volume 
!
! !REVISION HISTORY:
!     SJL  99.01.01:   Original SMP version
!     WS   99.04.13:   Added jfirst:jlast concept
!     SJL  99.07.15:   Merged c_core and d_core to this routine
!     WS   99.09.07:   Restructuring, cleaning, documentation
!     WS   99.10.18:   Walkthrough corrections; frozen for 1.0.7
!     WS   99.11.23:   Pruning of some 2-D arrays
!     SJL  99.12.23:   More comments; general optimization; reduction
!                      of redundant computation & communication
!     WS   00.05.14:   Modified ghost indices per Kevin's definition
!     WS   00.07.13:   Changed PILGRIM API
!     WS   00.08.28:   Cosmetic changes: removed old loop limit comments
!     AAM  00.08.30:   Introduced kfirst,klast
!     WS   00.12.01:   Replaced MPI_ON with SPMD; hs now distributed
!     WS   01.04.11:   PILGRIM optimizations for begin/endtransfer
!     WS   01.05.08:   Optimizations in the call of c_sw and d_sw
!     AAM  01.06.27:   Reinstituted 2D decomposition for use in ccm
!     WS   01.12.10:   Ghosted PT, code now uses mod_comm primitives
!     WS   01.12.31:   Removed vorticity damping, ghosted U,V,PT
!     WS   02.01.15:   Completed transition to mod_comm
!
!EOP
!---------------------------------------------------------------------
!BOC
! Local
      real(r8)   ub(im,jfirst:  jlast+1)
      real(r8)  wk1(im,jfirst-1:jlast+1)
      real(r8)  wk2(im,jfirst-ng_d:jlast+ng_d)
      real(r8)  wk3(im,jfirst-1:jlast+1)

      real(r8) p1d(im)

      real(r8) ptmp, pint, press
      real(r8) p1, p2
      real(r8)  rat, ycrit
      real(r8)  dt0, dt5
      real(r8)  dl, dp

      integer i, j, k
      integer ks
      integer js1g1, js2g0, js2g1, js2gc, js1gc, js1gcp1, js1gd
      integer jn2g0, jn1g1, jn1gc, jn1gcp1, jn1gd
      integer iord , jord
      integer ktot, ktotp

      real(r8)  zt_c   
      real(r8)  zt_d  
      real(r8)  tau, fac, pk4
      real(r8)  tiny
      parameter (tiny = 1.e-10)
      double precision pi

! Declare permanent local arrays
      integer  ifax(13)                      !ECMWF fft
      real(r8), allocatable, save :: trigs(:)
      real(r8), allocatable, save :: dc(:,:), de(:,:), sc(:), se(:)
      real(r8), allocatable, save :: fc(:), f0(:)

      real(r8), allocatable, save :: cdx(:,:), cdy(:,:)
      real(r8), allocatable, save :: dtdx(:), dtdxe(:), txe5(:), dtxe5(:)
      real(r8), allocatable, save :: dyce(:),   dx(:) ,  rdx(:),    cy(:)
      real(r8), allocatable, save :: dtdx2(:), dtdx4(:),  dxdt(:), dxe(:)
      real(r8), allocatable, save :: cye(:),    dycp(:),  rdxe(:)

      real(r8) :: uxx(im,jfirst-ng_d:jlast+ng_d+1)  ! temporary only
      real(r8) :: vxx(im,jfirst-ng_d-1:jlast+ng_d)  ! temporary only

      real(r8) rdy, dtdy, dydt, dtdy5, tdy5
      data  dt0 / 0./

      save ifax
      save dtdy, dydt, dtdy5, tdy5, rdy
      save dl, dp
      save zt_c, zt_d

#if defined( SPMD )
      integer :: incount, outcount
#endif

!******************************************************************
!******************************************************************
!
! IMPORTANT CODE OPTIONS - SEE BELOW
!
!******************************************************************
!******************************************************************

! Option for which version of geopk to use with yz decomposition.
! Version geopk16 (geopk16byte =true) uses 16-byte reals to preserve accuracy
!   through round-off (order of summation varies with z decomposition).
!   This version avoids transposes and instead involves semi-global
!   communication in Z.
! If geopk16byte=false, variables are transposed to/from xy decomposition
!   for use in geopk.
! On last small timestep (ipe=1) for D-grid, the version of geopk that uses transposes
!   is called regardless, as some transposed quantities are required for
!   the te_map phase.

      logical geopk16byte
!     data geopk16byte / .true. /
      data geopk16byte / .false. /
      logical geopkc16, geopkd16

      geopkc16 = .false.
      geopkd16 = .false.
      if (geopk16byte) then
        geopkc16 = .true.
        if (ipe .ne. 1) geopkd16 = .true.
      endif

!******************************************************************

      ktot = klast - kfirst + 1
      ktotp = ktot + 1

! Set general loop limits
! jfirst >= 1; jlast <= jm
      js1g1  = max(1,jfirst-1)
      js2g0  = max(2,jfirst)
      js2g1  = max(2,jfirst-1)
      jn2g0  = min(jm-1,jlast)
      jn1g1  = min(jm,jlast+1)

! Construct C-grid dependent loop limits
      js1gc  = max(1,jfirst-ng_c)     ! ng_c latitudes on S (starting at 1)
      js1gcp1= max(1,jfirst-ng_c-1)   ! ng_c+1 latitudes on S (starting at 1)
      jn1gc  = min(jm,jlast+ng_c)     ! ng_c latitudes on N (ending at jm)
      jn1gcp1= min(jm,jlast+ng_c+1)   ! ng_c+1 latitudes on N (ending at jm)
!
! WS: 00.04.13 : An ugly hack to handle JORD>1 JCD=1
!
      js2gc  = max(2,jfirst-ng_c)     ! NG latitudes on S (starting at 2)
      if ( ng_c == 1 .AND. ng_d > 1 ) THEN
        js2gc  = max(2,jfirst-2) 
      endif

! Construct D-grid dependent loop limits
      js1gd  = max(1,jfirst-ng_d)     ! ng_d latitudes on S (starting at 1)
      jn1gd  = min(jm,jlast+ng_d)     ! ng_d latitudes on N (ending at jm)

      if( abs(dt0-dt) > 0.1 ) then
        if ( .not. allocated( dtdx ) ) then
          allocate(dtdx(jm),dtdx2(jm), dtdx4(jm), dtdxe(jm), dxdt(jm), &
                   dxe(jm),  cye(jm),dycp(jm),rdxe(jm),                &
                   txe5(jm), dtxe5(jm),dyce(jm),                       &
                   dx(jm),rdx(jm),cy(jm) )

          allocate( trigs(3*im/2+1) )
          allocate( sc(js2g0:jn2g0),    se(js2g0:jn1g1)    )
          allocate( dc(im,js2g0:jn2g0), de(im,js2g0:jn1g1) )

          call fftfax(im, ifax, trigs)

! Determine ycrit such that effective DX >= DY
          pi  = 4.d0 * datan(1.d0)
          rat = float(im)/float(2*(jm-1))
          ycrit = acos( min(0.81, rat) ) * (180./pi)

          call pft_cf(im, jm, js2g0, jn2g0, jn1g1, sc, se, dc, de,  &
                      cosp, cose, ycrit)

          allocate( cdx(js2g0:jn1g1,kfirst:klast) )
          allocate( cdy(js2g0:jn1g1,kfirst:klast) )

          allocate( f0(jfirst-ng_s:jlast+ng_d) ) ! 000304 bug fix: ng_s not ng_d
          allocate( fc(js2gc:jn1gc) )

          do j=max(1,jfirst-ng_s),min(jm,jlast+ng_d)        ! 000304 bug fix
             f0(j) = (om+om)*sinp(j)
          enddo

! Compute coriolis parameter at cell corners.
          do j=js2gc, jn1gc                    ! Not the issue with ng_c = ng_d 
             fc(j) = 0.5*(f0(j) + f0(j-1))
          enddo

        endif

        dt0 = dt
        dt5 = 0.5*dt

        pi  = 4.d0 * datan(1.d0)
        dl  = (pi+pi)/im
        dp  = pi/(jm-1)

        rdy   = 1./(ae*dp)
        dtdy  = dt *rdy
        dtdy5 = dt5*rdy
        dydt  = (ae*dp) / dt
        tdy5  = 0.5/dtdy

        do j=2,jm-1
          dx(j)    = dl*ae*cosp(j)
          rdx(j)   = 1. / dx(j)
          dtdx(j)  = dt / dx(j)
          dxdt(j)  = dx(j) / dt
          dtdx2(j) = 0.5*dtdx(j)
          dtdx4(j) = 0.5*dtdx2(j)
          dycp(j)  = ae*dp/cosp(j)
          cy(j)    =  rdy * acosp(j)
        enddo

        do j=2,jm
          dxe(j)   = ae*dl*cose(j)
          rdxe(j)  = 1. / dxe(j)
          dtdxe(j) = dt / dxe(j)
          dtxe5(j) = 0.5*dtdxe(j)
          txe5(j)  = 0.5/dtdxe(j)
           cye(j)  =  1. / (ae*cose(j)*dp)
          dyce(j)  = ae*dp/cose(j)
        enddo

! C-grid
        zt_c = abs(umax*dt5) / (dl*ae)
!       CPP_PRT_PREFIX write(6,*) 'c-core: ', (180./pi)*acos(zt_c)
! D-grid
        zt_d = abs(umax*dt) / (dl*ae)
!       CPP_PRT_PREFIX write(6,*) 'd-coret: ', (180./pi)*acos(zt_d)

        if ( ptopin .ne. ptop) then
             write(6,*) 'PTOP as input to cd_core != ptop from dynamics_vars'
             stop
        endif

!-----------------------------------------
! Divergence damping coeff. dx(2)*dy/(TAU)
!-----------------------------------------
          do k=kfirst,klast
             press = 0.5 * ( ak(k)+ak(k+1) + (bk(k)+bk(k+1))*1.E5 )
             tau = 8. * (1.+ tanh(1.0*log(ptop/press)) )
             tau = max(1., tau) / (128.*abs(dt))
            do j=js2g0,jn1g1
               fac = tau * ae / cose(j)
               cdx(j,k) = fac*dp
               cdy(j,k) = fac*dl
            enddo
            CPP_PRT_PREFIX write(6,*)                   &
                 k, 0.01*press, cdx(js2g0,k)*ae*dl*cose(js2g0)/1.E6
          enddo
      endif

#if defined( SPMD )

      call t_startf('send_1')
      call mp_send3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_s, ng_d, u, 1)
      call mp_send3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_d, ng_s, v, 2)
      call t_stopf('send_1')
#endif

      if ( ipe == -1 .or. ns == 1 ) then          ! starting cd_core

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

         do k=kfirst,klast
            do j=jfirst,jlast
               do i=1,im
                  delpf(i,j,k) = delp(i,j,k)
               enddo
            enddo
            call pft2d( delpf(1,js2g0,k), sc(js2g0), dc(1,js2g0),  &
                        im, jn2g0-js2g0+1,   ifax, trigs )
         enddo
      endif

#if defined( SPMD )

      call t_startf('recv_1')
      call mp_barrier()
      call mp_recv3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_s, ng_d, u, 1)
      call mp_recv3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_d, ng_s, v, 2)

      call mp_send3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_d, ng_d, pt, 1)
      if ( ipe == -1 .or. ns == 1 ) then          ! starting cd_core
         call mp_send3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                           ng_d, ng_d, delpf, 2 )
      endif                          ! end if ipe = -1 check
      call mp_barrier()
      call mp_recv3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                        ng_d, ng_d, pt, 1)
      if ( ipe == -1 .or. ns == 1 ) then          ! starting cd_core
         call mp_recv3d_ns(im, jm, jfirst, jlast, kfirst, klast, &
                           ng_d, ng_d, delpf, 2)
      endif                          ! end if ipe = -1 check

      call t_stopf('recv_1')