inidat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

         do k = 1, plev
            do j = 1, plat
               do i = 1, plon
                  v3s_tmp(i,j,k) = splon_arr3d(i,k,j)
               enddo
            enddo
         enddo


         call wrap_get_vara_realx(ncid_ini, tid, strt3d, cnt3d, arr3d)
         t3_tmp(:plon,:plev,:plat) = arr3d(:,:,:)

         call wrap_get_vara_realx(ncid_ini, qid, strt3d, cnt3d, arr3d)
         q3_tmp(:plon,:plev,1,:plat) = arr3d(:,:,:)

! Initialize tracers if not read in from input data.
! Initialize all user tracers (advected and non-advectec to 0.)

         if (readtrace) then
            do m=2,pcnst+pnats
               call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arr3d)
               q3_tmp(:plon,:plev,m,:plat) = arr3d(:,:,:)
            end do
         endif
!         
! Add random perturbation to temperature if required
!
         if (pertlim.ne.0.0) then
            write(6,*)'INIDAT: Adding random perturbation bounded by +/-', &
                      pertlim,' to initial temperature field'
            do lat=1,plat
               do k=1,plev
                  do i=1,nlon(lat)
                     call random_number (pertval)
                     pertval = 2.*pertlim*(0.5 - pertval)
                     t3_tmp(i,k,lat) = t3_tmp(i,k,lat)*(1. + pertval)
                  end do
               end do
            end do
         endif
!
! Initialize tracers if not read in from input data.
! Initialize all user tracers (advected and non-advectec to 0.)
! Ensure sufficient constituent concentration at all gridpoints 
!
         if (.not. readtrace) then
            do lat=1,plat
               q3_tmp(:plon,:plev,ixcldw,lat) = 0.
               if (nusr_adv .gt. 0) then
                  q3_tmp(:,:plev,ixuadv:ixuadv+nusr_adv-1,lat) = 0.
               endif
               if (nusr_nad .gt. 0) then
                  q3_tmp(:,:plev,ixunad:ixunad+nusr_nad-1,lat) = 0.
               end if
               if (trace_gas) then
                  if (doRamp_ghg ) call ramp_ghg
                  call chem_init_mix(lat, ps_tmp(1,lat), q3_tmp(1,1,1,lat), nlon(lat))
               endif
               if (trace_test1 .or. trace_test2 .or. trace_test3) then
                  call initesttr( q3_tmp(1,1,1,lat) ,nlon(lat))
               endif
            end do
         endif

         do lat=1,plat
            call qneg3('INIDAT  ',lat   ,nlon(lat),plond   ,plev    , &
                       pcnst+pnats,qmin ,q3_tmp(1,1,1,lat))
         end do

         zgsint_tmp = 0.
!
! Compute integrals of mass, moisture, and geopotential height
!
!gg  Integrals of mass and moisture should be unnecessary in Lin-Rood dynamics
!gg  because they are conserved. What's left is the global geopotential...
!gg  Dunno if that's necessary or not, so I left it in.

         do lat = 1, plat
!              
! Accumulate average mass of atmosphere
!
            zgssum = 0.
            do i=1,nlon(lat)
               zgssum = zgssum + phis_tmp(i,lat)
            end do
            zgsint_tmp = zgsint_tmp + w(lat)*zgssum/nlon(lat)
         end do                  ! end of latitude loop
!
! Normalize average height
!
         zgsint_tmp = zgsint_tmp*.5/gravit
!
! Globally avgd sfc. partial pressure of dry air (i.e. global dry mass):
!
! SJL:
         tmass0 = 98222./gravit
!
! WS: Average pole information moved here: treat the global arrays
!
!-----------------------------------------------------------
! Average T, PS, PHIS and Q at the poles.       The initial
! conditions *should* already have these variables averaged,
! but do it here for safety -- no harm if it's already done.
!-----------------------------------------------------------

         call xpavg(phis_tmp(1, 1), plon)
         call xpavg(phis_tmp(1,plat), plon)
         call xpavg(ps_tmp(1,   1), plon)
         call xpavg(ps_tmp(1,plat), plon)

!$omp parallel do private(i, j, k, ic)
         do k = 1, plev
   
            call xpavg(t3_tmp(1,k,   1), plon)
            call xpavg(t3_tmp(1,k,plat), plon)

            do ic = 1, pcnst+pnats
               call xpavg(q3_tmp(1,k,ic,   1),plon)
               call xpavg(q3_tmp(1,k,ic,plat),plon)
            enddo

         enddo

         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('*')/)
      endif                     ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
      call copy_inidat
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
      deallocate ( ps_tmp )
      deallocate ( u3s_tmp )
      deallocate ( v3s_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 ( sicthk_tmp )
      deallocate ( snowhice_tmp )
!
      return
!EOC
   end subroutine read_inidat
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
!BOP
! !IROUTINE: copy_inidat --- Copy temporary arrays to model arrays 
!
! !INTERFACE: 
   subroutine copy_inidat

! !USES:
      use precision
      use pmgrid
      use prognostics
      use buffer
      use comsrf
      use phys_grid
      use tracers, only: ixcldw
      use phys_grid, only: get_ncols_p
#if ( defined SPMD )
      use mpishorthand
      use spmd_dyn, only : comm_y, comm_z
      use parutilitiesmodule, only: parcollective2d, BCSTOP
#endif

      implicit none
!------------------------------Commons----------------------------------
#include <comqfl.h>

! !DESCRIPTION:
!
! Copy temporary arrays to model arrays 
! note that the use statements below contain the definitions
! of the model arrays
!
! !REVISION HISTORY:
!
!   00.06.01   Grant      First attempt at modifying for LRDC
!   00.10.01   Lin        Various revisions
!   00.12.02   Sawyer     Use PILGRIM to scatter data sets
!   01.03.26   Sawyer     Added ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
      real(r8), allocatable :: tmpchunk3d(:,:,:)             
      real(r8), allocatable :: tmpchunk(:,:)             
      integer i, j, ic, k, m
      integer n
      integer ncol
      real(r8) :: pmx, pmn

!-----------------------------------------------------------------------

#if ( defined SPMD )
! dynamics variables
      if (myid_z .eq. 0) then
         call scatter( ps_tmp, strip2d, ps, comm_y )
         call scatter( phis_tmp, strip2d, phis, comm_y )
      endif
#if defined ( TWOD_YZ )
      call parcollective2d( comm_z, BCSTOP, plon, endlat-beglat+1, ps ) 
      call parcollective2d( comm_z, BCSTOP, plon, endlat-beglat+1, phis ) 
#endif

! Warning: beglat may be invalid for different sized staggered grids. - gg

! dynamics variables
      allocate( uv_local(plon,beglat:endlat,beglev:endlev) )
      call scatter( u3s_tmp, strip3dxyz, uv_local, mpicom )
!$omp parallel do private(i,j,k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               u3s(i,j,k) = uv_local(i,j,k)
            enddo
         enddo
      enddo

      call scatter( v3s_tmp, strip3dxyz, uv_local, mpicom )
!$omp parallel do private(i,j,k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               v3s(i,j,k) = uv_local(i,j,k)
            enddo
         enddo
      enddo
      deallocate(uv_local)

      call scatter( t3_tmp, strip3dxzy, t3, mpicom )
      allocate( q3_local(plon,beglev:endlev,pcnst+pnats,beglat:endlat) )

      call scatter( q3_tmp, strip3dq3old, q3_local, mpicom )
      do m=1,pcnst+pnats
!$omp parallel do private(i,j,k)
         do k=beglev,endlev
            do j=beglat,endlat
               do i=1,plon
                  q3(i,j,k,m) = q3_local(i,k,m,j)
               enddo
            enddo
         enddo
      enddo
      deallocate( q3_local )
#else

! dynamics variables
      ps(:,:) = ps_tmp(:,:)
      phis(:,:) = phis_tmp(:,:)

! dynamics variables

!$omp parallel do private(i, j, k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               u3s(i,j,k) = u3s_tmp(i,j,k)
               v3s(i,j,k) = v3s_tmp(i,j,k)
            enddo
         enddo
      enddo

!$omp parallel do private(i, j, k, ic)
      do j=beglat,endlat
         do k=beglev,endlev
            do i=1,plon
               t3(i,k,j) = t3_tmp(i,k,j)
            enddo
         enddo

         do ic=1,pcnst+pnats
            do k=beglev,endlev
               do i=1,plon
                  q3(i,j,k,ic) = q3_tmp(i,k,ic,j)
               enddo
            enddo
         enddo
      enddo

#endif
      allocate ( tmpchunk(pcols,begchunk:endchunk) )
      allocate ( tmpchunk3d(pcols,plevmx,begchunk:endchunk) )
! physics variables
      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
         ncol = get_ncols_p(i)
         srfflx_state2d(i)%ts(:ncol) = tmpchunk(:ncol,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
         ncol = get_ncols_p(i)
         surface_state2d(i)%tssub(:ncol,:) = tmpchunk3d(:ncol,:,i)
      end do
!
!JR cloud and cloud water initialization.  Does this belong somewhere else?
!
      if (masterproc) then
         qcwat_tmp(:,:,:) = q3_tmp(:,:,1,:)
         lcwat_tmp(:,:,:) = q3_tmp(:,:,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
         zgsint = zgsint_tmp
      endif
!
#if ( defined SPMD )
      call mpibcast (zgsint,1,mpir8,0,mpicom)
#endif
      deallocate ( tmpchunk )
      deallocate ( tmpchunk3d)
!EOC
   end subroutine copy_inidat
!-----------------------------------------------------------------------
end module inidat