inidat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

module inidat

!----------------------------------------------------------------------- 
! 
! Purpose: 
! 
! Method: 
! 
! Author: 
! 
!-----------------------------------------------------------------------

   use precision
   use chemistry, only: chem_init_mix

   real(r8), allocatable :: ps_tmp(:,:)                   
   real(r8), allocatable :: u3_tmp(:,:,:)             
   real(r8), allocatable :: v3_tmp(:,:,:)             
   real(r8), allocatable :: t3_tmp(:,:,:)             
   real(r8), allocatable :: q3_tmp(:,:,:,:) 
   real(r8), allocatable :: qcwat_tmp(:,:,:) 
   real(r8), allocatable :: lcwat_tmp(:,:,:) 
   real(r8), allocatable :: phis_tmp(:,:)        
   real(r8), allocatable :: landfrac_tmp(:,:)                  
   real(r8), allocatable :: landm_tmp(:,:)                  
   real(r8), allocatable :: sgh_tmp(:,:)                 
   real(r8), allocatable :: ts_tmp(:,:)                   
   real(r8), allocatable :: tsice_tmp(:,:)                   
   real(r8), allocatable :: tssub_tmp(:,:,:)         
   real(r8), allocatable :: dpsl_tmp(:,:)                
   real(r8), allocatable :: dpsm_tmp(:,:)                 
   real(r8), allocatable :: vort_tmp(:,:,:)            
   real(r8), allocatable :: div_tmp(:,:,:)             
   real(r8), allocatable :: sicthk_tmp(:,:)               
   real(r8), allocatable :: snowhice_tmp(:,:)                

   real(r8) tmassf_tmp
   real(r8) qmass1_tmp 
   real(r8) qmass2_tmp 
   real(r8) zgsint_tmp 
   real(r8) qmassf_tmp 

contains

   subroutine read_inidat
!-----------------------------------------------------------------------
!
! Purpose:
! Read initial dataset and spectrally truncate as appropriate.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.3 2002/05/02 04:19:39 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

      use pmgrid
      use rgrid, only: nlon
      use comsrf,    only: plevmx,srfflx_state
      use commap, only: w
      use physconst, only: gravit
      use constituents, only: pcnst, pnats, cnst_name, qmin
      use tracers, only: nusr_adv, nusr_nad, ixuadv, ixunad, ixcldw
      implicit none

      include 'netcdf.inc'

#include <comctl.h>
#include <comhyb.h>
#include <comqfl.h>
#include <comlun.h>
#include <perturb.h>
!
! Local workspace
!
      integer i,j,k,m,lat,irow  ! grid and constituent indices
      integer ihem              ! hemisphere index
      real(r8) pdelb(plond,plev) ! pressure diff between interfaces
! using "B" part of hybrid grid only
      real(r8) hyad (plev)      ! del (A)
      real(r8) pssum            ! surface pressure sum
      real(r8) dotproda         ! dot product
      real(r8) dotprodb         ! dot product
      real(r8) pertval          ! perturbation value
      real(r8) zgssum           ! partial sums of phis
!
! Netcdf related variables
!
      integer lonsiz, latsiz, levsiz ! Dimension sizes
      integer londimid, levdimid, latdimid ! Dimension ID's
      integer uid, vid, tid, qid ! Variable ID's
      integer tracid(pcnst+pnats) ! Variable ID's
      integer phisid, sghid, psid ! Variable ID's
      integer landmid
#if ( ! defined COUP_CSM )
      integer tsid, ts1id, ts2id, ts3id, ts4id ,tsiceid! Variable ID's
#endif
#if ( defined COUP_SOM )
      integer sicid
#endif
      integer snowhiceid           ! Variable ID's
      integer landfracid        ! Variable ID's

      integer strt2d(3)         ! start lon, lat, time indices for netcdf 2-d
      integer strt3d(4)         ! start lon, lev, lat, time for netcdf 3-d
      data strt2d/3*1/          ! Only index 2 will ever change
      data strt3d/4*1/          ! Only indices 2,3 will ever change

      integer cnt2d(3)          ! lon, lat, time counts for netcdf 2-d
      integer cnt3d(4)          ! lon, lat, lev, time counts for netcdf 2-d
      data cnt2d/plon,1,1/      ! 2-d arrs: Always grab only a "plon" slice
      data cnt3d/plon,plev,plat,1/ ! 3-d arrs: Always grab a full time slice

      integer ndims2d           ! number of dimensions
      integer dims2d(NF_MAX_VAR_DIMS) ! variable shape
      integer ndims3d           ! number of dimensions
      integer dims3d(NF_MAX_VAR_DIMS) ! variable shape
      integer tmptype
      integer natt, ret, attlen ! netcdf return values
      logical phis_hires        ! true => PHIS came from hi res topo
      real(r8) arr3d(plon,plev,plat)
      character*(NF_MAX_NAME) tmpname
      character*256 text
      character*80 trunits      ! tracer untis
!
!-----------------------------------------------------------------------
! Allocate memory for temporary arrays
!-----------------------------------------------------------------------
!
! Note if not masterproc still might need to allocate array for spmd case
! since each processor calls MPI_scatter 
!
      allocate ( ps_tmp(plond,plat) )                   
      allocate ( u3_tmp(plond,plev,plat) )             
      allocate ( v3_tmp(plond,plev,plat) )             
      allocate ( t3_tmp(plond,plev,plat) )             
      allocate ( q3_tmp(plond,plev,pcnst+pnats,plat) ) 
      allocate ( qcwat_tmp(plond,plev,plat) ) 
      allocate ( lcwat_tmp(plond,plev,plat) ) 
      allocate ( phis_tmp(plond,plat) )        
      allocate ( landm_tmp(plond,plat) )                  
      allocate ( sgh_tmp(plond,plat) )                 
      allocate ( ts_tmp(plond,plat) )                   
      allocate ( tsice_tmp(plond,plat) )                   
      allocate ( tssub_tmp(plond,plevmx,plat) )         
      allocate ( dpsl_tmp(plond,plat) )                
      allocate ( dpsm_tmp(plond,plat) )                 
      allocate ( vort_tmp(plond,plev,plat) )            
      allocate ( div_tmp(plond,plev,plat) )             
      allocate ( sicthk_tmp(plond,plat) )               
      allocate ( snowhice_tmp(plond,plat) )                
      allocate ( landfrac_tmp(plond,plat) )                
!
!-----------------------------------------------------------------------
! Read in input variables
!-----------------------------------------------------------------------

      if (masterproc) then
!
! Get dimension IDs and lengths 
!
         call wrap_inq_dimid  (ncid_ini, 'lat', latdimid)
         call wrap_inq_dimlen (ncid_ini, latdimid, latsiz)
         call wrap_inq_dimid  (ncid_ini, 'lev', levdimid)
         call wrap_inq_dimlen (ncid_ini, levdimid, levsiz)
         call wrap_inq_dimid  (ncid_ini, 'lon', londimid)
         call wrap_inq_dimlen (ncid_ini, londimid, lonsiz)
!
! Get variable id's 
! Check that all tracer units are in mass mixing ratios
!
         call wrap_inq_varid (ncid_ini, 'U'   , uid)
         call wrap_inq_varid (ncid_ini, 'V'   , vid)
         call wrap_inq_varid (ncid_ini, 'T'   , tid)
         call wrap_inq_varid (ncid_ini, 'Q'   , qid)
         call wrap_inq_varid (ncid_ini, 'PS'  , psid)
         call wrap_inq_varid (ncid_ini, 'PHIS', phisid)
         call wrap_inq_varid (ncid_ini, 'SGH' , sghid)
         call wrap_inq_varid (ncid_ini, 'LANDM', landmid)
#if ( ! defined COUP_CSM )
!
! For land-fraction check if the variable name LANDFRAC is on the dataset if not assume FLAND
!
         if ( nf_inq_varid(ncid_ini, 'LANDFRAC', landfracid ) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'LANDFRAC', landfracid)
         else
            call wrap_inq_varid (ncid_ini, 'FLAND', landfracid)
         end if
         call wrap_inq_varid (ncid_ini, 'TS', tsid)
         call wrap_inq_varid (ncid_ini, 'TSICE', tsiceid)
         call wrap_inq_varid (ncid_ini, 'TS1', ts1id)
         call wrap_inq_varid (ncid_ini, 'TS2', ts2id)
         call wrap_inq_varid (ncid_ini, 'TS3', ts3id)
         call wrap_inq_varid (ncid_ini, 'TS4', ts4id)
         call wrap_inq_varid (ncid_ini, 'SNOWHICE', snowhiceid)
#if ( defined COUP_SOM )
         call wrap_inq_varid (ncid_ini, 'SICTHK', sicid)
#endif
#endif
         if (readtrace) then
            do m=2,pcnst+pnats
               call wrap_inq_varid (NCID_INI,cnst_name(m), tracid(m))
               call wrap_get_att_text (NCID_INI,tracid(m),'units', &
                  trunits)
               if (trunits(1:5) .ne. 'KG/KG' .and. &
                  trunits(1:5) .ne. 'kg/kg') then
                  write(6,*)'INIDAT: tracer units for tracer = ', &
                     cnst_name(m),' must be in KG/KG'
                  call endrun
               endif
            end do
         end if
!
! Check dimension ordering for one 2-d and one 3-d field.
! Assume other arrays of like rank will have dimensions ordered the same.
!
         call wrap_inq_var (ncid_ini, uid, tmpname, tmptype, &
            ndims3d, dims3d, natt)
         if (dims3d(1).ne.londimid .or. dims3d(2).ne.levdimid .or. &
            dims3d(3).ne.latdimid .or. ndims3d.gt.4) then
            write(6,*)'INIDAT: Bad number of dims or ordering on 3d fld'
            call endrun
         end if
         call wrap_inq_var (ncid_ini, psid, tmpname, tmptype, &
            ndims2d, dims2d, natt)
         if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. &
            ndims2d.gt.3) then
            write(6,*)'INIDAT: Bad number of dims or ordering on 2d fld'
            call endrun
         end if
!
! Check for presence of 'from_hires' attribute to decide whether to filter
!
         ret = nf_inq_attlen (ncid_ini, phisid, 'from_hires', attlen)
         if (ret.eq.NF_NOERR .and. attlen.gt.256) then
            write(6,*)'INIDAT: from_hires attribute length is too long'
            call endrun
         end if
         ret = nf_get_att_text (ncid_ini, phisid, 'from_hires', text)
         if (ret.eq.NF_NOERR .and. text(1:4).eq.'true') then
            phis_hires = .true.
            write(6,*)'INIDAT: Will filter input PHIS: attribute ',  &
               'from_hires is true'
         else
            phis_hires = .false.
            write(6,*)'INIDAT: Will not filter input PHIS: attribute ', &
               'from_hires is either false or not present'
         end if
!
! Read in 2d fields.  
! For stand alone run: get surface temp and 4 (sub)surface temp fields
!
         do j=1,plat
            strt2d(2) = j
            if (ideal_phys .or. aqua_planet) then
               do i=1,nlon(j)
                  phis_tmp(i,j) = 0.
                  sgh_tmp (i,j) = 0.
               end do
            else
               call wrap_get_vara_realx (ncid_ini, phisid, strt2d, cnt2d, phis_tmp(1,j))
               call wrap_get_vara_realx (ncid_ini, sghid , strt2d, cnt2d, sgh_tmp(1,j))
            end if
            call wrap_get_vara_realx(ncid_ini, landmid, strt2d, cnt2d, landm_tmp(1,j))
            call wrap_get_vara_realx(ncid_ini, psid, strt2d, cnt2d, ps_tmp(1,j))
#if ( ! defined COUP_CSM )
            if (aqua_planet) then
               do i=1,nlon(j)
                  landfrac_tmp(i,j) = 0.
               end do
            else
               call wrap_get_vara_realx (ncid_ini, landfracid, strt2d, cnt2d, landfrac_tmp(1,j))
            endif
            call wrap_get_vara_realx (ncid_ini, tsid, strt2d, cnt2d, ts_tmp(1,j))
            call wrap_get_vara_realx (ncid_ini, tsiceid, strt2d, cnt2d, tsice_tmp(1,j))
            call wrap_get_vara_realx (ncid_ini, ts1id, strt2d, cnt2d, tssub_tmp(1,1,j))
            call wrap_get_vara_realx (ncid_ini, ts2id, strt2d, cnt2d, tssub_tmp(1,2,j))
            call wrap_get_vara_realx (ncid_ini, ts3id, strt2d, cnt2d, tssub_tmp(1,3,j))
            call wrap_get_vara_realx (ncid_ini, ts4id, strt2d, cnt2d, tssub_tmp(1,4,j))
!
! Set sea-ice thickness and snow cover:
!
#if ( defined COUP_SOM )
            call wrap_get_vara_realx(ncid_ini, sicid, strt2d, cnt2d, sicthk_tmp(1,j))
#endif
            call wrap_get_vara_realx(ncid_ini, snowhiceid, strt2d, cnt2d, snowhice_tmp(1,j))
#endif
         end do
!
! Read in 3d fields.  
! Copies are done instead of reading directly into 
! prognostic arrays to address netcdf slowness on Cray.
! Array syntax would be really nice here.
! Initialize tracers if not read in from input data.
! Initialize all user tracers (advected and non-advectec to 0.)
!
         call wrap_get_vara_realx(ncid_ini, uid, strt3d, cnt3d, arr3d)
         u3_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)

         call wrap_get_vara_realx(ncid_ini, vid, strt3d, cnt3d, arr3d)
         v3_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)

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

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

         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(:plon,:plev,:plat)
            end do
         else
            do m=2,pcnst+pnats
               q3_tmp(:plon,:plev,m,:plat) = 0.
            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
!
!-----------------------------------------------------------------------
! Spectrally truncate ps and its derivatives (dpsl and dpsm), phis, 
! u, v, t, vorticity (vort), divergence (div).
!-----------------------------------------------------------------------
!
         call spetru (ps_tmp, phis_tmp, u3_tmp, v3_tmp, t3_tmp, &
                      vort_tmp, div_tmp, dpsl_tmp, dpsm_tmp, phis_hires)
!
! 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 
! The following appears here for consistency with the SLD branch (mvertens).
!
         if (.not. readtrace) then
            do lat=1,plat
               q3_tmp(:plon,:plev,ixcldw,lat) = 0.
               if (nusr_adv .gt. 0) then
                  do m = ixuadv,ixuadv+nusr_adv-1
                     do k=1,plev
                        do i=1,nlon(lat)
                           q3_tmp(i,k,m,lat) = q3_tmp(i,k,1,lat)*10.**(m-ixuadv)
                        end do
                     end do
                  end do
               endif
               if (nusr_nad .gt. 0) then
                  do m = ixunad,ixunad+nusr_nad-1
                     do k=1,plev
                        do i=1,nlon(lat)
                           q3_tmp(i,k,m,lat) = q3_tmp(i,k,1,lat)*10.**(m-ixunad)
                        end do
                     end do
                  end do
               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    , &