sw_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

     do 2233 j=js2gc1,jn1gc          ! ffsl needed on N*ng S*(ng-1)
          ffsl(j) = .false.
          if( cose(j) < zt_c ) then
            do i=1,im
              if( abs(trm1(i,j)) > 1. ) then
                ffsl(j) = .true. 
                go to 2233
              endif
            enddo
          endif
2233    continue

   call tpcc( trm2, ymass, v2(1,jfirst-ng_d), trm1,          &
              cry(1,jfirst), im, jm, ng_d,                   &
              iord, jord, fx, fy(1,jfirst), ffsl, cose,      &
              jfirst, jlast, slope, qtmp, al, ar, a6 )

   do j=js2g0,jn2g0
         uc(1,j) = uc(1,j) + dtdx2(j)*(wk1(im,j)-wk1(1,j)) + dycp(j)*fy(1,j)
      do i=2,im
         uc(i,j) = uc(i,j) + dtdx2(j)*(wk1(i-1,j)-wk1(i,j)) + dycp(j)*fy(i,j)
      enddo
   enddo

   do j=js2g0,jlast
        do i=1,im-1
           vc(i,j) = vc(i,j) + dtdy5*(wk1(i,j-1)-wk1(i,j))-dxe(j)*fx(i+1,j)
        enddo
           vc(im,j) = vc(im,j) + dtdy5*(wk1(im,j-1)-wk1(im,j))-dxe(j)*fx(1,j)
   enddo

 end subroutine c_sw


!--------------------------------------------------------------------------
 subroutine d_sw( u,     v,      uc,       vc,                       &   
                  pt,    delp,   delpf,    cx3,                      &
                  cy3,    mfx,    mfy,     cdx,    cdy,              &
                  dtdx,  dtdxe,  dtxe5,  txe5,  dyce,   rdx,  cy,     &
                  dx,    f0,     js2g0,  jn1g1,                       &  
                  im,    jm,     jfirst, jlast, ng_d,   nq,           &
                  iord,  jord,   zt_d,  rcap,   tiny,                 &
                  dtdy,  dtdy5,  tdy5,   rdy,   cosp,   acosp,        &
                  cose, sinlon, coslon, sinl5, cosl5 )
!--------------------------------------------------------------------------
! Routine for shallow water dynamics on the D-grid

! !USES:

  use precision
  use tp_core

  implicit none

! INPUT:
  integer, intent(in):: im
  integer, intent(in):: jm
  integer, intent(in):: jfirst
  integer, intent(in):: jlast
  integer, intent(in):: iord
  integer, intent(in):: jord
  integer, intent(in):: js2g0
  integer, intent(in):: jn1g1
  integer, intent(in):: ng_d
  integer, intent(in):: nq

! Prognostic variables:
  real(r8), intent(in):: u(im,jfirst-ng_d:jlast+ng_d)          ! Wind in X
  real(r8), intent(in):: v(im,jfirst-ng_d:jlast+ng_d)          ! Wind in Y
  real(r8), intent(inout):: delp(im,jfirst:jlast)         ! Delta pressure
  real(r8), intent(inout):: pt(im,jfirst-ng_d:jlast+ng_d)! Potential temperature
  real(r8), intent(in):: delpf(im,jfirst-ng_d:jlast+ng_d)

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

  real(r8), intent(in):: dtdx(jm)
  real(r8), intent(in):: dtdxe(jm)
  real(r8), intent(in):: dx(jm)
  real(r8), intent(in):: rdx(jm)
  real(r8), intent(in):: cy(jm)
  real(r8), intent(in):: dyce(jm)
  real(r8), intent(in):: dtxe5(jm)
  real(r8), intent(in):: txe5(jm)

  real(r8), intent(in):: cdx(js2g0:jn1g1)
  real(r8), intent(in):: cdy(js2g0:jn1g1)
  real(r8), intent(in):: f0(jfirst-ng_d:jlast+ng_d)

  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_d
  real(r8), intent(in):: rcap
  real(r8), intent(in):: tiny
  real(r8), intent(in):: tdy5
  real(r8), intent(in):: rdy
  real(r8), intent(in):: dtdy
  real(r8), intent(in):: dtdy5

! INPUT/OUTPUT:
  real(r8), intent(inout):: uc(im,jfirst-ng_d:jlast+ng_d)
  real(r8), intent(inout):: vc(im,jfirst-2   :jlast+2 )
  real(r8), intent(inout):: cx3(im,jfirst-ng_d:jlast+ng_d)! Accumulated Courant number in X
  real(r8), intent(inout):: cy3(im,jfirst:jlast+1)        ! Accumulated Courant number in Y
  real(r8), intent(inout):: mfx(im,jfirst:jlast)          ! Mass flux in X  (unghosted)
  real(r8), intent(inout):: mfy(im,jfirst:jlast+1)        ! Mass flux in Y

! 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)   va(im,jfirst-1:jlast)
    real(r8)   ub(im,jfirst:  jlast+1)

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

    real(r8) fxj(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)  c1, c2
    real(r8)  cpt_sp, cpt_np

    logical ffsl(jm)
    logical sld

    integer i, j
    integer js2gd, jn2g0, jn2g1, jn2gd, jn1gd

! Set loop limits

  jn2g0 = min(jm-1,jlast)
  jn2g1 = min(jm-1,jlast+1)
  js2gd = max(2,jfirst-ng_d)     ! NG latitudes on S (starting at 1)
  jn2gd = min(jm-1,jlast+ng_d)   ! NG latitudes on S (ending at jm-1)
  jn1gd = min(jm,jlast+ng_d)     ! NG latitudes on N (ending at jm)

! get C-grid U-wind at poles.
  call upol5(uc(1,jfirst), vc(1,jfirst), im, jm, coslon, sinlon,        &
             cosl5, sinl5, jfirst, jlast)

  do j=js2gd,jn2gd                     ! crx needed on N*ng S*ng
     do i=1,im
        crx(i,j) = dtdx(j)*uc(i,j)
     enddo
  enddo

  do 2225 j=js2gd,jn2gd                ! ffsl needed on N*ng S*ng
     ffsl(j) = .false.
     if( cosp(j) < zt_d ) then
         do i=1,im
            if( abs(crx(i,j)) > 1. ) then
               ffsl(j) = .true. 
               go to 2225
            endif
         enddo
      endif
2225    continue

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

  do j=js2g0,jn2g0                         ! No ghosting
     do i=1,im
!       va(i,j) = 0.5*(CRY(i,j)+CRY(i,j+1))
        if( cry(i,j)*cry(i,j+1) > 0. ) then
           if( cry(i,j) > 0. ) then
              va(i,j) = cry(i,j)
           else
              va(i,j) = cry(i,j+1)         ! cry ghosted on N
           endif
        else
           va(i,j) = 0.
        endif
     enddo
  enddo

! transport polar filtered delp
      call tp2c(wk2(1,jfirst), va(1,jfirst), delpf(1,jfirst-ng_d),   &
                crx(1,jfirst-ng_d),cry(1,jfirst),im,jm,iord,jord,    &
                ng_d, xfx, yfx, ffsl,                                &
                rcap, acosp,crx(1,jfirst), ymass,                    &
                cosp, 0, jfirst, jlast)

! <<< Save necessary data for large time step tracer transport >>>
      if( nq > 0 ) then
          do j=js2g0,jn2g0                       ! No ghosting needed
            do i=1,im
              cx3(i,j) = cx3(i,j) + crx(i,j)
              mfx(i,j) = mfx(i,j) + xfx(i,j)
            enddo
          enddo

          do j=js2g0,jlast                      ! No ghosting needed
            do i=1,im
              cy3(i,j) = cy3(i,j) + cry(i,j)
              mfy(i,j) = mfy(i,j) + yfx(i,j)
            enddo
          enddo
      endif

     do j=js2g0,jn2g0                         ! No ghosting needed
        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

! Update delp
        do j=jfirst,jlast
          do i=1,im
! SAVE old delp: pressure thickness ~ "air density"
            wk1(i,j) = delp(i,j)
            trm2(i,j) = wk1(i,j) + wk2(i,j)
            delp(i,j) = trm2(i,j)
          enddo
        enddo

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

! Update pt.
      do j=jfirst,jlast
            do i=1,im
              pt(i,j) = (pt(i,j)*wk1(i,j)+wk2(i,j)) / trm2(i,j)
            enddo
      enddo

! Compute upwind biased kinetic energy at the four cell corners

! Start using ub as  v (CFL)  on B-grid (cell corners)
        do j=js2g0,jn1g1                           ! ub needed on N
            ub(1,j) = dtdy5*(vc(1,j) + vc(im,j))  
          do i=2,im
            ub(i,j) = dtdy5*(vc(i,j) + vc(i-1,j))
          enddo
        enddo

      call ytp(im, jm, fy(1,jfirst), v(1,jfirst-ng_d), ub(1,jfirst),  &
               ub(1,jfirst), ng_d, jord, 1, jfirst, jlast)
! end using ub as v (CFL) on B-grid

   do j=js2g0,jn1g1                 ! ub needed on N
       do i=1,im                
!             Need to have uc defined at the poles (j=1 and JNP)
          ub(i,j) = dtxe5(j)*(uc(i,j) + uc(i,j-1))
       enddo
   enddo

  do j=js2g0,jn1g1                       ! wk1 needed on N
          sld = .false.
          if( cose(j) < zt_d ) then
            do i=1,im
              if( abs(ub(i,j)) > 1. ) then    ! ub ghosted on N
                sld = .true. 
                go to 2235
              endif
            enddo
          endif
2235      continue

     call xtp(im,  sld, fxj, u(1,j), ub(1,j),         &
              iord, ub(1,j), cose(j), 0,               &
              slope, qtmp, al, ar, a6 )

     do i=1,im
        wk1(i,j) =  txe5(j)*fxj(i) + tdy5*fy(i,j)  ! fy ghosted on N
     enddo

  enddo

! Add divergence damping to vector invariant form of the momentum eqn
! (absolute vorticity is damped by ffsl scheme, therefore divergence damping
! provides more consistent dissipation to divergent part of the flow)

!--------------------------
! Perform divergence damping 
!--------------------------

        do j=max(2,jfirst-1), jn2g1                   ! fy need on NS (below)
            do i=1,im
              fy(i,j) = v(i,j)*cosp(j)      ! v ghosted on NS at least
            enddo
        enddo

        do j=js2g0,jn1g1                     ! wk2 needed on N (below)
! i=1
            wk2(1,j) = u(im,j) - u(1,j)    ! u ghosted on N at least
            do i=2,im
              wk2(i,j) = u(i-1,j) - u(i,j)
            enddo
        enddo

      if ( jfirst == 1 ) then
! j=2
           do i=1,im
              wk1(i,2) = wk1(i,2) - cdy(2)*fy(i, 2) + cdx(2)*wk2(i,2)
           enddo
      endif

        do j=max(3,jfirst),jn2g1         ! wk1 needed on N (after TP2D)
            do i=1,im            ! fy ghosted on NS, wk2 on N
              wk1(i,j) = wk1(i,j) + cdy(j)*(fy(i,j-1) - fy(i,j))  &
                                  + cdx(j)*wk2(i,j)
            enddo
        enddo

      if ( jlast == jm ) then
           do i=1,im
              wk1(i,jm) = wk1(i,jm) + cdy(jm)*fy(i,jm-1) + cdx(jm)*wk2(i,jm)
           enddo
      endif
!------------------------------------
! End divergence damping computation
!------------------------------------


! Compute Vorticity on the D grid

      do j=js2gd,jn1gd              ! wk3 needed on N*ng S*ng
          do i=1,im
            wk3(i,j) = u(i,j)*cose(j)   ! u ghosted on N*ng S*ng
          enddo
      enddo


      if ( jfirst ==  1 ) then
          c1 = 0.
          do i=1,im
            c1 = c1 + wk3(i,2)
          end do
          c1 = -c1*rdy*rcap

          do i=1,im
            wk2(i,1) = c1
          enddo
      endif

      if ( jlast == jm ) then
          c2 = 0.
          do i=1,im
            c2 = c2 + wk3(i,jm)
          end do
          c2 = c2*rdy*rcap

          do i=1,im
            wk2(i,jm) = c2
          enddo
      else

! This is an attempt to avoid ghosting u on N*(ng+1)
          do i=1,im
! DEBUG
!            wk2(i,jn2gd) = 0.0
! testing
             wk2(i,jn2gd) = 1.e30
          enddo
      endif

      do j=js2gd, min(jm-1,jlast+ng_d-1)
          do i=1,im-1
             wk2(i,j) = ( wk3(i,j) - wk3(i,j+1)) * cy(j)  +         &
                        (v(i+1,j) - v(i,j))    * rdx(j)
          enddo
           wk2(im,j) = (wk3(im,j) - wk3(im,j+1)) *  cy(j) +         &
                       (v(1,j) - v(im,j)) * rdx(j)
      enddo

! wk2 is relative vorticity at this point

      do j=max(1,jfirst-ng_d), jn1gd
          do i=1,im
             wk3(i,j) = wk2(i,j) + f0(j)
! wk3 is absolute vorticity
          enddo
      enddo

      call tp2d(va(1,jfirst), wk3(1,jfirst-ng_d), crx(1,jfirst-ng_d),  &
                cry(1,jfirst), im, jm, iord, jord, ng_d, fx,           &
                fy(1,jfirst), ffsl, crx(1,jfirst),                     &
                ymass, cosp, 0, jfirst, jlast)

      do j=js2g0,jlast
          do i=1,im-1
            uc(i,j) = dtdxe(j)*(wk1(i,j)-wk1(i+1,j)) + dyce(j)*fy(i,j)
          enddo
           uc(im,j) = dtdxe(j)*(wk1(im,j)-wk1(1,j)) + dyce(j)*fy(im,j)
      enddo

      do j=js2g0,jn2g0
          do i=1,im
            vc(i,j) = dtdy*(wk1(i,j)-wk1(i,j+1)) - dx(j)*fx(i,j)
          enddo
      enddo

 end subroutine d_sw

end module sw_core