inidat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                       pcnst+pnats,qmin ,q3_tmp(1,1,1,lat))
         end do
!
!-----------------------------------------------------------------------
! Integrals of mass, moisture and geopotential height
!-----------------------------------------------------------------------
!
! Compute pdel from "A" portion of hybrid vertical grid
!
         do k=1,plev
            hyad(k) = hyai(k+1) - hyai(k)
         end do
         do k=1,plev
            do i=1,plon
               pdela(i,k) = hyad(k)*ps0
            end do
         end do
!        
! Initialize mass and moisture integrals for summation
! in a third calculation loop (assures bit-for-bit compare
! with non-random history tape).
!
         tmassf_tmp = 0.
         qmass1_tmp = 0.
         qmass2_tmp = 0.
         zgsint_tmp = 0.
!
! Compute integrals of mass, moisture, and geopotential height
!
         do irow = 1,plat/2
            do ihem=1,2
               if (ihem.eq.1) then
                  lat = irow
               else
                  lat = plat - irow + 1
               end if
!              
! Accumulate average mass of atmosphere
!
               call pdelb0(ps_tmp(1,lat),pdelb,nlon(lat))
               pssum  = 0.
               zgssum = 0.
               do i=1,nlon(lat)
                  pssum  = pssum  + ps_tmp  (i,lat)
                  zgssum = zgssum + phis_tmp(i,lat)
               end do
               tmassf_tmp = tmassf_tmp + w(irow)*pssum/nlon(lat)
               zgsint_tmp = zgsint_tmp + w(irow)*zgssum/nlon(lat)
!
! Calculate global integrals needed for water vapor adjustment
!
               do k=1,plev
                  dotproda = 0.
                  dotprodb = 0.
                  do i=1,nlon(lat)
                     dotproda = dotproda + q3_tmp(i,k,1,lat)*pdela(i,k)
                     dotprodb = dotprodb + q3_tmp(i,k,1,lat)*pdelb(i,k)
                  end do
                  qmass1_tmp = qmass1_tmp + w(irow)*dotproda/nlon(lat)
                  qmass2_tmp = qmass2_tmp + w(irow)*dotprodb/nlon(lat)
               end do
            end do
         end do                  ! end of latitude loop
!
! Normalize average mass, height
!
         tmassf_tmp = tmassf_tmp*.5/gravit
         qmass1_tmp = qmass1_tmp*.5/gravit
         qmass2_tmp = qmass2_tmp*.5/gravit
         zgsint_tmp = zgsint_tmp*.5/gravit
         qmassf_tmp = qmass1_tmp + qmass2_tmp
!
! Globally avgd sfc. partial pressure of dry air (i.e. global dry mass):
!
         tmass0 = 98222./gravit
         if (ideal_phys) tmass0 = 100000./gravit
         write(6,800) tmassf_tmp,tmass0,qmassf_tmp
         write(6,810) zgsint_tmp
800      format('INIDAT: MASS OF INITIAL DATA BEFORE CORRECTION = ' &
                ,1p,e20.10,/,' DRY MASS WILL BE HELD = ',e20.10,/, &
                ' MASS OF MOISTURE AFTER REMOVAL OF NEGATIVES = ',e20.10) 
810      format(/69('*')/'INIDAT: Globally averaged geopotential ', &
                'height = ',f16.10,' meters'/69('*')/)
!
! Compute and apply an initial mass fix factor which preserves horizontal
! gradients of ln(ps).
!
         if (adiabatic .or. ideal_phys) then
            fixmas = tmass0/tmassf_tmp
         else
            fixmas = (tmass0 + qmass1_tmp)/(tmassf_tmp - qmass2_tmp)
         end if
         do lat=1,plat
            do i=1,nlon(lat)
               ps_tmp(i,lat) = ps_tmp(i,lat)*fixmas
            end do
         end do
      endif                     ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
      call copy_inidat
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
      deallocate ( ps_tmp )
      deallocate ( u3_tmp )
      deallocate ( v3_tmp )
      deallocate ( t3_tmp )
      deallocate ( q3_tmp )
      deallocate ( qcwat_tmp )
      deallocate ( lcwat_tmp )
      deallocate ( phis_tmp )
      deallocate ( landfrac_tmp )
      deallocate ( landm_tmp )
      deallocate ( sgh_tmp )
      deallocate ( ts_tmp )
      deallocate ( tsice_tmp )
      deallocate ( tssub_tmp )
      deallocate ( dpsl_tmp )
      deallocate ( dpsm_tmp )
      deallocate ( vort_tmp )
      deallocate ( div_tmp )
      deallocate ( sicthk_tmp )
      deallocate ( snowhice_tmp )
!
      return
   end subroutine read_inidat

!*********************************************************************C

   subroutine copy_inidat
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays 
! note that the use statements below contain the definitions
! of the model arrays
!-----------------------------------------------------------------------
      use pmgrid
      use prognostics
      use buffer
      use comsrf
      use phys_grid
      use tracers, only: ixcldw
#if ( defined SPMD )
      use mpishorthand
      use spmd_dyn, only: npes, compute_gsfactors
#endif
!-----------------------------------------------------------------------
      implicit none
!------------------------------Commons----------------------------------
#include <comqfl.h>
!
! Local workspace
!
      real(r8), allocatable :: tmpchunk3d(:,:,:)             
      real(r8), allocatable :: tmpchunk(:,:)             
      integer, parameter :: iend = i1+plon-1
      integer, parameter :: jend = j1+plat-1
      integer begj
      integer n,i,j

#if ( defined SPMD )
      integer :: numperlat         ! number of values per latitude band
      integer :: numsend(0:npes-1) ! number of items to be sent
      integer :: numrecv           ! number of items to be received
      integer :: displs(0:npes-1)  ! displacement array
#endif
!-----------------------------------------------------------------------
#ifdef HADVTEST
!
!JR Overwrite fields for no mountains solid-body rotation
!
      call hadvtest_init
#endif

      begj = beglatex + numbnd

!PW Dynamics fields
#if ( defined SPMD )
      numperlat = plond
      call compute_gsfactors (numperlat, numrecv, numsend, displs)

      call mpiscatterv (ps_tmp,     numsend, displs, mpir8, ps(1,beglat,1),   numrecv, mpir8, 0, mpicom)
      call mpiscatterv (phis_tmp,   numsend, displs, mpir8, phis(1,beglat),   numrecv, mpir8, 0, mpicom)
      call mpiscatterv (dpsl_tmp,   numsend, displs, mpir8, dpsl(1,beglat),   numrecv, mpir8, 0, mpicom)
      call mpiscatterv (dpsm_tmp,   numsend, displs, mpir8, dpsm(1,beglat),   numrecv, mpir8, 0, mpicom)

      numperlat = plndlv
      call compute_gsfactors (numperlat, numrecv, numsend, displs)

      call mpiscatterv (u3_tmp,   numsend, displs, mpir8, u3(i1,1,begj,1),    numrecv, mpir8, 0, mpicom)
      call mpiscatterv (v3_tmp,   numsend, displs, mpir8, v3(i1,1,begj,1),    numrecv, mpir8, 0, mpicom)
      call mpiscatterv (t3_tmp,   numsend, displs, mpir8, t3(i1,1,begj,1),    numrecv, mpir8, 0, mpicom)
      call mpiscatterv (vort_tmp, numsend, displs, mpir8, vort(1,1,beglat,1), numrecv, mpir8, 0, mpicom)
      call mpiscatterv (div_tmp,  numsend, displs, mpir8, div(1,1,beglat,1),  numrecv, mpir8, 0, mpicom)

      numperlat = plndlv*(pcnst+pnats)
      call compute_gsfactors (numperlat, numrecv, numsend, displs)

      call mpiscatterv (q3_tmp, numsend, displs, mpir8, q3(i1,1,1,begj,1), numrecv, mpir8, 0, mpicom)

#else

      u3(i1:iend,:,j1:jend,1) = u3_tmp(:plon,:,:)
      v3(i1:iend,:,j1:jend,1) = v3_tmp(:plon,:,:)
      t3(i1:iend,:,j1:jend,1) = t3_tmp(:plon,:,:)
      q3(i1:iend,:,:,j1:jend,1) = q3_tmp(:plon,:,:,:)

      ps(:plon,:,1)    = ps_tmp(:plon,:)
      phis(:plon,:)    = phis_tmp(:plon,:)
      dpsl(:plon,:)    = dpsl_tmp(:plon,:)
      dpsm(:plon,:)    = dpsm_tmp(:plon,:)
      vort(:plon,:,:,1)= vort_tmp(:plon,:,:)
      div(:plon,:,:,1) = div_tmp(:plon,:,:)
#endif

      allocate ( tmpchunk(pcols,begchunk:endchunk) )
      allocate ( tmpchunk3d(pcols,plevmx,begchunk:endchunk) )

!PW Physics fields
      call scatter_field_to_chunk(1,1,1,plond,landfrac_tmp,landfrac(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,landm_tmp,landm(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,sgh_tmp,sgh(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,tsice_tmp,tsice(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,ts_tmp,tmpchunk)
      do i =begchunk,endchunk
         srfflx_state2d(i)%ts(:) = tmpchunk(:,i)
      end do
#if ( defined COUP_SOM )
      call scatter_field_to_chunk(1,1,1,plond,sicthk_tmp,sicthk(1,begchunk))
#endif
      call scatter_field_to_chunk(1,1,1,plond,snowhice_tmp,snowhice(1,begchunk))

      call scatter_field_to_chunk(1,plevmx,1,plond,tssub_tmp,tmpchunk3d)
      do i =begchunk,endchunk
         surface_state2d(i)%tssub(:,:) = tmpchunk3d(:,:,i)
      end do
!
!JR cloud and cloud water initialization.  Does this belong somewhere else?
!
      if (masterproc) then
         qcwat_tmp(:plon,:,:) = q3_tmp(:plon,:,1,:)
         lcwat_tmp(:plon,:,:) = q3_tmp(:plon,:,ixcldw,:)
      endif
      call scatter_field_to_chunk(1,plev,1,plond,qcwat_tmp,qcwat(1,1,begchunk,1))
      call scatter_field_to_chunk(1,plev,1,plond,lcwat_tmp,lcwat(1,1,begchunk,1))
      call scatter_field_to_chunk(1,plev,1,plond,t3_tmp,tcwat(1,1,begchunk,1))
      cld   (:,:,:,1) = 0.
      do n=2,3
         cld(:,:,:,n) = 0.
         qcwat(:,:,:,n) = qcwat(:,:,:,1) 
         tcwat(:,:,:,n) = tcwat(:,:,:,1) 
         lcwat(:,:,:,n) = lcwat(:,:,:,1) 
      end do
!
! Global integerals
!
      if (masterproc) then
         tmassf = tmassf_tmp
         qmass1 = qmass1_tmp
         qmass2 = qmass2_tmp
         qmassf = qmassf_tmp
         zgsint = zgsint_tmp
      endif

#if ( defined SPMD )
      call mpibcast (tmass0,1,mpir8,0,mpicom)
      call mpibcast (tmassf,1,mpir8,0,mpicom)
      call mpibcast (qmass1,1,mpir8,0,mpicom)
      call mpibcast (qmass2,1,mpir8,0,mpicom)
      call mpibcast (qmassf,1,mpir8,0,mpicom)
      call mpibcast (zgsint,1,mpir8,0,mpicom)
#endif
      deallocate ( tmpchunk )
      deallocate ( tmpchunk3d)

   end subroutine copy_inidat

#ifdef HADVTEST
   subroutine hadvtest_init

      use pmgrid
      use rgrid
      use commap
      use physconst, only:

      implicit none

#include <comhyb.h>
#include <hadvtest.h>

      integer i,k,lat

      real(r8) h0, u0, small_r, big_r, theta, theta_c, lambda, lambda_c
      real(r8) alfa, dlam, pie
      real(r8) pie
!
! First: zero sgh and phis fields
!
      sgh_tmp(:,:) = 0.
      phis_tmp(:,:) = 0.
!
!JR Analytic IC and wind
!
      pie = acos(-1.)
!
!jr Define wind and constituent fields
!
      h0 = 1000.
      u0 = 2.*pie*rearth/(12.*86400.)
      big_r = rearth/3.
      theta_c = +60.*pie/180.   ! 60 deg north
      theta_c = -60.*pie/180.   ! 60 deg south
      theta_c = 0.              ! equator
      lambda_c = 0.             ! Greenwich
      lambda_c = 3.*pie/2.

      do lat=1,plat
         theta = clat(lat)
         do k=1,plev
            alfa = 0.
            if (k.eq.1) then
               alfa = 0.
            else if (k.eq.2) then
               alfa = 0.05
            else if (k.eq.plev-1) then
               alfa = 0.5*pie - 0.05
            else if (k.eq.plev) then
               alfa = 0.5*pie      ! blows north
            else
               alfa = (k-2)*pie/(2.*(plev-3))
            end if

            do i=1,nlon(lat)
               lambda = 2.*pie*(i-1)/nlon(lat)
!
!jr Use these settings in conjunction with theta_c to start the blob at 
!jr Greenwich
!
               usave(i,k,lat) = u0*(cos(theta)*cos(alfa) +  &
                  sin(theta)*cos(lambda-0.5*pie)*sin(alfa))
               vsave(i,k,lat) = -u0*sin(lambda-0.5*pie)*sin(alfa)
!
!jr Use these settings in conjunction with theta_c to start the blob at 270.
!
               usave(i,k,lat) = u0*(cos(theta)*cos(alfa) + &
                  sin(theta)*cos(lambda)*sin(alfa))
               vsave(i,k,lat) = -u0*sin(lambda)*sin(alfa)
               u3_tmp(i,k,lat) = usave(i,k,lat)
               v3_tmp(i,k,lat) = vsave(i,k,lat)
               dlam = lambda - lambda_c
               small_r = rearth*acos(sin(theta_c)*sin(theta) +  &
                  cos(theta_c)*cos(theta)*cos(dlam))
               q3_tmp(i,k,1,lat) = 0.
               if (small_r .lt. big_r) then
                  q3_tmp(i,k,1,lat) = h0/2.*(1. + cos(pie*small_r/big_r))
               end if
!
!jr Stick Q into T to test spectral advection (of what's in T)
!jr Or put 300 in T.
!
               t3_tmp(i,k,lat) = 300.
               t3_tmp(i,k,lat) = q3_tmp(i,k,1,lat)
            end do
         end do
!
!jr Save surface pressure for future timesteps.  Set to 1.e5 everywhere
!
         do i=1,nlon(lat)
            ps_tmp(i,lat) = ps0
            pssave(i,lat) = ps_tmp(i,lat)
         end do
      end do

      return
   end subroutine hadvtest_init
#endif

end module inidat