!### 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:
!
! Open and read meteo files. Perform some meteo dependend calculations.
!
!
! !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.
!
! !REMARKS: --------------------------------------
!
! * All subroutines have a parameter holding its name:
! character(len=*), parameter :: rname =...
! for the traceback macro to work
!
! * A large number of surface pressure fields is 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 occure.
! 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' .
!
!
! !INTERFACE: ------------------------------------
!
module meteo
! !USES: ---------------------------------------
!
use GO , only : gol, goErr, goPr, goLabel
use GO , only : TDate
use grid , only : TllGridInfo, TLevelInfo
use tmm , only : TtmMeteo
!
use dims , only : nregions, nregions_all
!
use meteodata , only : Set, Check
use meteodata , only : TimeInterpolation
!
use meteodata , only : lli, lli_z
use meteodata , only : levi, levi_ec
!
use meteodata , only : sp1_dat, sp2_dat, sp_dat, spm_dat
use meteodata , only : tsp_dat
use meteodata , only : phlb_dat, m_dat
use meteodata , only : mfu_dat, mfv_dat, mfw_dat
use meteodata , only : pu_dat, pv_dat, pw_dat
use meteodata , only : temper_dat, humid_dat, gph_dat, omega_dat
use meteodata , only : lwc_dat, iwc_dat, cc_dat, cco_dat, ccu_dat
use meteodata , only : entu_dat, entd_dat, detu_dat, detd_dat
use meteodata , only : oro_dat
use meteodata , only : lsmask_dat
use meteodata , only : albedo_dat
use meteodata , only : sr_ecm_dat
use meteodata , only : sr_ols_dat
use meteodata , only : ci_dat
use meteodata , only : sst_dat
use meteodata , only : u10m_dat, v10m_dat
use meteodata , only : sshf_dat, slhf_dat
use meteodata , only : ewss_dat, nsss_dat
use meteodata , only : lsp_dat
use meteodata , only : cp_dat
!
use meteodata , only : sd_dat
use meteodata , only : swvl1_dat
use meteodata , only : src_dat
use meteodata , only : d2m_dat, t2m_dat
use meteodata , only : ssr_dat, ssrd_dat
use meteodata , only : str_dat, strd_dat
use meteodata , only : skt_dat
!
use meteodata , only : nveg, tv_dat, cvl_dat, cvh_dat
!
use meteodata , only : blh_dat
use meteodata , only : sf_dat
use meteodata , only : g10m_dat
use meteodata , only : ch4fire_dat
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
public :: sp1_dat, sp2_dat, sp_dat, spm_dat
public :: tsp_dat
public :: phlb_dat, m_dat
public :: mfu_dat, mfv_dat, mfw_dat
public :: pu_dat, pv_dat, pw_dat
public :: temper_dat, humid_dat, gph_dat, omega_dat
public :: lwc_dat, iwc_dat, cc_dat, cco_dat, ccu_dat
public :: entu_dat, entd_dat, detu_dat, detd_dat
public :: oro_dat
public :: lsmask_dat
public :: albedo_dat
public :: sr_ecm_dat
public :: sr_ols_dat
public :: ci_dat
public :: sst_dat
public :: u10m_dat
public :: v10m_dat
public :: sshf_dat
public :: slhf_dat
public :: ewss_dat, nsss_dat
public :: lsp_dat, cp_dat
public :: sd_dat
public :: swvl1_dat
public :: src_dat
public :: d2m_dat, t2m_dat
public :: ssr_dat, ssrd_dat
public :: str_dat, strd_dat
public :: skt_dat
public :: nveg, tv_dat, cvl_dat, cvh_dat
public :: blh_dat
public :: sf_dat
public :: g10m_dat
public :: lli, lli_z
public :: levi, levi_ec
! !PRIVATE TYPES: ------------------------------
!
type TMeteoField ! storage for 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: ----------------------------------
!
interface Setup
module procedure Setup_2d
module procedure Setup_3d
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)
!EOP
contains
! ==========================================================
! ===
! === main routines
! ===
! ==========================================================
subroutine Meteo_Init_Grids( status )
use Grid , only : Init
use dims , only : nregions, region_name
use dims , only : xbeg, xend, dx, xref, im
use dims , only : ybeg, yend, dy, yref, jm
use dims , only : nlon360, nlat180
use dims , only : echlevs, lme, a_ec, b_ec
use geometry , only : geomtryv
! --- in/out -------------------------------
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Init_Grids'
! --- local -----------------------------
integer :: n
real :: dlon, dlat
! --- begin ----------------------------
call goLabel(rname)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! setup horizontal grids
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! define mid of first west/south cells, spacing, and size:
do n = 0, nregions_all
! grid spacing:
dlon = real(xend(n)-xbeg(n))/im(n)
dlat = real(yend(n)-ybeg(n))/jm(n)
! define grid:
call Init( lli(n), xbeg(n)+dlon/2.0, dlon, im(n), &
ybeg(n)+dlat/2.0, dlat, jm(n), status, &
name=trim(region_name(n)) )
IF_NOTOK_RETURN(status=1)
end do
! zonal grids:
do n = 0, nregions_all
! grid spacing:
dlat = real(yend(n)-ybeg(n))/jm(n)
! define grid:
call Init( lli_z(n), 0.0, 360.0, 1, &
ybeg(n)+dlat/2.0, dlat, jm(n), status, &
name=trim(region_name(n))//'_z' )
IF_NOTOK_RETURN(status=1)
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)
! old at, bt
do n = 1, nregions
call geomtryv( n )
end do
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! ok
status = 0
call goLabel()
end subroutine Meteo_Init_Grids
! ***
subroutine Meteo_Done_Grids( status )
use Grid , only : Done
use Dims , only : nregions
! --- in/out -------------------------------
integer, intent(out) :: status
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Done_Grids'
! --- local -----------------------------
integer :: n
! --- begin --------------------------------
call goLabel(rname)
!
! done grid definitions
!
! horizontal grids:
do n = 0, nregions_all
call Done( lli(n), status )
IF_NOTOK_RETURN(status=1)
end do
! zonal grids:
do n = 0, nregions_all
call Done( lli_z(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)
! ***
! ok
status = 0
call goLabel()
end subroutine Meteo_Done_Grids
! ***
subroutine Meteo_Init( status )
use GO , only : TrcFile, Init, Done, ReadRc
use Binas , only : p_global
use TMM , only : Init
use dims , only : nregions, region_name
use dims , only : im, jm
use dims , only : nlon360, nlat180
use dims , only : echlevs, lme, a_ec, b_ec
use dims , only : lm, lmax_conv
use meteodata , only : Init
use global_data, only : rcfile
! --- in/out -------------------------------
integer, intent(out) :: status
! --- const --------------------------------------
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
! --- begin ----------------------------
call goLabel(rname)
! open rcfile:
call Init( rcF, rcfile, status )
IF_NOTOK_RETURN(status=1)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 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
! local grid sizes
imr = im(region)
jmr = jm(region)
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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! current pressure:
call Init( sp_dat(region), 'sp', 'Pa', 'computed', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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, &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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, &
(/1,imr/), (/1,jmr/), 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, &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,lmr/), &
rcF, (/'* ','ml ','humid'/), region, status )
IF_NOTOK_RETURN(status=1)
!
! ** computed *************************************
!
! no extra cells
halo = 0
! geopotential height(m) (lm+1, halo=0)
call Init( gph_dat(region), 'gph', 'm', 'computed', &
(/1,imr/), (/1,jmr/), halo, (/1,lmr+1/), &
'no-sourcekey', .false., 'no-destkey', status )
IF_NOTOK_RETURN(status=1)
! vertical velocity (Pa/s) (lm+1, halo=0)
call Init( omega_dat(region), 'omega', 'Pa/s', 'computed', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), region, status )
IF_NOTOK_RETURN(status=1)
! entd entrainement downdraft (im,jm,lmax_conv)
call Init_MeteoData( entd_dat(region), 'ed', 'kg/m2/s', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,lmax_conv/), &
rcF, (/'* ','ml ','convec'/), 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', &
(/1,imr/), (/1,jmr/), 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', '%', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','skt '/), region, status )
IF_NOTOK_RETURN(status=1)
#ifndef with_tmm_ecearth
! boundary layer height (m) ; instant
call Init_MeteoData( blh_dat(region), 'blh', 'm', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','blh '/), region, status )
IF_NOTOK_RETURN(status=1)
#endif
! ~~~ average field (accumulated)
! surface sensible heat flux (W/m2) ; time aver
call Init_MeteoData( sshf_dat(region), 'sshf', 'W/m2', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','nsss '/), region, status )
IF_NOTOK_RETURN(status=1)
! convective precipitation (m/s) ; time aver
call Init_MeteoData( cp_dat(region), 'cp', 'm/s', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','lsp '/), region, status )
IF_NOTOK_RETURN(status=1)
! surface solar radiation ( W/m2 ) ; time aver
call Init_MeteoData( ssr_dat(region), 'ssr', 'W/m2', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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
#ifndef with_tmm_ecearth
! 10m wind gust (m/s)
call Init_MeteoData( g10m_dat(region), 'g10m', 'm/s', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','g10m '/), region, status )
IF_NOTOK_RETURN(status=1)
#endif
! ~~~ in tmpp daily averages
! sea ice:
call Init_MeteoData( ci_dat(region), 'ci', '1', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','swvl1 '/), 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, '%', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), 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', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','srols '/), region, status )
IF_NOTOK_RETURN(status=1)
! ~~~ macc emissions
! wildfire emissions of CH4
call Init_MeteoData( ch4fire_dat(region), 'ch4fire', 'kg m**-2 s**-1', &
(/1,imr/), (/1,jmr/), halo, (/1,1/), &
rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','ch4fire '/), region, status )
IF_NOTOK_RETURN(status=1)
end do ! regions
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! close rcfile:
call Done( rcF, status )
IF_NOTOK_RETURN(status=1)
! ok
status = 0
call goLabel()
end subroutine Meteo_Init
! ***
!
! 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.
!
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:
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 --------------------------------
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( 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)
#ifndef with_tmm_ecearth
call Done( g10m_dat(n), status )
IF_NOTOK_RETURN(status=1)
#endif
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)
#ifndef with_tmm_ecearth
call Done( blh_dat(n), status )
IF_NOTOK_RETURN(status=1)
#endif
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)
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)
! ***
call Done( ch4fire_dat(n), status )
IF_NOTOK_RETURN(status=1)
! ***
end do ! regions
! ok
status = 0
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 --------------------------------
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( 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)
#ifndef with_tmm_ecearth
call Alloc( blh_dat(region), status )
IF_NOTOK_RETURN(status=1)
#endif
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)
#ifndef with_tmm_ecearth
call Alloc( g10m_dat(region), status )
IF_NOTOK_RETURN(status=1)
#endif
call Alloc( swvl1_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)
! ***
call Alloc( ch4fire_dat(region), status )
IF_NOTOK_RETURN(status=1)
! ***
end do ! regions
! ok
status = 0
call goLabel()
end subroutine Meteo_Alloc
! **************************************************************
! ***
! *** setup
! ***
! **************************************************************
!------------------------------------------------------------------------------
! TM5 !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: Meteo_Setup_Mass
!
! !DESCRIPTION:
!
! !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.
!
! !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). %data
! of SP and SP1 are only updated and used.
!
! !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
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
!EOP
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Setup_Mass'
! --- local -----------------------------
logical :: do_check_pressure
logical :: do_isfirst
integer :: n, p
integer :: idater(6)
real, allocatable :: dm_dt(:,:,:)
real :: dt_sec
integer :: l
real :: tol_rms, tol_diff
#ifdef with_prism
integer :: hour1, dhour
#endif
! --- begin --------------------------------
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
!call wrtgol( 'setup mass meteo from: ', tr1, ' to: ', tr2 ); call goPr
!
! ** mass fluxes *************************************
!
! loop over regions:
do n = 1, nregions_all
!if (mfu_dat(n)%used) then
! write (gol,'("mfu,mfv ",a)') trim(lli(n)%name); call goPr
!end if
! update horizontal u flux (unbalanced)
call Setup_MFUV( mfu_dat(n), mfv_dat(n), (/tr1,tr2/), lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
end do ! regions
! **
! loop over regions:
do n = 1, nregions_all
!if (mfw_dat(n)%used) then
! write (gol,'("mfw ",a)') trim(lli(n)%name); call goPr
!end if
! update vertical flux;
! tendency of surface pressure is by-product of vertical flux from spectral fields
! or filled with zero's
call Setup_MFW( mfw_dat(n), tsp_dat(n), (/tr1,tr2/), lli(n), 'n', levi, 'w', status )
IF_NOTOK_RETURN(status=1)
end do ! regions
!
! ** surface pressures *****************************
!
! loop over regions:
do n = 1, nregions_all
! skip ?
if ( .not. sp1_dat(n)%used ) cycle
!write (gol,'("sp1 ",a)') trim(lli(n)%name); call goPr
! 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( sp2_dat(n), (/tr1,tr1/), lli(n), '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)
! global field (first region) ?
! then match with average global surface pressure to ensure
! global mass balance;
! otherwise, match with parent grid:
if ( n == 1 ) then
call Match( 'area-aver', 'n', lli(0), sp_region0, &
lli(n), sp1_dat(n)%data(1:im(n),1:jm(n),1), status )
IF_NOTOK_RETURN(status=1)
else
p = parent(n)
call Match( 'area-aver', 'n', lli(p), sp1_dat(p)%data(1:im(p),1:jm(p),1), &
lli(n), sp1_dat(n)%data(1:im(n),1:jm(n),1), status )
IF_NOTOK_RETURN(status=1)
end if
! 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
!write (gol,'(" Reading initial sp from restart file")'); call goPr
call Restart_Read( status, region=n, 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
! fill halo cells:
call FillHalo_pm( n, status )
IF_NOTOK_RETURN(status=1)
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 ! regions
! **
! loop over regions:
do n = 1, nregions_all
! skip ?
if ( .not. sp2_dat(n)%used ) cycle
!write (gol,'("sp2 ",a)') trim(lli(n)%name); call goPr
#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( sp2_dat(n), (/tr1,tr1/), lli(n), '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(1:im(n),1:jm(n),1) = &
sp2_dat(n)%data1(1:im(n),1:jm(n),1)
sp2_dat(n)%tr2 = tr1 + IncrDate(sec=nint(dt_sec))
sp2_dat(n)%data1(1:im(n),1:jm(n),1) = &
sp2_dat(n)%data2(1:im(n),1:jm(n),1) - tsp_dat(n)%data(1:im(n),1:jm(n),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:
!write (gol,'("PLSX - setup SP2_dat at ",i2)') tr2%hour; call goPr
call Setup( sp2_dat(n), (/tr2,tr2/), lli(n), '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( sp2_dat(n), (/tr2,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
#endif
! end "ifdef with_prism"
! global field (first region) ?
! then match with average global surface pressure to ensure global mass balance;
! otherwise, match with parent grid:
if ( n == 1 ) then
call Match( 'area-aver', 'n', lli(0), sp_region0, &
lli(n), sp2_dat(n)%data(1:im(n),1:jm(n),1), status )
IF_NOTOK_RETURN(status=1)
else
p = parent(n)
call Match( 'area-aver', 'n', lli(p), sp2_dat(p)%data(1:im(p),1:jm(p),1), &
lli(n), sp2_dat(n)%data(1:im(n),1:jm(n),1), status )
IF_NOTOK_RETURN(status=1)
end if
end do ! 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
! loop over regions:
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
!write (gol,'("balance ",a)') trim(lli(n)%name); call goPr
! 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(im(n),jm(n),lm(n)) )
! mass change (kg) :
call FillMassChange( dm_dt, lli(n), levi, &
sp1_dat(n)%data(1:im(n),1:jm(n),1), &
sp2_dat(n)%data(1:im(n),1:jm(n),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)
!
do l = 1, lm(n)
call Match( 'sum', 'u', lli(p), pu_dat(p)%data1(0:im(p),1:jm(p),l), &
lli(n), pu_dat(n)%data1(0:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
!
do l = 1, lm(n)
call Match( 'sum', 'v', lli(p), pv_dat(p)%data1(1:im(p),1:jm(p)+1,l), &
lli(n), pv_dat(n)%data1(1:im(n),1:jm(n)+1,l), status )
IF_NOTOK_RETURN(status=1)
end do
!
do l = 0, lm(n)
call Match( 'sum', 'n', lli(p), pw_dat(p)%data1(1:im(p),1:jm(p),l), &
lli(n), pw_dat(n)%data1(1:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
!
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 CheckMassBalance( lli(n), &
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)), &
sp1_dat(n)%data(1:im(n),1:jm(n),1), &
sp2_dat(n)%data(1:im(n),1:jm(n),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))
call BalanceMassFluxes( lli(n), &
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)), &
pw_dat(n)%data1, dm_dt, lli(parent(n)), dt_sec, status )
IF_NOTOK_RETURN(status=1)
! 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 CheckMassBalance( lli(n), &
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)), &
sp1_dat(n)%data(1:im(n),1:jm(n),1), &
sp2_dat(n)%data(1:im(n),1:jm(n),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
! periodic boundary
if ( xcyc(n) == 1 ) then
pu_dat(n)%data1(0 ,:,:) = pu_dat(n)%data1(im(n),:,:)
pu_dat(n)%data1(im(n)+1,:,:) = pu_dat(n)%data1(1 ,:,:)
end if
! >>> 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)
!
do l = 1, lm(n)
call Match( 'sum', 'u', lli(p), pu_dat(p)%data2(0:im(p),1:jm(p),l), &
lli(n), pu_dat(n)%data2(0:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
!
do l = 1, lm(n)
call Match( 'sum', 'v', lli(p), pv_dat(p)%data2(1:im(p),1:jm(p)+1,l), &
lli(n), pv_dat(n)%data2(1:im(n),1:jm(n)+1,l), status )
IF_NOTOK_RETURN(status=1)
end do
!
do l = 0, lm(n)
call Match( 'sum', 'n', lli(p), pw_dat(p)%data2(1:im(p),1:jm(p),l), &
lli(n), pw_dat(n)%data2(1:im(n),1:jm(n),l), status )
IF_NOTOK_RETURN(status=1)
end do
!
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 CheckMassBalance( lli(n), &
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)), &
sp1_dat(n)%data(1:im(n),1:jm(n),1), &
sp2_dat(n)%data(1:im(n),1:jm(n),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:
! 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 BalanceMassFluxes( lli(n), &
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)), &
pw_dat(n)%data2, dm_dt, lli(parent(n)), dt_sec, status )
IF_NOTOK_RETURN(status=1)
! 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 CheckMassBalance( lli(n), &
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)), &
sp1_dat(n)%data(1:im(n),1:jm(n),1), &
sp2_dat(n)%data(1:im(n),1:jm(n),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
! periodic boundary
if ( xcyc(n) == 1 ) then
pu_dat(n)%data2(0 ,:,:) = pu_dat(n)%data2(im(n),:,:)
pu_dat(n)%data2(im(n)+1,:,:) = pu_dat(n)%data2(1 ,:,:)
end if
end if ! filled2
! >>>
! clear
deallocate( dm_dt )
end do ! regions
#endif
!
! ** done **************************************************
!
! ok
status = 0
call goLabel()
end subroutine Meteo_Setup_Mass
! ***
subroutine Meteo_Setup_Other( tr1, tr2, status, isfirst )
use GO , only : TDate, NewDate, rTotal, wrtgol
use GO , only : operator(-), operator(+), operator(/)
use GO , only : InterpolFractions
use datetime , only : tau2date
use tracer_data, only : Par_Check_Domain
use dims , only : nregions, im, jm, lm, parent
use dims , only : lmax_conv
use dims, only : xcyc
use Dims, only : czeta
use meteodata , only : SetData
use global_data, only : region_dat
#ifndef without_convection
use global_data, only : conv_dat
#endif
use Phys , only : ConvCloudDim
#ifdef with_prism
use restart , only : Restart_Read
#endif
! --- in/out -------------------------------
type(TDate), intent(in) :: tr1, tr2
integer, intent(out) :: status
logical, intent(in), optional :: isfirst
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Meteo_Setup_Other'
! --- local -----------------------------
logical :: do_isfirst
integer :: n, p
integer :: idater(6)
real, allocatable :: dm_dt(:,:,:)
real :: dt_sec
integer :: i, j, l
integer :: lsave
real :: tote, totd, maxe
real, pointer :: dxyp(:)
type(TDate) :: tmid
real :: alfa1, alfa2
integer :: iveg
! --- begin --------------------------------
call goLabel(rname)
! initial call ?
if ( present(isfirst) ) then
do_isfirst = isfirst
else
do_isfirst = .false.
end if
!call wrtgol( 'setup other meteo from: ', tr1, ' to: ', tr2 ); call goPr
!
! ** orography *****************************
!
! read orographies (if necessary):
do n = 1, nregions_all
!if (oro_dat(n)%used) then; write (gol,'("oro ",a)') trim(lli(n)%name); call goPr; end if
call Setup( oro_dat(n), (/tr1,tr2/), 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)
! interpolate:
spm_dat(n)%data(1:im(n),1:jm(n),1) = &
alfa1 * sp1_dat(n)%data(1:im(n),1:jm(n),1) + &
alfa2 * sp2_dat(n)%data(1:im(n),1:jm(n),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
! 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):
call Setup_TQ( temper_dat(n), humid_dat(n), (/tr1,tr2/), lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
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( temper_dat(n), (/tr1,tr2/), 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( humid_dat(n), (/tr1,tr2/), 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
!if (gph_dat(n)%used) then; write (gol,'("gph ",a)') trim(lli(n)%name); call goPr; end if
! 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( lwc_dat(n), (/tr1,tr2/), lli(n), 'n', levi, 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup( iwc_dat(n), (/tr1,tr2/), lli(n), 'n', levi, 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup_CloudCovers( cc_dat(n), cco_dat(n), ccu_dat(n), (/tr1,tr2/), lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
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):
call Setup_Convec( entu_dat(n), entd_dat(n), detu_dat(n), detd_dat(n), &
omega_dat(n), gph_dat(n), (/tr1,tr2/), lli(n), levi, status )
IF_NOTOK_RETURN(status=1)
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
! loop over grid cells
do j = 1, jm(n)
do i = 1, im(n)
! 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
!cmk calculate the rates in kg/gridbox and scale with czeta
dxyp => region_dat(n)%dxyp
do j = 1, jm(n)
do i = 1, im(n)
! 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
!
! ** surface fields *****************************
!
do n = 1, nregions_all
!write (gol,'("surface fields ",a)') trim(lli(n)%name); call goPr
! * lsmask
call Setup( lsmask_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * albedo
call Setup( albedo_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * sr_ecm
call Setup( sr_ecm_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * sr_ols
call Setup( sr_ols_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * sea ice
call Setup( ci_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * sea surface temperature
call Setup( sst_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * u10m
call Setup( u10m_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * v10m
call Setup( v10m_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * skin reservoir content
call Setup( src_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * 2m dewpoint temperature
call Setup( d2m_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * 2m temperature
call Setup( t2m_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * slhf
call Setup( slhf_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * sshf
call Setup( sshf_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * surface stress
call Setup( ewss_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup( nsss_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * convective precipitation
call Setup( cp_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * large scale stratiform precipitation
call Setup( lsp_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * surface solar radiation
call Setup( ssr_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup( ssrd_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * surface thermal radiation
call Setup( str_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup( strd_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * skin temperature
call Setup( skt_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
#ifndef with_tmm_ecearth
! * boundary layer height
call Setup( blh_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
#endif
! * snow fall and depth
call Setup( sf_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
call Setup( sd_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
#ifndef with_tmm_ecearth
! * g10m
call Setup( g10m_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
#endif
! * soil water level 1
call Setup( swvl1_dat(n), (/tr1,tr2/), 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( tv_dat(n,iveg), (/tr1,tr2/), 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( cvl_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * high vegetation cover
call Setup( cvh_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
! * MACC emissions
call Setup( ch4fire_dat(n), (/tr1,tr2/), lli(n), 'n', status )
IF_NOTOK_RETURN(status=1)
end do ! regions
!
! ** done ********************************************
!
! ok
status = 0
call goLabel()
end subroutine Meteo_Setup_Other
!------------------------------------------------------------------------------
! 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 -
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------------
!BOC
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/SetupSetup'
! --- local ----------------------------------
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 -----------------------------
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
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:
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:
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:
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
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, if needed.
!
! Then set md%data according to its type of interpolation
! (see TimeInterpolation in meteodata.F90). For constant type,
! %data => %data1.
!\\
!\\
! !INTERFACE:
!
subroutine Setup_2d( 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 ParTools , only : myid, root, Par_Broadcast
!
! !INPUT/OUTPUT PARAMETERS:
!
type(TMeteoData), intent(inout) :: md
!
! !INPUT PARAMETERS:
!
type(TDate), intent(in) :: tr(2)
type(TllGridInfo), intent(in) :: lli
character(len=1), intent(in) :: nuv
!
! !OUTPUT PARAMETERS:
!
integer, intent(out) :: status
!
! !REVISION HISTORY:
! 29 Mar 2010 - P. Le Sager -
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------------
!BOC
! --- const --------------------------------------
character(len=*), parameter :: rname = mname//'/Setup_2d'
! --- 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
! --- begin -----------------------------
call goLabel(rname)
!write (gol,'(" ",a,": ",a,l2)') rname, trim(md%name), md%used; call goPr
! leave if not in use:
if ( .not. md%used ) then
call goLabel(); status=0; return
end if
! not changed by default
md%changed = .false.
!
! Get time interval valid for that met field, and read/copy flags
!
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 field
!
! read or copy primairy field ?
if ( data1_read ) then
if (myid==root) then ! only root does IO
! fill data:
select case ( md%name )
!
! special routine for surface pressure
!
case ( 'sp', 'sps' )
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),1), &
md%tmi1, status )
IF_NOTOK_RETURN(status=1)
!
! special routine for Olsson surface roughness:
!
case ( 'srols' )
call Read_SR_OLS( tmmd, md%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),1), &
md%tmi1, status )
IF_NOTOK_RETURN(status=1)
!
! general field
!
case default
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data1_tref, data1_t1, data1_t2, &
lli, nuv, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),1), &
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, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),1), &
status )
IF_NOTOK_RETURN(status=1)
end if
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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 or copy secondary field ?
if ( data2_read ) then
if (myid==root) then ! only root does IO
! fill data:
select case ( md%name )
!
! special routine for surface pressure
!
case ( 'sp', 'sps' )
call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, md%data2(md%is(1):md%is(2),md%js(1):md%js(2),1), &
md%tmi2, status )
IF_NOTOK_RETURN(status=1)
!
! general field
!
case default
call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), &
data2_tref, data2_t1, data2_t2, &
lli, nuv, md%data2(md%is(1):md%is(2),md%js(1):md%js(2),1), &
md%tmi2, status )
IF_NOTOK_RETURN(status=1)
end select
! write meteofiles ?
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, md%data2(md%is(1):md%is(2),md%js(1):md%js(2),1), &
status )
IF_NOTOK_RETURN(status=1)
end if
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply time interpolation:
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
!
! done
!
! ok
status = 0
call goLabel()
end subroutine Setup_2d
! ***
subroutine Setup_3d( md, tr, lli, nuv, levi, nw, status )
use GO , only : TDate, wrtgol, operator(/=)
use Grid , only : TllGridInfo, TLevelInfo
use TMM , only : TMeteoInfo, ReadField, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
type(TMeteoData), intent(inout) :: md
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_3d'
! --- 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(:,:)
! --- begin -----------------------------
call goLabel(rname)
! leave if not in use:
if ( .not. md%used ) then
call goLabel(); status=0; return
end if
! not changed by default
md%changed = .false.
!
! time stuff
!
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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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(md%is(1):md%is(2),md%js(1):md%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, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),md%ls(1):md%ls(2)), &
md%tmi1, status )
IF_NOTOK_RETURN(status=1)
! write meteofiles ?
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, md%data1(md%is(1):md%is(2),md%js(1):md%js(2),md%ls(1):md%ls(2)), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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(md%is(1):md%is(2),md%js(1):md%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, md%data2(md%is(1):md%is(2),md%js(1):md%js(2),md%ls(1):md%ls(2)), &
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, md%data2(md%is(1):md%is(2),md%js(1):md%js(2),md%ls(1):md%ls(2)), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply time interpolation:
call TimeInterpolation( md, tr, status )
IF_NOTOK_RETURN(status=1)
!
! done
!
! ok
status = 0
call goLabel()
end subroutine Setup_3d
! **************************************************************
! ***
! *** mass fluxes
! ***
! **************************************************************
subroutine Setup_MFUV( md_mfu, md_mfv, tr, lli, levi, status )
use GO , only : TDate, wrtgol, operator(/=)
use Grid , only : TllGridInfo, TLevelInfo
use TMM , only : TMeteoInfo, Read_MFUV, WriteField
use meteodata, only : TMeteoData, TimeInterpolation
use ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
type(TMeteoData), intent(inout) :: md_mfu
type(TMeteoData), intent(inout) :: md_mfv
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_MFUV'
! --- 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_spu(:,:)
real, allocatable :: tmp_spv(:,:)
! --- begin -----------------------------
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
call goLabel(); status=0; return
end if
! not changed by default
md_mfu%changed = .false.
md_mfv%changed = .false.
!
! time stuff
!
! (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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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(md_mfu%is(1)-1:md_mfu%is(2),md_mfu%js(1):md_mfu%js(2) ) )
allocate( tmp_spv(md_mfv%is(1) :md_mfv%is(2),md_mfv%js(1):md_mfv%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, &
md_mfu%data1(md_mfu%is(1)-1:md_mfu%is(2),&
md_mfu%js(1) :md_mfu%js(2),&
md_mfu%ls(1)+1:md_mfu%ls(2)), &
md_mfu%tmi1, &
tmp_spv, &
md_mfv%data1(md_mfv%is(1) :md_mfv%is(2) , &
md_mfv%js(1) :md_mfv%js(2)+1, &
md_mfv%ls(1)+1:md_mfv%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, md_mfu%data1(md_mfu%is(1)-1:md_mfu%is(2),&
md_mfu%js(1) :md_mfu%js(2),&
md_mfu%ls(1)+1:md_mfu%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, md_mfv%data1(md_mfv%is(1) :md_mfv%is(2) , &
md_mfv%js(1) :md_mfv%js(2)+1, &
md_mfv%ls(1)+1:md_mfv%ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_spu )
deallocate( tmp_spv )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_mfu%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( md_mfv%data1, root, status )
IF_NOTOK_RETURN(status=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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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(md_mfu%is(1)-1:md_mfu%is(2),md_mfu%js(1):md_mfu%js(2) ) )
allocate( tmp_spv(md_mfv%is(1) :md_mfv%is(2),md_mfv%js(1):md_mfv%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, &
md_mfu%data2(md_mfu%is(1)-1:md_mfu%is(2),&
md_mfu%js(1) :md_mfu%js(2),&
md_mfu%ls(1)+1:md_mfu%ls(2)), &
md_mfu%tmi2, &
tmp_spv, &
md_mfv%data2(md_mfv%is(1) :md_mfv%is(2) , &
md_mfv%js(1) :md_mfv%js(2)+1, &
md_mfv%ls(1)+1:md_mfv%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, md_mfu%data2(md_mfu%is(1)-1:md_mfu%is(2),&
md_mfu%js(1) :md_mfu%js(2),&
md_mfu%ls(1)+1:md_mfu%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, md_mfv%data2(md_mfv%is(1) :md_mfv%is(2) , &
md_mfv%js(1) :md_mfv%js(2)+1, &
md_mfv%ls(1)+1:md_mfv%ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_spu )
deallocate( tmp_spv )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_mfu%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( md_mfv%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply 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
!
! ok
status = 0
call goLabel()
end subroutine Setup_MFUV
! ***
subroutine Setup_MFW( md_mfw, md_tsp, tr, lli, nuv, levi, nw, status )
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 ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
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'
! --- 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(:,:)
! --- begin -----------------------------
call goLabel(rname)
! leave if not in use:
if ( .not. md_mfw%used ) then
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
call goLabel(); status=1; return
end if
! not changed by default
md_mfw%changed = .false.
md_tsp%changed = .false.
!
! time stuff
!
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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2)) )
! fill data:
call Read_MFW( tmmd, md_mfw%sourcekey, &
data1_tref, data1_t1, data1_t2, &
lli, levi, &
tmp_sp, md_mfw%data1(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2),md_mfw%ls(1):md_mfw%ls(2)), &
md_tsp%data1(md_tsp%is(1):md_tsp%is(2),md_tsp%js(1):md_tsp%js(2),1), &
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, md_mfw%data1(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2),md_mfw%ls(1):md_mfw%ls(2)), &
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, md_tsp%data1(md_tsp%is(1):md_tsp%is(2),md_tsp%js(1):md_tsp%js(2),1), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_mfw%data1, root, status )
IF_NOTOK_RETURN(status=1)
!
call Par_Broadcast( md_tsp%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2)) )
! fill data:
call Read_MFW( tmmd, md_mfw%sourcekey, &
data2_tref, data2_t1, data2_t2, &
lli, levi, &
tmp_sp, md_mfw%data2(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2),md_mfw%ls(1):md_mfw%ls(2)), &
md_tsp%data2(md_tsp%is(1):md_tsp%is(2),md_tsp%js(1):md_tsp%js(2),1), &
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, md_mfw%data2(md_mfw%is(1):md_mfw%is(2),md_mfw%js(1):md_mfw%js(2),md_mfw%ls(1):md_mfw%ls(2)), &
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, &
md_tsp%data2(md_tsp%is(1):md_tsp%is(2),md_tsp%js(1):md_tsp%js(2),1), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_mfw%data2, root, status )
IF_NOTOK_RETURN(status=1)
!
call Par_Broadcast( md_tsp%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply 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
!
! ok
status = 0
call goLabel()
end subroutine Setup_MFW
! **************************************************************
! ***
! *** temperature and humidity
! ***
! **************************************************************
subroutine Setup_TQ( 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 ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
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(:,:)
! --- begin -----------------------------
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
call goLabel(); status=0; return
end if
! not changed by default
md_T%changed = .false.
md_Q%changed = .false.
!
! time stuff
!
! (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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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(md_T%is(1):md_T%is(2),md_T%js(1):md_T%js(2) ) )
! fill data:
call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_sp, &
md_T%data1(md_T%is(1):md_T%is(2),&
md_T%js(1):md_T%js(2),&
md_T%ls(1):md_T%ls(2)), &
md_T%tmi1, &
md_Q%data1(md_Q%is(1):md_Q%is(2), &
md_Q%js(1):md_Q%js(2), &
md_Q%ls(1):md_Q%ls(2) ), &
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, md_T%data1(md_T%is(1):md_T%is(2),&
md_T%js(1):md_T%js(2),&
md_T%ls(1):md_T%ls(2)), &
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, md_Q%data1(md_Q%is(1):md_Q%is(2), &
md_Q%js(1):md_Q%js(2), &
md_Q%ls(1):md_Q%ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_T%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( md_Q%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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(md_Q%is(1):md_Q%is(2),md_Q%js(1):md_Q%js(2)) )
! fill data:
call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
tmp_sp, &
md_T%data2(md_T%is(1):md_T%is(2),&
md_T%js(1):md_T%js(2),&
md_T%ls(1):md_T%ls(2)), &
md_T%tmi2, &
md_Q%data2(md_Q%is(1):md_Q%is(2), &
md_Q%js(1):md_Q%js(2), &
md_Q%ls(1):md_Q%ls(2) ), &
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, md_T%data2(md_T%is(1):md_T%is(2),&
md_T%js(1):md_T%js(2),&
md_T%ls(1):md_T%ls(2)), &
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, md_Q%data2(md_Q%is(1):md_Q%is(2), &
md_Q%js(1):md_Q%js(2), &
md_Q%ls(1):md_Q%ls(2) ), &
status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( md_T%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( md_Q%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply 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
!
! ok
status = 0
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
! ***
subroutine Meteo_CheckPressure( n, status )
use ParTools , only : myid, root, Par_AllReduce, ntracetloc
use dims , only : idate, newsrun
use dims , only : xcyc, im, jm
use redgridZoom, only : calc_pdiff
use io_hdf , only : io_write2d_32d, DFACC_CREATE
! --- in/out -----------------------------
integer, intent(in) :: n ! region
integer, intent(out) :: status
! --- const --------------------------------------
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 ------------------------------
call goLabel(rname)
! compare 'advected' pressure with read pressure
if ( .not. newsrun ) then
! compute difference between 'advected' pressure sp and read pressure sp1;
! if no tracer have to be transported on this pe, use dummy:
if ( ntracetloc == 0 ) then
pdiffmax_l = 0.0
else
call calc_pdiff( n, sp1_dat(n)%data(:,:,1), &
sp_dat(n)%data(:,:,1), pdiffmax_l )
end if
! compute maximum over all pe's:
call Par_AllReduce( pdiffmax_l, pdiffmax, 'max', status )
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. : ",es8.2," [Pa]")') pdiffmax; call goErr
write (gol,'(" treshold : ",es8.2," [Pa]")') pdiffmax_treshold; call goErr
write (gol,'("pressure arrays saved to local `pressure.hdf`")'); call goErr
if ( myid == root ) 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
call goErr; status=1; return
end if ! max diff
end if ! no newsrun
! ok
status = 0
call goLabel()
end subroutine Meteo_CheckPressure
! **************************************************************
! ***
! *** 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
call goLabel(rname)
! Pa/s = kg/s / m2 * g
! init with mass flux; revert sign from upward to downard, devide by gravity accelaration:
omega%data = - mfw%data * grav ! Pa/s m2
! loop over levels:
do l = 1, size(omega%data,3)
! devide by cell area (m2) :
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
call goLabel()
end subroutine Compute_Omega
! **************************************************************
! ***
! *** convective fluxes
! ***
! **************************************************************
subroutine Setup_Convec( 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 ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
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'
! --- 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(:,:)
! --- begin -----------------------------
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
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
!
! (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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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
write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
end if
! surface pressure field:
allocate( tmp_sp(entu%is(1):entu%is(2),entu%js(1):entu%js(2)) )
! fill data:
call Read_Convec( tmmd, entu%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
omega%data, omega%tmi, &
gph%data, gph%tmi, &
tmp_sp, &
entu%data1, entu%tmi1, entd%data1, entd%tmi1, &
detu%data1, detu%tmi1, detd%data1, 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, entu%data1, 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, entd%data1, 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, detu%data1, 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, detd%data1, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( entu%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( entd%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( detu%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( detd%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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
write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr
end if
! surface pressure field:
allocate( tmp_sp(entu%is(1):entu%is(2),entu%js(1):entu%js(2)) )
! fill data2:
call Read_Convec( tmmd, entu%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
omega%data, omega%tmi, &
gph%data, gph%tmi, &
tmp_sp, &
entu%data2, entu%tmi2, entd%data2, entd%tmi2, &
detu%data2, detu%tmi2, detd%data2, 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, entu%data2, 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, entd%data2, 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, detu%data2, 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, detd%data2, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( entu%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( entd%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( detu%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( detd%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply 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
!
! ok
status = 0
call goLabel()
end subroutine Setup_Convec
! **************************************************************
! ***
! *** cloud cover
! ***
! **************************************************************
subroutine Setup_CloudCovers( 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 ParTools , only : myid, root, Par_Broadcast
! --- in/out ----------------------------------
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'
! --- 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(:,:)
! --- begin -----------------------------
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
call goLabel(); status=0; return
end if
! not changed by default
cc%changed = .false.
cco%changed = .false.
ccu%changed = .false.
!
! time stuff
!
! (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 field
!
! primairy field ?
if ( data1_read ) then
if (myid==root) 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(cc%is(1):cc%is(2),cc%js(1):cc%js(2)) )
! fill data:
call Read_CloudCovers( tmmd, cc%sourcekey, &
data1_tref, data1_t1, data1_t2, lli, levi, &
tmp_sp, cc%data1, cc%tmi1, &
cco%data1, cco%tmi1, ccu%data1, 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, cc%data1, 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, cco%data1, 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, ccu%data1, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( cc%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( cco%data1, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( ccu%data1, root, status )
IF_NOTOK_RETURN(status=1)
! 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
! secondary field ?
if ( data2_read ) then
if (myid==root) 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(cc%is(1):cc%is(2),cc%js(1):cc%js(2)) )
! fill data2:
call Read_CloudCovers( tmmd, cc%sourcekey, &
data2_tref, data2_t1, data2_t2, lli, levi, &
tmp_sp, cc%data2, cc%tmi2, &
cco%data2, cco%tmi2, ccu%data2, 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, cc%data2, 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, cco%data2, 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, ccu%data2, status )
IF_NOTOK_RETURN(status=1)
end if
! clear
deallocate( tmp_sp )
end if ! root ?
! send to other processors if necessary:
call Par_Broadcast( cc%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( cco%data2, root, status )
IF_NOTOK_RETURN(status=1)
call Par_Broadcast( ccu%data2, root, status )
IF_NOTOK_RETURN(status=1)
! 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
!
! apply 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
!
! ok
status = 0
call goLabel()
end subroutine Setup_CloudCovers
! ***
!
! sp -> phlb -> mass
!
! NOTE: assume that halo cells in sp have been filled correctly ...
!
subroutine Pressure_to_Mass( region, status )
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
! --- in/out ----------------------------------
integer, intent(in) :: region
integer, intent(out) :: status
! --- const ----------------------------------
character(len=*), parameter :: rname = mname//'/Pressure_to_Mass'
! --- local ----------------------------------
integer :: imr, jmr, lmr
integer :: l
! --- begin ----------------------------------
! local grid size:
imr = im(region) ; jmr = jm(region) ; lmr = lm(region)
! fill pressure boundaries (Pa):
if ( phlb_dat(region)%used ) then
call HPressure( levi, sp_dat(region)%data(1:imr,1:jmr,1), &
phlb_dat(region)%data(1:imr,1:jmr,:), 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(1:imr,1:jmr,:), '*', 'm2', status ) ! kg
IF_NOTOK_RETURN(status=0)
end if
! ok
status = 0
end subroutine Pressure_to_Mass
! ***
subroutine FillHalo_pm( region, status )
use dims , only : im, jm
use dims , only : xcyc
! --- in/out ----------------------------------
integer, intent(in) :: region
integer, intent(out) :: status
! --- const ----------------------------------
character(len=*), parameter :: rname = mname//'/FillHalo_pm'
! --- local ----------------------------------
integer :: imr, jmr
! --- begin ----------------------------------
! local grid size:
imr = im(region) ; jmr = jm(region)
! cyclic boundary conditions ?
if ( xcyc(region) == 1 ) then
! fill western halo cells with east-most cells:
sp_dat(region)%data( -1:0 ,:,1) = sp_dat(region)%data(imr-1:imr,:,1)
phlb_dat(region)%data( -1:0 ,:,:) = phlb_dat(region)%data(imr-1:imr,:,:)
m_dat(region)%data( -1:0 ,:,:) = m_dat(region)%data(imr-1:imr,:,:)
! fill eastern halo cells with west-most cells:
sp_dat(region)%data(imr+1:imr+2,:,1) = sp_dat(region)%data( 1:2 ,:,1)
phlb_dat(region)%data(imr+1:imr+2,:,:) = phlb_dat(region)%data( 1:2 ,:,:)
m_dat(region)%data(imr+1:imr+2,:,:) = m_dat(region)%data( 1:2 ,:,:)
end if
! ok
status = 0
end subroutine FillHalo_pm
! ***
! mass -> phlb -> sp
subroutine Mass_to_Pressure( region, status )
use Binas, only : grav
use Grid , only : AreaOper
use dims , only : im, jm, lm
! --- in/out ----------------------------------
integer, intent(in) :: region
integer, intent(out) :: status
! --- const ----------------------------------
character(len=*), parameter :: rname = mname//'/Mass_to_Pressure'
! --- local ----------------------------------
integer :: imr, jmr, lmr
integer :: l
! --- begin ----------------------------------
! local grid size:
imr = im(region) ; jmr = jm(region); 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(1:imr,1:jmr,l) = phlb_dat(region)%data(1:imr,1:jmr,l+1) &
+ m_dat(region)%data(1:imr,1:jmr,l)*grav ! kg m/s2 = Pa m2
end do
! devide by grid cell area:
call AreaOper( lli(region), phlb_dat(region)%data(1:imr,1:jmr,:), '/', 'm2', status ) ! Pa
IF_NOTOK_RETURN(status=0)
! copy surface pressure:
spm_dat(region)%data(1:imr,1:jmr,1) = phlb_dat(region)%data(1:imr,1:jmr,1) ! Pa
! ok
status = 0
end subroutine Mass_to_Pressure
! ***
subroutine compute_gph( region, status )
use Dims, only : cdebug, kdebug, okdebug, kmain
use Dims, only : idate, itau, revert, newsrun
use Dims, only : im, jm, lm
use Dims, only : at, bt
use binas, only : grav
use global_data, only : mass_dat
use datetime, only : tstamp
use ParTools, only : myid, root
! --- in/out ----------------------------------
integer, intent(in) :: region
integer, intent(out) :: status
! --- const ----------------------------------
character(len=*), parameter :: rname = mname//'/compute_gph'
! --- local ----------------------------------
real,dimension(:,:,:),pointer :: gph, t, q
real,dimension(:,:,:),pointer :: ps
integer :: i,j,l
real :: tv,pdown,pup
! --- begin -----------------------------
! leave if not in use:
if ( .not. gph_dat(region)%used ) then
!write (*,'(" gph_dat not in use")')
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
!write (gol,'(a,": not changed")') rname; call goErr
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
if ( myid==root .and. cdebug ) call tstamp(kdebug,itau,rname)
if ( myid==root .and. okdebug ) call tstamp(kmain ,itau,rname)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! compute geo potential height
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
gph(:,:,1) = oro_dat(region)%data(:,:,1)/grav ! oro is stored in g*m
do l=1,lm(region)-1
do j=1,jm(region)
do i=1,im(region)
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 !note_mk: now in lm+1
gph(:,:,lm(region)+1) = 200000.0
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! done
nullify( ps )
nullify( t )
nullify( q )
nullify( gph )
! ok
status = 0
end subroutine compute_gph
end module meteo