inidat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

module inidat

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

   use precision
   use comspe
   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 :: tl_tmp(:,:,:) 
   real(r8), allocatable :: tm_tmp(:,:,:) 
   real(r8), allocatable :: ql_tmp(:,:,:) 
   real(r8), allocatable :: qm_tmp(:,:,:) 
   real(r8), allocatable :: phis_tmp(:,:)        
   real(r8), allocatable :: phisl_tmp(:,:)        
   real(r8), allocatable :: phism_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 :: 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.22.4.3.4.1 2002/05/17 17:04:42 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

    use pmgrid
    use pspect
    use rgrid
    use comsrf,    only: plevmx
    use commap
    use physconst, only: rair, 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) phi(2,psp/2)                ! used in spectral truncation of phis
    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
    real(r8) tmp1                        ! tmp space
    integer ii                           ! index
!
! 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#endif
    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 ( tl_tmp(plond,plev,plat) )
    allocate ( tm_tmp(plond,plev,plat) )
    allocate ( ql_tmp(plond,plev,plat) )
    allocate ( qm_tmp(plond,plev,plat) )
    allocate ( phis_tmp(plond,plat) )        
    allocate ( phisl_tmp(plond,plat) )        
    allocate ( phism_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 ( 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, divergence (div).
!-----------------------------------------------------------------------
!
       call spetru (ps_tmp   ,phis_tmp  ,u3_tmp  ,v3_tmp  ,t3_tmp   , &
                    q3_tmp   ,div_tmp   ,dpsl_tmp,dpsm_tmp,tl_tmp   , &
                    tm_tmp   ,ql_tmp    ,qm_tmp  ,phi     ,phisl_tmp, &
                    phism_tmp,phis_hires)
!
! For sld do not use the spectrally truncated q3 reread the q field
!
       call wrap_get_vara_realx(ncid_ini, qid, strt3d, cnt3d, arr3d)
       q3_tmp(:plon,:plev,1,:plat) = arr3d(:plon,:plev,:plat)
!
! 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
                do m = ixuadv,ixuadv+nusr_adv-1
                   q3_tmp(:nlon(lat),:plev,m,lat) = &
                                            q3_tmp(:nlon(lat),:plev,1,lat)*10.**(m-ixuadv)
                end do
             endif
             if (nusr_nad .gt. 0) then
                do m = ixunad,ixunad+nusr_nad-1
                   q3_tmp(:nlon(lat),:plev,m,lat) = &
                                            q3_tmp(:nlon(lat),:plev,1,lat)*10.**(m-ixunad)
                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    , &
                     pcnst+pnats,qmin ,q3_tmp(1,1,1,lat))
       end do
!
! Compute ln(Ps*) (Ritchie & Tanguay, 1995) in spectral space
!
       tmp1 = 1./(rair*t0(plev))
       do ii = 1,psp/2
          i = 2*ii - 1
          lnpstar(i  ) = -phi(1,ii)*tmp1
          lnpstar(i+1) = -phi(2,ii)*tmp1
       end do
!
!-----------------------------------------------------------------------
! Integrals of mass, moisture and geopotential height
!-----------------------------------------------------------------------
!
! Compute pdel from "A" portion of hybrid vertical grid
!