#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!