scan2.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

subroutine scan2 (ztodt, cwava, etamid)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Second gaussian latitude scan, converts from spectral coefficients to 
! grid point values, from poles to equator, with read/calculate/write cycle.
! 
! Method: 
! The latitude pair loop in this routine is multitasked.
!
! The grid point values of ps, t, u, v, z (vorticity), and d (divergence)
! are calculated and stored for each latitude from the spectral coefficients.
! In addition, the pressure-surface corrections to the horizontal diffusion
! are applied and the global integrals of the constituent fields are 
! computed for the mass fixer.
!
! Author: 
! Original version:  CCM1
!
!-----------------------------------------------------------------------
!
! $Id: scan2.F90,v 1.10 2001/10/19 17:50:32 eaton Exp $
! $Author: eaton $
!
!-----------------------------------------------------------------------

   use precision
   use pmgrid
   use comslt
   use prognostics
   use rgrid
   use mpishorthand
!-----------------------------------------------------------------------
   implicit none
!------------------------------Commons----------------------------------
#include <comqfl.h>
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
!
! Input arguments
!
   real(r8), intent(in) :: ztodt                ! twice the timestep unless nstep = 0
   real(r8), intent(in) :: cwava(plat)          ! weight applied to global integrals
   real(r8), intent(in) :: etamid(plev)         ! vertical coords at midpoints 
!
!---------------------------Local workspace-----------------------------
!
   real(r8) hw2a(pcnst)   ! component of constituent global mass integral (mass weighting is 
                          ! based upon the "A" portion of the hybrid grid)
   real(r8) hw2b(pcnst)   ! component of constituent global mass integral (mass weighting is 
                          ! based upon the "B" portion of the hybrid grid)
   real(r8) hw3a(pcnst)   ! component of constituent global mass integral (mass weighting is 
                          ! based upon the "A" portion of the hybrid grid)
   real(r8) hw3b(pcnst)   ! component of constituent global mass integral (mass weighting is 
                          ! based upon the "B" portion of the hybrid grid)
   real(r8) hwxa(pcnst,4)
   real(r8) hwxb(pcnst,4)
   real(r8) hw2al(pcnst,plat)   ! |------------------------------------
   real(r8) hw2bl(pcnst,plat)   ! | latitudinal contributions to the
   real(r8) hw3al(pcnst,plat)   ! | components of global mass integrals
   real(r8) hw3bl(pcnst,plat)   ! |
   real(r8) hwxal(pcnst,4,plat) ! |
   real(r8) hwxbl(pcnst,4,plat) ! |-----------------------------------
!                                
! Symmetric fourier coefficient arrays for all variables transformed 
! from spherical harmonics (see subroutine grcalc)
!                                
   real(r8) grdpss(plond,begirow:endirow)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8) grzs(plond,plev,begirow:endirow)   ! sum(n) of z(n,m)*P(n,m) 
   real(r8) grds(plond,plev,begirow:endirow)   ! sum(n) of d(n,m)*P(n,m)
   real(r8) gruhs(plond,plev,begirow:endirow)  ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grvhs(plond,plev,begirow:endirow)  ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grths(plond,plev,begirow:endirow)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8) grpss(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8) grus(plond,plev,begirow:endirow)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grvs(plond,plev,begirow:endirow)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grts(plond,plev,begirow:endirow)   ! sum(n) of t(n,m)*P(n,m)
   real(r8) grpls(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8) grpms(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*H(n,m)
!
! Antisymmetric fourier coefficient arrays for all variables transformed
! from spherical harmonics (see grcalc)
!
   real(r8) grdpsa(plond,begirow:endirow)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8) grza(plond,plev,begirow:endirow)   ! sum(n) of z(n,m)*P(n,m)
   real(r8) grda(plond,plev,begirow:endirow)   ! sum(n) of d(n,m)*P(n,m)
   real(r8) gruha(plond,plev,begirow:endirow)  ! sum(n)K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grvha(plond,plev,begirow:endirow)  ! sum(n)K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grtha(plond,plev,begirow:endirow)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8) grpsa(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8) grua(plond,plev,begirow:endirow)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grva(plond,plev,begirow:endirow)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8) grta(plond,plev,begirow:endirow)   ! sum(n) of t(n,m)*P(n,m)
   real(r8) grpla(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8) grpma(plond,begirow:endirow)       ! sum(n) of lnps(n,m)*H(n,m)

   integer m,n, irow            ! indices
   integer lat,j                ! latitude indices
!
!-----------------------------------------------------------------------
!

   call t_startf ('grcalc')

#if ( defined SPMD )

!$OMP PARALLEL DO PRIVATE (J)

   do j=begirow,endirow
      call grcalcs (j, ztodt, grts(1,1,j), grths(1,1,j), grds(1,1,j), &
                    grzs(1,1,j), grus(1,1,j), gruhs(1,1,j), grvs(1,1,j), grvhs(1,1,j), &
                    grpss(1,j), grdpss(1,j), grpms(1,j), grpls(1,j))

      call grcalca (j, ztodt, grta(1,1,j), grtha(1,1,j), grda(1,1,j), &
                    grza(1,1,j), grua(1,1,j), gruha(1,1,j), grva(1,1,j), grvha(1,1,j), &
                    grpsa(1,j), grdpsa(1,j), grpma(1,j), grpla(1,j))
   end do

#else

!$OMP PARALLEL DO PRIVATE (LAT, J)

   do lat=beglat,endlat
      if (lat > plat/2) then
         j = plat - lat + 1
         call grcalcs (j, ztodt, grts(1,1,j), grths(1,1,j), grds(1,1,j), &
                       grzs(1,1,j), grus(1,1,j), gruhs(1,1,j), grvs(1,1,j), grvhs(1,1,j), &
                       grpss(1,j), grdpss(1,j), grpms(1,j), grpls(1,j))
      else
         j = lat
         call grcalca (j, ztodt, grta(1,1,j), grtha(1,1,j), grda(1,1,j), &
                       grza(1,1,j), grua(1,1,j), gruha(1,1,j), grva(1,1,j), grvha(1,1,j), &
                       grpsa(1,j), grdpsa(1,j), grpma(1,j), grpla(1,j))
      end if
   end do

#endif

   call t_stopf ('grcalc')

!#if ( defined SPMD )
!   call t_startf ('exchange')
!   call exchange (grdpss, grzs, grds, gruhs, grvhs, &
!                  grths, grpss, grus, grvs, grts, &
!                  grpls, grpms, &
!                  grdpsa, grza, grda, gruha, grvha, &
!                  grtha, grpsa, grua, grva, grta, &
!                  grpla, grpma)
!   call t_stopf ('exchange')
!#endif

   call t_startf('spegrd')

!$OMP PARALLEL DO PRIVATE (LAT, J, IROW)

   do lat=beglat,endlat
      j = j1 - 1 + lat
      irow = lat
      if (lat > plat/2) irow = plat - lat + 1
      call spegrd (ztodt, lat, cwava(lat), qfcst(i1,1,1,lat), &
                   etamid, ps(1,lat,n3), u3(i1,1,j,n3), v3(i1,1,j,n3), t3(i1,1,j,n3), &
                   qminus(i1,1,1,j), vort(1,1,lat,n3), div(1,1,lat,n3), hw2al(1,lat), hw2bl(1,lat), &
                   hw3al(1,lat), hw3bl(1,lat), hwxal(1,1,lat), hwxbl(1,1,lat), q3(i1,1,1,j,n3m1), &
                   grdpss(1,irow), grzs(1,1,irow), grds(1,1,irow), gruhs(1,1,irow), grvhs(1,1,irow), &
                   grths(1,1,irow), grpss(1,irow), grus(1,1,irow), grvs(1,1,irow), grts(1,1,irow), &
                   grpls(1,irow), grpms(1,irow), grdpsa(1,irow), grza(1,1,irow), grda(1,1,irow), &
                   gruha(1,1,irow), grvha(1,1,irow), grtha(1,1,irow), grpsa(1,irow), grua(1,1,irow), &
                   grva(1,1,irow), grta(1,1,irow), grpla(1,irow), grpma(1,irow), dps(1,lat), &
                   dpsl(1,lat), dpsm(1,lat), nlon(lat))
                   
   end do

   call t_stopf('spegrd')

#ifdef SPMD
#ifdef TIMING_BARRIERS
   call t_startf ('sync_realloc5')
   call mpibarrier (mpicom)
   call t_stopf ('sync_realloc5')
#endif
   call t_startf('realloc5')
   call realloc5 (hw2al   ,hw2bl   ,hw3al   ,hw3bl   ,tmass   , &
                  hw1lat  ,hwxal   ,hwxbl   )
   call t_stopf('realloc5')
#endif

!
! Accumulate and normalize global integrals for mass fixer (dry mass of
! atmosphere is held constant).
!
   call t_startf ('scan2_single')
   tmassf = 0.
   do lat=1,plat
      tmassf = tmassf + tmass(lat)
   end do
   tmassf = tmassf*.5
!
! Initialize moisture and mass integrals
!
   hw1(1) = 0.
   do m=1,pcnst
      hw2a(m) = 0.
      hw2b(m) = 0.
      hw3a(m) = 0.
      hw3b(m) = 0.
      do n=1,4
         hwxa(m,n) = 0.
         hwxb(m,n) = 0.
      end do
   end do
!
! Compute alpha for water ONLY
!
   do lat=1,plat
      hw1(1) = hw1(1) + hw1lat(1,lat)
      hw2a(1) = hw2a(1) + hw2al(1,lat)
      hw2b(1) = hw2b(1) + hw2bl(1,lat)
      hw3a(1) = hw3a(1) + hw3al(1,lat)
      hw3b(1) = hw3b(1) + hw3bl(1,lat)
   end do
!
   qmassf = hw1(1)
   if (adiabatic .or. ideal_phys) then
      fixmas = tmass0/tmassf
   else
      fixmas = (tmass0 + qmassf)/tmassf
   end if
!
   hw2(1)    = hw2a(1) + fixmas*hw2b(1)
   hw3(1)    = hw3a(1) + fixmas*hw3b(1)
   if(hw3(1) .ne. 0.) then
      alpha(1)  = ( hw1(1) - hw2(1) )/hw3(1)
   else
      alpha(1)  = 1.
   endif
!
! Compute alpha for non-water constituents
!
   do m = 2,pcnst
      hw1(m) = 0.
      do lat=1,plat
         hw1(m) = hw1(m) + hw1lat(m,lat)
      end do
      do n = 1,4
         do lat=1,plat
            hwxa(m,n) = hwxa(m,n) + hwxal(m,n,lat)
            hwxb(m,n) = hwxb(m,n) + hwxbl(m,n,lat)
         end do
      end do
      hw2a(m) = hwxa(m,1) - alpha(1)*hwxa(m,2)
      hw2b(m) = hwxb(m,1) - alpha(1)*hwxb(m,2)
      hw3a(m) = hwxa(m,3) - alpha(1)*hwxa(m,4)
      hw3b(m) = hwxb(m,3) - alpha(1)*hwxb(m,4)