#include <misc.h>
#include <params.h>
#if ( defined SCAM )
#include <max.h>
#endif
module inidat 1,4
!-----------------------------------------------------------------------
!
! Purpose: Read initial dataset and spectrally truncate as appropriate.
!
! Method: Initialize one or a few fields at a time, to minimize the
! memory requirements
!
! Author:
! Modified: P. Worley, to implement initialization of subsets
! of fields. (8/03)
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use comsrf
use abortutils, only: endrun
#if ( defined SCAM )
use pmgrid, only:plev
use getnetcdfdata
use scamMod, only :initlonidx,initlatidx, initTimeIdx,modelfile,analysisfile
#include <runtype.h>
!-----------------------------------------------------------------------
#include <comlun.h>
!-----------------------------------------------------------------------
#include <netcdf.inc>
!-----------------------------------------------------------------------
#endif
real(r8), allocatable :: ps_tmp(:,:)
real(r8), allocatable :: u3_tmp(:,:,:)
real(r8), allocatable :: v3_tmp(:,:,:)
real(r8), allocatable :: t3_tmp(:,:,:)
real(r8), allocatable :: q_tmp(:,:,:)
real(r8), allocatable :: qm_tmp(:,:,:)
real(r8), allocatable :: phis_tmp(:,:)
real(r8), allocatable :: landfrac_tmp(:,:)
real(r8), allocatable :: tsocn_tmp(:,:)
real(r8), allocatable :: icefrac_tmp(:,:)
real(r8), allocatable :: landm_tmp(:,:)
real(r8), allocatable :: pblht_tmp(:,:)
real(r8), allocatable :: tpert_tmp(:,:)
real(r8), allocatable :: qpert_tmp(:,:)
real(r8), allocatable :: cld_tmp(:,:,:)
real(r8), allocatable :: qcwat_tmp(:,:,:)
real(r8), allocatable :: tcwat_tmp(:,:,:)
real(r8), allocatable :: lcwat_tmp(:,:,:)
real(r8), allocatable :: sgh_tmp(:,:)
real(r8), allocatable :: tsice_tmp(:,:)
real(r8), allocatable :: tsice_rad_tmp(:,:)
real(r8), allocatable :: tbot_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(:,:)
#ifdef HADVTEST
real(r8), allocatable :: u_save(:,:,:)
real(r8), allocatable :: v_save(:,:,:)
real(r8), allocatable :: ps_save(:,:)
#endif
real(r8) tmassf_tmp
real(r8) qmass1_tmp
real(r8) qmass2_tmp
real(r8) zgsint_tmp
real(r8) qmassf_tmp
#if ( defined BFB_CAM_SCAM_IOP )
real(r8), allocatable :: ps_sav(:,:)
real(r8), allocatable :: u3_sav(:,:,:)
real(r8), allocatable :: v3_sav(:,:,:)
real(r8), allocatable :: t3_sav(:,:,:)
real(r8), allocatable :: phis_sav(:,:)
#endif
logical read_tsicerad
logical read_tbot
logical read_pblh
logical read_tpert
logical read_qpert
logical read_cloud
logical read_qcwat
logical read_tcwat
logical read_lcwat
#if ( defined SCAM )
logical fill_ends !
real(r8) dplevs( MAX_DATASET_LEVS +1 )
integer dplev
real(r8) coldat(plev )
real(r8) srfdat
integer STATUS
logical have_surfdat
real(r8) surfdat
#endif
contains
subroutine read_inidat(datatype) 1,10
!-----------------------------------------------------------------------
!
! Purpose:
! Read initial dataset and spectrally truncate as appropriate.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
implicit none
integer, optional, intent(in) :: datatype ! in: model or analysis dataset
integer type
#if (defined SCAM )
if (.not.present (datatype)) then
write(6,*)'must specify datatype when calling this routine as part of SCAM'
call endrun
else
type=datatype
endif
call read_inidat_scam(type)
#endif
! Dynamics initialization
call read_inidat_uv
call read_inidat_t
call read_inidat_qps
call read_inidat_phis
! Physics initialization
call read_inidat_phy2D
call read_inidat_cld
call read_inidat_qpert
#ifdef HADVTEST
! HADVTEST reinitialization
call hadvtest_init
#endif
#if ( defined SCAM )
STATUS = NF_CLOSE( NCID_ini )
#endif
return
end subroutine read_inidat
#if (defined SCAM )
!*********************************************************************C
subroutine read_inidat_scam(datatype) 1,2
!-----------------------------------------------------------------------
!
! Purpose:
! Read initial level information for single column version of cam
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use time_manager, only: get_curr_date
implicit none
character*(MAX_PATH_LEN) initdata
real(r8) time( MAX_RECORDS )
integer, intent(in) :: datatype ! in: model or analysis dataset
integer max_recs
parameter(max_recs=MAX_RECORDS)
integer i, delta
integer date( MAX_RECORDS ), datesec( MAX_RECORDS )
integer time_dimID, time_varID,tmpid
integer lev_dimID
integer date_varID, datesec_varID
integer ntime
integer curr_date
integer :: ncsec ! current time of day [seconds]
integer :: yr, mon, day ! year, month, and day component
integer iexcept, ieee_handler, myhandler
logical use_nf_real ! used for 4 byte mode
!
! Trap ieee exceptions on SUN for debugging purposes
!
#if ( defined sun )
iexcept = ieee_handler( 'set', 'common', myhandler )
if ( iexcept .ne. 0 ) write( 6,* ) 'ieee trapping not supported here'
#endif
#if USE_4BYTE_REAL
use_nf_real = .true.
#else
use_nf_real = .false.
#endif
have_surfdat = .false.
!
! Open initial dataset
!
if ( datatype .eq. MODEL ) then
initdata = modelfile
fill_ends = .TRUE.
endif
if ( datatype .eq. ANAL ) then
initdata = analysisfile
fill_ends = .FALSE.
endif
STATUS = NF_OPEN( initdata, NF_NOWRITE, ncid_ini )
if( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:','Cant open model data netcdf file', &
initdata
STATUS = NF_CLOSE( NCID_ini )
return
endif
! -----------------------------------------------------------------
! Get dimension IDs and sizes
! -----------------------------------------------------------------
STATUS = NF_INQ_DIMID ( NCID_INI, 'time', time_dimID )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:', 'Cant get variable dim for time'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_INQ_DIMLEN( NCID_INI, time_dimID, ntime )
if ( ntime .GT. MAX_RECORDS ) then
write( *, * ) 'ERROR: readinitdata.F:',initdata, 'has too many records', &
'It should contain at most', max_recs ,' records'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_INQ_DIMID ( NCID_INI, 'lev', lev_dimID )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:', 'Cant get variable dim for lev'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_INQ_DIMLEN( NCID_INI, lev_dimID, dplev )
! -----------------------------------------------------------------
! Get time variables. Variable called 'time' is current day of year
! -----------------------------------------------------------------
STATUS = NF_INQ_VARID( NCID_INI, 'time', time_varID )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:','Cant get variable ID for time'
STATUS = NF_CLOSE( NCID_INI )
return
endif
if (use_nf_real) then
STATUS = NF_GET_VAR_REAL( NCID_INI, time_varID, time )
else
STATUS = NF_GET_VAR_DOUBLE( NCID_INI, time_varID, time )
endif
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata:','Cant get variable time'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_INQ_VARID( NCID_INI, 'date', date_varID )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:','Cant get variable ID for date'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_GET_VAR_INT( NCID_INI, date_varID, date )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata:','Cant get variable date'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_INQ_VARID( NCID_INI, 'datesec', datesec_varID )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata.F:','Cant get variable ID for datasec'
STATUS = NF_CLOSE( NCID_INI )
return
endif
STATUS = NF_GET_VAR_INT( NCID_INI, datesec_varID, datesec )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata:','Cant get variable datesec'
STATUS = NF_CLOSE( NCID_INI )
return
endif
!
! Find the index into the date array that corresponds to the date that
! is closest to mcdate.
!
initTimeIdx = 1
call get_curr_date(yr,mon,day,ncsec)
curr_date = yr*10000 + mon*100 + day
delta = abs( curr_date - mod( date( 1 )/100, 100 ))
do i = 1, ntime
if ( delta .GT. abs( curr_date - mod( date( i )/100, 100 ))) then
initTimeIdx = i
delta = abs( curr_date - mod( date( i )/100, 100 ))
endif
end do
call getncdata( NCID_INI, initLatIdx, initLonIdx, initTimeIdx, &
'lev', dplevs, STATUS )
if ( STATUS .NE. NF_NOERR ) then
write( 6,* )'ERROR - readinitdata:','Cant get variable lev'
return
endif
return
end subroutine read_inidat_scam
#endif
!
!*********************************************************************C
subroutine read_inidat_uv 1,29
!-----------------------------------------------------------------------
!
! Purpose:
! Read U and V from initial dataset, spectrally truncate,
! and generate VORT and DIV.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use pmgrid
use rgrid, only: nlon
use prognostics
use spetru
#if ( defined SPMD )
use mpishorthand
use spmd_dyn, only: npes, compute_gsfactors
#endif
implicit none
include 'netcdf.inc'
#include <comctl.h>
#include <comhyb.h>
#include <comqfl.h>
#include <comlun.h>
!
! Local workspace
!
integer i,j,k,m,lat,irow ! grid and constituent indices
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer uid, vid ! Variable ID's
integer strt3d(4) ! start lon, lev, lat, time for netcdf 3-d
data strt3d/4*1/ ! Only indices 2,3 will ever change
integer cnt3d(4) ! lon, lat, lev, time counts for netcdf 2-d
data cnt3d/plon,plev,plat,1/ ! 3-d arrs: Always grab a full time slice
integer ndims3d ! number of dimensions
integer dims3d(NF_MAX_VAR_DIMS) ! variable shape
integer tmptype
integer natt, ret, attlen ! netcdf return values
real(r8), pointer, dimension(:,:,:) :: arr3d
character*(NF_MAX_NAME) tmpname
character*256 text
character*80 trunits ! tracer untis
!
! Temporary array copy related variables
!
!
integer, parameter :: iend = i1+plon-1
integer, parameter :: jend = j1+plat-1
integer begj
integer n
#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
!
if (masterproc) then
!
!-----------------------------------------------------------------------
! Allocate memory for temporary U and V arrays
!-----------------------------------------------------------------------
!
allocate ( u3_tmp(plond,plev,plat) )
allocate ( v3_tmp(plond,plev,plat) )
#if ( defined BFB_CAM_SCAM_IOP )
allocate ( u3_sav(plond,plev,plat) )
allocate ( v3_sav(plond,plev,plat) )
#endif
!
! 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
!
call wrap_inq_varid (ncid_ini, 'U' , uid)
call wrap_inq_varid (ncid_ini, 'V' , vid)
!
! Check dimension ordering for U field.
! Assume V field will have dimensions ordered the same.
!
call wrap_inq_var (ncid_ini, uid, tmpname, tmptype, ndims3d, dims3d, natt)
#if ( defined SCAM )
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'U', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, u3_tmp(1,:,1), STATUS )
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'V', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, v3_tmp(1,:,1), STATUS )
!
! Allocate memory for temporary DIV and VORT arrays
!
allocate ( vort_tmp(plond,plev,plat) )
allocate ( div_tmp(plond,plev,plat) )
vort_tmp(:,:,:)=0.
div_tmp(:,:,:)=0.
#else
if (dims3d(1).ne.londimid .or. dims3d(2).ne.levdimid .or. &
dims3d(3).ne.latdimid .or. ndims3d.gt.4) then
call endrun ('INIDAT: Bad number of dims or ordering on 3d fld')
end if
!
! Allocate work space for reading in fields
!
allocate ( arr3d(plon,plev,plat) )
!
! Read in fields.
! Copies are done instead of reading directly into
! prognostic arrays for netcdf performance.
!
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)
!
! Release work space
!
deallocate ( arr3d )
!
! Allocate memory for temporary DIV and VORT arrays
!
allocate ( vort_tmp(plond,plev,plat) )
allocate ( div_tmp(plond,plev,plat) )
!
!-----------------------------------------------------------------------
! Spectrally truncate u, v, vorticity (vort), divergence (div).
!-----------------------------------------------------------------------
!
#if ( defined BFB_CAM_SCAM_IOP )
u3_sav(:plon,:plev,:plat)=u3_tmp(:plon,:plev,:plat)
v3_sav(:plon,:plev,:plat)=v3_tmp(:plon,:plev,:plat)
call spetru_uv (u3_sav, v3_sav, vort_tmp, div_tmp)
#else
call spetru_uv (u3_tmp, v3_tmp, vort_tmp, div_tmp)
#endif
#endif
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
#if ( defined SPMD )
begj = beglatex + numbnd
numperlat = plndlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
!
! Scatter U, V, VORT, DIV
!
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 (vort_tmp, numsend, displs, mpir8, vort(1,1,beglat,1), numrecv, mpir8, 0, mpicom)
call mpiscatterv (div_tmp, numsend, displs, mpir8, div(1,1,beglat,1), numrecv, mpir8, 0, mpicom)
#else
u3(i1:iend,:,j1:jend,1) = u3_tmp(:plon,:,:)
v3(i1:iend,:,j1:jend,1) = v3_tmp(:plon,:,:)
vort(:plon,:,:,1)= vort_tmp(:plon,:,:)
div(:plon,:,:,1) = div_tmp(:plon,:,:)
#endif
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
deallocate ( u3_tmp )
deallocate ( v3_tmp )
#if ( defined BFB_CAM_SCAM_IOP )
deallocate ( u3_sav )
deallocate ( v3_sav )
#endif
deallocate ( vort_tmp )
deallocate ( div_tmp )
#if ( defined SPMD )
else ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays, allocating temporary arrays
! as necessary to make mpiscatterv happy
!-----------------------------------------------------------------------
!
begj = beglatex + numbnd
numperlat = plndlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
!
! Scatter U
!
allocate ( u3_tmp(plond,plev,plat) )
call mpiscatterv (u3_tmp, numsend, displs, mpir8, u3(i1,1,begj,1), numrecv, mpir8, 0, mpicom)
deallocate ( u3_tmp )
!
! Scatter V
!
allocate ( v3_tmp(plond,plev,plat) )
call mpiscatterv (v3_tmp, numsend, displs, mpir8, v3(i1,1,begj,1), numrecv, mpir8, 0, mpicom)
deallocate ( v3_tmp )
!
! Scatter VORT
!
allocate ( vort_tmp(plond,plev,plat) )
call mpiscatterv (vort_tmp, numsend, displs, mpir8, vort(1,1,beglat,1), numrecv, mpir8, 0, mpicom)
deallocate ( vort_tmp )
!
! Scatter DIV
!
allocate ( div_tmp(plond,plev,plat) )
call mpiscatterv (div_tmp, numsend, displs, mpir8, div(1,1,beglat,1), numrecv, mpir8, 0, mpicom)
deallocate ( div_tmp )
!
#endif
endif ! end of if-not-masterproc
!
return
end subroutine read_inidat_uv
!
!*********************************************************************C
subroutine read_inidat_t 1,36
!-----------------------------------------------------------------------
!
! Purpose:
! Read T, TCWAT, and TBOT from initial dataset.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use prognostics
use buffer
use phys_grid
use pmgrid
use rgrid, only: nlon
use phys_buffer, only: pbuf, pbuf_times, pbuf_get_fld_idx
use spetru
#if ( defined SPMD )
use mpishorthand
use spmd_dyn, only: npes, compute_gsfactors
#endif
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
real(r8) pertval ! perturbation value
real(r8), allocatable :: tmpchunk(:,:)
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer tid ! Variable ID's
integer tcwatid
#if ( ! defined COUP_CSM )
integer tbotid ! Variable ID's
#endif
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
real(r8), pointer, dimension(:,:,:) :: arr3d
character*(NF_MAX_NAME) tmpname
character*256 text
character*80 trunits ! tracer untis
!
! Temporary array copy related variables
!
!
integer, parameter :: iend = i1+plon-1
integer, parameter :: jend = j1+plat-1
integer begj
integer n
real(r8), pointer, dimension(:,:,:,:) :: tcwat
#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
!
! logical flags to track which "extra" fields are indeed in IC file
!
read_tbot = .false.
read_tcwat = .false.
!
if (masterproc) then
!
!-----------------------------------------------------------------------
! Allocate memory for temporary arrays
!-----------------------------------------------------------------------
!
allocate ( t3_tmp(plond,plev,plat) )
#if ( defined BFB_CAM_SCAM_IOP )
allocate ( t3_sav(plond,plev,plat) )
#endif
allocate ( tcwat_tmp(plond,plev,plat) )
allocate ( tbot_tmp (plond,plat) )
allocate ( tmpchunk(pcols,begchunk:endchunk) )
!
! 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 and check dimensioning for T and TBOT fields.
! Assume TCWAT field will have dimensions ordered the same as T
!
call wrap_inq_varid (ncid_ini, 'T' , tid)
call wrap_inq_var (ncid_ini, tid, tmpname, tmptype, ndims3d, dims3d, natt)
#if ( ! defined SCAM )
if (dims3d(1).ne.londimid .or. dims3d(2).ne.levdimid .or. &
dims3d(3).ne.latdimid .or. ndims3d.gt.4) then
call endrun ('INIDAT: Bad number of dims or ordering on 3d fld')
end if
#endif
#if ( ! defined COUP_CSM )
if ( nf_inq_varid(ncid_ini, 'TBOT', tbotid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'TBOT' , tbotid)
read_tbot = .true.
call wrap_inq_var (ncid_ini, tbotid, tmpname, tmptype, ndims2d, dims2d, natt)
#if ( ! defined SCAM )
if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. ndims2d.gt.3) then
call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
end if
#endif
end if
#endif
if ( nf_inq_varid(ncid_ini, 'TCWAT', tcwatid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'TCWAT', tcwatid)
read_tcwat = .true.
end if
!
! Allocate more local work space
!
allocate ( arr3d(plon,plev,plat) )
!
!-----------------------------------------------------------------------
! Read in T field.
! A copy is done instead of reading directly into
! t3_tmp for netcdf performance.
!-----------------------------------------------------------------------
!
#if ( defined SCAM )
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'T', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, t3_tmp(1,:,1), STATUS )
#else
call wrap_get_vara_realx(ncid_ini, tid, strt3d, cnt3d, arr3d)
t3_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
!
! 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 t.
!-----------------------------------------------------------------------
!
#if ( defined BFB_CAM_SCAM_IOP )
t3_sav(:plon,:plev,:plat)=t3_tmp(:plon,:plev,:plat)
call spetru_t( t3_sav)
#else
call spetru_t( t3_tmp)
#endif
#endif
!
!-----------------------------------------------------------------------
! Read in TCWAT field.
! A copy is done instead of reading directly into
! tcwat_tmp for netcdf performance.
!-----------------------------------------------------------------------
!
#if ( defined SCAM )
if (read_tcwat) then
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TCWAT', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, tcwat_tmp(1,:,1), STATUS )
else
tcwat_tmp(:plon,:,:) = t3_tmp(:plon,:,:)
write(6,*) 'Warning: TCWAT not found on IC file; initialized with T'
end if
#else
if (read_tcwat) then
call wrap_get_vara_realx(ncid_ini, tcwatid, strt3d, cnt3d, arr3d)
tcwat_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
else
tcwat_tmp(:plon,:,:) = t3_tmp(:plon,:,:)
write(6,*) 'Warning: TCWAT not found on IC file; initialized with T'
end if
#endif
!
! Release work space
!
deallocate ( arr3d )
!
!-----------------------------------------------------------------------
! Read in TBOT field.
!-----------------------------------------------------------------------
!
#if ( ! defined COUP_CSM )
if (aqua_planet) then
tbot_tmp(:plon,:) = 0.
else
if (read_tbot) then
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TBOT', tbot_tmp, STATUS )
#else
do j=1,plat
strt2d(2) = j
call wrap_get_vara_realx (ncid_ini, tbotid , strt2d, cnt2d, tbot_tmp (1,j))
end do
#endif
endif
end if
#endif
if (.not. read_tbot) then
tbot_tmp(:plon,:) = t3_tmp(:plon,plev,:)
write(6,*) 'Warning: TBOT not found on IC file; initialized with lowest level of T'
endif
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
!
! Scatter T
!
#if ( defined SPMD )
begj = beglatex + numbnd
numperlat = plndlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
call mpiscatterv (t3_tmp, numsend, displs, mpir8, t3(i1,1,begj,1), numrecv, mpir8, 0, mpicom)
#else
t3(i1:iend,:,j1:jend,1) = t3_tmp(:plon,:,:)
#endif
!
! Scatter TCWAT
!
! (cloud water initialization should be done in their own packages. Do it
! here for now since moving it will change answers.)
m = pbuf_get_fld_idx('TCWAT')
tcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
call scatter_field_to_chunk(1,plev,1,plond,tcwat_tmp,tcwat(:,:,:,1))
!
! Scatter TBOT
!
call scatter_field_to_chunk(1,1,1,plond,tbot_tmp,tmpchunk)
do i =begchunk,endchunk
surface_state2d(i)%tbot(:) = tmpchunk(:,i)
end do
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
deallocate ( t3_tmp )
#if ( defined BFB_CAM_SCAM_IOP )
deallocate ( t3_sav )
#endif
deallocate ( tcwat_tmp )
deallocate ( tbot_tmp )
deallocate ( tmpchunk )
!
#if ( defined SPMD )
else ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays, allocating temporary arrays
! as necessary to make mpiscatterv happy
!-----------------------------------------------------------------------
!
begj = beglatex + numbnd
numperlat = plndlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
!
! Scatter T
!
allocate ( t3_tmp(plond,plev,plat) )
call mpiscatterv (t3_tmp, numsend, displs, mpir8, t3(i1,1,begj,1), numrecv, mpir8, 0, mpicom)
deallocate ( t3_tmp )
!
! Scatter TCWAT
!
m = pbuf_get_fld_idx('TCWAT')
tcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
allocate ( tcwat_tmp(plond,plev,plat) )
call scatter_field_to_chunk(1,plev,1,plond,tcwat_tmp,tcwat(:,:,:,1))
deallocate ( tcwat_tmp )
!
! Scatter TBOT
!
allocate ( tbot_tmp (plond,plat) )
allocate ( tmpchunk(pcols,begchunk:endchunk) )
call scatter_field_to_chunk(1,1,1,plond,tbot_tmp,tmpchunk)
do i =begchunk,endchunk
surface_state2d(i)%tbot(:) = tmpchunk(:,i)
end do
deallocate ( tbot_tmp )
deallocate ( tmpchunk )
!
#endif
endif ! end of if-not-masterproc
!
! Finish initializing TCWAT
!
if (pbuf_times > 1) then
do n = 2, pbuf_times
tcwat(:,:,:,n) = tcwat(:,:,:,1)
end do
end if
!
return
!
end subroutine read_inidat_t
!
!*********************************************************************C
subroutine read_inidat_qps 1,71
!-----------------------------------------------------------------------
!
! Purpose:
! Read PS, Q, QCWAT, and LCWAT from initial dataset, spectrally
! truncate PS, generate DPSL and DPSM, and compute integrals of
! mass and mositure.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use prognostics
use buffer
use phys_grid
use pmgrid
use rgrid, only: nlon
use commap, only: w
use physconst, only: gravit
use phys_buffer, only: pbuf, pbuf_times, pbuf_get_fld_idx
use constituents, only: pcnst, pnats, cnst_name, qmin, cnst_read_iv, cnst_get_ind
use chemistry, only: chem_implements_cnst, chem_init_cnst
use aerosol_intr, only: aerosol_implements_cnst, aerosol_init_cnst
use tracers, only: tracers_implements_cnst, tracers_init_cnst
use cldcond, only: cldcond_implements_cnst, cldcond_init_cnst
use spetru
#if ( defined SPMD )
use mpishorthand
use spmd_dyn, only: npes, compute_gsfactors
#endif
implicit none
include 'netcdf.inc'
#include <comctl.h>
#include <comhyb.h>
#include <comqfl.h>
#include <comlun.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
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer tracid(pcnst+pnats) ! Variable ID's
integer psid ! Variable ID's
integer qcwatid
integer lcwatid
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
real(r8), pointer, dimension(:,:,:) :: arr4d
real(r8), pointer, dimension(:,:,:) :: arr3d
real(r8), pointer, dimension(:,:,:) :: arr2d
character*(NF_MAX_NAME) tmpname
character*256 text
character*80 trunits ! tracer untis
!
! Temporary array copy related variables
!
integer, parameter :: iend = i1+plon-1
integer, parameter :: jend = j1+plat-1
integer begj
integer n
integer :: ixcldice, ixcldliq
real(r8), pointer, dimension(:,:,:,:) :: qcwat, lcwat
#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
!
!-----------------------------------------------------------------------
! Initialize Q and tracers
!-----------------------------------------------------------------------
!
! logical flags to track which "extra" fields are indeed in IC file
!
read_qcwat = .false.
read_lcwat = .false.
!
! Allocate memory for temporary arrays
!
allocate ( qm_tmp(plon,plev,beglat:endlat) )
allocate ( arr3d(plon,plev,plat) )
!
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
!
! Guard: Check that "Q" is on IC file
!
if (cnst_read_iv(1) .and. &
nf_inq_varid(ncid_ini, cnst_name(1), tracid(1) ) /= NF_NOERR ) then
write(6,*) 'Error: ',cnst_name(1), ' not found on IC file'
call endrun ('READ_INIDAT_QPS')
end if
do m=1,pcnst+pnats
if (cnst_read_iv(m) .and. &
nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
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
call endrun ('INIDAT: tracer units for tracer = '//cnst_name(m)//' must be in KG/KG')
endif
end if
end do
if ( nf_inq_varid(ncid_ini, 'QCWAT', qcwatid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'QCWAT', qcwatid)
read_qcwat = .true.
end if
if ( nf_inq_varid(ncid_ini, 'LCWAT', lcwatid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'LCWAT', lcwatid)
read_lcwat = .true.
end if
!
! Check dimension ordering for one 3-d field.
! Assume other arrays of like rank will have dimensions ordered the same.
!
call wrap_inq_var (ncid_ini, tracid(1), tmpname, tmptype, ndims3d, dims3d, natt)
#if ( ! defined SCAM )
if (dims3d(1).ne.londimid .or. dims3d(2).ne.levdimid .or. &
dims3d(3).ne.latdimid .or. ndims3d.gt.4) then
call endrun ('INIDAT: Bad number of dims or ordering on 3d fld')
end if
#endif
!
! Allocate local (masterproc-only) work space
!
allocate ( q_tmp(plon,plev,plat) )
allocate ( qcwat_tmp(plond,plev,plat) )
allocate ( lcwat_tmp(plond,plev,plat) )
!
! Read in QCWAT
!
if (read_qcwat) then
#if ( defined SCAM )
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'QCWAT', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, qcwat_tmp(1,:,1), STATUS )
#else
call wrap_get_vara_realx(ncid_ini, qcwatid, strt3d, cnt3d, arr3d)
qcwat_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
#endif
end if
!
! Read in or initialize LCWAT
!
if (read_lcwat) then
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'LCWAT', coldat, STATUS )
lcwat_tmp(1,:,1) = coldat(:)
#else
call wrap_get_vara_realx(ncid_ini, lcwatid, strt3d, cnt3d, arr3d)
lcwat_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
#endif
else
lcwat_tmp = 0.0
end if
!
endif ! end of if-masterproc
!
! Read in other 3d fields.
! Copies are done instead of reading directly into
! prognostic arrays for netcdf performance.
! Constituents not read from initial file are initialized by the package
! that implements them.
!
#if ( defined SPMD )
numperlat = plnlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
begj = beglatex + numbnd
#endif
do m=1,pcnst+pnats
if (masterproc) then
if (cnst_read_iv(m) .and. &
nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
cnst_name(m), coldat, STATUS )
arr3d(1,:,1) = coldat(:)
#else
call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arr3d)
#endif
else
write(6,*) 'Warning: Not reading ',cnst_name(m), ' from IC file.'
arr3d = 0.
if (cldcond_implements_cnst(cnst_name(m))) then
call cldcond_init_cnst(cnst_name(m), arr3d)
write(6,*) ' ', cnst_name(m), ' initialized by "cldcond_init_cnst"'
else if (chem_implements_cnst(cnst_name(m))) then
call chem_init_cnst(cnst_name(m), arr3d)
write(6,*) ' ', cnst_name(m), ' initialized by "chem_init_cnst"'
else if (tracers_implements_cnst(cnst_name(m))) then
call tracers_init_cnst(cnst_name(m), arr3d)
write(6,*) ' ', cnst_name(m), ' initialized by "tracers_init_cnst"'
else if (aerosol_implements_cnst(cnst_name(m))) then
call aerosol_init_cnst(cnst_name(m), arr3d)
write(6,*) ' ', cnst_name(m), ' initialized by "aerosol_init_cnst"'
else
write(6,*) ' ', cnst_name(m), ' set to 0.'
end if
endif
do lat=1,plat
call qneg3('INIDAT ', lat, nlon(lat), plon, plev, 1, &
qmin(m), arr3d(:plon,:plev,lat))
end do
if ( cnst_name(m) == 'Q' ) then
q_tmp = arr3d
if ( .not. read_qcwat) then
qcwat_tmp(:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
write(6,*) 'Warning: QCWAT not found on IC file; initialized with ',cnst_name(m)
endif
else if ( (cnst_name(m) == 'CLDLIQ' .or. cnst_name(m) == 'CLDICE') .and. &
.not. read_lcwat ) then
lcwat_tmp (:plon,:plev,:plat) = lcwat_tmp (:plon,:plev,:plat) + arr3d(:plon,:plev,:plat)
write(6,*) 'Warning: LCWAT not found on IC file; ',cnst_name(m),' added to LCWAT instead'
end if
endif
#if ( defined SPMD )
!
! Scatter Q
!
call mpiscatterv (arr3d, numsend, displs, mpir8, qm_tmp, numrecv, mpir8, 0, mpicom)
do j=0,endlat-beglat
do k=1,plev
do i=1,plon
q3(i1+i-1,k,m,begj+j,1) = qm_tmp(i,k,beglat+j)
enddo
enddo
enddo
#else
do j=1,plat
do k=1,plev
do i=1,plon
q3(i1+i-1,k,m,j1+j-1,1) = arr3d(i,k,j)
enddo
enddo
enddo
#endif
end do
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
if (masterproc) then
!
! Scatter QCWAT
!
m = pbuf_get_fld_idx('QCWAT')
qcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
call scatter_field_to_chunk(1,plev,1,plond,qcwat_tmp,qcwat(:,:,:,1))
!
! Scatter LCWAT
!
m = pbuf_get_fld_idx('LCWAT')
lcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
call scatter_field_to_chunk(1,plev,1,plond,lcwat_tmp,lcwat(:,:,:,1))
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
! (Keep q_tmp for use in global integrals)
!-----------------------------------------------------------------------
!
deallocate ( arr3d )
deallocate ( qm_tmp )
deallocate ( qcwat_tmp )
deallocate ( lcwat_tmp )
!
#if ( defined SPMD )
else ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
deallocate ( arr3d )
deallocate ( qm_tmp )
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays, allocating temporary arrays
! as necessary to make mpiscatterv happy
!-----------------------------------------------------------------------
!
! Scatter QCWAT
!
m = pbuf_get_fld_idx('QCWAT')
qcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
allocate ( qcwat_tmp(plond,plev,plat) )
call scatter_field_to_chunk(1,plev,1,plond,qcwat_tmp,qcwat(:,:,:,1))
deallocate ( qcwat_tmp )
!
! Scatter LCWAT
!
m = pbuf_get_fld_idx('LCWAT')
lcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
allocate ( lcwat_tmp(plond,plev,plat) )
call scatter_field_to_chunk(1,plev,1,plond,lcwat_tmp,lcwat(:,:,:,1))
deallocate ( lcwat_tmp )
#endif
endif ! end of if-not-masterproc
!
! Finish initializing QCWAT and LCWAT
!
if (pbuf_times > 1) then
do n = 2, pbuf_times
qcwat(:,:,:,n) = qcwat(:,:,:,1)
lcwat(:,:,:,n) = lcwat(:,:,:,1)
end do
end if
!
!-----------------------------------------------------------------------
! Initialize PS
!-----------------------------------------------------------------------
!
if (masterproc) then
!
!-----------------------------------------------------------------------
! Allocate memory for temporary arrays
!-----------------------------------------------------------------------
!
#if ( defined BFB_CAM_SCAM_IOP )
allocate ( ps_sav(plond,plat) )
#endif
allocate ( ps_tmp(plond,plat) )
allocate ( dpsl_tmp(plond,plat) )
allocate ( dpsm_tmp(plond,plat) )
!
! 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, 'lon', londimid)
call wrap_inq_dimlen (ncid_ini, londimid, lonsiz)
!
! Get variable id
!
call wrap_inq_varid (ncid_ini, 'PS' , psid)
!
! Check dimension ordering for PS.
!
call wrap_inq_var (ncid_ini, psid, tmpname, tmptype, ndims2d, dims2d, natt)
!
! Read in PS.
!
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'PS', ps_tmp, STATUS )
dpsl_tmp(:,:) = 0.
dpsm_tmp(:,:) = 0.
#else
if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. ndims2d.gt.3) then
call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
end if
do j=1,plat
strt2d(2) = j
call wrap_get_vara_realx(ncid_ini, psid, strt2d, cnt2d, ps_tmp(1,j))
end do
!
!-----------------------------------------------------------------------
! Spectrally truncate ps and its derivatives (dpsl and dpsm).
!-----------------------------------------------------------------------
!
#if ( defined BFB_CAM_SCAM_IOP )
ps_sav(:plon,:plat)=ps_tmp(:plon,:plat)
call spetru_ps (ps_sav, dpsl_tmp, dpsm_tmp)
#else
call spetru_ps (ps_tmp, dpsl_tmp, dpsm_tmp)
#endif
#endif
!
!-----------------------------------------------------------------------
! Prepare for integrals of mass and moisture.
!-----------------------------------------------------------------------
!
! Compute pdel from "A" portion of hybrid vertical grid
!
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.
!
! Compute integrals of mass, moisture, and geopotential height
!
#if ( defined SCAM )
do irow = 1,plat
#else
do irow = 1,plat/2
#endif
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.
do i=1,nlon(lat)
pssum = pssum + ps_tmp (i,lat)
end do
tmassf_tmp = tmassf_tmp + w(irow)*pssum/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 + q_tmp(i,k,lat)*pdela(i,k)
dotprodb = dotprodb + q_tmp(i,k,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
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
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)
#if (! defined SCAM )
!
! 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
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
#if ( defined SPMD )
!
! Scatter PS, DPSL, and DPSM
!
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 (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)
#else
ps(:plon,:,1) = ps_tmp(:plon,:)
dpsl(:plon,:) = dpsl_tmp(:plon,:)
dpsm(:plon,:) = dpsm_tmp(:plon,:)
#endif
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
deallocate ( q_tmp )
deallocate ( ps_tmp )
#if ( defined BFB_CAM_SCAM_IOP )
deallocate ( ps_sav )
#endif
deallocate ( dpsl_tmp )
deallocate ( dpsm_tmp )
!
#if ( defined SPMD )
else ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays, allocating temporary arrays
! as necessary to make mpiscatterv happy
!-----------------------------------------------------------------------
!
numperlat = plond
call compute_gsfactors (numperlat, numrecv, numsend, displs)
!
! Scatter PS
!
allocate ( ps_tmp(plond,plat) )
call mpiscatterv (ps_tmp, numsend, displs, mpir8, ps(1,beglat,1), numrecv, mpir8, 0, mpicom)
deallocate ( ps_tmp )
!
! Scatter DPSL
!
allocate ( dpsl_tmp(plond,plat) )
call mpiscatterv (dpsl_tmp, numsend, displs, mpir8, dpsl(1,beglat), numrecv, mpir8, 0, mpicom)
deallocate ( dpsl_tmp )
!
! Scatter DPSM
!
allocate ( dpsm_tmp(plond,plat) )
call mpiscatterv (dpsm_tmp, numsend, displs, mpir8, dpsm(1,beglat), numrecv, mpir8, 0, mpicom)
deallocate ( dpsm_tmp )
#endif
endif ! end of if-not-masterproc
!
!-----------------------------------------------------------------------
! Broadcast global integrals
!-----------------------------------------------------------------------
!
if (masterproc) then
tmassf = tmassf_tmp
qmass1 = qmass1_tmp
qmass2 = qmass2_tmp
qmassf = qmassf_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)
#endif
!
return
!
end subroutine read_inidat_qps
!
!*********************************************************************C
subroutine read_inidat_phis 1,23
!-----------------------------------------------------------------------
!
! Purpose:
! Read PHIS from initial dynamics dataset and spectrally truncate.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use pmgrid
use rgrid, only: nlon
use commap, only: w
use physconst, only: gravit
use prognostics
use spetru
#if ( defined SPMD )
use mpishorthand
use spmd_dyn, only: npes, compute_gsfactors
#endif
implicit none
include 'netcdf.inc'
#include <comctl.h>
#include <comhyb.h>
#include <comqfl.h>
#include <comlun.h>
!
! Local workspace
!
integer i,j,k,m,lat,irow ! grid and constituent indices
integer ihem ! hemisphere index
real(r8) hyad (plev) ! del (A)
real(r8) zgssum ! partial sums of phis
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer phisid ! Variable ID's
integer strt2d(3) ! start lon, lat, time indices for netcdf 2-d
data strt2d/3*1/ ! Only index 2 will ever change
integer cnt2d(3) ! lon, lat, time counts for netcdf 2-d
data cnt2d/plon,1,1/ ! 2-d arrs: Always grab only a "plon" slice
integer ndims2d ! number of dimensions
integer dims2d(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
character*(NF_MAX_NAME) tmpname
character*256 text
character*80 trunits ! tracer untis
!
! Temporary array copy related variables
!
!
integer n
#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
!
!-----------------------------------------------------------------------
! Allocate memory for temporary arrays
!-----------------------------------------------------------------------
!
allocate ( phis_tmp(plond,plat) )
#if ( defined BFB_CAM_SCAM_IOP )
allocate ( phis_sav(plond,plat) )
#endif
!
if (masterproc) then
!
!-----------------------------------------------------------------------
! Prepare for reading input variables
!-----------------------------------------------------------------------
!
! 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
!
call wrap_inq_varid (ncid_ini, 'PHIS', phisid)
!
! Check dimension orderings for PHIS field.
!
call wrap_inq_var (ncid_ini, phisid, tmpname, tmptype, ndims2d, dims2d, natt)
#if ( ! defined SPMD )
if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. ndims2d.gt.3) then
call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
end if
#endif
!
! 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
call endrun ('INIDAT: from_hires attribute length is too long')
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 PHIS.
!
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'PHIS', phis_tmp , STATUS )
if (ideal_phys .or. aqua_planet) then
do j=1,plat
strt2d(2) = j
do i=1,nlon(j)
phis_tmp(i,j) = 0.
end do
end do
end if
#else
do j=1,plat
strt2d(2) = j
if (ideal_phys .or. aqua_planet) then
do i=1,nlon(j)
phis_tmp(i,j) = 0.
end do
else
call wrap_get_vara_realx (ncid_ini, phisid, strt2d, cnt2d, phis_tmp(1,j))
end if
end do
!
!-----------------------------------------------------------------------
! Spectrally truncate phis.
!-----------------------------------------------------------------------
!
#if ( defined BFB_CAM_SCAM_IOP )
phis_sav(:plon,:plat)=phis_tmp(:plon,:plat)
call spetru_phis (phis_sav, phis_hires)
#else
call spetru_phis (phis_tmp, phis_hires)
#endif
#endif
!
!-----------------------------------------------------------------------
! Prepare for integral of geopotential height
!-----------------------------------------------------------------------
!
! Compute pdel from "A" portion of hybrid vertical grid
!
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 height integral for summation
! in a third calculation loop (assures bit-for-bit compare
! with non-random history tape).
!
zgsint_tmp = 0.
!
! Compute integral of geopotential height
!
#if ( defined SCAM )
do irow = 1,plat
#else
do irow = 1,plat/2
#endif
do ihem=1,2
if (ihem.eq.1) then
lat = irow
else
lat = plat - irow + 1
end if
!
! Accumulate average height
!
zgssum = 0.
do i=1,nlon(lat)
zgssum = zgssum + phis_tmp(i,lat)
end do
zgsint_tmp = zgsint_tmp + w(irow)*zgssum/nlon(lat)
end do
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):
!
write(6,810) zgsint_tmp
810 format(/69('*')/'INIDAT: Globally averaged geopotential ', &
'height = ',f16.10,' meters'/69('*')/)
endif ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
#if ( defined SPMD )
numperlat = plond
call compute_gsfactors (numperlat, numrecv, numsend, displs)
call mpiscatterv (phis_tmp, numsend, displs, mpir8, phis(1,beglat), numrecv, mpir8, 0, mpicom)
#else
phis(:plon,:) = phis_tmp(:plon,:)
#endif
!
! Global integerals
!
if (masterproc) then
zgsint = zgsint_tmp
endif
#if ( defined SPMD )
call mpibcast (zgsint,1,mpir8,0,mpicom)
#endif
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
deallocate ( phis_tmp )
#if ( defined BFB_CAM_SCAM_IOP )
deallocate ( phis_sav )
#endif
!
return
end subroutine read_inidat_phis
!*********************************************************************C
subroutine read_inidat_phy2D 1,47
!-----------------------------------------------------------------------
!
! Purpose:
! Read initial 2D physics fields.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use pmgrid
use rgrid, only: nlon
use commap, only: clat
use physconst, only: gravit
use constituents, only: pcnst, pnats, cnst_name, qmin, cnst_read_iv
use chemistry, only: chem_implements_cnst, chem_init_cnst
use tracers, only: tracers_implements_cnst, tracers_init_cnst
implicit none
include 'netcdf.inc'
#include <comctl.h>
#include <comlun.h>
!
! Local workspace
!
integer i,j,k,m,lat,irow ! grid and constituent indices
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer sghid ! Variable ID's
integer landmid
integer pblhtid
integer tpertid
#if ( ! defined COUP_CSM )
integer ts1id, ts2id, ts3id, ts4id ,tsiceid,tsice_rad_id! Variable ID's
#endif
#if ( defined COUP_SOM )
integer sicid
integer icefracid
integer tsocnid
#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
data cnt2d/plon,1,1/ ! 2-d arrs: Always grab only a "plon" slice
integer ndims2d ! number of dimensions
integer dims2d(NF_MAX_VAR_DIMS) ! variable shape
integer tmptype
integer natt, ret, attlen ! netcdf return values
character*(NF_MAX_NAME) tmpname
character*256 text
!
!-----------------------------------------------------------------------
! 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 ( landm_tmp(plond,plat) )
allocate ( sgh_tmp(plond,plat) )
allocate ( tsice_tmp(plond,plat) )
allocate ( tsice_rad_tmp(plond,plat) )
allocate ( tssub_tmp(plond,plevmx,plat) )
allocate ( sicthk_tmp(plond,plat) )
allocate ( snowhice_tmp(plond,plat) )
allocate ( landfrac_tmp(plond,plat) )
allocate ( pblht_tmp(plond,plat) )
allocate ( tpert_tmp(plond,plat) )
#if ( defined COUP_SOM )
allocate ( icefrac_tmp(plond,plat) )
allocate ( tsocn_tmp(plond,plat) )
#endif
!
!-----------------------------------------------------------------------
! Read in input variables
!-----------------------------------------------------------------------
!
! logical flags to track which "extra" fields are indeed in IC file
!
read_tsicerad = .false.
read_pblh = .false.
read_tpert = .false.
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
!
call wrap_inq_varid (ncid_ini, 'SGH' , sghid)
if (nf_inq_varid (ncid_ini, 'LANDM_COSLAT', landmid) /= nf_noerr) then
write(6,*)'INIDAT: LANDM_COSLAT not found on initial dataset.'
write(6,*)' Need to run definesurf to create it.'
write(6,*)' This field became a requirement as of cam2_0_2_dev43'
call endrun ()
end if
if ( nf_inq_varid(ncid_ini, 'PBLH', pblhtid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'PBLH' , pblhtid)
read_pblh = .true.
end if
if ( nf_inq_varid(ncid_ini, 'TPERT', tpertid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'TPERT', tpertid)
read_tpert = .true.
end if
#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
if ( nf_inq_varid(ncid_ini, 'TSICERAD', tsice_rad_id) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'TSICERAD', tsice_rad_id)
read_tsicerad = .true.
end if
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)
call wrap_inq_varid (ncid_ini, 'ICEFRAC', icefracid)
call wrap_inq_varid (ncid_ini, 'TSOCN', tsocnid)
#endif
#endif
!
! Check dimension ordering for one 2D field.
! Assume other arrays of like rank will have dimensions ordered the same.
!
call wrap_inq_var (ncid_ini, landmid, tmpname, tmptype, ndims2d, dims2d, natt)
#if ( ! defined SCAM )
if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. ndims2d.gt.3) then
call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
end if
#endif
#if ( defined SCAM )
!
! Read in fields.
do j=1,plat
strt2d(2) = j
if (ideal_phys .or. aqua_planet) then
do i=1,nlon(j)
sgh_tmp (i,j) = 0.
end do
else
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'SGH', sgh_tmp, STATUS)
end if
landm_tmp(:,:) = 0.
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'LANDM_COSLAT', landm_tmp, STATUS)
#if ( ! defined COUP_CSM )
if (aqua_planet) then
do i=1,nlon(j)
landfrac_tmp(i,j) = 0.
end do
else
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'LANDFRAC', landfrac_tmp, STATUS)
endif
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TSICE', tsice_tmp, STATUS )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TS1', srfdat, STATUS )
tssub_tmp(1,1,1)=srfdat
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TS2', srfdat, STATUS )
tssub_tmp(1,2,1) = srfdat
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TS3', srfdat, STATUS )
tssub_tmp(1,3,1) = srfdat
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TS4', srfdat, STATUS )
tssub_tmp(1,4,1) = srfdat
if (read_tsicerad) then
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TSICERAD', tsice_rad_tmp, STATUS )
else
tsice_rad_tmp(:plon,j) = tsice_tmp(:plon,j)
if (j .eq. 1) &
write(6,*) 'Warning: TSICERAD not found on IC file; initialized with TSICE'
end if
if (read_pblh) then
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'PBLH', pblht_tmp, STATUS )
else
pblht_tmp (:plon,j) = 0.
if (j .eq. 1) &
write(6,*) 'Warning: PBLH not found on IC file; initialized to 0.'
end if
if (read_tpert) then
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TPERT', tpert_tmp, STATUS )
else
tpert_tmp (:plon,j) = 0.
if (j .eq. 1) &
write(6,*) 'Warning: TPERT not found on IC file; initialized to 0.'
end if
!
! Set sea-ice thickness and snow cover:
!
#if ( defined COUP_SOM )
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'SICTHK', sicthk_tmp, STATUS)
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'ICEFRAC', icefrac_tmp, STATUS)
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'TSOCN', tsocn_tmp, STATUS)
#endif
call getncdata (ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'SNOWHICE', snowhice_tmp, STATUS)
#endif
end do
#else ! if not defined SCAM
!
! Read in fields.
do j=1,plat
strt2d(2) = j
if (ideal_phys .or. aqua_planet) then
do i=1,nlon(j)
sgh_tmp (i,j) = 0.
end do
else
call wrap_get_vara_realx (ncid_ini, sghid , strt2d, cnt2d, sgh_tmp(1,j))
end if
landm_tmp(:plon,j) = 0.
call wrap_get_vara_realx(ncid_ini, landmid, strt2d, cnt2d, landm_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, 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))
if (read_tsicerad) then
call wrap_get_vara_realx (ncid_ini, tsice_rad_id, strt2d, cnt2d, tsice_rad_tmp(1,j))
else
tsice_rad_tmp(:plon,j) = tsice_tmp(:plon,j)
if (j .eq. 1) &
write(6,*) 'Warning: TSICERAD not found on IC file; initialized with TSICE'
end if
if (read_pblh) then
call wrap_get_vara_realx (ncid_ini, pblhtid, strt2d, cnt2d, pblht_tmp(1,j ))
else
pblht_tmp (:plon,j) = 0.
if (j .eq. 1) &
write(6,*) 'Warning: PBLH not found on IC file; initialized to 0.'
end if
if (read_tpert) then
call wrap_get_vara_realx (ncid_ini, tpertid, strt2d, cnt2d, tpert_tmp(1,j ))
else
tpert_tmp (:plon,j) = 0.
if (j .eq. 1) &
write(6,*) 'Warning: TPERT not found on IC file; initialized to 0.'
end if
!
! 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))
call wrap_get_vara_realx(ncid_ini, icefracid, strt2d, cnt2d, icefrac_tmp(1,j))
call wrap_get_vara_realx(ncid_ini, tsocnid, strt2d, cnt2d, tsocn_tmp(1,j))
#endif
call wrap_get_vara_realx(ncid_ini, snowhiceid, strt2d, cnt2d, snowhice_tmp(1,j))
#endif
end do
#endif ! end of ifdef SCAM
endif ! end of if-masterproc
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
call copy_inidat_phy2D
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
deallocate ( landm_tmp )
deallocate ( sgh_tmp )
deallocate ( tsice_tmp )
deallocate ( tsice_rad_tmp )
deallocate ( tssub_tmp )
deallocate ( sicthk_tmp )
deallocate ( snowhice_tmp )
deallocate ( pblht_tmp )
deallocate ( tpert_tmp )
deallocate ( landfrac_tmp )
#if ( defined COUP_SOM )
deallocate ( icefrac_tmp )
deallocate ( tsocn_tmp )
#endif
!
return
end subroutine read_inidat_phy2D
!*********************************************************************C
subroutine copy_inidat_phy2D 1,18
!-----------------------------------------------------------------------
! Copy temporary 2D arrays to model physics arrays
! note that the use statements below contain the definitions
! of the model arrays
!-----------------------------------------------------------------------
use prognostics
use buffer
use phys_grid
use phys_buffer, only: pbuf, pbuf_times, pbuf_get_fld_idx
#if ( defined SPMD )
use mpishorthand
use spmd_dyn, only: npes, compute_gsfactors
#endif
!-----------------------------------------------------------------------
implicit none
!------------------------------Commons----------------------------------
#include <comqfl.h>
!
! Local workspace
!
real(r8), allocatable :: tmpchunk3d(:,:,:)
real(r8), allocatable :: tmpchunk(:,:)
integer n,i,j, m, c
real(r8), pointer, dimension(:,:,:,:) :: qcwat, lcwat
#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
!-----------------------------------------------------------------------
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,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
call scatter_field_to_chunk(1, 1,1,plond,tsice_rad_tmp,tsice_rad(1,begchunk))
call scatter_field_to_chunk(1, 1,1,plond,pblht_tmp,pblht(1 ,begchunk ))
call scatter_field_to_chunk(1, 1,1,plond,tpert_tmp,tpert(1 ,begchunk ))
#if ( defined COUP_SOM )
call scatter_field_to_chunk(1,1,1,plond,sicthk_tmp,sicthk(1,begchunk))
call scatter_field_to_chunk(1,1,1,plond,icefrac_tmp,icefrac(1,begchunk))
call scatter_field_to_chunk(1,1,1,plond,tsocn_tmp,tsocn(1,begchunk))
! define an initial ocean fraction and non-land ice fraction
! The 1st "where" stmt used to be done in update_srf_fractions (dev45)
do c=begchunk,endchunk
where (icefrac(:pcols,c) + landfrac(:pcols,c) > 1.0)
icefrac(:pcols,c) = 1. - landfrac(:pcols,c)
end where
where (landfrac(:pcols,c) < 1.)
aice(:pcols,c) = icefrac(:pcols,c)/(1. - landfrac(:pcols,c))
elsewhere
aice(:pcols,c) = 0.
end where
ocnfrac(:pcols,c) = 1. - landfrac(:pcols,c) - icefrac(:pcols,c)
end do
write(6,*)'INIDAT: ocnfrac=',ocnfrac(1,begchunk)
!
! Master needs global landfrac
!
call gather_chunk_to_field(1,1,1,plon,landfrac,landfrac_field)
! write(6,*)'INIDAT iam=',iam,' landfrac=',landfrac
! write(6,*)'INIDAT iam=',iam,' landfrac_field=',landfrac_field
!
!JR Could read in Focn from initial dataset if available
Focn(:,:) = 0.
#else
Focn(:,:) = inf
frzmlt(:,:) = 0. ! needs to be 0, otherwise test in tstm always true
tsocn(:,:) = inf
#endif
deallocate ( tmpchunk )
deallocate ( tmpchunk3d)
return
end subroutine copy_inidat_phy2D
!
!*********************************************************************C
subroutine read_inidat_cld 1,15
!-----------------------------------------------------------------------
!
! Purpose:
! Read CLD from initial dataset.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use pmgrid
use phys_grid
use phys_buffer, only: pbuf, pbuf_times, pbuf_get_fld_idx
implicit none
include 'netcdf.inc'
#include <comlun.h>
!
! Local workspace
!
integer m,n
real(r8), pointer, dimension(:,:,:,:) :: cld
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer cldid ! variable ID
integer strt3d(4) ! start lon, lev, lat, time for netcdf 3-d
data strt3d/4*1/ ! Only indices 2,3 will ever change
integer cnt3d(4) ! lon, lat, lev, time counts for netcdf 2-d
data cnt3d/plon,plev,plat,1/ ! 3-d arrs: Always grab a full time slice
integer ndims3d ! number of dimensions
integer dims3d(NF_MAX_VAR_DIMS) ! variable shape
integer tmptype
integer natt ! netcdf return value
real(r8), pointer, dimension(:,:,:) :: arr3d
character*(NF_MAX_NAME) tmpname
!
!-----------------------------------------------------------------------
! Allocate memory for temporary array
!-----------------------------------------------------------------------
!
! Note if not masterproc still might need to allocate array for spmd case
! since each processor calls scatter routine
!
allocate ( cld_tmp (plond,plev,plat) )
!
!-----------------------------------------------------------------------
! Read in input variable
!-----------------------------------------------------------------------
!
! logical flags to track whether CLD field is indeed in IC file
!
read_cloud = .false.
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
!
if ( nf_inq_varid(ncid_ini, 'CLOUD', cldid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'CLOUD', cldid )
read_cloud = .true.
!
! Check dimension ordering for CLOUD field.
!
call wrap_inq_var (ncid_ini, cldid, tmpname, tmptype, ndims3d, dims3d, natt)
#if (! defined SCAM )
if (dims3d(1).ne.londimid .or. dims3d(2).ne.levdimid .or. &
dims3d(3).ne.latdimid .or. ndims3d.gt.4) then
call endrun ('INIDAT: Bad number of dims or ordering on 3d fld')
end if
#endif
end if
!
! Allocate local (masterproc-only) work space
!
allocate ( arr3d(plon,plev,plat) )
!
! Read in 3d field.
! Copy is done instead of reading directly into
! field array for netcdf performance.
!
if (read_cloud) then
#if ( defined SCAM )
call getinterpncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'CLOUD', have_surfdat, surfdat, fill_ends, &
dplevs, dplev, cld_tmp(1,:,1), STATUS )
#else
call wrap_get_vara_realx(ncid_ini, cldid , strt3d, cnt3d, arr3d)
cld_tmp (:plon,:plev,:plat) = arr3d(:plon,:plev,:plat)
#endif
else
cld_tmp (:plon,:,:) = 0.
write(6,*) 'Warning: CLOUD not found on IC file; initialized to 0.'
end if
!
! Release work space
!
deallocate ( arr3d )
endif ! end of if-masterproc
m = pbuf_get_fld_idx('CLD')
cld => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
call scatter_field_to_chunk(1,plev,1,plond,cld_tmp,cld(:,:,:,1))
if (pbuf_times > 1) then
do n = 2, pbuf_times
cld (:,:,:,n) = cld (:,:,:,1)
end do
end if
deallocate ( cld_tmp )
return
!
end subroutine read_inidat_cld
!
!*********************************************************************C
subroutine read_inidat_qpert 1,14
!-----------------------------------------------------------------------
!
! Purpose:
! Read QPERT from initial dataset.
!
!-----------------------------------------------------------------------
!
! $Id: inidat.F90,v 1.20.4.18.10.1 2004/05/20 18:36:04 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use pmgrid
use phys_grid
use buffer
implicit none
include 'netcdf.inc'
#include <comlun.h>
!
! Local workspace
!
integer j
real(r8), allocatable :: tmpchunk(:,:)
!
! Netcdf related variables
!
integer lonsiz, latsiz, levsiz ! Dimension sizes
integer londimid, levdimid, latdimid ! Dimension ID's
integer qpertid ! variable ID
integer strt2d(3) ! start lon, lat, time indices for netcdf 2-d
data strt2d/3*1/ ! Only index 2 will ever change
integer cnt2d(3) ! lon, lat, time counts for netcdf 2-d
data cnt2d/plon,1,1/ ! 2-d arrs: Always grab only a "plon" slice
integer ndims2d ! number of dimensions
integer dims2d(NF_MAX_VAR_DIMS) ! variable shape
integer tmptype
integer natt ! netcdf return value
character*(NF_MAX_NAME) tmpname
!
!-----------------------------------------------------------------------
! Initialize chunked data structure
!-----------------------------------------------------------------------
!
qpert(:,:,:) = 0.
#if ( ! defined COUP_CSM )
!
!-----------------------------------------------------------------------
! Allocate memory for temporary array
!-----------------------------------------------------------------------
!
! Note if not masterproc still might need to allocate array for spmd case
! since each processor calls scatter routine
!
allocate ( qpert_tmp(plond,plat) )
allocate ( tmpchunk(pcols,begchunk:endchunk) )
!
!-----------------------------------------------------------------------
! Read in input variable
!-----------------------------------------------------------------------
!
! logical flags to track whether QPERT field is indeed in IC file
!
read_qpert = .false.
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 and check dimension ordering for QPERT field.
!
if ( nf_inq_varid(ncid_ini, 'QPERT', qpertid) == NF_NOERR ) then
call wrap_inq_varid (ncid_ini, 'QPERT', qpertid)
read_qpert = .true.
!
call wrap_inq_var (ncid_ini, qpertid, tmpname, tmptype, ndims2d, dims2d, natt)
#if ( ! defined SCAM )
if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. ndims2d.gt.3) then
call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
end if
#endif
end if
!
! Read in QPERT.
if (read_qpert) then
#if ( defined SCAM )
call getncdata( ncid_ini, initLatIdx, initLonIdx, initTimeIdx, &
'QPERT', qpert_tmp, STATUS )
#else
do j=1,plat
strt2d(2) = j
call wrap_get_vara_realx (ncid_ini, qpertid, strt2d, cnt2d, qpert_tmp(1,j))
end do
#endif
else
qpert_tmp (:plon,:) = 0.
write(6,*) 'Warning: QPERT not found on IC file; initialized to 0.'
end if
end if
call scatter_field_to_chunk(1,1,1,plond,qpert_tmp,tmpchunk)
do j =begchunk,endchunk
qpert(:,1,j) = tmpchunk(:,j)
end do
deallocate ( qpert_tmp )
deallocate ( tmpchunk )
#endif
return
!
end subroutine read_inidat_qpert
!
!*********************************************************************C
#ifdef HADVTEST
!-----------------------------------------------------------------------
!
! Purpose:
! Initialize fields for solid body rotation test.
!
! Author: J. Rosinski
! Modified: P. Worley (for low memory version of inidat)
!
!-----------------------------------------------------------------------
subroutine hadvtest_init 1,5
use pmgrid
use rgrid
use commap
use prognostics
use dynconst, only: rearth
implicit none
#include <comhyb.h>
#include <hadvtest.h>
!
!---------------------------Local workspace-----------------------------
!
integer i,j,k,lon,lat
real(r8) h0, u0, small_r, big_r, theta, theta_c, lambda, lambda_c
real(r8) alfa, dlam, pie
!
!-----------------------------------------------------------------------
!
!
! First: allocate save fields
!
allocate( u_save(plond,plev,beglat:endlat) )
allocate( v_save(plond,plev,beglat:endlat) )
allocate( ps_save(plond,beglat:endlat) )
!
! Next: zero sgh and phis fields
!
sgh(:,:) = 0.
phis(:,:) = 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=beglat,endlat
j = j1 - 1 + lat
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 lon=1,nlon(lat)
i = i1 -1 + lon
lambda = 2.*pie*(lon-1)/nlon(lat)
!
!jr Use these settings in conjunction with theta_c to start the blob at
!jr Greenwich
!
u_save(lon,k,lat) = u0*(cos(theta)*cos(alfa) + &
sin(theta)*cos(lambda-0.5*pie)*sin(alfa))
v_save(lon,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.
!
u_save(lon,k,lat) = u0*(cos(theta)*cos(alfa) + &
sin(theta)*cos(lambda)*sin(alfa))
v_save(lon,k,lat) = -u0*sin(lambda)*sin(alfa)
u3(i,k,j,1) = u_save(lon,k,lat)
v3(i,k,j,1) = v_save(lon,k,lat)
dlam = lambda - lambda_c
small_r = rearth*acos(sin(theta_c)*sin(theta) + &
cos(theta_c)*cos(theta)*cos(dlam))
q3(i,k,1,j,1) = 0.
if (small_r .lt. big_r) then
q3(i,k,1,j,1) = 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(i,k,j,1) = 300.
t3(i,k,j,1) = q3(i,k,1,j,1)
end do
end do
!
!jr Set surface pressure to 1.e5 everywhere
!
do lon=1,nlon(lat)
ps(lon,lat,1) = ps0
ps_save(lon,lat) = ps(lon,lat,1)
end do
end do
return
end subroutine hadvtest_init
!
!*********************************************************************C
#endif
end module inidat