inidat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

       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 ( tl_tmp )
    deallocate ( tm_tmp )
    deallocate ( ql_tmp )
    deallocate ( qm_tmp )
    deallocate ( phis_tmp )        
    deallocate ( phisl_tmp )        
    deallocate ( phism_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 ( div_tmp )
    deallocate ( sicthk_tmp )
    deallocate ( snowhice_tmp )
!
    return
  end subroutine read_inidat

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

  subroutine copy_inidat
!-----------------------------------------------------------------------
!
! Purpose:
! Copy temporary arrays to model arrays 
! note that the use statements below contain the definitions
! of the model arrays
!
!-----------------------------------------------------------------------

    use precision
    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

#include <comqfl.h>
!
!---------------------------Local workspace-----------------------------
!
    real(r8), allocatable :: tmpchunk3d(:,:,:)             
    real(r8), allocatable :: tmpchunk(:,:)             
    integer, parameter :: iend = i1+plon-1 ! last "real model" i-index in extended grid
    integer, parameter :: jend = j1+plat-1 ! last "real model" j-index in extended grid
    integer begj, endj
    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 flat-earth solid-body rotation
!
    call hadvtest_init
#endif

    begj = beglatex + numbnd
    endj = begj + numlats - 1

!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 (phisl_tmp ,numsend, displs, mpir8,phisl (1,beglat)   ,numrecv, mpir8,0,mpicom)
    call mpiscatterv (phism_tmp ,numsend, displs, mpir8,phism (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 (div_tmp   ,numsend, displs, mpir8,div   (1,1,beglat,1) ,numrecv, mpir8,0,mpicom)
    call mpiscatterv (tl_tmp    ,numsend, displs, mpir8,tl    (1,1,beglat) ,numrecv, mpir8,0,mpicom)
    call mpiscatterv (tm_tmp    ,numsend, displs, mpir8,tm    (1,1,beglat) ,numrecv, mpir8,0,mpicom)
    call mpiscatterv (ql_tmp    ,numsend, displs, mpir8,ql    (1,1,beglat) ,numrecv, mpir8,0,mpicom)
    call mpiscatterv (qm_tmp    ,numsend, displs, mpir8,qm    (1,1,beglat) ,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)

    call mpibcast   (lnpstar   ,psp,mpir8,0 , mpicom)

#else

    ps    (:,:,1) = ps_tmp    (:,:)
    phis  (:,:)   = phis_tmp  (:,:)
    phisl (:,:)   = phisl_tmp (:,:)
    phism (:,:)   = phism_tmp (:,:)
    dpsl  (:,:)   = dpsl_tmp  (:,:)
    dpsm  (:,:)   = dpsm_tmp  (:,:)

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

    div   (:,:,:,1) = div_tmp   (:,:,:)
    tl    (:,:,:) = tl_tmp    (:,:,:)
    tm    (:,:,:) = tm_tmp    (:,:,:)
    ql    (:,:,:) = ql_tmp    (:,:,:)
    qm    (:,:,:) = qm_tmp    (:,:,:)

    q3    (i1:iend,:,:,j1:jend,1) = q3_tmp(:,:,:,:)
#endif
    dpsmm1(:,:)   = dpsm (:,:)
    dpsmp1(:,:)   = dpsm (:,:)
    dpslm1(:,:)   = dpsl (:,:)
    dpslp1(:,:)   = dpsl (:,:)
    tlm1  (:,:,:) = tl   (:,:,:)
    tmm1  (:,:,:) = tm   (:,:,:)
    ed1   (:,:,:) = 0.

!PW Physics fields
    allocate ( tmpchunk(pcols,begchunk:endchunk) )
    allocate ( tmpchunk3d(pcols,plevmx,begchunk:endchunk) )
    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,2
         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 physconst, only:
    use commap

    implicit none

#include <comhyb.h>
#include <hadvtest.h>
!
!---------------------------Local workspace-----------------------------
!
    integer i   !
    integer k   ! - indices
    integer 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