#include <misc.h>module sw_core!! This module contains vertical independent part of the Lagrangian dynamicscontains!-------------------------------------------------------------------------- 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          endif2222    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          endif3333      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