sw_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>

module sw_core
!
! This module contains vertical independent part of the Lagrangian dynamics

contains

!--------------------------------------------------------------------------
 subroutine c_sw(u,    v,  pt, delp,                       &
                 uc,   vc,   ptc,   delpf,   ptk,                    &
                cosp,   acosp,   cose,   coslon,   sinlon,           &
                cosl5,  sinl5,   dxdt,   dxe,      dtdx2,            &
                dtdx4,  dtxe5,   rdxe,   dycp,     dydt,             & 
                dtdy5,  cye,     fc,     ifax,     trigs,            &
                dc,     sc,      zt_c,   tiny,     rcap,             &
                im,     jm,      jfirst, jlast,    ng_c,             &
                ng_d,   ng_s,    js2g0,  jn2g0,    js2gc,            &
                jn1gc,   iord,   jord        )
!--------------------------------------------------------------------------
! Routine for shallow water dynamics on the C-grid

! !USES:

  use precision
  use tp_core
  use pft_module

  implicit none

! INPUT:
  integer, intent(in):: im
  integer, intent(in):: jm
  integer, intent(in):: jfirst
  integer, intent(in):: jlast
  integer, intent(in):: js2g0
  integer, intent(in):: jn2g0
  integer, intent(in):: js2gc
  integer, intent(in):: jn1gc
  integer, intent(in):: iord
  integer, intent(in):: jord
  integer, intent(in):: ng_c
  integer, intent(in):: ng_s
  integer, intent(in):: ng_d

! polar filter related input arrays:
  integer, intent(in)::  ifax(13)                      !ECMWF fft
  real(r8), intent(in):: trigs(3*im/2+1)
  real(r8), intent(in):: dc(im,js2g0:jn2g0)
  real(r8), intent(in):: sc(js2g0:jn2g0)

! Prognostic variables:
  real(r8), intent(in):: u(im,jfirst-ng_c:jlast+ng_c+1)
  real(r8), intent(in):: v(im,jfirst-ng_c-1:jlast+ng_c)  ! Wind in Y
  real(r8), intent(in):: pt(im,jfirst-ng_c:jlast+ng_c)   ! Wind in Y
  real(r8), intent(in):: delp(im,jfirst:jlast)         ! Delta pressure
  real(r8), intent(in):: delpf(im,jfirst-ng_c:jlast+ng_c)

  real(r8), intent(in):: fc(js2gc:jn1gc)

  real(r8), intent(in):: cosp(jm)
  real(r8), intent(in):: acosp(jm)
  real(r8), intent(in):: cose(jm)

  real(r8), intent(in):: dxdt(jm)
  real(r8), intent(in):: dxe(jm)
  real(r8), intent(in):: rdxe(jm)
  real(r8), intent(in):: dtdx2(jm)
  real(r8), intent(in):: dtdx4(jm)
  real(r8), intent(in):: dtxe5(jm)
  real(r8), intent(in):: dycp(jm)
  real(r8), intent(in)::  cye(jm)

  real(r8), intent(in):: sinlon(im)
  real(r8), intent(in):: coslon(im)
  real(r8), intent(in):: sinl5(im)
  real(r8), intent(in):: cosl5(im)

  real(r8), intent(in):: zt_c
  real(r8), intent(in):: rcap
  real(r8), intent(in):: tiny
  real(r8), intent(in):: dydt
  real(r8), intent(in):: dtdy5

! Output:
  real(r8), intent(out):: uc(im,jfirst-ng_c:jlast+ng_c)
  real(r8), intent(out):: vc(im,jfirst-2   :jlast+2 )
  real(r8), intent(out):: ptc(im,jfirst:jlast)
  real(r8), intent(out):: ptk(im,jfirst:jlast)

!--------------------------------------------------------------
! Local 
    real(r8)   fx(im,jfirst:jlast)
    real(r8)  xfx(im,jfirst:jlast)
    real(r8) trm2(im,jfirst:jlast)

    real(r8)  wk1(im,jfirst-1:jlast+1)
    real(r8)  cry(im,jfirst-1:jlast+1)
    real(r8)   fy(im,jfirst-1:jlast+1)
    real(r8) ymass(im,jfirst:  jlast+1) 
    real(r8)  yfx(im,jfirst:  jlast+1)
    real(r8) trm1(im,jfirst-1:jlast+2)

    real(r8)   va(im,jfirst-1:jlast)
    real(r8)   ub(im,jfirst:  jlast+1)

    real(r8)  wk2(im,jfirst-ng_d:jlast+ng_d)
    real(r8)   u2(im,jfirst-ng_d:jlast+ng_d)
    real(r8)   v2(im,jfirst-ng_d:jlast+ng_d)
    real(r8)  crx(im,jfirst-ng_d:jlast+ng_d)

    real(r8) fxj(im)
    real(r8) p1d(im)
    real(r8) cx1(im)

    real(r8) qtmp(-im/3:im+im/3)
    real(r8) slope(-im/3:im+im/3)
    real(r8) al(-im/3:im+im/3)
    real(r8) ar(-im/3:im+im/3)
    real(r8) a6(-im/3:im+im/3)

    real(r8) us, vs, un, vn
    real(r8) p1ke, p2ke

    logical ffsl(jm)
    logical sld

    integer i, j, im2
    integer js1g1
    integer js2g1
    integer js2gc1
    integer js2gcp1
    integer jn2gc
    integer jn1g1

    im2 = im/2

! Set loop limits

    js1g1 = max(1,jfirst-1)
    js2g1 = max(2,jfirst-1)
    js2gcp1 = max(2,jfirst-ng_c-1)   ! NG-1 latitudes on S (starting at 2)
    jn1g1 = min(jm,jlast+1)
    jn2gc = min(jm-1,jlast+ng_c)   ! NG latitudes on N (ending at jm-1)
 
!
! Treat the special case of ng_c = 1
!
    if ( ng_c == 1 .AND. ng_d > 1 ) THEN
        js2gc1 = js2gc
    else
        js2gc1 = max(2,jfirst-ng_c+1)   ! NG-1 latitudes on S (starting at 2)
    endif

    do j=js2gcp1,jn2gc                ! u2 local partition only, not jlast
       do i=1,im-1
          v2(i,j) = v(i,j) + v(i+1,j)
       enddo
       v2(im,j) = v(im,j) + v(1,j)
    enddo

    do j=js2gc,jn2gc
       do i=1,im
          u2(i,j) = u(i,j) + u(i,j+1)
       enddo
    enddo

        if ( jfirst == 1 ) then
! Projection at SP
          us = 0.
          vs = 0.

          do i=1,im2
            us = us + (u2(i+im2,2)-u2(i,2))*sinlon(i)         &
                    + (v2(i,2)-v2(i+im2,2))*coslon(i)
            vs = vs + (u2(i+im2,2)-u2(i,2))*coslon(i)         &
                    + (v2(i+im2,2)-v2(i,2))*sinlon(i)
          enddo

          us = us/im
          vs = vs/im

! SP
          do i=1,im2
            u2(i,1)  = -us*sinlon(i) - vs*coslon(i)
            v2(i,1)  =  us*coslon(i) - vs*sinlon(i)
            u2(i+im2,1)  = -u2(i,1)
            v2(i+im2,1)  = -v2(i,1)
          enddo

          p1ke = 0.125*(u2(1, 1)**2 + v2(1, 1)**2)

        endif

        if ( jlast == jm ) then

! Projection at NP
          un = 0.
          vn = 0.

          j = jm-1
          do i=1,im2
            un = un + (u2(i+im2,j)-u2(i,j))*sinlon(i)        &
                    + (v2(i+im2,j)-v2(i,j))*coslon(i)
            vn = vn + (u2(i,j)-u2(i+im2,j))*coslon(i)        &
                    + (v2(i+im2,j)-v2(i,j))*sinlon(i)
          enddo

          un = un/im
          vn = vn/im

! NP
          do i=1,im2
            u2(i,jm) = -un*sinlon(i) + vn*coslon(i)
            v2(i,jm) = -un*coslon(i) - vn*sinlon(i)
            u2(i+im2,jm) = -u2(i,jm)
            v2(i+im2,jm) = -v2(i,jm)
          enddo

          p2ke = 0.125*(u2(1,jm)**2 + v2(1,jm)**2)

        endif

! A -> C
        do j=js2gc,jn2gc                ! uc needed N*ng S*ng
! i=1
            uc(1,j) = 0.25*(u2(1,j)+u2(im,j))
          do i=2,im
            uc(i,j) = 0.25*(u2(i,j)+u2(i-1,j))
          enddo
        enddo

        do j=js2gc,jn1gc                ! vc needed N*ng, S*ng (for ycc)
          do i=1,im
            vc(i,j) = 0.25*(v2(i,j)+v2(i,j-1))  ! v2 needed N*ng S*(ng+1)
          enddo
        enddo

        do j=js2g1,jn1g1                   ! cry needed on NS
          do i=1,im
            cry(i,j) = dtdy5*vc(i,j)
          enddo
        enddo

        do j=js2g0,jn1g1                     ! ymass needed on NS
          do i=1,im
            ymass(i,j) = cry(i,j)*cose(j)
          enddo
        enddo

! New va definition
        do j=js2g1,jn2g0                     ! va needed on S (for YCC, iv==1)
          do i=1,im
            va(i,j) = 0.5*(cry(i,j)+cry(i,j+1))
          enddo
        enddo

! SJL: Check if FFSL integer fluxes need to be computed

        do 2222 j=js2gc,jn2gc                ! ffsl needed on N*sg S*sg
          do i=1,im
            crx(i,j) = uc(i,j)*dtdx2(j)
          enddo
          ffsl(j) = .false.
          if( cosp(j) < zt_c ) then
            do i=1,im
              if( abs(crx(i,j)) > 1. ) then
                ffsl(j) = .true. 
                go to 2222
              endif
            enddo
          endif
2222    continue

! 2D transport of polar filtered delp (for computing fluxes!)
! Update is done on the unfiltered delp

   call tp2c( trm2,  va(1,jfirst),  delpf(1,jfirst-ng_c),    &
              crx(1,jfirst-ng_c), cry(1,jfirst),             &
              im, jm, iord, jord, ng_c, xfx,                 &
              yfx, ffsl, rcap, acosp,                        &
              crx(1,jfirst), ymass, cosp,                    &
              0, jfirst, jlast)

   do j=js2g0,jn2g0                      ! xfx not ghosted
      if( ffsl(j) ) then
         do i=1,im
           xfx(i,j) = xfx(i,j)/sign(max(abs(crx(i,j)),tiny),crx(i,j))
         enddo
      endif
   enddo

! pt-advection using pre-computed mass fluxes
! use ub below as the storage for pt (wk2) increment
! WS 99.09.20 : pt, crx need on N*ng S*ng, yfx on N

    call tp2c(ub ,va(1,jfirst), pt(1,jfirst-ng_c),        &
              crx(1,jfirst-ng_c), cry(1,jfirst),          &
              im, jm,  iord, jord, ng_c, fx,              &
              fy(1,jfirst), ffsl, rcap, acosp,            &
              xfx, yfx, cosp, 1, jfirst, jlast)

#if defined ( ALT_PFT )
! Alternative way of polar filter; nothing to be done here
#else
     call pft2d(trm2(1,js2g0), sc(js2g0), dc(1,js2g0), im,  &
                jn2g0-js2g0+1,  ifax, trigs )
     call pft2d( ub(1,js2g0), sc(js2g0), dc(1,js2g0),  im,  &
                 jn2g0-js2g0+1,  ifax, trigs )
#endif

    do j=jfirst,jlast
       do i=1,im
          ptk(i,j) = delp(i,j) + trm2(i,j)
          ptc(i,j) = (pt(i,j)*delp(i,j) + ub(i,j))/ptk(i,j)
       enddo
    enddo

!------------------
! Momentum equation
!------------------

     call ycc(im, jm, fy, vc(1,jfirst-2), va(1,jfirst-1),   &
              va(1,jfirst-1), jord, 1, jfirst, jlast)

     do j=js2g1,jn2g0

          do i=1,im
            cx1(i) = dtdx4(j)*u2(i,j)
          enddo

          sld = .false.
          if( cosp(j) < zt_c ) then
            do i=1,im
              if( abs(cx1(i)) > 1. ) then
                sld = .true. 
                go to 3333
              endif
            enddo
          endif
3333      continue

          p1d(im) = uc(1,j)
          do i=1,im-1
            p1d(i) = uc(i+1,j)
          enddo

          call xtp(im,   sld, fxj, p1d, cx1, iord,   &
                   cx1,  cosp(j),  0,   slope,       &
                   qtmp, al,       ar,  a6)
 
          do i=1,im
            wk1(i,j) = dxdt(j)*fxj(i) + dydt*fy(i,j)
          enddo
     enddo


     if ( jfirst == 1 ) then
          do i=1,im
            wk1(i,1) = p1ke
          enddo
     endif

     if ( jlast == jm ) then
          do i=1,im
            wk1(i,jm) = p2ke
          enddo
     endif

! crx redefined as trm1
     do j=js2gc1,jn1gc
! i=1
            trm1(1,j) = dtxe5(j)*u(im,j)
          do i=2,im
            trm1(i,j) = dtxe5(j)*u(i-1,j)
          enddo
     enddo

     do j=js1g1,jlast                       ! cry needed on S
          do i=1,im
            cry(i,j) = dtdy5*v(i,j)
          enddo
     enddo

     do j=jfirst,jlast
          do i=1,im
            ymass(i,j) = cry(i,j)*cosp(j)       ! ymass actually unghosted
          enddo
     enddo

     do j=js2g0,jlast
          do i=1,im
            trm2(i,j) = 0.5*(cry(i,j)+cry(i,j-1)) ! cry ghosted on S 
          enddo
     enddo

!    Compute absolute vorticity on the C-grid.

     if ( jfirst == 1 ) then
          do i=1,im
            wk2(i,1) = 0.
          enddo
     endif

     do j=js2gc,jn2gc
          do i=1,im
            wk2(i,j) = uc(i,j)*cosp(j)
          enddo
     enddo

     if ( jlast == jm ) then
          do i=1,im
            wk2(i,jm) = 0.
          enddo
     endif

     do j=js2gc1,jn1gc
! The computed absolute vorticity on C-Grid is assigned to v2
          v2(1,j) = fc(j) + (wk2(1,j-1)-wk2(1,j))*cye(j) +     &
                    (vc(1,j) - vc(im,j))*rdxe(j)

          do i=2,im
             v2(i,j) = fc(j) + (wk2(i,j-1)-wk2(i,j))*cye(j) +  &
                       (vc(i,j) - vc(i-1,j))*rdxe(j)
          enddo
     enddo