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_tmp800    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#endifend module inidat