#include <misc.h>
#include <params.h>
#if ( defined SCAM )
#include <max.h>
#endif


subroutine initcom 2,21

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Initialize Model commons, including COMCON, COMHYB, COMMAP, COMSPE,
! and COMTRCNM
! 
! Method: 
! 
! Author: 
! Original version:  CCM1
! Standardized:      L. Bath, Jun 1992
!                    L. Buja, Feb 1996
!
!-----------------------------------------------------------------------
!
! $Id: initcom.F90,v 1.16.2.7.4.1 2004/05/20 18:36:05 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use pspect
   use comspe
   use rgrid
   use gauaw_mod, only: gauaw
   use commap
   use dynconst, only: rearth, ra, dynconsti
   use physconst, only: rair
   use time_manager, only: get_step_size
   use abortutils, only: endrun
#if ( defined SCAM )
   use scamMod, only :columnlat,columnlon
#endif
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comfft.h>
!-----------------------------------------------------------------------
#include <comhyb.h>
!-----------------------------------------------------------------------

!
! Local workspace
!
   real(r8) zsi(plat)      ! sine of latitudes
   real(r8) zw(plat)       ! Gaussian weights
   real(r8) zra2           ! ra squared
   real(r8) zalp(2*pspt)   ! Legendre function array
   real(r8) zdalp(2*pspt)  ! Derivative array
   real(r8) zslat          ! sin of lat  and cosine of colatitude

   integer i           ! longitude index
   integer j           ! Latitude index
   integer k           ! Level index
   integer kk          ! Level index
   integer kkk         ! Level index
   integer m,lm,mr,lmr ! Indices for legendre array
   integer n           ! Index for legendre array
   integer nkk         ! Print control variables
   integer ik1         ! Print index temporary variable
   integer ik2         ! Print index temporary variable
   integer itmp        ! Dimension of polynomial arrays temporary.
   integer iter        ! Iteration index
   real(r8)    zdt         ! Time step for settau

   logical lprint      ! Debug print flag
   integer irow        ! Latitude pair index
   integer lat         ! Latitude index

   real(r8) xlat           ! Latitude (radians)
   real(r8) pi             ! Mathematical pi (3.14...)
   real(r8) dtime          ! timestep size [seconds]
!
!-----------------------------------------------------------------------
#if ( !defined SCAM )
   call dynconsti
#endif
!
   lprint = masterproc .and. .FALSE.

   dtime = get_step_size()

   call hdinti  (rearth  ,dtime   )
!
! Initialize commap.  Set hybrid level dependent arrays
!
   call hycoef

#if ( !defined SCAM )
!
! NMAX dependent arrays
!
   if (pmmax.gt.plon/2) then
      call endrun ('INITCOM:mmax=ptrm+1 .gt. plon/2')
   end if
#endif

   zra2 = ra*ra
   do j=2,pnmax
      sq(j)  = j*(j-1)*zra2
      rsq(j) = 1./sq(j)
   end do
   sq(1)  = 0.
   rsq(1) = 0.
!
! MMAX dependent arrays
!
   do j=1,pmmax
      xm(j) = j-1
   end do
#if ( !defined SCAM )
!
! Gaussian latitude dependent arrays
!
   call gauaw(zsi     ,zw      ,plat    )
   do irow=1,plat/2
      slat(irow) = zsi(irow)
      w(irow)              = zw(irow)
      w(plat - irow + 1)   = zw(irow)
      cs(irow)  = 1. - zsi(irow)*zsi(irow)
      xlat = asin(slat(irow))
      clat(irow) = -xlat
      clat(plat - irow + 1) = xlat
   end do

   do lat=1,plat
      latdeg(lat) = clat(lat)*45./atan(1._r8)
   end do
#endif
!
! Integration matrices of hydrostatic equation(href) and conversion
! term(a).  href computed as in ccm0 but isothermal bottom ecref
! calculated to conserve energy
!
   do k=1,plev
      do kk=1,plev
         href(kk,k) = 0.
         ecref(kk,k) = 0.
      end do
   end do
!
! Mean atmosphere energy conversion term is consistent with continiuty
! Eq.  In ecref, 1st index = column; 2nd index = row of matrix.
! Mean atmosphere energy conversion term is energy conserving
!
   do k=1,plev
      ecref(k,k) = 0.5/hypm(k) * hypd(k)
      do kk=1,k-1
         ecref(kk,k) = 1./hypm(k) * hypd(kk)
      end do
   end do
!
! Reference hydrostatic integration matrix consistent with conversion
! term for energy conservation.  In href, 1st index = column; 
! 2nd index = row of matrix.
!
   do k = 1,plev
      do kk = k,plev
         href(kk,k) = ecref(k,kk)*hypd(kk)/hypd(k)
      end do
   end do
!
! Print statements
!
   if (lprint) then
      nkk = plev/13
      if (mod(plev,13).ne.0) nkk = nkk + 1
      write(6,*)' '
      write(6,*)'INITCOM: Hydrostatic matrix href'
      do kk=1,nkk
         ik1 = 1 + (kk-1)*13
         ik2 = min0( ik1+12, plev )
         write(6,9920) (k,k=ik1,ik2)
         do kkk=1,plev
            write(6,9910) kkk,(href(kkk,k),k=ik1,ik2)
         end do
      end do
      write(6,*)' '
      write(6,*)'INITCOM: Thermodynamic matrix ecref'
      do kk=1,nkk
         ik1 = 1 + (kk-1)*13
         ik2 = min0( ik1+12, plev )
         write(6,9920) (k,k=ik1,ik2)
         do kkk=1,plev
            write(6,9910) kkk,(ecref(kkk,k),k=ik1,ik2)
         end do
      end do
   end if
!
! Multiply href by r
!
   do k=1,plev
      do kk=1,plev
         href(kk,k) = href(kk,k)*rair
      end do
   end do
!
! Compute truncation parameters
!
   if (masterproc) then
      write(6,9950) ptrm,ptrn,ptrk
   end if
!
! Compute semi-implicit timestep constants (COMSPE)
!
   zdt = dtime
   if (.not.nlres) zdt = 0.5*zdt
!
   call settau(zdt)
#if ( defined SCAM )

   do j=1,plat
      slat(j) = 1.0_r8 * sin(4.0_r8*atan(1.0_r8)*columnLat/180._r8)
      w(j)   = 2.0_r8/plat
      cs(j)  = 10._r8 - slat(j)*slat(j)
!         itmp = 2*pspt - 1
!         call phcs  (zalp    ,zdalp   ,itmp    ,zslat    )
!         call reordp(j       ,itmp    ,zalp    ,zdalp   )
   end do

!
! Latitude array (list of latitudes in radians)
!
   xlat = asin(slat(1))
   clat(1) = xlat

   clat(1)=columnLat*atan(1._r8)/45.
   latdeg(1) = clat(1)*45._r8/atan(1._r8)
   clon(1,1)   = 4.0_r8*atan(1._r8)*mod((columnLon+360._r8),360._r8)/180._r8
   londeg(1,1) = mod((columnLon+360._r8),360._r8)

#else
!
! Compute constants related to Legendre transforms
! Compute and reorder ALP and DALP
!
   allocate( alp  (pspt,plat/2) )
   allocate( dalp (pspt,plat/2) )
   do j=1,plat/2
      zslat = slat(j)
      itmp = 2*pspt - 1
      call phcs  (zalp    ,zdalp   ,itmp    ,zslat    )
      call reordp(j       ,itmp    ,zalp    ,zdalp   )
   end do
!
! Copy and save local ALP and DALP
!
   allocate( lalp  (lpspt,plat/2) )
   allocate( ldalp (lpspt,plat/2) )
   do j=1,plat/2
      do lm=1,numm(iam)
         m = locm(lm,iam)
         mr = nstart(m)
         lmr = lnstart(lm)
         do n=1,nlen(m)
            lalp(lmr+n,j) = alp(mr+n,j)
            ldalp(lmr+n,j) = dalp(mr+n,j)
         end do
      end do
   end do
!
! Determine whether full or reduced grid
!
   fullgrid = .true.
   do j=1,plat
      if (masterproc) then
         write(6,*)'nlon(',j,')=',nlon(j),' wnummax(',j,')=',wnummax(j)
      end if
      if (nlon(j).lt.plon) fullgrid = .false.
   end do
!
! Mirror latitudes south of south pole
!
   lat = 1
   do j=j1-2,1,-1
      nlonex(j) = nlon(lat)
      lat = lat + 1
   end do
   nlonex(j1-1) = nlon(1)     ! south pole
!
! Real latitudes
!
   j = j1
   do lat=1,plat
      nlonex(j) = nlon(lat)
      j = j + 1
   end do
   nlonex(j1+plat) = nlon(plat)  ! north pole
!
! Mirror latitudes north of north pole
!
   lat = plat
   do j=j1+plat+1,platd
      nlonex(j) = nlon(lat)
      lat = lat - 1
   end do
!
! Longitude array
!
   pi = 4.0*atan(1.0)
   do lat=1,plat
      do i=1,nlon(lat)
         londeg(i,lat) = (i-1)*360./nlon(lat)
         clon(i,lat)   = (i-1)*2.0*pi/nlon(lat)
      end do
   end do

   do j=1,plat/2
      nmmax(j) = wnummax(j) + 1
   end do
   do m=1,pmmax
      do irow=1,plat/2
         if (nmmax(irow) .ge. m) then
            beglatpair(m) = irow
            goto 10
         end if
      end do
      call endrun ('INITCOM: Should not ever get here')
10    continue
   end do
!
! Set up trigonometric tables for fft
!
   do j=1,plat
      call set99(trig(1,j),ifax(1,j),nlon(j))
   end do
!
! Set flag indicating dynamics grid is now defined.
! NOTE: this ASSUMES initcom is called after spmdinit.  The setting of nlon done here completes
! the definition of the dynamics grid.
!
   dyngrid_set = .true.
#endif

   return

9910 format( 1x,i3,13f9.5)
9920 format(/,      13i9)
9950 format(/,'     Truncation Parameters',/,'     NTRM = ',i4,/, &
      '     NTRN = ',i4,/,'     NTRK = ',i4,/)

end subroutine initcom