!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "tmm.inc"
!
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: METEO
!
! !DESCRIPTION: Routines to initialize/finalize meteo grids and data, allocate
! datasets, and fill them. Include wrappers around read/write
! meteo files.
! Perform some meteo dependend calculations (omega, gph,
! mass <=> pressure, ...)
!
!
! !REVISION HISTORY:
!
! 09 Jun 2010 - P. Le Sager
! - Fix in METEO_SETUP_MASS when reading restart files.
! - Added some (protex) doc.
! - Merge updates from EC-Earth project.
!
! 10 Aug 2010 - Arjo Segers
! - Reset previous fix since it gives differences after a restart.
! - Use 'pw_dat' instead of 'mfw_dat' for massflux balancing;
! otherwise 'mfw_dat' is changed by matching its values in a zoom
! region with the parent, and this would give tiny differences
! between a long run and two smaller runs with a restart in between.
! - Reformatted protex comments.
!
! 10 Nov 2011 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
! 24 Oct 2013 - Ph. Le Sager - 6 new routines for parallel I/O
!
! !REMARKS:
!
! (1) Several surface pressure fields are used:
!
! sp1,sp2 : Surface pressures at begin and end of dynamic time step.
! Their values are interpolated between surface pressures
! read from the meteorological archive (in real(4) !)
! or received from the meteorological model.
! Fields from the meteorological archive are stored into
! the 'sp2' structure, and copied from there into 'sp1'.
!
! sp : Actual surface pressure due to advection.
! In theory this field is equal to 'sp1' at the begin of a timestep,
! but due to numerical inacuracies ( real(4) vs real(8) )
! tiny differeces occur.
! Therefore, this field is stored and restored in case of restart.
!
! spm Surface pressure for the mid of the time interval,
! thus in between 'sp1' and 'sp2' .
!
! (2) FIXME ZOOM : already coded, just need to check if it works as expected
!
! !INTERFACE:
!
MODULE METEO
!
! !USES:
!
use GO, only : gol, goErr, goPr, goLabel
use GO, only : TDate
use partools, only : isRoot
use grid, only : TllGridInfo, TLevelInfo
use tmm, only : TtmMeteo
!
use dims, only : nregions, nregions_all, okdebug
use tm5_distgrid, only : dgrid, Get_DistGrid, GATHER, SCATTER, UPDATE_HALO
use tm5_distgrid, only : SCATTER_J_BAND, SCATTER_I_BAND
USE METEODATA
IMPLICIT NONE
PRIVATE
!
! !PUBLIC MEMBER FUNCTIONS:
!
public :: Meteo_Init_Grids, Meteo_Done_Grids
public :: Meteo_Init, Meteo_Done, Meteo_Alloc
public :: Meteo_Setup_Mass, Meteo_Setup_Other
public :: Set, Check
public :: TimeInterpolation
!
! !PRIVATE TYPES:
!
type TMeteoField ! storage for a single meteo field:
logical :: used
character(len=16) :: name
character(len=16) :: unit
integer :: is(2), js(2), ls(2) ! shapes
real, pointer :: data(:,:,:)
end type TMeteoField
!
! !INTERFACE:
!
#ifdef with_parallel_io_meteo
interface Setup
module procedure Setup_2d_parallel_io
module procedure Setup_3d_parallel_io
end interface
#else
interface Setup
module procedure Setup_2d
module procedure Setup_3d
end interface
#endif
! Following are not striclty needed, since called here only once each
interface Setup_MFUV
module procedure Setup_MFUV_parallel_io
module procedure Setup_MFUV_serial_io
end interface
interface Setup_MFW
module procedure Setup_MFW_parallel_io
module procedure Setup_MFW_serial_io
end interface
interface Setup_CONVEC
module procedure Setup_CONVEC_parallel_io
module procedure Setup_CONVEC_serial_io
end interface
interface Setup_CLOUDCOVERS
module procedure Setup_CLOUDCOVERS_parallel_io
module procedure Setup_CLOUDCOVERS_serial_io
end interface
interface Setup_DIFFUS
module procedure Setup_DIFFUS_parallel_io
module procedure Setup_DIFFUS_serial_io
end interface
!
! !PRIVATE DATA MEMBERS:
!
character(len=*), parameter :: mname = 'Meteo'
type(TtmMeteo), save :: tmmd ! interface to TM meteo data
real :: sp_region0(1,1) ! single cell global surface pressure (region 0)
#ifdef with_tmm_tm5
logical, save :: use_tiedtke
#endif
!
!EOP
!------------------------------------------------------------------------
CONTAINS
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_INIT_GRIDS
!
! !DESCRIPTION: initialize grids and levels for each region. Grids on the
! local domain are simply copied from the DistGrid object.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_INIT_GRIDS( status )
!
! !USES:
!
use Grid, only : Init
use dims, only : region_name
use dims, only : xbeg, xend, dx, xref, im
use dims, only : ybeg, yend, dy, yref, jm
use dims, only : echlevs, lme, a_ec, b_ec
use geometry, only : geomtryv
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 19 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Meteo_Init_Grids'
integer :: n
real :: dlon, dlat
! --- begin ----------------------------
if (okdebug) call goLabel(rname)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! setup horizontal grids for the 0th one-cell grid
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! grid spacing:
dlon = real(xend(0)-xbeg(0))/im(0)
dlat = real(yend(0)-ybeg(0))/jm(0)
! define grid:
call Init( lli(0), xbeg(0)+dlon/2.0, dlon, im(0), &
ybeg(0)+dlat/2.0, dlat, jm(0), status, &
name=trim(region_name(0)) )
IF_NOTOK_RETURN(status=1)
! zonal grids
dlat = real(yend(0)-ybeg(0))/jm(0)
call Init( lli_z(0), 0.0, 360.0, 1, &
ybeg(0)+dlat/2.0, dlat, jm(0), status, &
name=trim(region_name(0))//'_z' )
IF_NOTOK_RETURN(status=1)
global_lli(0) = lli(0)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! local horizontal grid : get info from Distributed Grid
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
do n=1, nregions_all
call Get_DistGrid( dgrid(n), lli=lli(n), lli_z=lli_z(n), global_lli=global_lli(n) )
! correct name (it defines file to read data)
lli(n)%name = trim(region_name(n))
lli_z(n)%name = trim(region_name(n))//'_z'
global_lli(n)%name = trim(region_name(n))
end do
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! level definition
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! setup parent level definition:
call Init( levi_ec, lme, a_ec, b_ec, status ) ! ecmwf levels
IF_NOTOK_RETURN(status=1)
! setup level definition:
call Init( levi, levi_ec, echlevs, status ) ! tm half level selection
IF_NOTOK_RETURN(status=1)
! determine "old" at, bt of dims module
call geomtryv( )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_INIT_GRIDS
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_DONE_GRIDS
!
! !DESCRIPTION: finalize all grids and levels used for met fields.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_DONE_GRIDS( status )
!
! !USES:
!
use Grid, only : Done
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 19 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Meteo_Done_Grids'
integer :: n
! --- begin --------------------------------
if (okdebug) call goLabel(rname)
! horizontal (local) and zonal grids
do n = 0, nregions_all
call Done( lli(n), status )
IF_NOTOK_RETURN(status=1)
call Done( lli_z(n), status )
IF_NOTOK_RETURN(status=1)
end do
! horizontal (global) grids
do n = 1, nregions_all
call Done( global_lli(n), status )
IF_NOTOK_RETURN(status=1)
end do
! done parent level definition:
call Done( levi_ec, status )
IF_NOTOK_RETURN(status=1)
! level definition:
call Done( levi, status )
IF_NOTOK_RETURN(status=1)
! done
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_DONE_GRIDS
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_INIT
!
! !DESCRIPTION: Init met fields, i.e. nullify pointers, store shape, and set
! if needed (ie used) according to meteo.rc.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_INIT( status )
!
! !USES:
!
use GO, only : TrcFile, Init, Done, ReadRc
use Binas, only : p_global
use TMM, only : Init
use dims, only : im, jm, lm, lmax_conv
use meteodata, only : Init
use global_data, only : rcfile
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Meteo_Init'
! --- local -----------------------------
integer :: region, n
integer :: imr, jmr, lmr
integer :: halo
type(TrcFile) :: rcF
integer :: iveg
character(len=4) :: sveg
integer :: i01, i02, j01, j02
! --- begin ----------------------------
if (okdebug) call goLabel(rname)
! open rcfile:
call Init( rcF, rcfile, status )
IF_NOTOK_RETURN(status=1)
#ifdef with_tmm_tm5
! are convection fluxes derived (Tiedtke, sub files) or from IFS (convec files)?
call ReadRc( rcF, 'tiedtke', use_tiedtke, status )
IF_NOTOK_RETURN(status=1)
#endif
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! meteo database
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! setup interface to TM meteo:
call Init( tmmd, rcF, status )
IF_NOTOK_RETURN(status=1)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! define meteo data
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! global mean surface pressure
sp_region0 = p_global
! setup meteo fields: not in use, not allocated:
do region = 1, nregions_all
call Get_DistGrid( dgrid(region), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 )
lmr = lm(region)
!
! ** surface pressure *************************************
!
! two extra horizontal cells
halo = 2
! end of interval; also reads for sp1 and spm :
call Init_MeteoData( sp2_dat(region), 'sp', 'Pa', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','ml','sp'/), region, status )
IF_NOTOK_RETURN(status=1)
! check time interpolation:
if ( sp2_dat(region)%tinterp(1:6) /= 'interp' ) then
write (gol,'("surface pressure should be interpolated:")'); call goErr
write (gol,'(" requested tinterp : ",a)') trim(sp2_dat(region)%tinterp); call goErr
call goErr; status=1; return
end if
! start of interval (copied from sp2_dat):
call Init( sp1_dat(region), 'sp', 'Pa', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! current pressure:
call Init( sp_dat(region), 'sp', 'Pa', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! surface pressure at mid of dynamic time interval:
call Init( spm_dat(region), 'sp', 'Pa', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
!
! ** 3D pressure and mass **************************
!
! two extra horizontal cells (same as surface pressures)
halo = 2
! pressure at half levels (lm+1):
call Init( phlb_dat(region), 'phlb', 'Pa', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! air mass:
call Init( m_dat(region), 'm', 'kg', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
!
! ** massfluxes *************************************
!
! ~~ vertical
! no extra cells
halo = 0
! vertical flux (kg/s)
call Init_MeteoData( mfw_dat(region), 'mfw', 'kg/s', &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
rcF, (/'* ','ml ','mflux_w'/), region, status )
IF_NOTOK_RETURN(status=1)
! vertical flux (kg/s) : BALANCED
! NOTE: data is copied from mfw, thus use same tinterp
! for correct allocation of data arrays
call Init( pw_dat(region), 'pw', 'kg/s', mfw_dat(region)%tinterp, &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! tendency of surface pressure:
call Init_MeteoData( tsp_dat(region), 'tsp', 'Pa/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','ml ','mflux_w'/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~ horizontal
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:imr,1:jmr ,1:lmr) in pu_t(0:imr+1,0:jmr+1,0:lmr)
! ^ ^ ^ ^ too large !
! pv_tmpp --> pv(1:imr,1:jmr+1,1:lmr) in pv_t(0:imr+1,0:jmr+1,0:lmr)
! ^ ^ ^ ^ too large !
! The extra cells are implemented below as halo cells.
! one extra cell
halo = 1
!! east/west flux (kg/s)
!call Init( mfu_dat(region), 'mfu', 'kg/s', tinterp_curr, &
! (/0,imr/), (/1,jmr/), halo, (/1,lmr/), &
! sourcekey_curr, write_meteo, status, default=destkey_curr )
!IF_NOTOK_RETURN(status=1)
!! south/north flux (kg/s)
!call Init( mfv_dat(region), 'mfv', 'kg/s', tinterp_curr, &
! (/1,imr/), (/0,jmr/), halo, (/1,lmr/), &
! sourcekey_curr, write_meteo, status, default=destkey_curr )
!IF_NOTOK_RETURN(status=1)
! east/west flux (kg/s)
call Init_MeteoData( mfu_dat(region), 'mfu', 'kg/s', &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
rcF, (/'* ','ml ','mflux_uv'/), region, status )
IF_NOTOK_RETURN(status=1)
! south/north flux (kg/s)
call Init_MeteoData( mfv_dat(region), 'mfv', 'kg/s', &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
rcF, (/'* ','ml ','mflux_uv'/), region, status )
IF_NOTOK_RETURN(status=1)
!! east/west flux (kg/s) : BALANCED
!call Init( pu_dat(region), 'pu', 'kg/s', 'computed', &
! (/0,imr/), (/1,jmr/), halo, (/1,lmr/), &
! 'no-sourcekey', .false., 'no-destkey', status )
!IF_NOTOK_RETURN(status=1)
!
!! south/north flux (kg/s) : BALANCED
!call Init( pv_dat(region), 'pv', 'kg/s', 'computed', &
! (/1,imr/), (/0,jmr/), halo, (/1,lmr/), &
! 'no-sourcekey', .false., 'no-destkey', status )
!IF_NOTOK_RETURN(status=1)
halo = 1
! east/west flux (kg/s) : BALANCED
! NOTE: data is copied from mfu, thus use same tinterp
! for correct allocation of data arrays
call Init( pu_dat(region), 'pu', 'kg/s', mfu_dat(region)%tinterp, &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! south/north flux (kg/s) : BALANCED
! NOTE: data is copied from mfv, thus use same tinterp
! for correct allocation of data arrays
call Init( pv_dat(region), 'pv', 'kg/s', mfv_dat(region)%tinterp, &
(/i01,i02/), (/j01,j02/), halo, (/0,lmr/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
!
! ** temperature *************************************
!
! no extra cells
halo = 0
! temperature (K) (halo=0)
call Init_MeteoData( temper_dat(region), 'T', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','temper'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! ** humidity *************************************
!
! no extra cells
halo = 0
! humidity (kg/kg) (halo = 0)
call Init_MeteoData( humid_dat(region), 'Q', 'kg/kg', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','humid'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! ** computed *************************************
!
halo = 1 ! halo needed for station output in USER_OUTPUT_AEROCOM
! geopotential height(m) (lm+1, halo=0)
call Init( gph_dat(region), 'gph', 'm', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! no extra cells
halo = 0
! vertical velocity (Pa/s) (lm+1, halo=0)
call Init( omega_dat(region), 'omega', 'Pa/s', 'computed', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
!
! ** clouds *************************************
!
! no extra cells
halo = 0
! lwc liquid water content (kg/kg) (halo=0)
call Init_MeteoData( lwc_dat(region), 'CLWC', 'kg/kg', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','cloud'/), region, status )
IF_NOTOK_RETURN(status=1)
! iwc ice water content (kg/kg) (halo=0)
call Init_MeteoData( iwc_dat(region), 'CIWC', 'kg/kg', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','cloud'/), region, status )
IF_NOTOK_RETURN(status=1)
! cc cloud cover (fraction) (halo=0)
call Init_MeteoData( cc_dat(region), 'CC', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','cloud'/), region, status )
IF_NOTOK_RETURN(status=1)
! cco cloud cover overhead = above bottom of box (fraction) (halo=0)
call Init_MeteoData( cco_dat(region), 'CCO', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','cloud'/), region, status )
IF_NOTOK_RETURN(status=1)
! ccu cloud cover underfeet = below top of box (fraction) (halo=0)
call Init_MeteoData( ccu_dat(region), 'CCU', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','cloud'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! ** convection *************************************
!
! no extra cells
halo = 0
! entu entrainement updraft
call Init_MeteoData( entu_dat(region), 'eu', 'kg/m2/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), region, status )
IF_NOTOK_RETURN(status=1)
#ifdef with_tmm_tm5
if (.not. use_tiedtke .and.(entu_dat(region)%tinterp(1:4) /= 'aver')) then
write(gol,'("unexpected time interpolation ''",a,"''")') trim(entu_dat(region)%tinterp)
call goErr
TRACEBACK; status=1; return
endif
#endif
! entd entrainement downdraft (im,jm,lmax_conv)
call Init_MeteoData( entd_dat(region), 'ed', 'kg/m2/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), region, status )
IF_NOTOK_RETURN(status=1)
! detu detrainement updraft
call Init_MeteoData( detu_dat(region), 'du', 'kg/m2/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), region, status )
IF_NOTOK_RETURN(status=1)
! detd detrainement downdraft
call Init_MeteoData( detd_dat(region), 'dd', 'kg/m2/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! ** diffusion *************************************
!
! no extra cells
halo = 0
! turbulent diffusion coefficient for heat:
! dT/dt = d/dz ( Kz dT/dz )
! defined half levels,
! probably top of layer is archived (0,..,L-1), implicit zero for surface:
call Init_MeteoData( kzz_dat(region), 'K', 'm2/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), &
rcF, (/'* ','ml ','diffus'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! *** surface fields
!
! no extra cells
halo = 0
! orography (m*[g])
call Init_MeteoData( oro_dat(region), 'oro', 'm m/s2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.const','sfc.an ','oro '/), region, status )
IF_NOTOK_RETURN(status=1)
! land/sea mask (%)
call Init_MeteoData( lsmask_dat(region), 'lsm', '%', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.const','sfc.an ','lsm '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ instantaneous fields
! sea surface temperatue:
call Init_MeteoData( sst_dat(region), 'sst', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','sst '/), region, status )
IF_NOTOK_RETURN(status=1)
! 10m u wind (m/s)
call Init_MeteoData( u10m_dat(region), 'u10m', 'm/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','u10m '/), region, status )
IF_NOTOK_RETURN(status=1)
! 10m v wind (m/s)
call Init_MeteoData( v10m_dat(region), 'v10m', 'm/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','v10m '/), region, status )
IF_NOTOK_RETURN(status=1)
! skin reservoir content (m water) ; instant
call Init_MeteoData( src_dat(region), 'src', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','src '/), region, status )
IF_NOTOK_RETURN(status=1)
! 2 meter dewpoint temperature (K) ; instant
call Init_MeteoData( d2m_dat(region), 'd2m', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','d2m '/), region, status )
IF_NOTOK_RETURN(status=1)
! 2 meter temperature (K) ; instant
call Init_MeteoData( t2m_dat(region), 't2m', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','t2m '/), region, status )
IF_NOTOK_RETURN(status=1)
! skin temperature (K) ; instant
call Init_MeteoData( skt_dat(region), 'skt', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','skt '/), region, status )
IF_NOTOK_RETURN(status=1)
! boundary layer height (m) ; instant
call Init_MeteoData( blh_dat(region), 'blh', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','blh '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ average field (accumulated)
! surface sensible heat flux (W/m2) ; time aver
call Init_MeteoData( sshf_dat(region), 'sshf', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','sshf '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface latent heat flux (W/m2) ; time aver
call Init_MeteoData( slhf_dat(region), 'slhf', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','slhf '/), region, status )
IF_NOTOK_RETURN(status=1)
! east-west surface stress (N/m2); time aver
call Init_MeteoData( ewss_dat(region), 'ewss', 'N/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ewss '/), region, status )
IF_NOTOK_RETURN(status=1)
! north-south surface stress (N/m2) ; time aver
call Init_MeteoData( nsss_dat(region), 'nsss', 'N/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','nsss '/), region, status )
IF_NOTOK_RETURN(status=1)
halo = 1 ! halo needed for station output in USER_OUTPUT_AEROCOM
! convective precipitation (m/s) ; time aver
call Init_MeteoData( cp_dat(region), 'cp', 'm/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','cp '/), region, status )
IF_NOTOK_RETURN(status=1)
! large scale stratiform precipitation (m/s) ; time aver
call Init_MeteoData( lsp_dat(region), 'lsp', 'm/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','lsp '/), region, status )
IF_NOTOK_RETURN(status=1)
! no extra cells
halo = 0
! surface solar radiation ( W/m2 ) ; time aver
call Init_MeteoData( ssr_dat(region), 'ssr', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ssr '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface solar radiation downwards ( W/m2 ) ; time aver
call Init_MeteoData( ssrd_dat(region), 'ssrd', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ssrd '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface thermal radiation ( W/m2 ) ; time aver
call Init_MeteoData( str_dat(region), 'str', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','str '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface thermal radiation downwards ( W/m2 ) ; time aver
call Init_MeteoData( strd_dat(region), 'strd', 'W/m2', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','strd '/), region, status )
IF_NOTOK_RETURN(status=1)
! snow fall (m water eqv); time aver
call Init_MeteoData( sf_dat(region), 'sf', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','sf '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ time averages in grib files tfc+[12,15] etc
! 10m wind gust (m/s)
call Init_MeteoData( g10m_dat(region), 'g10m', 'm/s', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','g10m '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ in tmpp daily averages
! sea ice:
call Init_MeteoData( ci_dat(region), 'ci', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','ci '/), region, status )
IF_NOTOK_RETURN(status=1)
! snow depth (m water eqv); day aver ?
call Init_MeteoData( sd_dat(region), 'sd', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','sd '/), region, status )
IF_NOTOK_RETURN(status=1)
! volumetric soil water layer 1 ( m3 water / m3 soil) ; day aver ?
call Init_MeteoData( swvl1_dat(region), 'swvl1', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','swvl1 '/), region, status )
IF_NOTOK_RETURN(status=1)
! soil temperature layer 1 (K)
call Init_MeteoData( stl1_dat(region), 'stl1', 'K', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','stl1 '/), region, status )
IF_NOTOK_RETURN(status=1)
! vegetation type (%) ; day aver
do iveg = 1, nveg
write (sveg,'("tv",i2.2)') iveg
call Init_MeteoData( tv_dat(region,iveg), sveg, '%', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status )
IF_NOTOK_RETURN(status=1)
end do
! low vegetation cover (0-1) ; day aver
call Init_MeteoData( cvl_dat(region), 'cvl', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status )
IF_NOTOK_RETURN(status=1)
! high vegetation cover (0-1) ; day aver
call Init_MeteoData( cvh_dat(region), 'cvh', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status )
IF_NOTOK_RETURN(status=1)
! albedo ; daily average
call Init_MeteoData( albedo_dat(region), 'albedo', '1', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','albedo '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface roughness (ecmwf,ncep)
call Init_MeteoData( sr_ecm_dat(region), 'sr', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','sr '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ monthly data
! surface roughness (olsson) ; monthly
call Init_MeteoData( sr_ols_dat(region), 'srols', 'm', &
(/i01,i02/), (/j01,j02/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','srols '/), region, status )
IF_NOTOK_RETURN(status=1)
end do ! regions
! allocate work arrays for gather/scatter for I/O on grid #1
! Note : COULD BE that large on ROOT only, and simply (1,1,1) on other
! processors. Just have to be careful with the setup routines. Commented for zoom.
!PLS allocate( rwork1( im(1), jm(1), 0:lmr+1 ) )
!PLS allocate( rwork2( im(1), jm(1), 0:lmr+1 ) )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! close rcfile:
call Done( rcF, status )
IF_NOTOK_RETURN(status=1)
! ok
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_INIT
!EOC
!
! Read multiple keys in rcfile to setup meteodata structure.
! The following keys are read:
!
! meteo.tinterp. <-- time interpolation
! tmm.sourcekey.. <-- input file name description
! tmm.output.. <-- write meteo ?
! tmm.destkey.. <-- output file name description
!
! where is first '*' and then set to the grid name,
! and is set to each of the provided keys.
!
! Called for region=1..nregions_all
SUBROUTINE INIT_METEODATA( md, name, unit, is, js, halo, ls, &
rcF, rcs, region, status )
use GO, only : TRcFile, ReadRc
use Dims, only : nregions, nregions_max
use MeteoData, only : TMeteoData, Init, Set
! --- in/out -------------------------------------
type(TMeteoData), intent(out) :: md
character(len=*), intent(in) :: name, unit
integer, intent(in) :: is(2), js(2)
integer, intent(in) :: halo
integer, intent(in) :: ls(2)
type(TRcFile), intent(inout) :: rcF
character(len=*), intent(in) :: rcs(:)
integer, intent(in) :: region
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Init_MeteoData'
! --- local -------------------------------------
character(len=10) :: tinterp
character(len=256) :: sourcekey
logical :: write_meteo
character(len=256) :: destkey
logical :: used
! --- begin -------------------------------------
! time interpolation :
call ReadRc( rcF, 'meteo.tinterp', rcs, tinterp, status )
IF_NOTOK_RETURN(status=1)
! source filenames:
call ReadRc( rcF, 'tmm.sourcekey.*', rcs, sourcekey, status, default='no-sourcekey' )
IF_ERROR_RETURN(status=1)
call ReadRc( rcF, 'tmm.sourcekey.'//trim(lli(region)%name), rcs, sourcekey, status, default=sourcekey )
IF_ERROR_RETURN(status=1)
! write flag:
call ReadRc( rcF, 'tmm.output.*', rcs, write_meteo, status, default=.false. )
IF_ERROR_RETURN(status=1)
call ReadRc( rcF, 'tmm.output.'//trim(lli(region)%name), rcs, write_meteo, status, default=write_meteo )
IF_ERROR_RETURN(status=1)
! destination filenames:
if ( write_meteo ) then
call ReadRc( rcF, 'tmm.destkey.*', rcs, destkey, status, default='no-destkey' )
IF_ERROR_RETURN(status=1)
call ReadRc( rcF, 'tmm.destkey.'//trim(lli(region)%name), rcs, destkey, status, default=destkey )
IF_ERROR_RETURN(status=1)
else
destkey = 'no-destkey'
end if
! define meteo data,
! but should be marked as 'used' to be allocated and filled:
call Init( md, name, unit, tinterp, is, js, halo, ls, &
sourcekey, write_meteo, destkey, status )
IF_NOTOK_RETURN(status=1)
! read this type of meteo ?
! only regions 1..nregions or the extra fiels above nregions_max
! could be in use:
! [all regions, but do "if test", because nregions may be different from nregions_max]
if ( (region <= nregions) .or. (region > nregions_max) ) then
call ReadRc( rcF, 'meteo.read.*', rcs, used, status, default=.false. )
IF_ERROR_RETURN(status=1)
call ReadRc( rcF, 'meteo.read.'//trim(lli(region)%name), rcs, used, status, default=used )
IF_ERROR_RETURN(status=1)
else
used = .false.
end if
! in use ?
call Set( md, status, used=used )
IF_NOTOK_RETURN(status=1)
! ok
status = 0
END SUBROUTINE INIT_METEODATA
! ***
SUBROUTINE METEO_DONE( status )
use TMM, only : Done
use Dims, only : nregions_all
use meteodata, only : Done
! --- in/out -------------------------------
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Done'
! --- local -----------------------------
integer :: n
integer :: iveg
! --- begin --------------------------------
if (okdebug) call goLabel(rname)
! interface to TM meteo:
call Done( tmmd, status )
IF_NOTOK_RETURN(status=1)
!
! done meteo data
!
! destroy meteo fields:
do n = 1, nregions_all
! ***
call Done( sp1_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sp2_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sp_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( spm_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( phlb_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( m_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( mfu_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( mfv_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( mfw_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( tsp_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( pu_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( pv_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( pw_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( temper_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( humid_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( gph_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( omega_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( lwc_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( iwc_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( cc_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( cco_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( ccu_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( entu_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( entd_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( detu_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( detd_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( kzz_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
call Done( oro_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( lsmask_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( albedo_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sr_ecm_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sr_ols_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( ci_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sst_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( u10m_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( v10m_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( g10m_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( src_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( d2m_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( t2m_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( blh_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sshf_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( slhf_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( ewss_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( nsss_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( cp_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( lsp_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( ssr_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( ssrd_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( str_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( strd_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( skt_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sd_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( sf_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( swvl1_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( stl1_dat(n), status )
IF_NOTOK_RETURN(status=1)
do iveg = 1, nveg
call Done( tv_dat(n,iveg), status )
IF_NOTOK_RETURN(status=1)
end do
call Done( cvl_dat(n), status )
IF_NOTOK_RETURN(status=1)
call Done( cvh_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
end do ! regions
! work arrays
!PLS deallocate(rwork1, rwork2)
! ok
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_DONE
! ***
SUBROUTINE METEO_ALLOC( status )
use dims, only : nregions_all
use meteodata, only : Alloc
! --- in/out -------------------------------
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Alloc'
! --- local -----------------------------
integer :: region
integer :: iveg
! --- begin --------------------------------
if (okdebug) call goLabel(rname)
! allocate meteo fields if in use:
do region = 1, nregions_all
! ***
call Alloc( sp1_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sp2_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sp_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( spm_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( phlb_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( m_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( mfu_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( mfv_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( mfw_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( tsp_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( pu_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( pv_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( pw_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( temper_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( humid_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( gph_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( omega_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( lwc_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( iwc_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( cc_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( cco_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( ccu_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( entu_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( entd_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( detu_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( detd_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( kzz_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
call Alloc( oro_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( lsmask_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( albedo_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sr_ecm_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sr_ols_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( ci_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sst_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( u10m_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( v10m_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( src_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( d2m_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( t2m_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( skt_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( blh_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sshf_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( slhf_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( ewss_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( nsss_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( cp_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( lsp_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( ssr_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( ssrd_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( str_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( strd_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sd_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( sf_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( g10m_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( swvl1_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( stl1_dat(region), status )
IF_NOTOK_RETURN(status=1)
do iveg = 1, nveg
call Alloc( tv_dat(region,iveg), status )
IF_NOTOK_RETURN(status=1)
end do
call Alloc( cvl_dat(region), status )
IF_NOTOK_RETURN(status=1)
call Alloc( cvh_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
end do ! regions
! ok
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_ALLOC
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_SETUP_MASS
!
! !DESCRIPTION: Set up Mass FLuxes and Surface Pressures
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_SETUP_MASS( tr1, tr2, status, isfirst, check_pressure )
!
! !USES:
!
use go, only : TDate, rTotal, operator(-), wrtgol
use go, only : IncrDate, operator(+), Get
use grid, only : Match, TllGridInfo, assignment(=), Done
use Grid, only : FillMassChange, BalanceMassFluxes, CheckMassBalance
use dims, only : nregions, im, jm, lm, parent
use dims, only : xcyc
use meteodata, only : SetData ! to copy %data and %tr from one MD to another
#ifdef with_prism
use meteodata, only : TimeInterpolation
#endif
use restart, only : Restart_Read
!
! !INPUT PARAMETERS:
!
type(TDate), intent(in) :: tr1, tr2
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
logical, intent(in), optional :: check_pressure
logical, intent(in), optional :: isfirst
!
! !REVISION HISTORY:
!
! 12 Mar 2010 - P. Le Sager - Fix when reading restart files. Added
! protex doc. Added comments.
! 9 Jun 2010 - P. Le Sager - Merged with updates for EC-Earth project.
!
! 10 Aug 2010, Arjo Segers
! Reset previous fix since it makes a restart different from a long run.
! Use 'pw_dat' instead of 'mfw_dat' since otherwise the later changed
! while matching a zoom region with its parent, and this would give
! tiny differences during a restart of a zoomed run.
!
! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
! !REMARKS:
!
! push of Surf Press is done with sp2 (the only one on which we call
! setup -ie the only one for which %data1 and %data2 matter). Only %data
! of SP and SP1 are updated and used, and not their %data1 and %data2.
!
!------------------------------------------------------------------------------
!EOP
character(len=*), parameter :: rname = mname//'/Meteo_Setup_Mass'
logical :: do_check_pressure, WestBorder, NorthBorder
logical :: do_isfirst
integer :: n, p
integer :: idater(6)
real, allocatable :: dm_dt(:,:,:)
real :: dt_sec
integer :: l, i0, i1, j0, j1, is, js, ie, je
real :: tol_rms, tol_diff
type(TllGridInfo) :: L_lli
real, pointer :: field(:,:), field_parent(:,:) ! work arrays
real, allocatable :: islice(:,:), jslice(:,:), bigIslice(:,:), bigJslice(:,:)
real, allocatable :: full_pu(:,:,:), full_pv(:,:,:), full_pw(:,:,:), full_dm_dt(:,:,:)
#ifdef with_prism
integer :: hour1, dhour
#endif
! only for zoom regions (and only if with advection), when matching wih parents:
real, allocatable :: field3D_cur(:,:,:), field3D_par(:,:,:),wrkarr1(:,:,:),wrkarr2(:,:,:)
! --- begin --------------------------------
if (okdebug) call goLabel(rname)
! check pressure ?
if ( present(check_pressure) ) then
do_check_pressure = check_pressure
else
do_check_pressure = .false.
end if
! initial call ?
if ( present(isfirst) ) then
do_isfirst = isfirst
else
do_isfirst = .false.
end if
!
! ** HORIZONTAL MASS FLUXES *************************************
!
do n = 1, nregions_all
L_lli = global_lli(n)
! update horizontal u flux (unbalanced) -
#ifdef with_parallel_io_meteo
call Setup_MFUV( n, mfu_dat(n), mfv_dat(n), (/tr1,tr2/), levi, status)
IF_NOTOK_RETURN(status=1)
#else
call Setup_MFUV( n, mfu_dat(n), mfv_dat(n), (/tr1,tr2/), L_lli, levi, status)
IF_NOTOK_RETURN(status=1)
#endif
end do
!
! ** VERTICAL MASS FLUX *************************************
!
do n = 1, nregions_all
L_lli = global_lli(n)
! update vertical flux;
! tendency of surface pressure is by-product of vertical flux from spectral fields
! or filled with zero's
#ifdef with_parallel_io_meteo
call Setup_MFW( n, mfw_dat(n), tsp_dat(n), (/tr1,tr2/), 'n', levi, 'w', status)
IF_NOTOK_RETURN(status=1)
#else
call Setup_MFW( n, mfw_dat(n), tsp_dat(n), (/tr1,tr2/), L_lli, 'n', levi, 'w', status)
IF_NOTOK_RETURN(status=1)
#endif
end do
!
! ** SURFACE PRESSURES : SP1, SP *****************************
!
REG: do n = 1, nregions_all
! skip ?
if ( .not. sp1_dat(n)%used ) cycle
L_lli = global_lli(n)
! Advance 'next' surface pressure (a/k/a sp2%data) to start of
! new interval tr1. If start of a new meteo interval, then data
! is automatically read from file, or recieved from coupler
! with OASIS/prism
call Setup( n, sp2_dat(n), (/tr1,tr1/), L_lli, 'n', status )
IF_NOTOK_RETURN(status=1)
! copy SP2 into SP1 (%data and %tr)
call SetData( sp1_dat(n), sp2_dat(n), status )
IF_NOTOK_RETURN(status=1)
! GATHER sp1 array (dummy if not root)
!-----------------
! ...of parent region, if any:
if ( n /= 1 ) then
p = parent(n)
if (isRoot) then
allocate( field_parent(im(p), jm(p)) )
else
allocate( field_parent(1,1) )
end if
call GATHER( dgrid(p), sp1_dat(p)%data(:,:,1), field_parent, sp1_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
end if
! ...of current region:
if (isRoot) then
allocate( field(im(n), jm(n)) )
else
allocate( field(1,1) )
end if
call GATHER( dgrid(n), sp1_dat(n)%data(:,:,1), field, sp1_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! MATCH surface pressures to ensure mass balance
!-----------------
if (isRoot) then
! IF global field (i.e first region) : match global region with one-cell
! world value (average global surface pressure), ELSE match with parent
if ( n == 1 ) then
call Match( 'area-aver', 'n', lli(0), sp_region0, &
global_lli(n), field, status )
IF_NOTOK_RETURN(status=1)
else
call Match( 'area-aver', 'n', global_lli(p), field_parent, &
global_lli(n), field, status )
IF_NOTOK_RETURN(status=1)
endif
end if
! SCATTER sp1 array, and clean
!-----------------
call SCATTER( dgrid(n), sp1_dat(n)%data(:,:,1), field, sp1_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
deallocate(field)
if (n/=1) deallocate(field_parent)
! Set SP
!-----------------
! Initial call ? then set current surface pressure to just
! read/advanced sp1.
! otherwise, sp remains filled with the advected pressure.
if ( do_isfirst ) then
!write (gol,'(" copy SP1 to SP ...")'); call goPr
!pls ! PLS - (not working in the current OASIS/IFS setup. Kept for reference.
!pls ! If no_restart and is_first then sp2 is filled with
!pls ! t+dhour sp. Get SP from SP2 then:
!pls #ifdef with_prism
!pls
!pls select case ( sp2_dat(n)%tinterp )
!pls case ( 'interp6' ) ; dhour = 6
!pls case ( 'interp3' ) ; dhour = 3
!pls case ( 'interp2' ) ; dhour = 2
!pls case ( 'interp1' ) ; dhour = 1
!pls case default
!pls write (gol,'("unsupported time interpolation:")'); call goErr
!pls write (gol,'(" md%tinterp : ",a)') trim(sp2_dat(n)%tinterp); call goErr
!pls TRACEBACK; status=1; return
!pls end select
!pls dt_sec = dhour * 3600.0 ! sec
!pls
!pls sp_dat(n)%data(1:im(n),1:jm(n),1) = &
!pls sp2_dat(n)%data(1:im(n),1:jm(n),1) - &
!pls tsp_dat(n)%data(1:im(n),1:jm(n),1) * dt_sec
!pls
!pls sp_dat(n)%tr = tr1
!pls
!pls ! copy sp into sp1 :
!pls call SetData( sp1_dat(n), sp_dat(n), status )
!pls IF_NOTOK_RETURN(status=1)
!pls
!pls #else
! copy sp1 into sp :
call SetData( sp_dat(n), sp1_dat(n), status )
IF_NOTOK_RETURN(status=1)
!pls #endif
! fill pressure and mass from sp
call Pressure_to_Mass( n, status )
IF_NOTOK_RETURN(status=1)
! eventually replace by fields in restart file, since meteo
! from hdf meteo files is in real(4) while computed
! pressures and mass are probably in real(8) ;
! not for coupled run, since this receives pressures from
! ifs.
#ifndef oasis4
!#ifndef with_prism
call Restart_Read( status, surface_pressure=.true., pressure=.true., air_mass=.true. )
IF_NOTOK_RETURN(status=1)
!AJS>>> don't do this! sp1 contains data interpolated between
! fields received from the archive or the coupled model,
! while sp contains the actual pressure after advection.
!! copy sp into sp1 (PLS, 29-3-2010)
!call SetData( sp1_dat(n), sp_dat(n), status )
!IF_NOTOK_RETURN(status=1)
!<<<
#endif
end if ! end first
!! fill initial pressure and mass arrays,
!! eventually apply cyclic boundaries to mass
!call Meteo_SetupMass( n, status )
!IF_NOTOK_RETURN(status=1)
! check 'advected' pressure ?
if ( do_check_pressure) then
! compare 'advected' pressure still in sp with just read
! pressure sp1 : diff b/w sp%data and sp1%data
call Meteo_CheckPressure( n, status )
IF_NOTOK_RETURN(status=1)
end if
END DO REG ! regions
!
! ** SURFACE PRESSURES : SP2 *****************************
!
REG2: do n = 1, nregions_all
! skip ?
if ( .not. sp2_dat(n)%used ) cycle
!write (gol,'("sp2 ",a)') trim(lli(n)%name); call goPr
! grid and bounds
L_lli = global_lli(n)
i0 = sp2_dat(n)%is(1)
i1 = sp2_dat(n)%is(2)
j0 = sp2_dat(n)%js(1)
j1 = sp2_dat(n)%js(2)
#ifdef with_prism
! sp2 for prism coupler is computed from : sp(t2) = sp(t1) + tsp*(t2-t1)
if ( sp2_dat(n)%sourcekey(1:6) == 'prism:' ) then
select case ( sp2_dat(n)%tinterp )
case ( 'interp6' ) ; dhour = 6
case ( 'interp3' ) ; dhour = 3
case ( 'interp2' ) ; dhour = 2
case ( 'interp1' ) ; dhour = 1
case default
write (gol,'("unsupported time interpolation:")'); call goErr
write (gol,'(" md%tinterp : ",a)') trim(sp2_dat(n)%tinterp); call goErr
TRACEBACK; status=1; return
end select
! current interval [tr1,tr2] at begin of dhour interval ?
call Get( tr1, hour=hour1 )
if ( modulo(hour1,dhour) == 0 ) then
! reset sp2_dat%data1 and sp2_dat%data2:
!PLS ---- original code -----
!PLS
!PLS ! o data1 : surface pressure received for tr1
!PLS write (gol,'(" fill sp2%data1 with prism received field ...")'); call goPr
!PLS ! set filled flags to false to force re-reading if necessary;
!PLS ! prism received lnsp fields are stored in cache
!PLS ! thus re-reading is fast and error-free
!PLS sp2_dat(n)%filled1 = .false.
!PLS sp2_dat(n)%filled2 = .false.
!PLS ! now (re)read :
!PLS write (gol,'("PLS - (re)setup SP2_dat at ",i2)') tr1%hour; call goPr
!PLS call Setup( sp2_dat(n), (/tr1,tr1/), lli(n), 'n', status )
!PLS IF_NOTOK_RETURN(status=1)
!PLS
!PLS ! o data2 : data1 + tsp * dhour*3600.0 with dhour 3 or 1 hour
!PLS write (gol,'(" compute sp2%data2 from sp2%data1 and sp tendency ...")'); call goPr
!PLS dt_sec = dhour * 3600.0 ! sec
!PLS sp2_dat(n)%data2(1:im(n),1:jm(n),1) = &
!PLS sp2_dat(n)%data1(1:im(n),1:jm(n),1) + tsp_dat(n)%data(1:im(n),1:jm(n),1) * dt_sec
!PLS sp2_dat(n)%tr2 = tr1 + IncrDate(sec=nint(dt_sec))
!PLS
!PLS ! o data : interpolation between data1 and data2
!PLS call wrtgol( ' interpolate sp2%data to : ', tr2 ); call goPr
!PLS call TimeInterpolation( sp2_dat(n), (/tr2,tr2/), status )
!PLS IF_NOTOK_RETURN(status=1)
! Read into sp2%data1 : surface pressure received for
! tr1 (truly tr1+dhour)
!write (gol,'(" fill sp2%data1 with prism received field ...")'); call goPr
! set filled flags to false to force re-reading if necessary;
! prism received lnsp fields are stored in cache
! thus re-reading is fast and error-free
!PLS: read sp from prism for t=tr1 into SP2%DATA1. This
! is truly sp at t = tr1 + dhour, according to the
! coupler settings.
sp2_dat(n)%filled1 = .false.
sp2_dat(n)%filled2 = .false.
call Setup( n, sp2_dat(n), (/tr1,tr1/), L_lli, 'n', status )
IF_NOTOK_RETURN(status=1)
! move %data1 to %data2, and get %data1 from %data2:
! data1 = data2 - tsp * dhour*3600.0
!write (gol,'(" compute sp2%data1 from sp2%data2 and sp tendency ...")'); call goPr
dt_sec = dhour * 3600.0 ! sec
sp2_dat(n)%data2(i0:i1,j0:j1,1) = sp2_dat(n)%data1(i0:i1,j0:j1,1)
sp2_dat(n)%tr2 = tr1 + IncrDate(sec=nint(dt_sec))
sp2_dat(n)%data1(i0:i1,j0:j1,1) = &
sp2_dat(n)%data2(i0:i1,j0:j1,1) - tsp_dat(n)%data(i0:i1,j0:j1,1) * dt_sec
endif ! endif "it is beginning of coupling interval"
! Once SP2_DAT contains data1 and data2 valid for a dhour
! interval, %data is simply interpolated between %data1 and
! %data2:
!call wrtgol( ' interpolate sp2%data to : ', tr2 ); call goPr
call TimeInterpolation( sp2_dat(n), (/tr2,tr2/), status )
IF_NOTOK_RETURN(status=1)
!pls else
!pls
!pls ! sp2_dat contains data1 and data2 valid for a dhour interval;
!pls ! set %data to interpolation between %data1 and %data2:
!pls call wrtgol( ' interpolate sp2%data to : ', tr2 ); call goPr
!pls call TimeInterpolation( sp2_dat(n), (/tr2,tr2/), status )
!pls IF_NOTOK_RETURN(status=1)
!pls
!pls end if
else
! PLS: this one is never used apparently...
! AJS: it might be used in a partial coupling with only some fields
! exchanged and others read; this was often the case during the
! first coupling experiments, and might be useful for testing
! advance 'next' surface pressure to end of interval:
call Setup( n, sp2_dat(n), (/tr2,tr2/), L_lli, 'n', status )
IF_NOTOK_RETURN(status=1)
end if ! endif "it is prism sourcekey"
#else
! advance 'next' surface pressure to end of interval:
call Setup( n, sp2_dat(n), (/tr2,tr2/), L_lli, 'n', status )
IF_NOTOK_RETURN(status=1)
#endif /* WITH_PRISM */
! GATHER sp2 array (dummy if not root)
!-----------------
! ...of parent region, if any:
if ( n /= 1 ) then
p = parent(n)
if (isRoot) then
allocate( field_parent(im(p), jm(p)) )
else
allocate( field_parent(1,1) )
end if
call GATHER( dgrid(p), sp2_dat(p)%data(:,:,1), field_parent, sp2_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
end if
! ...of current region:
if (isRoot) then
allocate( field(im(n), jm(n)) )
else
allocate( field(1,1) )
end if
call GATHER( dgrid(n), sp2_dat(n)%data(:,:,1), field, sp2_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! MATCH surface pressures to ensure mass balance
!-----------------
if (isRoot) then
! IF global field (i.e first region) : match global region with one-cell
! world value (average global surface pressure), ELSE match with parent
if ( n == 1 ) then
call Match( 'area-aver', 'n', lli(0), sp_region0, &
global_lli(n), field, status )
IF_NOTOK_RETURN(status=1)
else
call Match( 'area-aver', 'n', global_lli(p), field_parent, &
global_lli(n), field, status )
IF_NOTOK_RETURN(status=1)
endif
end if
! SCATTER sp2 array, and clean
!-----------------
call SCATTER( dgrid(n), sp2_dat(n)%data(:,:,1), field, sp2_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
deallocate(field)
if (n/=1) deallocate(field_parent)
END DO REG2 ! regions
#ifndef without_advection
!
! ** MASS BALANCE *****************************
!
! NOTE: since only the surface pressure gradient is used,
! it is not necessary to use the data1 and data2 arrays
do n = 1, nregions_all
! skip ?
if ( .not. pu_dat(n)%used ) cycle
if ( .not. pv_dat(n)%used ) cycle
if ( .not. pw_dat(n)%used ) cycle
L_lli = global_lli(n)
i0 = sp2_dat(n)%is(1)
i1 = sp2_dat(n)%is(2)
j0 = sp2_dat(n)%js(1)
j1 = sp2_dat(n)%js(2)
! local indices and tile location (is, ie, js, je must be equal to i0, i1, j0, j1 BTW)
CALL GET_DISTGRID( dgrid(n), &
I_STRT=is, I_STOP=ie, &
J_STRT=js, J_STOP=je, &
hasWestBorder=WestBorder, hasNorthBorder=NorthBorder)
! length of time step between sp1 and sp2:
dt_sec = rTotal( sp2_dat(n)%tr(1) - sp1_dat(n)%tr(1), 'sec' )
! allocate temporary array:
allocate(dm_dt(i0:i1,j0:j1,lm(n)))
! mass change (kg) :
call FillMassChange( dm_dt, lli(n), levi, &
sp1_dat(n)%data(i0:i1,j0:j1,1), &
sp2_dat(n)%data(i0:i1,j0:j1,1), &
status )
IF_NOTOK_RETURN(status=1)
! mass tendency (kg/s) :
dm_dt = dm_dt / dt_sec ! kg/s
! >>> data1 >>>
! initial guess for balanced fluxes are unbalanced fluxes:
pu_dat(n)%data1 = mfu_dat(n)%data1
pu_dat(n)%filled1 = mfu_dat(n)%filled1
pu_dat(n)%tr1 = mfu_dat(n)%tr1
pv_dat(n)%data1 = mfv_dat(n)%data1
pv_dat(n)%filled1 = mfv_dat(n)%filled1
pv_dat(n)%tr1 = mfv_dat(n)%tr1
pw_dat(n)%data1 = mfw_dat(n)%data1
pw_dat(n)%filled1 = mfw_dat(n)%filled1
pw_dat(n)%tr1 = mfw_dat(n)%tr1
! MATCH WITH PARENT GRID IF NECESSARY
! -----------------------------------
! note strange indexing:
! pu_dat(n)%data1( 0:im(n), 1:jm(n) , 1:lm(n) )
! pv_dat(n)%data1( 1:im(n), 1:jm(n)+1, 1:lm(n) )
if ( n >1 ) then
p = parent(n)
! gather whole-region arrays
if (isRoot) then
allocate(field3D_cur(0:im(n),1:jm(n)+1,lm(n)))
allocate(field3D_par(0:im(p),1:jm(p)+1,lm(p)))
allocate(wrkarr1(im(n),jm(n),lm(n)))
allocate(wrkarr2(im(p),jm(p),lm(p)))
else
allocate( field3D_cur(1,1,1) )
allocate( field3D_par(1,1,1) )
allocate( wrkarr1(1,1,1) )
allocate( wrkarr2(1,1,1) )
end if
!for slice scattering
allocate(islice(j0:j1,lm(n)))
allocate(jslice(i0:i1,lm(n)))
if (isRoot) then
allocate(bigIslice(1:jm(n),lm(n)))
allocate(bigJslice(1:im(n),lm(n)))
else
allocate(bigIslice(1,1))
allocate(bigJslice(1,1))
end if
!------- U ----------------
call GATHER( dgrid(n), pu_dat(n)%data1, wrkarr1, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pu_dat(p)%data1, wrkarr2, pu_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
field3D_cur(1:im(n),1:jm(n),:) = wrkarr1
field3D_cur( 0,1:jm(n),:) = field3D_cur(im(n),1:jm(n),:) ! E-W periodicity
field3D_par(1:im(p),1:jm(p),:) = wrkarr2
field3D_par( 0,1:jm(p),:) = field3D_par(im(p),1:jm(p),:) ! E-W periodicity
do l = 1, lm(n)
call Match( 'sum', 'u', global_lli(p), field3D_par(0:im(p),1:jm(p),l), &
global_lli(n), field3D_cur(0:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
if(isRoot) wrkarr1 = field3D_cur(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pu_dat(n)%data1, wrkarr1, pu_dat(n)%halo, status )
! scatter extra column field3D_cur(0,:,:) - needed only for noncyclic
! zoom-region, for others update_halo takes care of it [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = field3D_cur(0,1:jm(n),:)
CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1)
if(WestBorder)pu_dat(n)%data1(0,j0:j1,:) = islice
!------- V ----------------
call GATHER( dgrid(n), pv_dat(n)%data1, wrkarr1, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pv_dat(p)%data1, wrkarr2, pv_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
field3D_cur(1:im(n),1:jm(n),:) = wrkarr1
field3D_cur(1:im(n),jm(n)+1,:) = field3D_cur(1:im(n),1,:) ! donnut periodicity
field3D_par(1:im(p),1:jm(p),:) = wrkarr2
field3D_par(1:im(p),jm(p)+1,:) = field3D_par(1:im(p),1,:) ! donnut periodicity
do l = 1, lm(n)
call Match( 'sum', 'v', global_lli(p), field3D_par(1:im(p),1:jm(p)+1,l), &
global_lli(n), field3D_cur(1:im(n),1:jm(n)+1,l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
if(isRoot) wrkarr1 = field3D_cur(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pv_dat(n)%data1, wrkarr1, pv_dat(n)%halo, status )
! scatter extra north row field3D_cur(:,j1+1,:)
if(isRoot)bigJslice = field3D_cur(:,jm(n)+1,:)
CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n))
if(NorthBorder)pv_dat(n)%data1(i0:i1,jm(n)+1,:)= jslice
!------- W ----------------
call GATHER( dgrid(n), pw_dat(n)%data1, wrkarr1, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pw_dat(p)%data1, wrkarr2, pw_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
do l = 0, lm(n)
call Match( 'sum', 'v', global_lli(p), wrkarr2(1:im(p),1:jm(p),l), &
global_lli(n), wrkarr1(1:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
call SCATTER( dgrid(n), pw_dat(n)%data1, wrkarr1, pw_dat(n)%halo, status )
!----- Done
deallocate(field3D_cur, field3d_par, wrkarr1, wrkarr2, bigJslice,&
bigIslice, jslice, islice)
end if
!#ifdef with_prism
! skip initial mass balance; relative large differences might exist
! between pressure imposed by mass fluxes and pressure according to
! surface pressure tendencies since the later is based on:
!
! sp(t-1)+tsp(t-1) _ *
! _ - o-------* sp(t), sp(t)+tsp(t)
! sp(t-1) o
!
! PLS : I do not understand that diagram... tsp is for an
! interval, and sp for a point in time. This may be
! wrong then. What we had at the first time step was:
!
! sp(t+1)+tsp(t:t+1) _ *
! _ - => sp(t) to sp(t)+tsp(t:t+1)
! sp(t+1) o
!
! AJS : This describes what the CTM received before the above
! described update. The 'tsp' was *not* for an interval but
! an instantaneous field describing the 'direction' of the surface
! pressure in time (you might call this 'tendency', but that is a
! dangerous word in GEMS IFS-CTM coupling context).
! Thus, at time 't-1' the only estimate of 'sp(t)' we could make was:
! sp(t-1)+tsp(t-1)
! At time 't' we received the actual 'sp(t)' and this was of course
! different from the initial guess.
!
! PLS : Just need to be sure that we have the correct sp to start
! with. Code above has been modified, so that we have:
!
! sp(t)+tsp(t:t+1) _ *
! _ - => sp(t) to sp(t)+tsp(t:t+1)
! sp(t) o
!
!#else
! CHECK INITIAL MASS BALANCE:
! -----------------------------------
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! tolerance for difference between sp from mass fluxes and sp from tendency:
tol_rms = 1.0e-4 ! max rms
tol_diff = 1.0e-3 ! max absolute difference
CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data1, pu_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data1, pv_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
call CheckMassBalance( lli(n), &
pu_dat(n)%data1(i0-1:i1, j0:j1 , 1:lm(n) ), &
pv_dat(n)%data1( i0:i1, j0:j1+1, 1:lm(n) ), &
sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
dt_sec, tol_rms, tol_diff, status )
if (status/=0) then
write (gol,'("initial mass imbalance too large for region ",i2)') n; call goErr
call goErr; status=1; return
end if
!#endif
! BALANCE HORIZONTAL FLUXES
! -----------------------------------
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! needs to be done globally... so gather data
if (isRoot) then
allocate(full_pu( 0:im(n), 1:jm(n), 0:lm(n)) ) ! must have same number of levels as mfu
allocate(full_pv( 1:im(n), 1:jm(n)+1, 0:lm(n)) )
allocate(full_pw( 1:im(n), 1:jm(n), 0:lm(n)) ) ! used also as temp arr in comm
allocate(full_dm_dt(im(n), jm(n), lm(n)) )
else
allocate( full_pu(1,1,1) )
allocate( full_pv(1,1,1) )
allocate( full_pw(1,1,1) )
allocate( full_dm_dt(1,1,1))
end if
!for slice scattering
allocate(islice(j0:j1,0:lm(n)))
allocate(jslice(i0:i1,0:lm(n)))
if (isRoot) then
allocate(bigIslice(1:jm(n),0:lm(n)))
allocate(bigJslice(1:im(n),0:lm(n)))
else
allocate(bigIslice(1,1))
allocate(bigJslice(1,1))
end if
call GATHER( dgrid(n), pu_dat(n)%data1, full_pw, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
full_pu(1:im(n),1:jm(n),:) = full_pw
full_pu(0,:,:) = full_pu(im(n),:,:) ! East-West periodicity
end if
call GATHER( dgrid(n), pv_dat(n)%data1, full_pw, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
full_pv(1:im(n),1:jm(n),:) = full_pw
full_pv(:,jm(n)+1,:) = full_pv(:,1,:) ! donut periodicity
end if
call GATHER( dgrid(n), dm_dt, full_dm_dt, 0, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(n), pw_dat(n)%data1, full_pw, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
! PRINT*, "BEFORE BMF:"
! print*, minval(full_pu(0:im(n),1:jm(n) ,1:lm(n))), maxval(full_pu(0:im(n),1:jm(n) ,1:lm(n)))
! print*, minval(full_pv(1:im(n),1:jm(n)+1,1:lm(n))), maxval(full_pv(1:im(n),1:jm(n)+1,1:lm(n)))
! print*, minval(full_pw), maxval(full_pw)
call BalanceMassFluxes( global_lli(n), &
full_pu(0:im(n),1:jm(n) ,1:lm(n)), &
full_pv(1:im(n),1:jm(n)+1,1:lm(n)), &
full_pw, full_dm_dt, global_lli(parent(n)), dt_sec, status )
IF_NOTOK_RETURN(status=1)
! PRINT*, "AFTER BMF:"
! print*, minval(full_pu(0:im(n),1:jm(n) ,1:lm(n))), maxval(full_pu(0:im(n),1:jm(n) ,1:lm(n)))
! print*, minval(full_pv(1:im(n),1:jm(n)+1,1:lm(n))), maxval(full_pv(1:im(n),1:jm(n)+1,1:lm(n)))
! print*, minval(full_pw), maxval(full_pw)
end if
call SCATTER( dgrid(n), pw_dat(n)%data1, full_pw, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if(isRoot) full_pw = full_pu(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pu_dat(n)%data1, full_pw, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! scatter extra column full_pu(0,:,:) - needed only for noncyclic zoom
! region, for others update_halo takes care of it [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = full_pu(0,1:jm(n),:)
CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1)
if(WestBorder)pu_dat(n)%data1(0,j0:j1,0:lm(n)) = islice
if(isRoot) full_pw = full_pv(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pv_dat(n)%data1, full_pw, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! Scatter PV(:,jm+1,:)
if(isroot) bigJslice=full_pv(1:im(n),jm(n)+1,:)
CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n))
if(NorthBorder)pv_dat(n)%data1(i0:i1,jm(n)+1,0:lm(n))=jslice
! CHECK FINAL MASS BALANCE:
! -----------------------------------
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! tolerance for difference between sp from mass fluxes and sp from tendency:
tol_rms = 1.0e-7 ! max rms
tol_diff = 1.0e-6 ! max absolute difference
CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data1, pu_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data1, pv_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
! print*, "sum before 4", sum(pv_dat(n)%data1(i0:i1, j0+1:j1+1, 1:lm(n)))
call CheckMassBalance( lli(n), &
pu_dat(n)%data1(i0-1:i1, j0:j1 , 1:lm(n) ), &
pv_dat(n)%data1( i0:i1, j0:j1+1, 1:lm(n) ), &
sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
dt_sec, tol_rms, tol_diff, status )
if (status/=0) then
write (gol,'("final mass imbalance too large for region ",i2)') n; call goErr
call goErr; status=1; return
end if
!done
deallocate(full_pw, full_pu, full_pv, full_dm_dt, bigJslice, bigIslice,&
jslice, islice)
! >>> data2 >>>
if ( any((/mfu_dat%filled2,mfv_dat%filled2,mfw_dat%filled2/)) ) then
! check ...
if ( .not. all((/mfu_dat(n)%filled2,mfv_dat(n)%filled2,mfw_dat(n)%filled2/)) ) then
write (gol,'("either none or all secondary data should be in use:")'); call goErr
write (gol,'(" mfu_dat%filled2 : ",l1)') mfu_dat(n)%filled2; call goErr
write (gol,'(" mfv_dat%filled2 : ",l1)') mfv_dat(n)%filled2; call goErr
write (gol,'(" mfw_dat%filled2 : ",l1)') mfw_dat(n)%filled2; call goErr
call goErr; status=1; return
end if
! initial guess for balanced fluxes are unbalanced fluxes:
pu_dat(n)%data2 = mfu_dat(n)%data2
pu_dat(n)%filled2 = .true.
pu_dat(n)%tr2 = mfu_dat(n)%tr2
pv_dat(n)%data2 = mfv_dat(n)%data2
pv_dat(n)%filled2 = .true.
pv_dat(n)%tr2 = mfv_dat(n)%tr2
pw_dat(n)%data2 = mfw_dat(n)%data2
pw_dat(n)%filled2 = .true.
pw_dat(n)%tr2 = mfw_dat(n)%tr2
! MATCH WITH PARENT GRID IF NECESSARY
! -----------------------------------
! note strange indexing:
! pu_dat(n)%data2( 0:im(n), 1:jm(n) , 1:lm(n) )
! pv_dat(n)%data2( 1:im(n), 1:jm(n)+1, 1:lm(n) )
if ( n >1 ) then
p = parent(n)
! gather whole-region arrays
if (isRoot) then
allocate(field3D_cur(0:im(n),1:jm(n)+1,lm(n)))
allocate(field3D_par(0:im(p),1:jm(p)+1,lm(p)))
allocate(wrkarr1(im(n),jm(n),lm(n)))
allocate(wrkarr2(im(p),jm(p),lm(p)))
else
allocate( field3D_cur(1,1,1) )
allocate( field3D_par(1,1,1) )
allocate( wrkarr1(1,1,1) )
allocate( wrkarr2(1,1,1) )
end if
!for slice scattering
allocate(islice(j0:j1,lm(n)))
allocate(jslice(i0:i1,lm(n)))
if (isRoot) then
allocate(bigIslice(1:jm(n),lm(n)))
allocate(bigJslice(1:im(n),lm(n)))
else
allocate(bigIslice(1,1))
allocate(bigJslice(1,1))
end if
!------- U ----------------
call GATHER( dgrid(n), pu_dat(n)%data2, wrkarr1, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pu_dat(p)%data2, wrkarr2, pu_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
field3D_cur(1:im(n),1:jm(n),:) = wrkarr1
field3D_cur( 0,1:jm(n),:) = field3D_cur(im(n),1:jm(n),:) ! E-W periodicity
field3D_par(1:im(p),1:jm(p),:) = wrkarr2
field3D_par( 0,1:jm(p),:) = field3D_par(im(p),1:jm(p),:) ! E-W periodicity
do l = 1, lm(n)
call Match( 'sum', 'u', global_lli(p), field3D_par(0:im(p),1:jm(p),l), &
global_lli(n), field3D_cur(0:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
if(isRoot) wrkarr1 = field3D_cur(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pu_dat(n)%data2, wrkarr1, pu_dat(n)%halo, status )
! scatter extra column field3D_cur(0,:,:) - needed only for non-cyclic
! zoom-region, for others update_halo takes care of it [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = field3D_cur(0,1:jm(n),:)
CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1)
if(WestBorder)pu_dat(n)%data2(0,j0:j1,:) = islice
!------- V ----------------
call GATHER( dgrid(n), pv_dat(n)%data2, wrkarr1, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pv_dat(p)%data2, wrkarr2, pv_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
field3D_cur(1:im(n),1:jm(n),:) = wrkarr1
field3D_cur(1:im(n),jm(n)+1,:) = field3D_cur(1:im(n),1,:) ! donnut periodicity
field3D_par(1:im(p),1:jm(p),:) = wrkarr2
field3D_par(1:im(p),jm(p)+1,:) = field3D_par(1:im(p),1,:) ! donnut periodicity
do l = 1, lm(n)
call Match( 'sum', 'v', global_lli(p), field3D_par(1:im(p),1:jm(p)+1,l), &
global_lli(n), field3D_cur(1:im(n),1:jm(n)+1,l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
if(isRoot) wrkarr1 = field3D_cur(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pv_dat(n)%data2, wrkarr1, pv_dat(n)%halo, status )
! scatter extra north row field3D_cur(:,j1+1,:)
if(isRoot)bigJslice = field3D_cur(:,jm(n)+1,:)
CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n))
if(NorthBorder)pv_dat(n)%data2(i0:i1,jm(n)+1,:)= jslice
!------- W ----------------
call GATHER( dgrid(n), pw_dat(n)%data2, wrkarr1, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(p), pw_dat(p)%data2, wrkarr2, pw_dat(p)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
do l = 0, lm(n)
call Match( 'sum', 'v', global_lli(p), wrkarr2(1:im(p),1:jm(p),l), &
global_lli(n), wrkarr1(1:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
end if
call SCATTER( dgrid(n), pw_dat(n)%data2, wrkarr1, pw_dat(n)%halo, status )
!----- Done
deallocate(field3D_cur, field3d_par, wrkarr1, wrkarr2, bigJslice,&
bigIslice, jslice, islice)
end if
! CHECK INITIAL MASS BALANCE:
! -----------------------------------
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data2, pu_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data2, pv_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
call CheckMassBalance( lli(n), &
pu_dat(n)%data2(i0-1:i1, j0:j1 , 1:lm(n) ), &
pv_dat(n)%data2( i0:i1, j0:j1+1, 1:lm(n) ), &
sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
dt_sec, 1.0e-4, 1.0e-3, status )
if (status/=0) then
write (gol,'("initial mass imbalance too large for region ",i2)') n; call goErr
call goErr; status=1; return
end if
! BALANCE HORIZONTAL FLUXES
! -----------------------------------
! balance horizontal fluxes:
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
if (isRoot) then
allocate(full_pu( 0:im(n), 1:jm(n), 0:lm(n)) ) ! must have same number of levels as mfu
allocate(full_pv( 1:im(n), 1:jm(n)+1, 0:lm(n)) )
allocate(full_pw( 1:im(n), 1:jm(n), 0:lm(n)) ) ! used also as temp arr in comm
allocate(full_dm_dt(im(n), jm(n), lm(n)) )
else
allocate( full_pu(1,1,1) )
allocate( full_pv(1,1,1) )
allocate( full_pw(1,1,1) )
allocate( full_dm_dt(1,1,1))
end if
!for slice scattering
allocate(islice(j0:j1,0:lm(n)))
allocate(jslice(i0:i1,0:lm(n)))
if (isRoot) then
allocate(bigIslice(1:jm(n),0:lm(n)))
allocate(bigJslice(1:im(n),0:lm(n)))
else
allocate(bigIslice(1,1))
allocate(bigJslice(1,1))
end if
call GATHER( dgrid(n), pu_dat(n)%data2, full_pw, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
full_pu(1:im(n),1:jm(n),:) = full_pw
full_pu(0,:,:) = full_pu(im(n),:,:) ! East-West periodicity
end if
call GATHER( dgrid(n), pv_dat(n)%data2, full_pw, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
full_pv(1:im(n),1:jm(n),:) = full_pw
full_pv(:,jm(n)+1,:) = full_pv(:,1,:) ! donut periodicity
end if
call GATHER( dgrid(n), dm_dt, full_dm_dt, 0, status )
IF_NOTOK_RETURN(status=1)
call GATHER( dgrid(n), pw_dat(n)%data2, full_pw, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if (isRoot) then
call BalanceMassFluxes( global_lli(n), &
full_pu(0:im(n),1:jm(n) ,1:lm(n)), &
full_pv(1:im(n),1:jm(n)+1,1:lm(n)), &
full_pw, full_dm_dt, global_lli(parent(n)), dt_sec, status )
IF_NOTOK_RETURN(status=1)
end if
call SCATTER( dgrid(n), pw_dat(n)%data2, full_pw, pw_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
if(isRoot) full_pw = full_pu(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pu_dat(n)%data2, full_pw, pu_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! scatter extra column full_pu(0,:,:) - needed only for noncyclic zoom
! regions, for others update_halo takes care of it [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = full_pu(0,1:jm(n),:)
CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1)
if(WestBorder) pu_dat(n)%data2(0,j0:j1,:) = islice
if(isRoot) full_pw = full_pv(1:im(n),1:jm(n),:)
call SCATTER( dgrid(n), pv_dat(n)%data2, full_pw, pv_dat(n)%halo, status )
IF_NOTOK_RETURN(status=1)
! Scatter PV(:,jm+1,:)
if(isroot) bigJslice=full_pv(1:im(n),jm(n)+1,:)
CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n))
if(NorthBorder)pv_dat(n)%data2(i0:i1,jm(n)+1,:)=jslice
! CHECK FINAL MASS BALANCE:
! -----------------------------------
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n))
CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data2, pu_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data2, pv_dat(n)%halo, status)
IF_NOTOK_RETURN(status=1)
call CheckMassBalance( lli(n), &
pu_dat(n)%data2(i0-1:i1, j0:j1 , 1:lm(n) ), &
pv_dat(n)%data2( i0:i1, j0:j1+1, 1:lm(n) ), &
sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), &
dt_sec, 1.0e-7, 1.0e-6, status )
if (status/=0) then
write (gol,'("final mass imbalance too large for region ",i2)') n; call goErr
call goErr; status=1; return
end if
deallocate(full_pw, full_pu, full_pv, full_dm_dt, bigJslice, bigIslice,&
jslice, islice)
end if ! filled2
! >>>
! clear
deallocate( dm_dt )
end do ! regions
#endif /* ADVECTION */
!------------
! Done
!------------
call done(l_lli, status)
IF_NOTOK_RETURN(status=1)
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_SETUP_MASS
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_SETUP_OTHER
!
! !DESCRIPTION:
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_SETUP_OTHER( tr1, tr2, status, isfirst )
!
! !USES:
!
use GO, only : TDate, NewDate, rTotal, wrtgol
use GO, only : operator(-), operator(+), operator(/)
use GO, only : InterpolFractions
use dims, only : nregions, im, jm, lm
use dims, only : lmax_conv
use dims, only : xcyc
use Dims, only : czeta
use global_data, only : region_dat
#ifndef without_convection
use global_data, only : conv_dat
#endif
use Phys, only : ConvCloudDim
!
! !INPUT PARAMETERS:
!
type(TDate), intent(in) :: tr1, tr2
logical, intent(in), optional :: isfirst
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Meteo_Setup_Other'
logical :: do_isfirst
integer :: n, p
integer :: i, j, l
integer :: lsave, i0, i1, j0, j1
real :: tote, totd, maxe
real, pointer :: dxyp(:)
type(TDate) :: tmid
real :: alfa1, alfa2
integer :: iveg
! --- begin --------------------------------
if (okdebug) call goLabel(rname)
! initial call ?
if ( present(isfirst) ) then
do_isfirst = isfirst
else
do_isfirst = .false.
end if
!
! ** orography *****************************
!
! read orographies (if necessary):
do n = 1, nregions_all
call setup( n, oro_dat(n), (/tr1,tr2/), global_lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
end do
!
! ** spm **************************************
!
! loop over regions:
do n = 1, nregions
! skip ?
if ( .not. spm_dat(n)%used ) cycle
!write (gol,'("spm ",a)') trim(lli(n)%name); call goPr
! mid time:
tmid = tr1 + ( tr2 - tr1 )/2
! deterimine weights to sp1 and sp2 :
call InterpolFractions( tmid, sp1_dat(n)%tr(1), sp2_dat(n)%tr(1), alfa1, alfa2, status )
IF_NOTOK_RETURN(status=1)
call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 )
! interpolate:
spm_dat(n)%data(i0:i1,j0:j1,1) = alfa1 * sp1_dat(n)%data(i0:i1,j0:j1,1) + &
alfa2 * sp2_dat(n)%data(i0:i1,j0:j1,1)
! store time:
spm_dat(n)%tr = (/tr1,tr2/)
end do ! regions
!
! ** omega **************************************
!
! loop over regions:
do n = 1, nregions_all
!if (omega_dat(n)%used) then; write (gol,'("omega ",a)') trim(lli(n)%name); call goPr; end if
! re-compute omega from vertical mass flux:
call Compute_Omega( omega_dat(n), lli(n), mfw_dat(n), status )
IF_NOTOK_RETURN(status=1)
end do ! regions
!
! ** temperature and humid **************************************
!
! loop over regions:
do n = 1, nregions_all
! ncep meteo requires conversion of virtual temperature using humidity ...
if ( (temper_dat(n)%sourcekey(1:4) == 'ncep') .or. (humid_dat(n)%sourcekey(1:4) == 'ncep') ) then
!write (gol,'("temper and humid ",a)') trim(lli(n)%name); call goPr
! read temperature and humidity (if necessary):
! #ifdef with_parallel_io_meteo [COMMENTED SINCE NOT TM5-NC SOURCEKEY]
! call setup_TQ( n, temper_dat(n), humid_dat(n), (/tr1,tr2/), levi, status)
! IF_NOTOK_RETURN(status=1)
! #else
call setup_TQ( n, temper_dat(n), humid_dat(n), (/tr1,tr2/), global_lli(n), levi, status)
IF_NOTOK_RETURN(status=1)
! #endif
else
!if (temper_dat(n)%used) then; write (gol,'("temper ",a)') trim(lli(n)%name); call goPr; end if
! read temperature (if necessary):
call setup( n, temper_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status)
IF_NOTOK_RETURN(status=1)
!if (humid_dat(n)%used) then; write (gol,'("humid ",a)') trim(lli(n)%name); call goPr; end if
! read humidity (if necessary):
call setup( n, humid_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status)
IF_NOTOK_RETURN(status=1)
end if
end do ! regions
!
! ** gph **************************************
!
! loop over regions:
do n = 1, nregions_all
!write (gol,'("gph for region ",a," (#",i1,") ", l2)') trim(lli(n)%name), n, gph_dat(n)%used; call goPr
! re-compute gph from pressure, temperature, and humidity:
call compute_gph( n, status )
IF_NOTOK_RETURN(status=1)
end do ! regions
!
! ** clouds **************************************
!
! loop over regions:
do n = 1, nregions
!if (any((/lwc_dat(n)%used,iwc_dat(n)%used,cc_dat(n)%used,cco_dat(n)%used,ccu_dat(n)%used/))) then
! write (gol,'("clouds ",a)') trim(lli(n)%name); call goPr
!end if
call setup( n, lwc_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status)
IF_NOTOK_RETURN(status=1)
call setup( n, iwc_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status)
IF_NOTOK_RETURN(status=1)
#ifdef with_parallel_io_meteo
call setup_CloudCovers( n, cc_dat(n), cco_dat(n), ccu_dat(n), (/tr1,tr2/), levi, status)
IF_NOTOK_RETURN(status=1)
#else
call setup_CloudCovers( n, cc_dat(n), cco_dat(n), ccu_dat(n), (/tr1,tr2/), global_lli(n), levi, status)
IF_NOTOK_RETURN(status=1)
#endif
end do
!
! ** convection **************************************
!
! loop over regions:
do n = 1, nregions
!if (entu_dat(n)%used) then; write (gol,'("convection ",a)') trim(lli(n)%name); call goPr; end if
! read (if necessary):
#ifdef with_parallel_io_meteo
call setup_Convec( n, entu_dat(n), entd_dat(n), detu_dat(n), detd_dat(n), &
omega_dat(n), gph_dat(n), (/tr1,tr2/), levi, status )
IF_NOTOK_RETURN(status=1)
#else
call setup_Convec( n, entu_dat(n), entd_dat(n), detu_dat(n), detd_dat(n), &
omega_dat(n), gph_dat(n), (/tr1,tr2/), global_lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
#endif
end do
#ifndef without_convection
! ~~ convective clouds
! loop over regions:
do n = 1, nregions
! skip ?
if ( .not. entu_dat(n)%used ) cycle
if ( .not. entd_dat(n)%used ) cycle
! update necessary ?
if ( any((/entu_dat(n)%changed,entd_dat(n)%changed/)) ) then
call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 )
! loop over grid cells
do j = j0, j1
do i = i0, i1
! compute convective cloud dimensions for this column:
call ConvCloudDim( 'u', size(detu_dat(n)%data,3), &
detu_dat(n)%data(i,j,:), entd_dat(n)%data(i,j,:),&
conv_dat(n)%cloud_base(i,j), &
conv_dat(n)%cloud_top (i,j), &
conv_dat(n)%cloud_lfs (i,j), &
status )
IF_NOTOK_RETURN(status=1)
end do ! i
end do ! j
end if ! changed
end do ! regions
#endif
! ~~ unit conversion
! loop over regions:
do n = 1, nregions
! skip ?
if ( .not. entu_dat(n)%used ) cycle
if ( .not. entd_dat(n)%used ) cycle
if ( .not. detu_dat(n)%used ) cycle
if ( .not. detd_dat(n)%used ) cycle
! update necessary ?
if ( any((/ entu_dat(n)%changed, entd_dat(n)%changed, &
detu_dat(n)%changed, detd_dat(n)%changed /)) ) then
call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 )
!cmk calculate the rates in kg/gridbox and scale with czeta
dxyp => region_dat(n)%dxyp
do j = j0, j1
do i = i0, i1
! kg/m2/s -> kg/gridbox/s * scale_factor
entu_dat(n)%data(i,j,:) = entu_dat(n)%data(i,j,:)*dxyp(j)*czeta
detu_dat(n)%data(i,j,:) = detu_dat(n)%data(i,j,:)*dxyp(j)*czeta
! ensure netto zero tracer transport by updraught in column
! (add difference between total entrement and detrement
! to level where entrement reaches maximum):
tote = sum( entu_dat(n)%data(i,j,:) )
totd = sum( detu_dat(n)%data(i,j,:) )
maxe = entu_dat(n)%data(i,j,1) ! changed: reported by PB feb 2003
lsave = 1
do l = 2, lmax_conv
if ( entu_dat(n)%data(i,j,l) > maxe ) then
maxe = entu_dat(n)%data(i,j,l)
lsave = l
end if
end do
entu_dat(n)%data(i,j,lsave) = entu_dat(n)%data(i,j,lsave) - tote + totd
! kg/m2/s -> kg/gridbox/s * scale_factor
entd_dat(n)%data(i,j,:) = entd_dat(n)%data(i,j,:)*dxyp(j)*czeta
detd_dat(n)%data(i,j,:) = detd_dat(n)%data(i,j,:)*dxyp(j)*czeta
! ensure netto zero tracer transport by downdraught in column
! (add difference between total entrement and detrement
! to level where entrement reaches maximum):
tote = sum( entd_dat(n)%data(i,j,:) ) ! total entrainement
totd = sum( detd_dat(n)%data(i,j,:) ) ! total detrainement
maxe = 0.0
lsave = lmax_conv
do l = 1, lmax_conv
if ( entd_dat(n)%data(i,j,l) > maxe ) then
maxe = entd_dat(n)%data(i,j,l)
lsave = l
end if
end do
entd_dat(n)%data(i,j,lsave) = entd_dat(n)%data(i,j,lsave) - tote + totd
end do
end do
end if ! changed ?
end do ! regions
!
! ** diffusion **************************************
!
! loop over regions:
do n = 1, nregions
#ifdef with_parallel_io_meteo
call Setup_Diffus( n, kzz_dat(n), (/tr1,tr2/), levi, status )
IF_NOTOK_RETURN(status=1)
#else
call Setup_Diffus( n, kzz_dat(n), (/tr1,tr2/), global_lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
#endif
end do ! regions
!
! ** surface fields *****************************
!
do n = 1, nregions_all
!write (gol,'("surface fields ",a)') trim(lli(n)%name); call goPr
! * lsmask
call setup( n, lsmask_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * albedo
call setup( n, albedo_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * sr_ecm
call setup( n, sr_ecm_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * sr_ols
call setup( n, sr_ols_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * sea ice
call setup( n, ci_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * sea surface temperature
call setup( n, sst_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * u10m
call setup( n, u10m_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * v10m
call setup( n, v10m_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * skin reservoir content
call setup( n, src_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * 2m dewpoint temperature
call setup( n, d2m_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * 2m temperature
call setup( n, t2m_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * slhf
call setup( n, slhf_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * sshf
call setup( n, sshf_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * surface stress
call setup( n, ewss_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
call setup( n, nsss_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * convective precipitation
call setup( n, cp_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * large scale stratiform precipitation
call setup( n, lsp_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * surface solar radiation
call setup( n, ssr_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
call setup( n, ssrd_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * surface thermal radiation
call setup( n, str_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
call setup( n, strd_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * skin temperature
call setup( n, skt_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * boundary layer height
call setup( n, blh_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * snow fall and depth
call setup( n, sf_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
call setup( n, sd_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * g10m
call setup( n, g10m_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * soil water level 1
call setup( n, swvl1_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * soil temperature level 1
call setup( n, stl1_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * vegetation types
do iveg = 1, nveg
select case ( iveg )
case ( 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 13, 16, 17, 18, 19 )
call setup( n, tv_dat(n,iveg), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
case ( 8, 12, 14, 15, 20 )
if ( tv_dat(n,iveg)%used ) tv_dat(n,iveg)%data = 0.0
case default
write (gol,'("do not know how to setup vegetation type ",i2)') iveg
call goErr; status=1; return
end select
end do
! * low vegetation cover
call setup( n, cvl_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
! * high vegetation cover
call setup( n, cvh_dat(n), (/tr1,tr2/), global_lli(n), 'n', status)
IF_NOTOK_RETURN(status=1)
end do ! regions
!
! ** done ********************************************
!
! ok
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_SETUP_OTHER
!EOC
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUPSETUP
!
! !DESCRIPTION: for one met data MD and one time range TR, returns the dates
! at begining and end of the met field interval that
! encompasses TR, and if the data for these dates (%data1 and
! %data2, resp.) must be read or copied.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUPSETUP( md, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
!
! !USES:
!
use GO, only : TDate, NewDate, IncrDate, AnyDate, IsAnyDate, Get, Set, wrtgol
use GO, only : rTotal, iTotal
use GO, only : operator(+), operator(-), operator(/)
use GO, only : operator(==), operator(/=), operator(<), operator(<=)
use meteodata, only : TMeteoData
use global_data, only : fcmode, tfcday0
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md
!
! !INPUT PARAMETERS:
!
type(TDate), intent(in) :: tr(2)
!
! !OUTPUT PARAMETERS:
!
logical, intent(out) :: data1_read, data1_copy
type(TDate), intent(out) :: data1_tref, data1_t1, data1_t2
logical, intent(out) :: data2_read, data2_copy
type(TDate), intent(out) :: data2_tref, data2_t1, data2_t2
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 29 Mar 2010 - P. Le Sager -
!
!EOP
!------------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/SetupSetup'
integer :: dth, baseh
integer :: year, month, day, hour, minu
type(TDate) :: tmid
type(TDate) :: tc(2)
integer :: dth_int
type(TDate) :: tprev, tnext
real :: dhr
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! default output:
data1_read = .false.
data1_copy = .false.
data2_read = .false.
data2_copy = .false.
!
! trap constant fields ...
!
! constant and already filled ? then leave
if ( (md%tinterp == 'const') .and. md%filled1 ) then
if (okdebug) call goLabel()
status = 0; return
end if
!
! fc stuff
!
! 3 hourly data only available up to 72h, then 6 hourly
if ( fcmode ) then
! number of hours from fcday 00:00 to end of requested interval:
dhr = rTotal( tr(2) - tfcday0, 'hour' )
! lower time resolution after a while ...
if ( tfcday0 < NewDate(year=2006,month=03,day=14) ) then
! after 12+72 hour ?
if ( dhr > 12.0 + 72.0 ) then
! convert time interpolation:
select case ( md%tinterp )
case ( 'aver3' )
write (gol,'("WARNING - convert time interpolation from `aver3` to `aver6`")'); call goPr
md%tinterp = 'aver6'
case ( 'interp3' )
write (gol,'("WARNING - convert time interpolation from `interp3` to `interp6`")'); call goPr
md%tinterp = 'interp6'
end select
end if ! > 72 hour
else
! after 12+96 hour ?
if ( dhr > 12.0 + 96.0 ) then
! convert time interpolation:
select case ( md%tinterp )
case ( 'aver3' )
write (gol,'("WARNING - convert time interpolation from `aver3` to `aver6`")'); call goPr
md%tinterp = 'aver6'
case ( 'interp3' )
write (gol,'("WARNING - convert time interpolation from `interp3` to `interp6`")'); call goPr
md%tinterp = 'interp6'
end select
end if ! > 96 hour
end if ! change in fc resolution
end if ! fcmode
!
! time stuff
!
! basic time resolution in hours
select case ( md%tinterp )
case ( 'const', 'month' )
! nothing to be set here ...
case ( 'aver24' )
! constant fields produced valid for [00,24]
dth = 24
baseh = 00
case ( 'aver24_3' )
! constant fields produced by tmpp valid for [21,21] = [09-12,09+12]
dth = 24
baseh = -3
case ( 'const3', 'interp3', 'aver3', 'cpl3' )
dth = 3
baseh = 0
case ( 'interp2', 'cpl2' )
dth = 2
baseh = 0
case ( 'const1', 'interp1', 'aver1', 'cpl1' )
dth = 1
baseh = 0
case ( 'const6', 'interp6', 'aver6', 'cpl6' )
dth = 6
baseh = 0
case ( 'interp6_3' )
dth = 6
baseh = 3
case default
write (gol,'("unsupported time interpolation : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end select
! set time parameters for field to be read:
select case ( md%tinterp )
!
! ** constant fields
!
case ( 'const' )
! read main field ?
data1_read = .not. md%filled1
! read or leave ?
if ( data1_read ) then
data1_tref = tr(1) ! <--- used for file names
data1_t1 = AnyDate()
data1_t2 = AnyDate()
else
! field valid around requested interval, thus leave:
if (okdebug) call goLabel()
status=0; return
end if
!
! ** constant fields, valid for complete month
!
case ( 'month' )
! extract time values for begin of current interval:
call Get( tr(1), year=year, month=month )
! interval for this month:
tc(1) = NewDate( year=year, month=month, day=01, hour=00 )
month = month + 1
if ( month > 12 ) then
month = 1
year = year + 1
end if
tc(2) = NewDate( year=year, month=month, day=01, hour=00 )
! check for strange values:
if ( (tr(1) < tc(1)) .or. (tc(2) < tr(2)) ) then
write (gol,'("determined invalid constant interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tc(1), ' - ', tc(2) ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
!write (gol,'(" WARNING - requested interval exceeds meteo interval; should be improved")')
end if
! read main field ?
if ( md%filled1 ) then
data1_read = md%tr1(1) /= tc(1)
else
data1_read = .true.
end if
! read or leave ?
if ( data1_read ) then
data1_tref = tr(1)
data1_t1 = tc(1)
data1_t2 = tc(2)
else
! field valid around requested interval, thus leave:
if (okdebug) call goLabel()
status=0; return
end if
!
! ** constant fields, valid for 24hr intervals [21:00,21:00]
! constant fields, valid for 6hr intervals [21:00,03:00] etc
! constant fields, valid for 3hr intervals [22:30,01:30] etc
!
case ( 'const6', 'const3' )
! extract time values for begin of current interval:
call Get( tr(1), year, month, day, hour, minu )
! round hour to 00/06/12/18 or 00/03/06/09/12/15/18/21 or 09
hour = dth * nint(real(hour+minu/60.0-baseh)/real(dth)) + baseh
! set mid of 3 or 6 hour interval:
tmid = NewDate( year, month, day, hour )
! interval with constant field
tc(1) = tmid - IncrDate(hour=dth)/2
tc(2) = tmid + IncrDate(hour=dth)/2
! check for strange values:
if ( (tr(1) < tc(1)) .or. (tc(2) < tr(2)) ) then
write (gol,'("determined invalid constant interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tc(1), ' - ', tc(2) ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end if
! read main field ?
if ( md%filled1 ) then
data1_read = md%tr1(1) /= tmid
else
data1_read = .true.
end if
! read or leave ?
if ( data1_read ) then
data1_tref = tmid
data1_t1 = tmid
data1_t2 = tmid
else
! field valid around requested interval, thus leave:
if (okdebug) call goLabel()
status=0; return
end if
!
! ** couple fields, valid for 3hr intervals [00:00,03:00] etc
! input filed valid for BEGIN of interval !
!
case ( 'cpl6', 'cpl3', 'cpl2', 'cpl1' )
! extract time values for begin of current interval:
call Get( tr(1), year, month, day, hour, minu )
! round hour to previous baseh + 00/03/06/09/12/15/18/21
hour = dth * floor(real(hour-baseh)/real(dth)) + baseh
! interval with constant field
tc(1) = NewDate( year, month, day, hour )
tc(2) = tc(1) + IncrDate(hour=dth)
! check for strange values:
if ( (tr(1) < tc(1)) .or. (tc(2) < tr(1)) ) then
write (gol,'("determined invalid first interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tc(1), ' - ', tc(2) ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end if
! read primary field ?
if ( md%filled1 ) then
! read new field if times are different:
data1_read = (md%tr1(1) /= tc(1)) .or. (md%tr1(2) /= tc(1))
else
! not filled yet, thus must read:
data1_read = .true.
end if
! read or leave ?
if ( data1_read ) then
data1_tref = tc(1) ! begin of time interval
data1_t1 = tc(1)
data1_t2 = tc(1)
end if
!
! ** average fields, valid for 3hr intervals [00:00,03:00] etc
! average fields, valid for 3hr intervals [00:00,06:00] etc
!
case ( 'aver1', 'aver3', 'aver6', 'aver24', 'aver24_3' )
! extract time values for begin of current interval:
call Get( tr(1), year, month, day, hour, minu )
! round hour to previous baseh + 00/03/06/09/12/15/18/21
hour = dth * floor(real(hour-baseh)/real(dth)) + baseh
! interval with constant field
tc(1) = NewDate( year, month, day, hour )
tc(2) = tc(1) + IncrDate(hour=dth)
! check for strange values:
if ( (tr(1) < tc(1)) .or. (tc(2) < tr(1)) ) then
write (gol,'("determined invalid first interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tc(1), ' - ', tc(2) ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end if
! read primary field ?
if ( md%filled1 ) then
! read new field if times are different:
data1_read = (md%tr1(1) /= tc(1)) .or. (md%tr1(2) /= tc(2))
else
! not filled yet, thus must read:
data1_read = .true.
end if
if ( data1_read ) then
data1_tref = tc(1)
data1_t1 = tc(1)
data1_t2 = tc(2)
end if
! setup reading of secondary data only if end of requested
! interval is later than primary interval:
if ( tc(2) < tr(2) ) then
! extract time values for end of requested interval:
call Get( tr(2), year, month, day, hour, minu )
! round hour to next baseh + 00/03/06/09/12/15/18/21
hour = dth * floor(real(hour+minu/60.0-baseh)/real(dth)) + baseh
! interval with constant field
tc(1) = NewDate( year, month, day ) + IncrDate(hour=hour)
tc(2) = tc(1) + IncrDate(hour=dth)
! check for strange values:
if ( (tr(2) < tc(1)) .or. (tc(2) < tr(2)) ) then
write (gol,'("determined invalid second interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tc(1), ' - ', tc(2) ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end if
! read secondary field ?
if ( md%filled2 ) then
! read new field if times are different;
data2_read = (md%tr2(1) /= tc(1)) .or. (md%tr2(2) /= tc(2))
else
! not filled yet, thus must read:
data2_read = .true.
end if
if ( data2_read ) then
data2_tref = tc(1)
data2_t1 = tc(1)
data2_t2 = tc(2)
end if
end if ! tr partly after primary interval
!
! ** interpolated between 6 hourly times 00/06/12/18
! interpolated between 6 hourly times 03/09/15/21
! interpolated between 3 hourly times 00/03/06/09/12/15/18/21
!
case ( 'interp6', 'interp6_3', 'interp3', 'interp2', 'interp1' )
! extract time values for begin of current interval:
call Get( tr(1), year, month, day, hour, minu )
! truncate hour to previous 00/06/12/18, 03/09/15/21,
! or 00/03/06/09/12/15/18/21
hour = dth * floor(real(hour+minu/60.0-baseh)/real(dth)) + baseh
! set begin of 3 or 6 hour interval:
tprev = NewDate( year, month, day, hour )
! extract time values for end of current interval:
call Get( tr(2), year, month, day, hour, minu )
! truncate hour to previous 00/06/12/18
hour = dth * ceiling(real(hour+minu/60.0-baseh)/real(dth)) + baseh
! set end of 3 or 6 hour interval:
tnext = NewDate( year, month, day, hour )
! checks:
! [tprev,tmax] should be dth hours
! [tprev,tmax] should contain [tr(1),tr(2)]
dth_int = iTotal(tnext-tprev,'hour')
if ( (tr(1) < tprev) .or. (tnext < tr(2)) .or. &
( (dth_int /= 0) .and. (dth_int /= dth) ) ) then
write (gol,'("determined invalid interpolation interval:")'); call goErr
call wrtgol( ' requested : ', tr(1), ' - ', tr(2) ); call goErr
call wrtgol( ' guessed : ', tprev, ' - ', tnext ); call goErr
write (gol,'(" for tinterp : ",a)') md%tinterp; call goErr
call goErr; status=1; return
end if
!
! . <-- previous field at dth hours
! o <-- latest interpolated field
! x <-- target
! o <-- next field at dth hours
! tr1 tr tr2
! --+--------------+------
! tprev tnext
!
! read main field ?
if ( md%filled1 ) then
! md%data should be defined in [tprev,tr]
data1_read = (md%tr1(1) < tprev) .or. (tr(2) < md%tr1(1))
else
data1_read = .true.
end if
if ( data1_read ) then
data1_tref = tprev
data1_t1 = tprev
data1_t2 = tprev
end if
! read second field ?
if ( md%filled2 ) then
! md%data should be defined for tnext
data2_read = md%tr2(1) /= tnext
else
data2_read = .true.
end if
if ( data2_read ) then
data2_tref = tnext
data2_t1 = tnext
data2_t2 = tnext
end if
!
! ** error ...
!
case default
write (gol,'("unsupported time interpolation : ",a)') md%tinterp ; call goErr
call goErr; status=1; return
end select
!
! set ref times
!
if ( fcmode ) then
! in forecast mode, tfcday0 is 00:00 at the day the forecast starts;
data1_tref = tfcday0
data2_tref = tfcday0
else
! dummy tref's : begin of day in which [data?_t1,data?_t2] starts:
data1_tref = data1_t1
if ( IsAnyDate(data1_tref) ) data1_tref = tr(1)
call Set( data1_tref, hour=0, min=0, sec=0, mili=0 )
data2_tref = data2_t1
if ( IsAnyDate(data2_tref) ) data2_tref = tr(1)
call Set( data2_tref, hour=0, min=0, sec=0, mili=0 )
end if
!
! trap double reading
!
! data already in data2 ?
if ( data1_read .and. md%filled2 ) then
if ( (data1_t1 == md%tr2(1)) .and. (data1_t2 == md%tr2(2)) ) then
data1_read = .false.
data1_copy = .true.
end if
end if
! data2 just read ?
if ( data2_read .and. data1_read ) then
! data2 is same as data ?
if ( (data2_tref == data1_tref) .and. &
(data2_t1 == data1_t1) .and. (data2_t2 == data1_t2) ) then
data2_read = .false.
data2_copy = .true.
end if
end if
!write (gol,'("SetupSetup:")'); call goPr
!write (gol,'(" fcmode : ",l1)') fcmode; call goPr
!call wrtgol( ' tfcday0 : ', tfcday0 ); call goPr
!write (gol,'(" md%tinterp : ",a)') trim(md%tinterp); call goPr
!call wrtgol( ' tr(1) : ', tr(1) ); call goPr
!call wrtgol( ' tr(2) : ', tr(2) ); call goPr
!write (gol,'(" 1 read,copy : ",2l2)') data1_read, data1_copy; call goPr
!call wrtgol( ' 1 tref : ', data1_tref ); call goPr
!call wrtgol( ' 1 t1 : ', data1_t1 ); call goPr
!call wrtgol( ' 1 t2 : ', data1_t2 ); call goPr
!write (gol,'(" 2 read,copy : ",2l2)') data2_read, data2_copy; call goPr
!call wrtgol( ' 2 tref : ', data2_tref ); call goPr
!call wrtgol( ' 2 t1 : ', data2_t1 ); call goPr
!call wrtgol( ' 2 t2 : ', data2_t2 ); call goPr
! ok
status = 0
if (okdebug) call goLabel()
end subroutine SetupSetup
!EOC
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_2D
!
! !DESCRIPTION: Fill md%data1 and md%data2 of a 2D met field type (md), with
! data for date tr(1) and tr(2) respectively (and if needed)
! through reading or copying. Also write to disk the met field
! if requested.
!
! Then set md%data according to its type of interpolation (see
! TimeInterpolation in meteodata.F90).
! For constant type, %data => %data1.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_2D( region, md, tr, lli, nuv, status )
!
! !USES:
!
use GO, only : TDate, wrtgol
use Grid, only : TllGridInfo
use TMM, only : ReadField, Read_SP, Read_SR_OLS, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md ! met field
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
type(TDate), intent(in) :: tr(2) ! dates
type(TllGridInfo), intent(in) :: lli ! grid (GLOBAL)
character(len=1), intent(in) :: nuv ! staggering
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status ! return code
!
! !REVISION HISTORY:
! 4 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Setup_2d'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, pointer :: field(:,:) ! work array
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! debug
!write (gol,'(a," @ ",a)') trim(md%name),trim(lli%name); call goPr
! not changed by default
md%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read or copy
call SetupSetup( md, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! Read/write primary field
! -------------------------
if ( data1_read ) then
! test
if ( md%ls(1) /= md%ls(2) ) then
write (gol,'("SETUP_2D called instead of SETUP_3D, field is 3D:")'); call goErr
write (gol, '(" md%ls(1:2) : ",2i3)') md%ls; call goErr
status=1; IF_NOTOK_RETURN(status=1)
end if
! Need whole region for I/O on root. Dummy else.
IF (isRoot) THEN
ALLOCATE( field( im(region), jm(region)) )
ELSE
ALLOCATE( field(1,1) )
END IF
! Read/write
IOroot : IF (isRoot) THEN
select case ( md%name )
case ( 'sp', 'sps' )
! special routine for surface pressure
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
case ( 'srols' )
! special routine for Olsson surface roughness:
call Read_SR_OLS( tmmd, md%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
case default
! general field
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, lli, &
nuv, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
end select
! write meteofiles
if ( md%putout ) then
call WriteField( tmmd, md%destkey, &
md%tmi1, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, FIELD, status )
IF_NOTOK_RETURN(status=1)
end if
END IF IOroot
CALL SCATTER( dgrid(region), md%data1(:,:,1), FIELD, md%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE( FIELD )
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md%data1 = md%data2
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
end if
! -------------------------
! Read/write (or copy or nothing) secondary field
! -------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else.
IF (isRoot) THEN
ALLOCATE( field( im(region), jm(region)) )
ELSE
ALLOCATE( field(1,1) )
END IF
! Read/write
IOroot2: IF (isRoot) THEN
select case ( md%name )
case ( 'sp', 'sps' )
! special routine for surface pressure
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
case default
! general field
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, lli, &
nuv, FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
end select
! write meteo files
if ( md%putout ) then
call WriteField( tmmd, md%destkey, &
md%tmi2, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, FIELD, status )
IF_NOTOK_RETURN(status=1)
end if
END IF IOroot2
CALL SCATTER( dgrid(region), md%data2(:,:,1), FIELD, md%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE( FIELD )
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array
md%data2 = md%data1
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
end if
! -------------------------
! time interpolation
! -------------------------
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! done
! -------------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_2D
!EOC
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_2D_PARALLEL_IO
!
! !DESCRIPTION: Same as SETUP_2D, except reading is done by every processes.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_2D_PARALLEL_IO( region, md, tr, tdlli, nuv, status )
!
! !USES:
!
use GO, only : TDate, wrtgol
use Grid, only : TllGridInfo
use TMM, only : ReadField, Read_SP, Read_SR_OLS, WriteField
! use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md ! met field
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
type(TDate), intent(in) :: tr(2) ! dates
type(TllGridInfo), intent(in) :: tdlli ! dummy.. grid is already determined by the region
character(len=1), intent(in) :: nuv ! staggering
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status ! return code
!
! !REVISION HISTORY:
! 18 Oct 2013 - Ph. Le Sager - v0
!
!EOP
!------------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Setup_2d_parallel_io'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
integer :: i1, i2, j1, j2
real, pointer :: field(:,:) ! work array
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
md%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read or copy
call SetupSetup( md, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! Read/write primary field
! -------------------------
if ( data1_read ) then
! test
if ( md%ls(1) /= md%ls(2) ) then
write (gol,'("SETUP_2D called instead of SETUP_3D, field is 3D:")'); call goErr
write (gol, '(" md%ls(1:2) : ",2i3)') md%ls; call goErr
status=1; IF_NOTOK_RETURN(status=1)
end if
! could get those bounds from md% directly
call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
allocate( field( i1:i2, j1:j2) ) !! bonds are not strictly required, could as well do (i2-i1+1, ..)
! Read/write
select case ( md%name )
case ( 'sp', 'sps' )
! special routine for surface pressure
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli(region), FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
case ( 'srols' )
! special routine for Olsson surface roughness:
call Read_SR_OLS( tmmd, md%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli(region), FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
case default
! general field
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, lli(region), &
nuv, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
end select
md%data1(i1:i2, j1:j2, 1) = field
deallocate( field )
! write meteofiles
if ( md%putout ) then
write(gol,*)"writing of remapped met field not tested yet.. SKIPPED." ; call goErr
TRACEBACK; status=1; return
!!!!! NEED SOMETHING SIMILAR FOR DATA2 BELOW
! ! Need whole region for I/O on root. Dummy else.
! IF (isRoot) THEN
! ALLOCATE( field( im(region), jm(region)) )
! ELSE
! ALLOCATE( field(1,1) )
! END IF
!
! CALL GATHER( dgrid(region), md%data1(:,:,1), FIELD_GLOBAL, md%halo, status)
!
! IF (isRoot) THEN
!
! call WriteField( tmmd, md%destkey, &
! md%tmi1, trim(md%name), trim(md%unit), &
! data1_tref, data1_t1, data1_t2, &
! GLOBAL_lli(region), nuv, FIELD, status )
! IF_NOTOK_RETURN(status=1)
!
! END IF
! DEALLOCATE( FIELD )
end if
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md%data1 = md%data2
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
end if
! -------------------------
! Read/write (or copy or nothing) secondary field
! -------------------------
if ( data2_read ) then
call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
allocate( field( i1:i2, j1:j2) )
select case ( md%name )
case ( 'sp', 'sps' )
! special routine for surface pressure
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli(region), FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
case default
! general field
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, lli(region), &
nuv, FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
end select
md%data2(i1:i2, j1:j2, 1) = FIELD
deallocate( field )
! write meteo files
!TODO if ( md%putout ) then
!TODO
!TODO end if
!TODO
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array
md%data2 = md%data1
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
end if
! -------------------------
! time interpolation
! -------------------------
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! done
! -------------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_2D_PARALLEL_IO
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_3D
!
! !DESCRIPTION: same as SETUP_2D, but for 3D fields by accounting for levels
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_3D( region, md, tr, lli, nuv, levi, nw, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, ReadField, WriteField
use dims, only : im, jm
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md ! met field
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
type(TDate), intent(in) :: tr(2) ! dates
type(TllGridInfo), intent(in) :: lli ! grid
character(len=1), intent(in) :: nuv ! horiz. staggering
type(TLevelInfo), intent(in) :: levi ! levels
character(len=1), intent(in) :: nw ! vertical staggering
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status ! return code
!
! !REVISION HISTORY:
! 4 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Setup_3d'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer :: field(:,:,:) ! work array (data)
integer :: is(2), js(2) ! work arrays (bounds)
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! debug
!write (gol,'(a," @ ",a)') trim(md%name),trim(lli%name); call goPr
! not changed by default
md%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read or copy
call SetupSetup( md, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! Read/write primary field
! -------------------------
if ( data1_read ) then
! Need whole region for I/O on root. Dummy else. Allocate global array for I/O
is = (/1,im(region)/)
js = (/1,jm(region)/)
IF (isRoot) THEN
ALLOCATE( FIELD( is(1):is(2), js(1):js(2), md%ls(1):md%ls(2) ))
ELSE
ALLOCATE( FIELD(1,1,1) )
END IF
! Read/write on root
IOroot : IF (isRoot) THEN
! safety check
if ( data1_t2 /= data1_t1 ) then
! write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
! call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
! call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
! write (gol,'("please decide what to do with surface pressures ... ")'); call goErr
! call goErr; status=1; return
write (gol,'("WARNING - using instant surface pressure for regridding temporal averaged 3D field ...")'); call goPr
end if
! surface pressure
allocate( tmp_sp( is(1):is(2), js(1):js(2) ) )
! fill data
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, levi, nw, &
tmp_sp, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
! write meteo file
if ( md%putout ) then
call WriteField( tmmd, md%destkey, &
md%tmi1, 'sp', trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, levi, nw, &
tmp_sp, FIELD, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
END IF IOroot
CALL SCATTER( dgrid(region), md%data1, FIELD, md%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE( FIELD )
! data array is filled now
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md%data1 = md%data2
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
IF (isRoot) THEN
ALLOCATE(field(im(region), jm(region), md%ls(1):md%ls(2)))
ELSE
ALLOCATE(field(1,1,1))
END IF
! Read/write
IOroot2 : IF (isRoot) THEN
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, levi, nw, &
tmp_sp, FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( md%putout ) then
call WriteField( tmmd, md%destkey, &
md%tmi2, 'sp', trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, levi, nw, &
tmp_sp, FIELD, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
END IF IOroot2
CALL SCATTER( dgrid(region), md%data2, FIELD, md%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE( FIELD )
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array
md%data2 = md%data1
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
end if
! -------------------------
! time interpolation
! -------------------------
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! done
! -------------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_3D
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_3D_PARALLEL_IO
!
! !DESCRIPTION: same as SETUP_3D, except reading is done by every processes.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_3D_PARALLEL_IO( region, md, tr, tdlli, nuv, levi, nw, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, ReadField, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md ! met field
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
type(TDate), intent(in) :: tr(2) ! dates
type(TllGridInfo), intent(in) :: tdlli ! dummy.. grid is already determined by the region
character(len=1), intent(in) :: nuv ! horiz. staggering
type(TLevelInfo), intent(in) :: levi ! levels
character(len=1), intent(in) :: nw ! vertical staggering
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status ! return code
!
! !REVISION HISTORY:
! 18 Oct 2013 - Ph. Le Sager - v0
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Setup_3d_parallel_io'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
integer :: i1, i2, j1, j2
real, allocatable :: tmp_sp(:,:)
real, pointer :: field(:,:,:) ! work array
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
md%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read or copy
call SetupSetup( md, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! Read/write primary field
! -------------------------
if ( data1_read ) then
! could get those bounds from md% directly
call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
allocate( field( i1:i2, j1:j2, md%ls(1):md%ls(2)))
! safety check
if ( data1_t2 /= data1_t1 ) then
! write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr
! call wrtgol( ' data1_t1 : ', data1_t1 ) ; call goErr
! call wrtgol( ' data1_t2 : ', data1_t2 ) ; call goErr
! write (gol,'("please decide what to do with surface pressures ... ")') ; call goErr
! TRACEBACK; status=1; return
write (gol,'("WARNING - using instant surface pressure for regridding temporal averaged 3D field ...")'); call goPr
end if
! surface pressure
allocate( tmp_sp( i1:i2, j1:j2 ) )
! read data
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli(region), nuv, levi, nw, &
tmp_sp, FIELD, md%tmi1, status )
IF_NOTOK_RETURN(status=1)
md%data1(i1:i2, j1:j2, md%ls(1):md%ls(2)) = field
! write meteo file
if ( md%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! IF (isRoot) THEN
! ALLOCATE(field(im(region), jm(region), md%ls(1):md%ls(2)))
! ELSE
! ALLOCATE(field(1,1,1))
! END IF
!
!
! CALL gather( dgrid(region), md%data1, FIELD, md%halo, status)
! IF_NOTOK_RETURN(status=1)
!
! !! NEED global_lli and to also gather SP
!
! if (isRoot) then
! call WriteField( tmmd, md%destkey, &
! md%tmi1, 'sp', trim(md%name), trim(md%unit), &
! data1_tref, data1_t1, data1_t2, &
! GLOBAL_lli(region), nuv, levi, nw, &
! tmp_sp, FIELD, status )
! IF_NOTOK_RETURN(status=1)
! end if
!
endif
DEALLOCATE( TMP_SP )
DEALLOCATE( FIELD )
! data array is filled now
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md%data1 = md%data2
! data array is filled now:
md%filled1 = .true.
md%tr1(1) = data1_t1
md%tr1(2) = data1_t2
md%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! could get those bounds from md% directly
call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
allocate( field( i1:i2, j1:j2, md%ls(1):md%ls(2)))
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ) ; call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ) ; call goErr
write (gol,'("please deceide what to do with surface pressures ... ")') ; call goErr
TRACEBACK; status=1; return
end if
! surface pressure
allocate( tmp_sp( i1:i2, j1:j2 ) )
! read data
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli(region), nuv, levi, nw, &
tmp_sp, FIELD, md%tmi2, status )
IF_NOTOK_RETURN(status=1)
md%data2(i1:i2, j1:j2, md%ls(1):md%ls(2)) = field
! write meteofiles
if ( md%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry. SKIPPED." ; call goErr
TRACEBACK; status=1; return
! IF (isRoot) THEN
! ALLOCATE(field(im(region), jm(region), md%ls(1):md%ls(2)))
! ELSE
! ALLOCATE(field(1,1,1))
! END IF
!
! CALL gather( dgrid(region), md%data2, FIELD, md%halo, status)
! IF_NOTOK_RETURN(status=1)
!
! !! NEED to also gather SP
!
! if(isroot)then
! call WriteField( tmmd, md%destkey, &
! md%tmi2, 'sp', trim(md%name), trim(md%unit), &
! data2_tref, data2_t1, data2_t2, &
! global_lli(region), nuv, levi, nw, &
! tmp_sp, FIELD, status )
! IF_NOTOK_RETURN(status=1)
! END IF
!
end if
! clear
DEALLOCATE( TMP_SP )
DEALLOCATE( FIELD )
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array
md%data2 = md%data1
! data array is filled now
md%filled2 = .true.
md%tr2(1) = data2_t1
md%tr2(2) = data2_t2
end if
! -------------------------
! time interpolation
! -------------------------
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
! -------------------------
! done
! -------------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_3D_PARALLEL_IO
!EOC
! **************************************************************
! ***
! *** Specific SETUP routines for MASS FLUXES
! ***
! **************************************************************
SUBROUTINE SETUP_MFUV_SERIAL_IO( region, md_mfu, md_mfv, tr, lli, levi, status )
! Set up MFU and MFV (horizontal fluxes)
! Read or copy %data1 and %data2, and get %data through time interpolation
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_MFUV, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: md_mfu
type(TMeteoData), intent(inout) :: md_mfv
type(TDate), intent(in) :: tr(2) ! time range
type(TllGridInfo), intent(in) :: lli
type(TLevelInfo), intent(in) :: levi
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/SETUP_MFUV_SERIAL_IO'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
logical :: NorthBorder, WestBorder ! tile location
real, allocatable :: tmp_spu(:,:)
real, allocatable :: tmp_spv(:,:)
! to read the entire region
real, pointer :: wrld_u(:,:,:), wrld_v(:,:,:), wrkarr(:,:,:)
integer, dimension(2) :: is, js, ls
integer :: halo, i1, i2, j1, j2
real, allocatable :: bigIslice(:,:), bigJslice(:,:), Islice(:,:), Jslice(:,:)
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( md_mfu%used .neqv. md_mfv%used ) then
write (gol,'("either none or both mfu and mfv should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. md_mfu%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
md_mfu%changed = .false.
md_mfv%changed = .false.
! local indices and tile location
CALL GET_DISTGRID( dgrid(region), &
I_STRT=i1, I_STOP=i2, &
J_STRT=j1, J_STOP=j2, &
hasWestBorder=WestBorder, hasNorthBorder=NorthBorder)
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from mfu only)
call SetupSetup( md_mfu, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Use fact that mfu and mfv have been allocated with the same bounds and halo
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = md_mfu%ls
halo = md_mfu%halo
IF (isRoot) THEN
ALLOCATE( wrld_u( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) )
ALLOCATE( wrld_v( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) )
ALLOCATE( wrkarr( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
wrld_v = 0.
wrld_u = 0.
allocate( bigIslice(jm(region),ls(1):ls(2)))
allocate( bigJslice(im(region),ls(1):ls(2)))
ELSE
ALLOCATE(wrld_u(1,1,1), wrld_v(1,1,1), wrkarr(1,1,1))
ALLOCATE( bigIslice(1,1), bigJslice(1,1) )
END IF
ALLOCATE( Islice(j1:j2,ls(1):ls(2)) )
ALLOCATE( Jslice(i1:i2,ls(1):ls(2)) )
if (isRoot) then ! only root does IO
! safety check ...
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_spu(is(1)-1:is(2),js(1):js(2) ) )
allocate( tmp_spv(is(1) :is(2),js(1):js(2)+1) )
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:imr,1:jmr ,1:lmr) in pu_t(0:imr+1,0:jmr+1,0:lmr)
! pv_tmpp --> pv(1:imr,1:jmr+1,1:lmr) in pv_t(0:imr+1,0:jmr+1,0:lmr)
! fill data:
call Read_MFUV( tmmd, md_mfu%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_spu, &
wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
md_mfu%tmi1, &
tmp_spv, &
wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
md_mfv%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( md_mfu%putout ) then
call WriteField( tmmd, md_mfu%destkey, &
md_mfu%tmi1, 'spu', trim(md_mfu%name), trim(md_mfu%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'u', levi, 'n', &
tmp_spu, wrld_u(is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
if ( md_mfv%putout ) then
call WriteField( tmmd, md_mfv%destkey, &
md_mfv%tmi1, 'spv', trim(md_mfv%name), trim(md_mfv%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'v', levi, 'n', &
tmp_spv, wrld_v(is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_spu )
deallocate( tmp_spv )
end if ! root ?
! Scatter U
if(isRoot) wrkarr = wrld_u(is(1):is(2),js(1):js(2),:)
CALL SCATTER( dgrid(region), md_mfu%data1, wrkarr, md_mfu%halo, status)
IF_NOTOK_RETURN(status=1)
! manually scatter wrld_u(is(1)-1,:,:). This is needed only with non-cyclic
! zoom regions, since any update_halo will overwrite is(1)-1. [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = wrld_u(0,js(1):js(2),:)
CALL SCATTER_I_BAND( dgrid(region), islice, bigIslice, status, iref=1)
if (WestBorder) md_mfu%data1(0,j1:j2,:) = islice
! Scatter V
if(isRoot) wrkarr = wrld_v(is(1):is(2),js(1):js(2),:)
CALL SCATTER( dgrid(region), md_mfv%data1, wrkarr, md_mfv%halo, status)
IF_NOTOK_RETURN(status=1)
! manually SCATTER wrld_v( :, js(2)+1 , :) : NORTH POLE HALO
if(isroot) bigJslice=wrld_v(is(1):is(2),jm(region)+1,:)
CALL SCATTER_J_BAND( dgrid(region), jslice, bigJslice, status, jref=jm(region))
if (NorthBorder) md_mfv%data1(i1:i2,jm(region)+1,:)=jslice
deallocate(wrld_u, wrld_v, wrkarr, bigIslice, bigJslice, Islice, Jslice)
! data array is filled now:
md_mfu%filled1 = .true.
md_mfu%tr1(1) = data1_t1
md_mfu%tr1(2) = data1_t2
md_mfu%changed = .true.
md_mfv%filled1 = .true.
md_mfv%tr1(1) = data1_t1
md_mfv%tr1(2) = data1_t2
md_mfv%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md_mfu%data1 = md_mfu%data2
md_mfv%data1 = md_mfv%data2
! data array is filled now:
md_mfu%filled1 = .true.
md_mfu%tr1(1) = data1_t1
md_mfu%tr1(2) = data1_t2
md_mfu%changed = .true.
md_mfv%filled1 = .true.
md_mfv%tr1(1) = data1_t1
md_mfv%tr1(2) = data1_t2
md_mfv%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = md_mfu%ls
halo = md_mfu%halo
IF (isRoot) THEN
allocate( wrld_u( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) )
allocate( wrld_v( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) )
allocate( wrkarr( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
wrld_v = 0.
wrld_u = 0.
allocate( bigIslice(jm(region),ls(1):ls(2)))
allocate( bigJslice(im(region),ls(1):ls(2)))
ELSE
ALLOCATE(wrld_u(1,1,1), wrld_v(1,1,1), wrkarr(1,1,1))
ALLOCATE( bigIslice(1,1), bigJslice(1,1) )
END IF
ALLOCATE( Islice(j1:j2,ls(1):ls(2)) )
ALLOCATE( Jslice(i1:i2,ls(1):ls(2)) )
if (isRoot) then ! only root does IO
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_spu(is(1)-1:is(2),js(1):js(2) ) )
allocate( tmp_spv(is(1) :is(2),js(1):js(2)+1) )
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:imr,1:jmr ,1:lmr) in pu_t(0:imr+1,0:jmr+1,0:lmr)
! pv_tmpp --> pv(1:imr,1:jmr+1,1:lmr) in pv_t(0:imr+1,0:jmr+1,0:lmr)
! fill data:
call Read_MFUV( tmmd, md_mfu%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
tmp_spu, &
wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
md_mfu%tmi2, &
tmp_spv, &
wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
md_mfv%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( md_mfu%putout ) then
call WriteField( tmmd, md_mfu%destkey, &
md_mfu%tmi2, 'spu', trim(md_mfu%name), trim(md_mfu%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'u', levi, 'n', &
tmp_spu, wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
endif
if ( md_mfv%putout ) then
call WriteField( tmmd, md_mfv%destkey, &
md_mfv%tmi2, 'spv', trim(md_mfv%name), trim(md_mfv%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'v', levi, 'n', &
tmp_spv, wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_spu )
deallocate( tmp_spv )
end if ! root
! Scatter U
if(isRoot) wrkarr = wrld_u(is(1):is(2),js(1):js(2),:)
CALL SCATTER( dgrid(region), md_mfu%data2, wrkarr, md_mfu%halo, status)
IF_NOTOK_RETURN(status=1)
! important for zoom regions only, since any update_halo will overwrite is(1)-1. [FIXME: could had a
! test around these 3 lines ]
if(isRoot) bigIslice = wrld_u(0,js(1):js(2),:)
CALL SCATTER_I_BAND( dgrid(region), islice, bigIslice, status, iref=1)
if (WestBorder) md_mfu%data2(0,j1:j2,:) = islice
! Scatter V
if(isRoot) wrkarr = wrld_v(is(1):is(2),js(1):js(2),:)
CALL SCATTER( dgrid(region), md_mfv%data2, wrkarr, md_mfv%halo, status)
IF_NOTOK_RETURN(status=1)
! manually SCATTER wrld_v( :, js(2)+1 , :) : NORTH POLE HALO
if(isroot) bigJslice=wrld_v(is(1):is(2),jm(region)+1,:)
CALL SCATTER_J_BAND( dgrid(region), jslice, bigJslice, status, jref=jm(region))
if (NorthBorder) md_mfv%data2(i1:i2,jm(region)+1,:)=jslice
DEALLOCATE(wrld_u, wrld_v, wrkarr, bigIslice, bigJslice, Islice, Jslice)
! data array is filled now:
md_mfu%filled2 = .true.
md_mfu%tr2(1) = data2_t1
md_mfu%tr2(2) = data2_t2
md_mfv%filled2 = .true.
md_mfv%tr2(1) = data2_t1
md_mfv%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from primary array:
md_mfu%data2 = md_mfu%data
md_mfv%data2 = md_mfv%data
! data array is filled now:
md_mfu%filled2 = .true.
md_mfu%tr2(1) = data2_t1
md_mfu%tr2(2) = data2_t2
md_mfv%filled2 = .true.
md_mfv%tr2(1) = data2_t1
md_mfv%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( md_mfu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( md_mfv, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_MFUV_SERIAL_IO
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_MFUV_PARALLEL_IO
!
! !DESCRIPTION: same as setup_mfuv, but with parallel I/O.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_MFUV_PARALLEL_IO( region, md_mfu, md_mfv, tr, levi, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_MFUV, WriteField
use dims, only : im, jm
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md_mfu
type(TMeteoData), intent(inout) :: md_mfv
type(TDate), intent(in) :: tr(2) ! time range
type(TLevelInfo), intent(in) :: levi
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 24 Oct 2013 - Ph. Le Sager - v0
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/SETUP_MFUV_PARALLEL_IO'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
logical :: NorthBorder, WestBorder ! tile location
real, allocatable :: tmp_spu(:,:)
real, allocatable :: tmp_spv(:,:)
real, pointer :: wrld_u(:,:,:), wrld_v(:,:,:)
integer, dimension(2) :: ls
integer :: halo, i1, i2, j1, j2
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( md_mfu%used .neqv. md_mfv%used ) then
write (gol,'("either none or both mfu and mfv should be in use")'); call goErr
TRACEBACK; status=1; return
end if
if ( .not. md_mfu%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
md_mfu%changed = .false.
md_mfv%changed = .false.
! local indices and tile location
CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, &
hasWestBorder=WestBorder, hasNorthBorder=NorthBorder)
!-------------------------
! time stuff & work arrays
!-------------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from mfu only)
call SetupSetup( md_mfu, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
if (data2_read .or. data1_read) then
ls = md_mfu%ls
halo = md_mfu%halo
ALLOCATE( wrld_u( i1-halo:i2+halo, j1-halo:j2+halo, ls(1):ls(2)) )
ALLOCATE( wrld_v( i1-halo:i2+halo, j1-halo:j2+halo, ls(1):ls(2)) )
wrld_v = 0.
wrld_u = 0.
! surface pressure field:
allocate( tmp_spu(i1-1:i2, j1:j2 ) )
allocate( tmp_spv(i1 :i2, j1:j2+1) )
end if
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! safety check
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
TRACEBACK; status=1; return
end if
! NOTE: strange old indexing:
! pu_tmpp --> pu(0:imr,1:jmr ,1:lmr) in pu_t(0:imr+1,0:jmr+1,0:lmr)
! pv_tmpp --> pv(1:imr,1:jmr+1,1:lmr) in pv_t(0:imr+1,0:jmr+1,0:lmr)
call Read_MFUV( tmmd, md_mfu%sourcekey, &
data1_tref, data1_t1, data1_t2, lli(region), levi, &
tmp_spu, &
wrld_u( i1-1:i2, j1:j2, ls(1)+1:ls(2) ), &
md_mfu%tmi1, &
tmp_spv, &
wrld_v( i1:i2, j1:j2+1, ls(1)+1:ls(2) ), &
md_mfv%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( md_mfu%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfu%destkey, &
! md_mfu%tmi1, 'spu', trim(md_mfu%name), trim(md_mfu%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'u', levi, 'n', &
! tmp_spu, wrld_u(is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
! status )
! IF_NOTOK_RETURN(status=1)
end if
if ( md_mfv%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfv%destkey, &
! md_mfv%tmi1, 'spv', trim(md_mfv%name), trim(md_mfv%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'v', levi, 'n', &
! tmp_spv, wrld_v(is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
! status )
! IF_NOTOK_RETURN(status=1)
end if
md_mfu%data1( i1-1:i2, j1:j2, ls(1):ls(2)) = wrld_u(i1-1:i2, j1:j2, :)
md_mfv%data1( i1 :i2, j1:j2+1, ls(1):ls(2)) = wrld_v(i1 :i2, j1:j2+1, :)
! data array is filled now:
md_mfu%filled1 = .true.
md_mfu%tr1(1) = data1_t1
md_mfu%tr1(2) = data1_t2
md_mfu%changed = .true.
md_mfv%filled1 = .true.
md_mfv%tr1(1) = data1_t1
md_mfv%tr1(2) = data1_t2
md_mfv%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md_mfu%data1 = md_mfu%data2
md_mfv%data1 = md_mfv%data2
! data array is filled now:
md_mfu%filled1 = .true.
md_mfu%tr1(1) = data1_t1
md_mfu%tr1(2) = data1_t2
md_mfu%changed = .true.
md_mfv%filled1 = .true.
md_mfv%tr1(1) = data1_t1
md_mfv%tr1(2) = data1_t2
md_mfv%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! safety check
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
TRACEBACK; status=1; return
end if
! fill data
call Read_MFUV( tmmd, md_mfu%sourcekey, &
data2_tref, data2_t1, data2_t2, lli(region), levi, &
tmp_spu, &
wrld_u( i1-1:i2, j1:j2, ls(1)+1:ls(2) ), &
md_mfu%tmi2, &
tmp_spv, &
wrld_v( i1:i2, j1:j2+1, ls(1)+1:ls(2) ), &
md_mfv%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
md_mfu%data2( i1-1:i2, j1:j2, ls(1):ls(2)) = wrld_u(i1-1:i2, j1:j2, :)
md_mfv%data2( i1 :i2, j1:j2+1, ls(1):ls(2)) = wrld_v(i1 :i2, j1:j2+1, :)
! write meteofiles
if ( md_mfu%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfu%destkey, &
! md_mfu%tmi2, 'spu', trim(md_mfu%name), trim(md_mfu%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'u', levi, 'n', &
! tmp_spu, wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), &
! status )
! IF_NOTOK_RETURN(status=1)
endif
if ( md_mfv%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfv%destkey, &
! md_mfv%tmi2, 'spv', trim(md_mfv%name), trim(md_mfv%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'v', levi, 'n', &
! tmp_spv, wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), &
! status )
! IF_NOTOK_RETURN(status=1)
end if
! data array is filled now:
md_mfu%filled2 = .true.
md_mfu%tr2(1) = data2_t1
md_mfu%tr2(2) = data2_t2
md_mfv%filled2 = .true.
md_mfv%tr2(1) = data2_t1
md_mfv%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from primary array:
md_mfu%data2 = md_mfu%data
md_mfv%data2 = md_mfv%data
! data array is filled now:
md_mfu%filled2 = .true.
md_mfu%tr2(1) = data2_t1
md_mfu%tr2(2) = data2_t2
md_mfv%filled2 = .true.
md_mfv%tr2(1) = data2_t1
md_mfv%tr2(2) = data2_t2
end if
! Clean
if (data2_read .or. data1_read) then
deallocate( tmp_spu )
deallocate( tmp_spv )
deallocate( wrld_u, wrld_v )
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( md_mfu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( md_mfv, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_MFUV_PARALLEL_IO
!EOC
! ***
subroutine Setup_MFW_serial_io( region, md_mfw, md_tsp, tr, lli, nuv, levi, nw, status )
! Set up MFW (vertical flux) and TSP (tendency surface pressure)
! Read or copy %data1 and %data2, and get %data through time interpolation
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, ReadField, Read_MFW, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: md_mfw
type(TMeteoData), intent(inout) :: md_tsp
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
character(len=1), intent(in) :: nuv
type(TLevelInfo), intent(in) :: levi
character(len=1), intent(in) :: nw
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Setup_MFW_serial_io'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer :: mfw(:,:,:), tsp(:,:) ! work arrays (data)
integer :: is(2), js(2), ls(2), halo ! work arrays (bounds)
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md_mfw%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! error if tsp is not in use ...
if ( .not. md_tsp%used ) then
write (gol,'("mfw is in use but tsp not ..")'); call goErr
if (okdebug) call goLabel()
status=1; return
end if
! not changed by default
md_mfw%changed = .false.
md_tsp%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read
call SetupSetup( md_mfw, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = md_mfw%ls
IF (isRoot) THEN
allocate( mfw(is(1):is(2), js(1):js(2), ls(1):ls(2) ))
allocate( tsp(is(1):is(2), js(1):js(2)) )
ELSE
allocate( mfw(1,1,1), tsp(1,1) )
END IF
if (isRoot) then ! only root does I/O
! safety check ...
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data:
call Read_MFW( tmmd, md_mfw%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli, levi, &
tmp_sp, mfw, &
tsp, &
md_mfw%tmi1, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( md_mfw%putout ) then
call WriteField( tmmd, md_mfw%destkey, &
md_mfw%tmi1, 'sp', trim(md_mfw%name), trim(md_mfw%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, levi, nw, &
tmp_sp, mfw, status )
IF_NOTOK_RETURN(status=1)
end if
if ( md_tsp%putout ) then
! use history from mfw ...
call WriteField( tmmd, md_tsp%destkey, &
md_mfw%tmi1, trim(md_tsp%name), trim(md_tsp%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, tsp, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root
CALL SCATTER( dgrid(region), md_mfw%data1, MFW, md_mfw%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), md_tsp%data1(:,:,1), TSP, md_tsp%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(MFW, TSP)
! data array is filled now:
md_mfw%filled1 = .true.
md_mfw%tr1(1) = data1_t1
md_mfw%tr1(2) = data1_t2
md_mfw%changed = .true.
!
md_tsp%filled1 = .true.
md_tsp%tr1(1) = data1_t1
md_tsp%tr1(2) = data1_t2
md_tsp%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md_mfw%data1 = md_mfw%data2
! data array is filled now:
md_mfw%filled1 = .true.
md_mfw%tr1(1) = data1_t1
md_mfw%tr1(2) = data1_t2
md_mfw%changed = .true.
!
md_tsp%filled1 = .true.
md_tsp%tr1(1) = data1_t1
md_tsp%tr1(2) = data1_t2
md_tsp%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = md_mfw%ls
IF (isRoot) THEN
allocate( mfw(is(1):is(2), js(1):js(2), ls(1):ls(2) ))
allocate( tsp(is(1):is(2), js(1):js(2)) )
ELSE
allocate( mfw(1,1,1), tsp(1,1) )
END IF
if (isRoot) then ! only root does IO
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data:
call Read_MFW( tmmd, md_mfw%sourcekey, &
data2_tref, data2_t1, data2_t2, &
lli, levi, tmp_sp, MFW, TSP, md_mfw%tmi2, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( md_mfw%putout ) then
call WriteField( tmmd, md_mfw%destkey, &
md_mfw%tmi2, 'sp', trim(md_mfw%name), trim(md_mfw%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, levi, nw, &
tmp_sp, MFW, status )
IF_NOTOK_RETURN(status=1)
end if
if ( md_tsp%putout ) then
! use history from mfw ...
call WriteField( tmmd, md_tsp%destkey, &
md_mfw%tmi2, trim(md_tsp%name), trim(md_tsp%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, TSP, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root
CALL SCATTER( dgrid(region), md_mfw%data2, MFW, md_mfw%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), md_tsp%data2(:,:,1), TSP, md_tsp%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(MFW, TSP)
! data array is filled now:
md_mfw%filled2 = .true.
md_mfw%tr2(1) = data2_t1
md_mfw%tr2(2) = data2_t2
!
md_tsp%filled2 = .true.
md_tsp%tr2(1) = data2_t1
md_tsp%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array:
md_mfw%data2 = md_mfw%data1
! data array is filled now:
md_mfw%filled2 = .true.
md_mfw%tr2(1) = data2_t1
md_mfw%tr2(2) = data2_t2
!
md_tsp%filled2 = .true.
md_tsp%tr2(1) = data2_t1
md_tsp%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( md_mfw, tr, status )
IF_NOTOK_RETURN(status=1)
!
call TimeInterpolation( md_tsp, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
end subroutine Setup_MFW_serial_io
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_MFW_PARALLEL_IO
!
! !DESCRIPTION: Same as SETUP_MFW_SERIAL_IO, but with parallel I/O :
!
! Set up MFW (vertical flux) and TSP (tendency surface pressure)
! Read or copy %data1 and %data2, and get %data through time interpolation
!
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_MFW_PARALLEL_IO( region, md_mfw, md_tsp, tr, nuv, levi, nw, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, ReadField, Read_MFW, WriteField
use dims, only : im, jm
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md_mfw
type(TMeteoData), intent(inout) :: md_tsp
type(TDate), intent(in) :: tr(2)
character(len=1), intent(in) :: nuv
type(TLevelInfo), intent(in) :: levi
character(len=1), intent(in) :: nw
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 24 Oct 2013 - Ph. Le Sager - v0
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Setup_MFW_PARALLEL_io'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer :: mfw(:,:,:), tsp(:,:) ! work arrays (data)
integer :: is(2), js(2), ls(2), halo ! work arrays (bounds)
integer :: i1, i2, j1, j2
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( .not. md_mfw%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! error if tsp is not in use ...
if ( .not. md_tsp%used ) then
write (gol,'("mfw is in use but tsp not ..")'); call goErr
if (okdebug) call goLabel()
status=1; return
end if
! not changed by default
md_mfw%changed = .false.
md_tsp%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read
call SetupSetup( md_mfw, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! work arrays
if (data1_read .or. data2_read) then
CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
is = (/i1,i2/)
js = (/j1,j2/)
ls = md_mfw%ls
allocate( mfw (is(1):is(2), js(1):js(2), ls(1):ls(2) ))
allocate( tsp (is(1):is(2), js(1):js(2) ))
allocate( tmp_sp(is(1):is(2), js(1):js(2) ))
end if
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! safety check
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ) ; call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ) ; call goErr
write (gol,'("please deceide what to do with surface pressures ... ")') ; call goErr
TRACEBACK; status=1; return
end if
! fill data
call Read_MFW( tmmd, md_mfw%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli(region), levi, &
tmp_sp, mfw, &
tsp, &
md_mfw%tmi1, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( md_mfw%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfw%destkey, &
! md_mfw%tmi1, 'sp', trim(md_mfw%name), trim(md_mfw%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli(region), nuv, levi, nw, &
! tmp_sp, mfw, status )
! IF_NOTOK_RETURN(status=1)
end if
if ( md_tsp%putout ) then
! use history from mfw ...
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
! TRACEBACK; status=1; return
! call WriteField( tmmd, md_tsp%destkey, &
! md_mfw%tmi1, trim(md_tsp%name), trim(md_tsp%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli(region), nuv, tsp, status )
! IF_NOTOK_RETURN(status=1)
end if
md_mfw%data1( i1:i2, j1:j2, ls(1):ls(2)) = mfw
md_tsp%data1( i1:i2, j1:j2, 1 ) = tsp
! data array is filled now:
md_mfw%filled1 = .true.
md_mfw%tr1(1) = data1_t1
md_mfw%tr1(2) = data1_t2
md_mfw%changed = .true.
!
md_tsp%filled1 = .true.
md_tsp%tr1(1) = data1_t1
md_tsp%tr1(2) = data1_t2
md_tsp%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md_mfw%data1 = md_mfw%data2
! data array is filled now:
md_mfw%filled1 = .true.
md_mfw%tr1(1) = data1_t1
md_mfw%tr1(2) = data1_t2
md_mfw%changed = .true.
!
md_tsp%filled1 = .true.
md_tsp%tr1(1) = data1_t1
md_tsp%tr1(2) = data1_t2
md_tsp%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! safety check
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! fill data:
call Read_MFW( tmmd, md_mfw%sourcekey, &
data2_tref, data2_t1, data2_t2, &
lli(region), levi, tmp_sp, MFW, TSP, md_mfw%tmi2, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( md_mfw%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! call WriteField( tmmd, md_mfw%destkey, &
! md_mfw%tmi2, 'sp', trim(md_mfw%name), trim(md_mfw%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, nuv, levi, nw, &
! tmp_sp, MFW, status )
! IF_NOTOK_RETURN(status=1)
end if
if ( md_tsp%putout ) then
write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr
TRACEBACK; status=1; return
! ! use history from mfw ...
! call WriteField( tmmd, md_tsp%destkey, &
! md_mfw%tmi2, trim(md_tsp%name), trim(md_tsp%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, nuv, TSP, status )
! IF_NOTOK_RETURN(status=1)
end if
md_mfw%data2( i1:i2, j1:j2, ls(1):ls(2)) = mfw
md_tsp%data2( i1:i2, j1:j2, 1 ) = tsp
! data array is filled now:
md_mfw%filled2 = .true.
md_mfw%tr2(1) = data2_t1
md_mfw%tr2(2) = data2_t2
!
md_tsp%filled2 = .true.
md_tsp%tr2(1) = data2_t1
md_tsp%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from secondary array:
md_mfw%data2 = md_mfw%data1
! data array is filled now:
md_mfw%filled2 = .true.
md_mfw%tr2(1) = data2_t1
md_mfw%tr2(2) = data2_t2
!
md_tsp%filled2 = .true.
md_tsp%tr2(1) = data2_t1
md_tsp%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( md_mfw, tr, status )
IF_NOTOK_RETURN(status=1)
!
call TimeInterpolation( md_tsp, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
if (data1_read .or. data2_read) then
deallocate(mfw, tsp)
deallocate( tmp_sp )
end if
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_MFW_PARALLEL_IO
!EOC
! **************************************************************
! ***
! *** temperature and humidity
! ***
! **************************************************************
subroutine Setup_TQ( region, md_T, md_Q, tr, lli, levi, status)
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_TQ, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: md_T
type(TMeteoData), intent(inout) :: md_Q
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
type(TLevelInfo), intent(in) :: levi
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Setup_TQ'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer :: T(:,:,:), Q(:,:,:) ! work array
integer :: is(2), js(2) ! work arrays (bounds)
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( md_T%used .neqv. md_Q%used ) then
write (gol,'("either none or both T and Q should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. md_T%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
md_T%changed = .false.
md_Q%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from T only)
call SetupSetup( md_T, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
IF (isRoot) THEN
ALLOCATE( T(is(1):is(2), js(1):js(2), md_T%ls(1):md_T%ls(2) ))
ALLOCATE( Q(is(1):is(2), js(1):js(2), md_Q%ls(1):md_Q%ls(2) ))
ELSE
ALLOCATE( T(1,1,1), Q(1,1,1) )
END IF
if (isRoot) then ! only root does IO
! safety check ...
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp( is(1):is(2), js(1):js(2) ) )
! fill data:
call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_sp, &
T, md_T%tmi1, &
Q, md_Q%tmi1, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( md_T%putout ) then
call WriteField( tmmd, md_T%destkey, &
md_T%tmi1, 'sp', trim(md_T%name), trim(md_T%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'n', &
tmp_sp, T, status )
IF_NOTOK_RETURN(status=1)
end if
if ( md_Q%putout ) then
call WriteField( tmmd, md_Q%destkey, &
md_Q%tmi1, 'sp', trim(md_Q%name), trim(md_Q%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'n', &
tmp_sp, Q, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! Distribute
CALL SCATTER( dgrid(region), md_T%data1, T, md_T%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), md_Q%data1, Q, md_Q%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(T, Q)
! data array is filled now:
md_T%filled1 = .true.
md_T%tr1(1) = data1_t1
md_T%tr1(2) = data1_t2
md_T%changed = .true.
md_Q%filled1 = .true.
md_Q%tr1(1) = data1_t1
md_Q%tr1(2) = data1_t2
md_Q%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
md_T%data1 = md_T%data2
md_Q%data1 = md_Q%data2
! data array is filled now:
md_T%filled1 = .true.
md_T%tr1(1) = data1_t1
md_T%tr1(2) = data1_t2
md_T%changed = .true.
md_Q%filled1 = .true.
md_Q%tr1(1) = data1_t1
md_Q%tr1(2) = data1_t2
md_Q%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
IF (isRoot) THEN
allocate( T(is(1):is(2), js(1):js(2), md_T%ls(1):md_T%ls(2) ))
allocate( Q(is(1):is(2), js(1):js(2), md_Q%ls(1):md_Q%ls(2) ))
ELSE
allocate( T(1,1,1), Q(1,1,1) )
END IF
if (isRoot) then ! only root does IO
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp( is(1):is(2), js(1):js(2)) )
! fill data:
call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
tmp_sp, &
T, md_T%tmi2, &
Q, md_Q%tmi2, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( md_T%putout ) then
call WriteField( tmmd, md_T%destkey, &
md_T%tmi2, 'sp', trim(md_T%name), trim(md_T%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'n', &
tmp_sp, T, status )
IF_NOTOK_RETURN(status=1)
endif
if ( md_Q%putout ) then
call WriteField( tmmd, md_Q%destkey, &
md_Q%tmi2, 'sp', trim(md_Q%name), trim(md_Q%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'n', &
tmp_sp, Q, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root
CALL SCATTER( dgrid(region), md_T%data2, T, md_T%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), md_Q%data2, Q, md_Q%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(T, Q)
! data array is filled now:
md_T%filled2 = .true.
md_T%tr2(1) = data2_t1
md_T%tr2(2) = data2_t2
md_Q%filled2 = .true.
md_Q%tr2(1) = data2_t1
md_Q%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data from primary array:
md_T%data2 = md_T%data1
md_Q%data2 = md_Q%data1
! data array is filled now:
md_T%filled2 = .true.
md_T%tr2(1) = data2_t1
md_T%tr2(2) = data2_t2
md_Q%filled2 = .true.
md_Q%tr2(1) = data2_t1
md_Q%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( md_T, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( md_Q, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
end subroutine Setup_TQ
! ***
! subroutine Meteo_SetupMass( n, status )
!
! use global_data, only : mass_dat
! use dims, only : newsrun
! use dims, only : xcyc, im, jm
! use geometry, only : geomtryv
!
! ! --- in/out -----------------------------
!
! integer, intent(in) :: n ! region
! integer, intent(out) :: status
!
! ! --- const --------------------------------------
!
! character(len=*), parameter :: rname = mname//'/Meteo_SetupMass'
!
! ! --- begin ------------------------------
!
! call goLabel(rname)
!
! ! compute initial pressure levels and mass ?
! if ( newsrun ) then
! call geomtryv( n )
! end if
!
! ! periodic boundary for m
! ! NOTE: m has been advected or created by geomtryv
! if ( xcyc(n) == 1 ) then
! mass_dat(n)%m_t(0 ,:,:) = mass_dat(n)%m_t(im(n),:,:)
! mass_dat(n)%m_t(im(n)+1,:,:) = mass_dat(n)%m_t(1 ,:,:)
! end if
!
! ! ok
! status = 0
! call goLabel()
!
! end subroutine Meteo_SetupMass
! ***
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: METEO_CHECKPRESSURE
!
! !DESCRIPTION: Compute difference b/w sp1_dat (read) and sp_dat (advected),
! and compare to threshold.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE METEO_CHECKPRESSURE( n, status )
!
! !USES:
!
use ParTools, only : Par_Reduce
use dims, only : idate, newsrun
use dims, only : xcyc, im, jm
use redgridZoom, only : calc_pdiff
#ifdef with_hdf4
use io_hdf, only : io_write2d_32d, DFACC_CREATE
#endif
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: n ! region
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Meteo_CheckPressure'
! maximum accepted pressure difference:
real, parameter :: pdiffmax_treshold = 1.0e2 ! Pa
! --- external -------------------------
integer(4), external :: sfStart, sfEnd
! --- local -----------------------------
real :: pdiffmax, pdiffmax_l
integer(4) :: io
! --- begin ------------------------------
if (okdebug) call goLabel(rname)
! compare 'advected' pressure with read pressure
if ( .not. newsrun ) then
! compute difference between 'advected' pressure sp and read pressure
! sp1, accounting for reduce grid if any
call calc_pdiff( n, pdiffmax_l )
! compute maximum over all pe's
call Par_Reduce( pdiffmax_l, 'max', pdiffmax, status, all=.true. )
IF_NOTOK_RETURN(status=1)
! check ...
if ( pdiffmax > pdiffmax_treshold ) then
write (gol,'("difference between advected and read-in pressure exceeds treshold :")'); call goErr
write (gol,'(" max diff. : ",es9.2," [Pa]")') pdiffmax; call goErr
write (gol,'(" treshold : ",es9.2," [Pa]")') pdiffmax_treshold; call goErr
write (gol,'("pressure arrays saved to local `pressure.hdf`")'); call goErr
#ifdef with_hdf4
if (isRoot) then
io = sfStart( 'pressure.hdf', DFACC_CREATE )
if ( io > 0 ) then
call io_write2d_32d( io, im(n)+4, 'LON', jm(n)+4, 'LAT', sp1_dat(n)%data(:,:,1), 'p' , idate )
call io_write2d_32d( io, im(n)+4, 'LON', jm(n)+4, 'LAT', sp_dat(n)%data(:,:,1), 'pold', idate )
status = sfend(io)
else
write (gol,'("writing pressures")'); call goErr
end if
end if ! root
#endif
TRACEBACK; status=1; return
end if ! max diff
end if ! no newsrun
! ok
status = 0
if (okdebug) call goLabel()
END SUBROUTINE METEO_CHECKPRESSURE
!EOC
! **************************************************************
! ***
! *** vertical velocity
! ***
! **************************************************************
subroutine Compute_Omega( omega, lli, mfw, status )
use binas, only : grav
use grid, only : TllGridInfo, AreaOper
use meteodata, only : TMeteoData
use tmm, only : SetHistory, AddHistory
! --- in/out ----------------------------------
type(TMeteoData), intent(inout) :: omega ! Pa/s downward
type(TllGridInfo), intent(in) :: lli
type(TMeteoData), intent(in) :: mfw ! kg/s upward
integer, intent(out) :: status
! --- const -----------------------------------
character(len=*), parameter :: rname = mname//'/Compute_Omega'
! --- local ----------------------------------
integer :: l
! --- begin ----------------------------------
! not in use ?
if ( .not. omega%used ) return
! leave if not in use:
if ( .not. mfw%used ) then
write (gol,'("omega (Pa/s) requires mfw (kg/s)")'); call goErr
call goErr; status=1; return
end if
if (okdebug) call goLabel(rname)
! Pa/s = kg/s / m2 * g
! init with mass flux; revert sign from upward to downard, divide by
! gravity accelaration
omega%data = - mfw%data * grav ! Pa/s m2
! loop over levels and divide by cell area (m2)
do l = 1, size(omega%data,3)
call AreaOper( lli, omega%data(:,:,l), '/', 'm2', status )
IF_NOTOK_RETURN(status=1)
end do
! info ..
!call SetHistory( omega%tmi, mfw%tmi, status )
!call AddHistory( omega%tmi, 'convert to Pa/s', status )
! ok
status = 0
if (okdebug) call goLabel()
end subroutine Compute_Omega
! **************************************************************
! ***
! *** Specific SETUP routine for CONVECTIVE FLUXES
! ***
! **************************************************************
subroutine Setup_Convec_serial_io( region, entu, entd, detu, detd, omega, gph, &
tr, lli, levi, status)
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_Convec, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: entu, entd, detu, detd
type(TMeteoData), intent(in) :: omega, gph
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
type(TLevelInfo), intent(in) :: levi
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Setup_Convec_serial_io'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
! to differentiate b/w local and global data set
real, pointer, dimension(:,:,:) :: L_entu, L_entd, L_detu, L_detd
real, pointer :: L_omega(:,:,:), L_gph(:,:,:)
integer, dimension(2) :: is, js, ls, auxls
integer :: halo
! --- begin -----------------------------
! leave if not in use:
if ( (.not. all((/entu%used,entd%used,detu%used,detd%used/)) ) &
.and. any((/entu%used,entd%used,detu%used,detd%used/)) ) then
write (gol,'("either none or all of entu/entd/detu/detd should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. entu%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! gph is required as input:
if (.not. gph%used) then
write (gol,'("gph should be in use to compute convective stuff from EC convective fluxes")'); call goErr
call goErr; status=1; return
end if
! omega is required as input:
if ( .not. omega%used ) then
write (gol,'("omega should be in use to compute convective stuff")'); call goErr
call goErr; status=1; return
end if
! not changed by default
entu%changed = .false.
entd%changed = .false.
detu%changed = .false.
detd%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from entu only)
call SetupSetup( entu, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Need whole region for I/O on root. Dummy else.
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = entu%ls
auxls = gph%ls
IF (isRoot) THEN
! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo=0
allocate( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate(L_gph (im(region), jm(region), auxls(1):auxls(2)) )
allocate(L_omega(im(region), jm(region), auxls(1):auxls(2)) )
ELSE
allocate( L_entu(1,1,1), L_entd(1,1,1), L_detu(1,1,1), L_detd(1,1,1))
allocate(L_gph (1,1,1))
allocate(L_omega(1,1,1))
END IF
CALL GATHER( dgrid(region), gph%data, L_gph, gph%halo, status)
IF_NOTOK_RETURN(status=1)
CALL GATHER( dgrid(region), omega%data, L_omega, omega%halo, status)
IF_NOTOK_RETURN(status=1)
! Read/write on root
IOroot : if (isRoot) then
!AJS ! safety check ...
!AJS if ( data1_t2 /= data1_t1 ) then
!AJS !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
!AJS !call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
!AJS !call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
!AJS !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
!AJS !call goErr; status=1; return
!AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!AJS end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data
call Read_Convec( tmmd, entu%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
L_omega, omega%tmi, &
L_gph, gph%tmi, &
tmp_sp, &
L_entu, entu%tmi1, L_entd, entd%tmi1, &
L_detu, detu%tmi1, L_detd, detd%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( entu%putout ) then
call WriteField( tmmd, entu%destkey, &
entu%tmi1, 'sp', trim(entu%name), trim(entu%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_entu, status )
IF_NOTOK_RETURN(status=1)
end if
if ( entd%putout ) then
call WriteField( tmmd, entd%destkey, &
entd%tmi1, 'sp', trim(entd%name), trim(entd%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_entd, status )
IF_NOTOK_RETURN(status=1)
end if
if ( detu%putout ) then
call WriteField( tmmd, detu%destkey, &
detu%tmi1, 'sp', trim(detu%name), trim(detu%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_detu, status )
IF_NOTOK_RETURN(status=1)
end if
if ( detd%putout ) then
call WriteField( tmmd, detd%destkey, &
detd%tmi1, 'sp', trim(detd%name), trim(detd%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_detd, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot
! Scatter & clean up
CALL SCATTER( dgrid(region), entu%data1, L_entu, entu%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), entd%data1, L_entd, entd%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), detu%data1, L_detu, detu%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), detd%data1, L_detd, detd%halo, status)
IF_NOTOK_RETURN(status=1)
deallocate(L_entu, L_entd, L_detu, L_detd, L_gph, L_omega)
! data array is filled now:
entu%filled1 = .true.
entu%tr1(1) = data1_t1
entu%tr1(2) = data1_t2
entu%changed = .true.
entd%filled1 = .true.
entd%tr1(1) = data1_t1
entd%tr1(2) = data1_t2
entd%changed = .true.
detu%filled1 = .true.
detu%tr1(1) = data1_t1
detu%tr1(2) = data1_t2
detu%changed = .true.
detd%filled1 = .true.
detd%tr1(1) = data1_t1
detd%tr1(2) = data1_t2
detd%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
entu%data1 = entu%data2
entd%data1 = entd%data2
detu%data1 = detu%data2
detd%data1 = detd%data2
! data array is filled now:
entu%filled1 = .true.
entu%tr1(1) = data1_t1
entu%tr1(2) = data1_t2
entu%changed = .true.
entd%filled1 = .true.
entd%tr1(1) = data1_t1
entd%tr1(2) = data1_t2
entd%changed = .true.
detu%filled1 = .true.
detu%tr1(1) = data1_t1
detu%tr1(2) = data1_t2
detu%changed = .true.
detd%filled1 = .true.
detd%tr1(1) = data1_t1
detd%tr1(2) = data1_t2
detd%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Need whole region for I/O on root. Dummy else
is = (/1,im(1)/)
js = (/1,jm(1)/)
ls = entu%ls
auxls = gph%ls
IF (isRoot) THEN
! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo
ALLOCATE( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE(L_gph (im(region),jm(region),auxls(1):auxls(2)))
ALLOCATE(L_omega (im(region),jm(region),auxls(1):auxls(2)))
ELSE
ALLOCATE( L_entu(1,1,1), L_entd(1,1,1), L_detu(1,1,1), L_detd(1,1,1))
ALLOCATE( L_gph(1,1,1), L_omega(1,1,1) )
END IF
CALL GATHER( dgrid(region), gph%data, L_gph, gph%halo, status)
IF_NOTOK_RETURN(status=1)
CALL GATHER( dgrid(region), omega%data, L_omega, omega%halo, status)
IF_NOTOK_RETURN(status=1)
! Read/write on root
IOroot2 : if (isRoot) then
!AJS ! safety check ...
!AJS if ( data2_t2 /= data2_t1 ) then
!AJS !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
!AJS !call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
!AJS !call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
!AJS !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
!AJS !call goErr; status=1; return
!AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!AJS end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data2
call Read_Convec( tmmd, entu%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
L_omega, omega%tmi, &
L_gph, gph%tmi, &
tmp_sp, &
L_entu, entu%tmi2, L_entd, entd%tmi2, &
L_detu, detu%tmi2, L_detd, detd%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( entu%putout ) then
call WriteField( tmmd, entu%destkey, &
entu%tmi2, 'sp', trim(entu%name), trim(entu%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_entu, status )
IF_NOTOK_RETURN(status=1)
end if
if ( entd%putout ) then
call WriteField( tmmd, entd%destkey, &
entd%tmi2, 'sp', trim(entd%name), trim(entd%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_entd, status )
IF_NOTOK_RETURN(status=1)
end if
if ( detu%putout ) then
call WriteField( tmmd, detu%destkey, &
detu%tmi2, 'sp', trim(detu%name), trim(detu%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_detu, status )
IF_NOTOK_RETURN(status=1)
end if
if ( detd%putout ) then
call WriteField( tmmd, detd%destkey, &
detd%tmi2, 'sp', trim(detd%name), trim(detd%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, '*', &
tmp_sp, L_detd, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot2
CALL SCATTER( dgrid(region), entu%data2, L_entu, entu%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), entd%data2, L_entd, entd%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), detu%data2, L_detu, detu%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), detd%data2, L_detd, detd%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE( L_entu, L_entd, L_detu, L_detd, L_gph, L_omega )
! data2 array is filled now:
entu%filled2 = .true.
entu%tr2(1) = data2_t1
entu%tr2(2) = data2_t2
entd%filled2 = .true.
entd%tr2(1) = data2_t1
entd%tr2(2) = data2_t2
detu%filled2 = .true.
detu%tr2(1) = data2_t1
detu%tr2(2) = data2_t2
detd%filled2 = .true.
detd%tr2(1) = data2_t1
detd%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
entu%data2 = entu%data1
entd%data2 = entd%data1
detu%data2 = detu%data1
detd%data2 = detd%data1
! data2 array is filled now:
entu%filled2 = .true.
entu%tr2(1) = data2_t1
entu%tr2(2) = data2_t2
entd%filled2 = .true.
entd%tr2(1) = data2_t1
entd%tr2(2) = data2_t2
detu%filled2 = .true.
detu%tr2(1) = data2_t1
detu%tr2(2) = data2_t2
detd%filled2 = .true.
detd%tr2(1) = data2_t1
detd%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( entu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( entd, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( detu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( detd, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_CONVEC_SERIAL_IO
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_CONVEC_PARALLEL_IO
!
! !DESCRIPTION: same as setup_convec_serial_io but with parallel i/o
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_CONVEC_PARALLEL_IO( region, entu, entd, detu, detd, omega, gph, &
tr, levi, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_Convec, WriteField
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: entu, entd, detu, detd
type(TMeteoData), intent(in) :: omega, gph
type(TDate), intent(in) :: tr(2)
type(TLevelInfo), intent(in) :: levi
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 24 Oct 2013 - Ph. Le Sager - v0
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/SETUP_CONVEC_PARALLEL_IO'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
! to differentiate b/w local and global data set
real, pointer, dimension(:,:,:) :: L_entu, L_entd, L_detu, L_detd
real, pointer :: L_omega(:,:,:), L_gph(:,:,:)
integer, dimension(2) :: is, js, ls, auxls
integer :: halo
integer :: i1, i2, j1, j2
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( (.not. all((/entu%used,entd%used,detu%used,detd%used/)) ) &
.and. any((/entu%used,entd%used,detu%used,detd%used/)) ) then
write (gol,'("either none or all of entu/entd/detu/detd should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. entu%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! gph is required as input:
if ( .not. gph%used ) then
write (gol,'("gph should be in use to compute convective stuff from EC convective fluxes")'); call goErr
call goErr; status=1; return
end if
! omega is required as input:
if ( .not. omega%used ) then
write (gol,'("omega should be in use to compute convective stuff")'); call goErr
call goErr; status=1; return
end if
! not changed by default
entu%changed = .false.
entd%changed = .false.
detu%changed = .false.
detd%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from entu only)
call SetupSetup( entu, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! work arrays
if (data1_read .or. data2_read) then
CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
is = (/i1,i2/)
js = (/j1,j2/)
ls = entu%ls
auxls = gph%ls
! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo=0
allocate( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
L_gph => gph%data
L_omega => omega%data
end if
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
!AJS if ( data1_t2 /= data1_t1 ) then
!AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!AJS end if
call Read_Convec( tmmd, entu%sourcekey, &
data1_tref, data1_t1, data1_t2, lli(region), levi, &
L_omega, omega%tmi, &
L_gph, gph%tmi, &
tmp_sp, &
L_entu, entu%tmi1, L_entd, entd%tmi1, &
L_detu, detu%tmi1, L_detd, detd%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( entu%putout ) then
! call WriteField( tmmd, entu%destkey, &
! entu%tmi1, 'sp', trim(entu%name), trim(entu%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_entu, status )
! IF_NOTOK_RETURN(status=1)
end if
if ( entd%putout ) then
! call WriteField( tmmd, entd%destkey, &
! entd%tmi1, 'sp', trim(entd%name), trim(entd%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_entd, status )
! IF_NOTOK_RETURN(status=1)
end if
if ( detu%putout ) then
! call WriteField( tmmd, detu%destkey, &
! detu%tmi1, 'sp', trim(detu%name), trim(detu%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_detu, status )
! IF_NOTOK_RETURN(status=1)
end if
if ( detd%putout ) then
! call WriteField( tmmd, detd%destkey, &
! detd%tmi1, 'sp', trim(detd%name), trim(detd%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_detd, status )
! IF_NOTOK_RETURN(status=1)
end if
entu%data1(i1:i2,j1:j2,:) = L_entu
entd%data1(i1:i2,j1:j2,:) = L_entd
detu%data1(i1:i2,j1:j2,:) = L_detu
detd%data1(i1:i2,j1:j2,:) = L_detd
! data array is filled now:
entu%filled1 = .true.
entu%tr1(1) = data1_t1
entu%tr1(2) = data1_t2
entu%changed = .true.
entd%filled1 = .true.
entd%tr1(1) = data1_t1
entd%tr1(2) = data1_t2
entd%changed = .true.
detu%filled1 = .true.
detu%tr1(1) = data1_t1
detu%tr1(2) = data1_t2
detu%changed = .true.
detd%filled1 = .true.
detd%tr1(1) = data1_t1
detd%tr1(2) = data1_t2
detd%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
entu%data1 = entu%data2
entd%data1 = entd%data2
detu%data1 = detu%data2
detd%data1 = detd%data2
! data array is filled now:
entu%filled1 = .true.
entu%tr1(1) = data1_t1
entu%tr1(2) = data1_t2
entu%changed = .true.
entd%filled1 = .true.
entd%tr1(1) = data1_t1
entd%tr1(2) = data1_t2
entd%changed = .true.
detu%filled1 = .true.
detu%tr1(1) = data1_t1
detu%tr1(2) = data1_t2
detu%changed = .true.
detd%filled1 = .true.
detd%tr1(1) = data1_t1
detd%tr1(2) = data1_t2
detd%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
!AJS if ( data2_t2 /= data2_t1 ) then
!AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!AJS end if
call Read_Convec( tmmd, entu%sourcekey, &
data2_tref, data2_t1, data2_t2, lli(region), levi, &
L_omega, omega%tmi, &
L_gph, gph%tmi, &
tmp_sp, &
L_entu, entu%tmi2, L_entd, entd%tmi2, &
L_detu, detu%tmi2, L_detd, detd%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
! if ( entu%putout ) then
! call WriteField( tmmd, entu%destkey, &
! entu%tmi2, 'sp', trim(entu%name), trim(entu%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_entu, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( entd%putout ) then
! call WriteField( tmmd, entd%destkey, &
! entd%tmi2, 'sp', trim(entd%name), trim(entd%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_entd, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( detu%putout ) then
! call WriteField( tmmd, detu%destkey, &
! detu%tmi2, 'sp', trim(detu%name), trim(detu%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_detu, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( detd%putout ) then
! call WriteField( tmmd, detd%destkey, &
! detd%tmi2, 'sp', trim(detd%name), trim(detd%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, '*', &
! tmp_sp, L_detd, status )
! IF_NOTOK_RETURN(status=1)
! end if
entu%data2(i1:i2,j1:j2,:) = L_entu
entd%data2(i1:i2,j1:j2,:) = L_entd
detu%data2(i1:i2,j1:j2,:) = L_detu
detd%data2(i1:i2,j1:j2,:) = L_detd
! data2 array is filled now:
entu%filled2 = .true.
entu%tr2(1) = data2_t1
entu%tr2(2) = data2_t2
entd%filled2 = .true.
entd%tr2(1) = data2_t1
entd%tr2(2) = data2_t2
detu%filled2 = .true.
detu%tr2(1) = data2_t1
detu%tr2(2) = data2_t2
detd%filled2 = .true.
detd%tr2(1) = data2_t1
detd%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
entu%data2 = entu%data1
entd%data2 = entd%data1
detu%data2 = detu%data1
detd%data2 = detd%data1
! data2 array is filled now:
entu%filled2 = .true.
entu%tr2(1) = data2_t1
entu%tr2(2) = data2_t2
entd%filled2 = .true.
entd%tr2(1) = data2_t1
entd%tr2(2) = data2_t2
detu%filled2 = .true.
detu%tr2(1) = data2_t1
detu%tr2(2) = data2_t2
detd%filled2 = .true.
detd%tr2(1) = data2_t1
detd%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( entu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( entd, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( detu, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( detd, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
if (data1_read .or. data2_read) then
deallocate(L_entu, L_entd, L_detu, L_detd)
deallocate( tmp_sp )
nullify(L_gph, L_omega)
end if
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_CONVEC_PARALLEL_IO
!EOC
! **************************************************************
! ***
! *** diffusive fluxes
! ***
! **************************************************************
SUBROUTINE SETUP_DIFFUS_SERIAL_IO( region, Kzz, tr, lli, levi, status )
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_Diffus, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: Kzz
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
type(TLevelInfo), intent(in) :: levi
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/SETUP_DIFFUS_SERIAL_IO'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer, dimension(:,:,:) :: L_kzz ! work arrays (data)
integer :: is(2), js(2), ls(2) ! work arrays (bounds)
! --- begin -----------------------------
! not in use ?
if ( .not. Kzz%used ) then
status=0; return
end if
if (okdebug) then
call goLabel(rname)
write(gol,'(" ",a,": ",a,l2)') rname, "Diffus", Kzz%used; call goPr
end if
! not changed by default
Kzz%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read
call SetupSetup( Kzz, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Allocate global arrays for I/O
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = Kzz%ls
IF (isRoot) THEN
ALLOCATE( L_kzz( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ELSE
ALLOCATE( L_kzz(1,1,1) )
END IF
IOroot : if (isRoot) then
!! safety check ...
!if ( data1_t2 /= data1_t1 ) then
! !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
! !call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
! !call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
! !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
! !call goErr; status=1; return
! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data:
call Read_Diffus( tmmd, Kzz%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_sp, &
L_Kzz, Kzz%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( Kzz%putout ) then
call WriteField( tmmd, Kzz%destkey, &
Kzz%tmi1, 'sp', trim(Kzz%name), trim(Kzz%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'w', &
tmp_sp, L_Kzz, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot
! Wrap up
CALL SCATTER( dgrid(region), Kzz%data1, L_kzz, kzz%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(L_kzz)
! data array is filled now:
Kzz%filled1 = .true.
Kzz%tr1(1) = data1_t1
Kzz%tr1(2) = data1_t2
Kzz%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
Kzz%data1 = Kzz%data2
! data array is filled now:
Kzz%filled1 = .true.
Kzz%tr1(1) = data1_t1
Kzz%tr1(2) = data1_t2
Kzz%changed = .true.
end if
! secondary field ?
if ( data2_read ) then
! Allocate global arrays for I/O
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = kzz%ls
IF (isRoot) THEN
ALLOCATE( L_kzz( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ELSE
ALLOCATE( L_kzz(1,1,1) )
END IF
! Read/write
IOroot2 : IF (isRoot) THEN
!! safety check ...
!if ( data2_t2 /= data2_t1 ) then
! !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
! !call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
! !call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
! !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
! !call goErr; status=1; return
! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data2:
call Read_Diffus( tmmd, Kzz%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
tmp_sp, &
L_Kzz, Kzz%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( Kzz%putout ) then
call WriteField( tmmd, Kzz%destkey, &
Kzz%tmi2, 'sp', trim(Kzz%name), trim(Kzz%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'w', &
tmp_sp, L_Kzz, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot2
CALL SCATTER( dgrid(region), Kzz%data2, L_kzz, kzz%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(L_kzz)
! data2 array is filled now:
Kzz%filled2 = .true.
Kzz%tr2(1) = data2_t1
Kzz%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
Kzz%data2 = Kzz%data1
! data2 array is filled now:
Kzz%filled2 = .true.
Kzz%tr2(1) = data2_t1
Kzz%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( Kzz, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_DIFFUS_SERIAL_IO
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_DIFFUS_PARALLEL_IO
!
! !DESCRIPTION:
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_DIFFUS_PARALLEL_IO( region, Kzz, tr, levi, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_Diffus, WriteField
!
! !INPUT/OUTPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: Kzz
type(TDate), intent(in) :: tr(2)
type(TLevelInfo), intent(in) :: levi
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 3 Dec 2013 - Ph. Le Sager -
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/SETUP_DIFFUS_PARALLEL_IO'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:)
real, pointer, dimension(:,:,:) :: L_kzz ! work arrays (data)
integer :: is(2), js(2), ls(2) ! work arrays (bounds)
integer :: i1, i2, j1, j2
! --- begin -----------------------------
! not in use ?
if ( .not. Kzz%used ) then
status=0; return
end if
if (okdebug) then
call goLabel(rname)
write(gol,'(" ",a,": ",a,l2)') rname, "Diffus", Kzz%used; call goPr
end if
! not changed by default
Kzz%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read
call SetupSetup( Kzz, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! work arrays
IF (data1_read .OR. data2_read) THEN
CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
is = (/i1,i2/)
js = (/j1,j2/)
ls = kzz%ls
ALLOCATE( L_kzz( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( tmp_sp(is(1):is(2), js(1):js(2) ) )
ENDIF
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
!! safety check ...
!if ( data1_t2 /= data1_t1 ) then
! !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
! !call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
! !call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
! !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
! !call goErr; status=1; return
! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!end if
call Read_Diffus( tmmd, Kzz%sourcekey, &
data1_tref, data1_t1, data1_t2, lli(region), levi, &
tmp_sp, &
L_Kzz, Kzz%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
!TODO ! write meteofiles ?
!TODO if ( Kzz%putout ) then
!TODO call WriteField( tmmd, Kzz%destkey, &
!TODO Kzz%tmi1, 'sp', trim(Kzz%name), trim(Kzz%unit), &
!TODO data1_tref, data1_t1, data1_t2, &
!TODO lli, 'n', levi, 'w', &
!TODO tmp_sp, L_Kzz, status )
!TODO IF_NOTOK_RETURN(status=1)
!TODO end if
!TODO
kzz%data1(i1:i2,j1:j2,:) = L_Kzz
! data array is filled now:
Kzz%filled1 = .true.
Kzz%tr1(1) = data1_t1
Kzz%tr1(2) = data1_t2
Kzz%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
Kzz%data1 = Kzz%data2
! data array is filled now:
Kzz%filled1 = .true.
Kzz%tr1(1) = data1_t1
Kzz%tr1(2) = data1_t2
Kzz%changed = .true.
end if
! secondary field ?
if ( data2_read ) then
!! safety check ...
!if ( data2_t2 /= data2_t1 ) then
! !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
! !call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
! !call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
! !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
! !call goErr; status=1; return
! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
!end if
! fill data2:
call Read_Diffus( tmmd, Kzz%sourcekey, &
data2_tref, data2_t1, data2_t2, lli(region), levi, &
tmp_sp, &
L_Kzz, Kzz%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
!TODO ! write meteofiles ?
!TODO if ( Kzz%putout ) then
!TODO call WriteField( tmmd, Kzz%destkey, &
!TODO Kzz%tmi2, 'sp', trim(Kzz%name), trim(Kzz%unit), &
!TODO data2_tref, data2_t1, data2_t2, &
!TODO lli, 'n', levi, 'w', &
!TODO tmp_sp, L_Kzz, status )
!TODO IF_NOTOK_RETURN(status=1)
!TODO end if
!TODO
Kzz%data1(i1:i2,j1:j2,:) = L_Kzz
! data2 array is filled now:
Kzz%filled2 = .true.
Kzz%tr2(1) = data2_t1
Kzz%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
Kzz%data2 = Kzz%data1
! data2 array is filled now:
Kzz%filled2 = .true.
Kzz%tr2(1) = data2_t1
Kzz%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( Kzz, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
if (data1_read .or. data2_read) then
deallocate( tmp_sp, L_Kzz)
end if
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_DIFFUS_PARALLEL_IO
!EOC
! **************************************************************
! ***
! *** Specific SETUP routine for CLOUD COVER
! ***
! **************************************************************
SUBROUTINE SETUP_CLOUDCOVERS_serial_io( region, cc, cco, ccu, tr, lli, levi, status )
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_CloudCovers, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use dims, only : im, jm
! --- in/out ----------------------------------
integer, intent(in) :: region ! region number
type(TMeteoData), intent(inout) :: cc, cco, ccu
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
type(TLevelInfo), intent(in) :: levi
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/SETUP_CLOUDCOVERS_serial_io'
! --- local ----------------------------------
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:) ! surface pressure
real, pointer, dimension(:,:,:) :: L_cc, L_cco, L_ccu ! work arrays (data)
integer :: is(2), js(2), ls(2) ! work arrays (bounds)
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( (.not. all((/cc%used,cco%used,ccu%used/)) ) .and. any((/cc%used,cco%used,ccu%used/)) ) then
write (gol,'("either none or all of cc/cco/ccu should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. cc%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
cc%changed = .false.
cco%changed = .false.
ccu%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from cc only)
call SetupSetup( cc, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! Allocate global arrays for I/O
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = cc%ls
IF (isRoot) THEN
ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ELSE
ALLOCATE( L_cc(1,1,1), L_cco(1,1,1), L_ccu(1,1,1) )
END IF
! Read/write on root
IOroot : if (isRoot) then
! safety check ...
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data:
call Read_CloudCovers( tmmd, cc%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_sp, &
L_cc, cc%tmi1, &
L_cco, cco%tmi1, &
L_ccu, ccu%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles
if ( cc%putout ) then
call WriteField( tmmd, cc%destkey, &
cc%tmi1, 'sp', trim(cc%name), trim(cc%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_cc, status )
IF_NOTOK_RETURN(status=1)
end if
if ( cco%putout ) then
call WriteField( tmmd, cco%destkey, &
cco%tmi1, 'sp', trim(cco%name), trim(cco%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_cco, status )
IF_NOTOK_RETURN(status=1)
end if
if ( ccu%putout ) then
call WriteField( tmmd, ccu%destkey, &
ccu%tmi1, 'sp', trim(ccu%name), trim(ccu%unit), &
data1_tref, data1_t1, data1_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_ccu, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot
! Wrap up
CALL SCATTER( dgrid(region), cc%data1, L_cc, cc%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), cco%data1, L_cco, cco%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), ccu%data1, L_ccu, ccu%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(L_cc, L_cco, L_ccu)
! data array is filled now:
cc%filled1 = .true.
cc%tr1(1) = data1_t1
cc%tr1(2) = data1_t2
cc%changed = .true.
cco%filled1 = .true.
cco%tr1(1) = data1_t1
cco%tr1(2) = data1_t2
cco%changed = .true.
ccu%filled1 = .true.
ccu%tr1(1) = data1_t1
ccu%tr1(2) = data1_t2
ccu%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
cc%data1 = cc%data2
cco%data1 = cco%data2
ccu%data1 = ccu%data2
! data array is filled now:
cc%filled1 = .true.
cc%tr1(1) = data1_t1
cc%tr1(2) = data1_t2
cc%changed = .true.
cco%filled1 = .true.
cco%tr1(1) = data1_t1
cco%tr1(2) = data1_t2
cco%changed = .true.
ccu%filled1 = .true.
ccu%tr1(1) = data1_t1
ccu%tr1(2) = data1_t2
ccu%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! Allocate global arrays for I/O
is = (/1,im(region)/)
js = (/1,jm(region)/)
ls = cc%ls
IF (isRoot) THEN
! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo
ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ELSE
ALLOCATE( L_cc(1,1,1), L_cco(1,1,1), L_ccu(1,1,1) )
END IF
! Read/write
IOroot2 : IF (isRoot) THEN
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
call goErr; status=1; return
end if
! surface pressure field:
allocate( tmp_sp(is(1):is(2),js(1):js(2)) )
! fill data2:
call Read_CloudCovers( tmmd, cc%sourcekey, data2_tref, &
data2_t1, data2_t2, lli, levi, &
tmp_sp, &
L_cc, cc%tmi2, &
L_cco, cco%tmi2, &
L_ccu, ccu%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
if ( cc%putout ) then
call WriteField( tmmd, cc%destkey, &
cc%tmi2, 'sp', trim( cc%name), trim( cc%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_cc, status )
IF_NOTOK_RETURN(status=1)
end if
if ( cco%putout ) then
call WriteField( tmmd, cco%destkey, &
cco%tmi2, 'sp', trim(cco%name), trim(cco%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_cco, status )
IF_NOTOK_RETURN(status=1)
end if
if ( ccu%putout ) then
call WriteField( tmmd, ccu%destkey, &
ccu%tmi2, 'sp', trim(ccu%name), trim(ccu%unit), &
data2_tref, data2_t1, data2_t2, &
lli, 'n', levi, 'n', &
tmp_sp, L_ccu, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if IOroot2
! Wrap up
CALL SCATTER( dgrid(region), cc%data2, L_cc, cc%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), cco%data2, L_cco, cco%halo, status)
IF_NOTOK_RETURN(status=1)
CALL SCATTER( dgrid(region), ccu%data2, L_ccu, ccu%halo, status)
IF_NOTOK_RETURN(status=1)
DEALLOCATE(L_cc, L_cco, L_ccu)
! data2 array is filled now:
cc%filled2 = .true.
cc%tr2(1) = data2_t1
cc%tr2(2) = data2_t2
cco%filled2 = .true.
cco%tr2(1) = data2_t1
cco%tr2(2) = data2_t2
ccu%filled2 = .true.
ccu%tr2(1) = data2_t1
ccu%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
cc%data2 = cc%data1
cco%data2 = cco%data1
ccu%data2 = ccu%data1
! data2 array is filled now:
cc%filled2 = .true.
cc%tr2(1) = data2_t1
cc%tr2(2) = data2_t2
cco%filled2 = .true.
cco%tr2(1) = data2_t1
cco%tr2(2) = data2_t2
ccu%filled2 = .true.
ccu%tr2(1) = data2_t1
ccu%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( cc, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( cco, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( ccu, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_CLOUDCOVERS_SERIAL_IO
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SETUP_CLOUDCOVERS_PARALLEL_IO
!
! !DESCRIPTION:
!\\
!\\
! !INTERFACE:
!
SUBROUTINE SETUP_CLOUDCOVERS_PARALLEL_IO( region, cc, cco, ccu, tr, levi, status )
!
! !USES:
!
use GO, only : TDate, wrtgol, operator(/=)
use Grid, only : TllGridInfo, TLevelInfo
use TMM, only : TMeteoInfo, Read_CloudCovers, WriteField
use dims, only : im, jm
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region ! region number
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: cc, cco, ccu
type(TDate), intent(in) :: tr(2)
type(TLevelInfo), intent(in) :: levi
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 24 Oct 2013 - Ph. Le Sager - v0
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/SETUP_CLOUDCOVERS_PARALLEL_IO'
logical :: data1_read, data1_copy
type(TDate) :: data1_tref, data1_t1, data1_t2
logical :: data2_read, data2_copy
type(TDate) :: data2_tref, data2_t1, data2_t2
real, allocatable :: tmp_sp(:,:) ! surface pressure
real, pointer, dimension(:,:,:) :: L_cc, L_cco, L_ccu ! work arrays (data)
integer :: is(2), js(2), ls(2) ! work arrays (bounds)
integer :: i1, i2, j1, j2
! --- begin -----------------------------
if (okdebug) call goLabel(rname)
! leave if not in use:
if ( (.not. all((/cc%used,cco%used,ccu%used/)) ) .and. any((/cc%used,cco%used,ccu%used/)) ) then
write (gol,'("either none or all of cc/cco/ccu should be in use")'); call goErr
call goErr; status=1; return
end if
if ( .not. cc%used ) then
if (okdebug) call goLabel()
status=0; return
end if
! not changed by default
cc%changed = .false.
cco%changed = .false.
ccu%changed = .false.
!------------------
! time stuff
!------------------
! get time interval of met field and check if data from start and/or end
! of interval must be read (sufficient to setup from cc only)
call SetupSetup( cc, tr, &
data1_read, data1_copy, data1_tref, data1_t1, data1_t2, &
data2_read, data2_copy, data2_tref, data2_t1, data2_t2, &
status )
IF_NOTOK_RETURN(status=1)
! work arrays
IF (data1_read .OR. data2_read) THEN
CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
is = (/i1,i2/)
js = (/j1,j2/)
ls = cc%ls
ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) )
ALLOCATE( tmp_sp(is(1):is(2),js(1):js(2)) )
ENDIF
!--------------------------
! read/write primary field
!--------------------------
if ( data1_read ) then
! safety check
if ( data1_t2 /= data1_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr
call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr
call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr
write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr
TRACEBACK; status=1; return
end if
call Read_CloudCovers( tmmd, cc%sourcekey, &
data1_tref, data1_t1, data1_t2, lli(region), levi, &
tmp_sp, &
L_cc, cc%tmi1, &
L_cco, cco%tmi1, &
L_ccu, ccu%tmi1, &
status )
IF_NOTOK_RETURN(status=1)
! ! write meteofiles
! if ( cc%putout ) then
! call WriteField( tmmd, cc%destkey, &
! cc%tmi1, 'sp', trim(cc%name), trim(cc%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_cc, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( cco%putout ) then
! call WriteField( tmmd, cco%destkey, &
! cco%tmi1, 'sp', trim(cco%name), trim(cco%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_cco, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( ccu%putout ) then
! call WriteField( tmmd, ccu%destkey, &
! ccu%tmi1, 'sp', trim(ccu%name), trim(ccu%unit), &
! data1_tref, data1_t1, data1_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_ccu, status )
! IF_NOTOK_RETURN(status=1)
! end if
cc%data1(i1:i2,j1:j2,:) = L_cc
cco%data1(i1:i2,j1:j2,:) = L_cco
ccu%data1(i1:i2,j1:j2,:) = L_ccu
! data array is filled now:
cc%filled1 = .true.
cc%tr1(1) = data1_t1
cc%tr1(2) = data1_t2
cc%changed = .true.
cco%filled1 = .true.
cco%tr1(1) = data1_t1
cco%tr1(2) = data1_t2
cco%changed = .true.
ccu%filled1 = .true.
ccu%tr1(1) = data1_t1
ccu%tr1(2) = data1_t2
ccu%changed = .true.
else if ( data1_copy ) then
! copy data from secondary array:
cc%data1 = cc%data2
cco%data1 = cco%data2
ccu%data1 = ccu%data2
! data array is filled now:
cc%filled1 = .true.
cc%tr1(1) = data1_t1
cc%tr1(2) = data1_t2
cc%changed = .true.
cco%filled1 = .true.
cco%tr1(1) = data1_t1
cco%tr1(2) = data1_t2
cco%changed = .true.
ccu%filled1 = .true.
ccu%tr1(1) = data1_t1
ccu%tr1(2) = data1_t2
ccu%changed = .true.
end if
!--------------------------
! read/write secondary field
!--------------------------
if ( data2_read ) then
! safety check ...
if ( data2_t2 /= data2_t1 ) then
write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr
call wrtgol( ' data2_t1 : ', data2_t1 ) ; call goErr
call wrtgol( ' data2_t2 : ', data2_t2 ) ; call goErr
write (gol,'("please deceide what to do with surface pressures ... ")') ; call goErr
call goErr; status=1; return
end if
call Read_CloudCovers( tmmd, cc%sourcekey, data2_tref, &
data2_t1, data2_t2, lli(region), levi, &
tmp_sp, &
L_cc, cc%tmi2, &
L_cco, cco%tmi2, &
L_ccu, ccu%tmi2, &
status )
IF_NOTOK_RETURN(status=1)
! ! write meteofiles ?
! if ( cc%putout ) then
! call WriteField( tmmd, cc%destkey, &
! cc%tmi2, 'sp', trim( cc%name), trim( cc%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_cc, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( cco%putout ) then
! call WriteField( tmmd, cco%destkey, &
! cco%tmi2, 'sp', trim(cco%name), trim(cco%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_cco, status )
! IF_NOTOK_RETURN(status=1)
! end if
! if ( ccu%putout ) then
! call WriteField( tmmd, ccu%destkey, &
! ccu%tmi2, 'sp', trim(ccu%name), trim(ccu%unit), &
! data2_tref, data2_t1, data2_t2, &
! lli, 'n', levi, 'n', &
! tmp_sp, L_ccu, status )
! IF_NOTOK_RETURN(status=1)
! end if
cc%data2(i1:i2,j1:j2,:) = L_cc
cco%data2(i1:i2,j1:j2,:) = L_cco
ccu%data2(i1:i2,j1:j2,:) = L_ccu
! data2 array is filled now:
cc%filled2 = .true.
cc%tr2(1) = data2_t1
cc%tr2(2) = data2_t2
cco%filled2 = .true.
cco%tr2(1) = data2_t1
cco%tr2(2) = data2_t2
ccu%filled2 = .true.
ccu%tr2(1) = data2_t1
ccu%tr2(2) = data2_t2
else if ( data2_copy ) then
! copy data2 from primary array:
cc%data2 = cc%data1
cco%data2 = cco%data1
ccu%data2 = ccu%data1
! data2 array is filled now:
cc%filled2 = .true.
cc%tr2(1) = data2_t1
cc%tr2(2) = data2_t2
cco%filled2 = .true.
cco%tr2(1) = data2_t1
cco%tr2(2) = data2_t2
ccu%filled2 = .true.
ccu%tr2(1) = data2_t1
ccu%tr2(2) = data2_t2
end if
!------------------
! time interpolation
!------------------
call TimeInterpolation( cc, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( cco, tr, status )
IF_NOTOK_RETURN(status=1)
call TimeInterpolation( ccu, tr, status )
IF_NOTOK_RETURN(status=1)
!------------------
! done
!------------------
if (data1_read .or. data2_read) then
deallocate( tmp_sp )
deallocate(L_cc, L_cco, L_ccu)
end if
status = 0
if (okdebug) call goLabel()
END SUBROUTINE SETUP_CLOUDCOVERS_PARALLEL_IO
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: PRESSURE_TO_MASS
!
! !DESCRIPTION: Get Air Mass: from surface pressure (sp), get pressure at
! box boundaries (so-called half-levels, phlb), and then air
! mass (m_dat).
!\\
!\\
! !INTERFACE:
!
SUBROUTINE PRESSURE_TO_MASS( region, status )
!
! !USES:
!
use Binas, only : grav
use Grid, only : HPressure
!use Grid, only : FillMass
use Grid, only : AreaOper
use dims, only : im, jm, lm
use dims, only : xcyc
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
! !REMARKS: (old remark: "assume that halo cells in sp have been filled
! correctly..." still valid???)
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Pressure_to_Mass'
integer :: l, i0, i1, j0, j1, lmr
! --- begin ----------------------------------
! Local grid size
i0 = sp_dat(region)%is(1)
i1 = sp_dat(region)%is(2)
j0 = sp_dat(region)%js(1)
j1 = sp_dat(region)%js(2)
lmr = lm(region)
! Fill pressure boundaries (Pa)
if ( phlb_dat(region)%used ) then
call HPressure( levi, sp_dat(region)%data(i0:i1, j0:j1, 1), &
phlb_dat(region)%data(i0:i1, j0:j1, :), status )
IF_NOTOK_RETURN(status=0)
end if
! Fill air mass (kg)
if ( m_dat(region)%used ) then
!call FillMass( m_dat(region)%data(1:imr,1:jmr,:), lli(region), levi, &
! sp_dat(region)%data(1:imr,1:jmr,1), status )
!IF_NOTOK_RETURN(status=0)
! Pressure difference between top and bottom of layer
do l = 1, lmr
m_dat(region)%data(:,:,l) = phlb_dat(region)%data(:,:,l) - phlb_dat(region)%data(:,:,l+1) ! Pa
end do
! Convert to kg/m2
m_dat(region)%data = m_dat(region)%data / grav ! Pa/g = kg/m2
! Convert to kg
call AreaOper( lli(region), m_dat(region)%data(i0:i1, j0:j1, :), '*', 'm2', status ) ! kg
IF_NOTOK_RETURN(status=0)
end if
! ok
status = 0
END SUBROUTINE PRESSURE_TO_MASS
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: MASS_TO_PRESSURE
!
! !DESCRIPTION: get 3D and surface (spm) pressures from 3D Air Mass.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE MASS_TO_PRESSURE( region, status )
!
! !USES:
!
use Binas, only : grav
use Grid, only : AreaOper
use dims, only : lm
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/Mass_to_Pressure'
integer :: l, i0, i1, j0, j1, lmr
! --- begin ----------------------------------
! Local grid size
i0 = sp_dat(region)%is(1)
i1 = sp_dat(region)%is(2)
j0 = sp_dat(region)%js(1)
j1 = sp_dat(region)%js(2)
lmr = lm(region)
! Fill pressure at half level boundaries:
! o zero in space:
phlb_dat(region)%data(:,:,lmr+1) = 0.0 ! kg m/s2 = Pa m2
! o add for each level pressure gradient:
do l = lmr, 1, -1
phlb_dat(region)%data(i0:i1, j0:j1, l) = phlb_dat(region)%data(i0:i1, j0:j1, l+1) &
+ m_dat(region)%data(i0:i1, j0:j1, l ) * grav ! kg m/s2 = Pa m2
end do
! Divide by grid cell area
call AreaOper( lli(region), phlb_dat(region)%data(i0:i1, j0:j1, :), '/', 'm2', status ) ! Pa
IF_NOTOK_RETURN(status=0)
! copy surface pressure
spm_dat(region)%data(i0:i1, j0:j1, 1) = phlb_dat(region)%data(i0:i1, j0:j1, 1) ! Pa
! ok
status = 0
END SUBROUTINE MASS_TO_PRESSURE
!EOC
!--------------------------------------------------------------------------
! TM5 !
!--------------------------------------------------------------------------
!BOP
!
! !IROUTINE: COMPUTE_GPH
!
! !DESCRIPTION: compute geopotential height
!\\
!\\
! !INTERFACE:
!
SUBROUTINE COMPUTE_GPH( region, status )
!
! !USES:
!
use Dims, only : itau, lm
use Dims, only : at, bt
use binas, only : grav
use datetime, only : tstamp
!
! !INPUT PARAMETERS:
!
integer, intent(in) :: region
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!
!EOP
!------------------------------------------------------------------------
!BOC
character(len=*), parameter :: rname = mname//'/compute_gph'
! --- local ----------------------------------
real,dimension(:,:,:),pointer :: gph, t, q
real,dimension(:,:,:),pointer :: ps
integer :: i,j,l,i0,i1,j0,j1
real :: tv,pdown,pup
! --- begin -----------------------------
! leave if not in use:
if ( .not. gph_dat(region)%used ) then
if (okdebug) then
write (gol,'(a," not used on : ",i2)') trim(gph_dat(region)%name),region; call goPr
endif
status=0; return
end if
! other meteo required:
if ( (.not. temper_dat(region)%used) .or. (.not. humid_dat(region)%used) &
.or. (.not. sp_dat(region)%used) .or. (.not. oro_dat(region)%used)) then
write (gol,'("computation of gph requires temper, humid, sp, and oro")'); call goErr
TRACEBACK; status=1; return
end if
! leave if input did not change:
if ( (.not. sp_dat(region)%changed) .and. &
(.not. temper_dat(region)%changed) .and. &
(.not. humid_dat(region)%changed) ) then
if (okdebug) then
write (gol,'(a,": not changed for region ",i2)') rname, region; call goErr
endif
status=0
return
end if
! field will be changed ...
gph_dat(region)%changed = .true.
! pointers to meteo field
ps => sp_dat(region)%data
t => temper_dat(region)%data
q => humid_dat(region)%data
gph => gph_dat(region)%data
! bounds w/o halo (same as: call Get_DistGrid( dgrid(region), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 )
i0 = gph_dat(region)%is(1)
i1 = gph_dat(region)%is(2)
j0 = gph_dat(region)%js(1)
j1 = gph_dat(region)%js(2)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! compute geo potential height
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
gph(i0:i1,j0:j1,1) = oro_dat(region)%data(i0:i1,j0:j1,1)/grav ! oro is stored in g*m
do l=1,lm(region)-1
do j=j0,j1
do i=i0,i1
tv = t(i,j,l)*(1. + 0.608*q(i,j,l))
pdown = at(l) + bt(l)*ps(i,j,1)
pup = at(l+1) + bt(l+1)*ps(i,j,1)
! rgas in different units!
gph(i,j,l+1) = gph(i,j,l) + tv*287.05*alog(pdown/pup)/grav
! note dec 2002 (MK) gph now from 1--->lm+1
end do
end do
end do
!set top of atmosphere at 200 km
gph(:,:,lm(region)+1) = 200000.0
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nullify( ps )
nullify( t )
nullify( q )
nullify( gph )
status = 0
END SUBROUTINE COMPUTE_GPH
!EOC
END MODULE METEO