!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" #include "tmm.inc" ! !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: METEO ! ! !DESCRIPTION: Routines to initialize/finalize meteo grids and data, allocate ! datasets, and fill them. Include wrappers around read/write ! meteo files. ! Perform some meteo dependend calculations (omega, gph, ! mass <=> pressure, ...) ! ! ! !REVISION HISTORY: ! ! 09 Jun 2010 - P. Le Sager ! - Fix in METEO_SETUP_MASS when reading restart files. ! - Added some (protex) doc. ! - Merge updates from EC-Earth project. ! ! 10 Aug 2010 - Arjo Segers ! - Reset previous fix since it gives differences after a restart. ! - Use 'pw_dat' instead of 'mfw_dat' for massflux balancing; ! otherwise 'mfw_dat' is changed by matching its values in a zoom ! region with the parent, and this would give tiny differences ! between a long run and two smaller runs with a restart in between. ! - Reformatted protex comments. ! ! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! ! (1) Several surface pressure fields are used: ! ! sp1,sp2 : Surface pressures at begin and end of dynamic time step. ! Their values are interpolated between surface pressures ! read from the meteorological archive (in real(4) !) ! or received from the meteorological model. ! Fields from the meteorological archive are stored into ! the 'sp2' structure, and copied from there into 'sp1'. ! ! sp : Actual surface pressure due to advection. ! In theory this field is equal to 'sp1' at the begin of a timestep, ! but due to numerical inacuracies ( real(4) vs real(8) ) ! tiny differeces occur. ! Therefore, this field is stored and restored in case of restart. ! ! spm Surface pressure for the mid of the time interval, ! thus in between 'sp1' and 'sp2' . ! ! !INTERFACE: ! MODULE METEO ! ! !USES: ! use GO, only : gol, goErr, goPr, goLabel use GO, only : TDate use partools, only : isRoot use grid, only : TllGridInfo, TLevelInfo use tmm, only : TtmMeteo ! use dims, only : nregions, nregions_all, okdebug use tm5_distgrid, only : dgrid, Get_DistGrid, GATHER, SCATTER, UPDATE_HALO use tm5_distgrid, only : SCATTER_J_BAND, SCATTER_I_BAND USE METEODATA IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Meteo_Init_Grids, Meteo_Done_Grids public :: Meteo_Init, Meteo_Done, Meteo_Alloc public :: Meteo_Setup_Mass, Meteo_Setup_Other public :: Set, Check public :: TimeInterpolation ! ! !PRIVATE TYPES: ! type TMeteoField ! storage for a single meteo field: logical :: used character(len=16) :: name character(len=16) :: unit integer :: is(2), js(2), ls(2) ! shapes real, pointer :: data(:,:,:) end type TMeteoField ! ! !INTERFACE: ! ! diff b/w para & serial: serial has LLI argument ! diff b/w 2d and 2d_n : size of 1st argument (TMeteoData) ! diff b/w 2d and 3d : two extra arguments for levels info interface Setup module procedure Setup_2d_para module procedure Setup_2d_n_para module procedure Setup_3d_para module procedure Setup_2d_serial module procedure Setup_2d_n_serial module procedure Setup_3d_serial end interface interface Setup_CONVEC module procedure Setup_CONVEC_para module procedure Setup_CONVEC_serial end interface interface Setup_CLOUDCOVERS module procedure Setup_CLOUDCOVERS_para module procedure Setup_CLOUDCOVERS_serial end interface ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'Meteo' type(TtmMeteo), save :: tmmd ! interface to TM meteo data real :: sp_region0(1,1) ! single cell global surface pressure (region 0) #ifdef with_tmm_tm5 logical, save :: use_tiedtke #endif ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: METEO_INIT_GRIDS ! ! !DESCRIPTION: initialize grids and levels for each region. Grids on the ! local domain are simply copied from the DistGrid object. !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_INIT_GRIDS( status ) ! ! !USES: ! use Grid, only : Init use dims, only : region_name use dims, only : xbeg, xend, dx, xref, im use dims, only : ybeg, yend, dy, yref, jm use dims, only : echlevs, lme, a_ec, b_ec use geometry, only : geomtryv ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Meteo_Init_Grids' integer :: n real :: dlon, dlat ! --- begin ---------------------------- call goLabel(rname) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! setup horizontal grids for the 0th one-cell grid ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! grid spacing: dlon = real(xend(0)-xbeg(0))/im(0) dlat = real(yend(0)-ybeg(0))/jm(0) ! define grid: call Init( lli(0), xbeg(0)+dlon/2.0, dlon, im(0), & ybeg(0)+dlat/2.0, dlat, jm(0), status, & name=trim(region_name(0)) ) IF_NOTOK_RETURN(status=1) ! zonal grids dlat = real(yend(0)-ybeg(0))/jm(0) call Init( lli_z(0), 0.0, 360.0, 1, & ybeg(0)+dlat/2.0, dlat, jm(0), status, & name=trim(region_name(0))//'_z' ) IF_NOTOK_RETURN(status=1) global_lli(0) = lli(0) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! local horizontal grid : get info from Distributed Grid ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do n=1, nregions_all call Get_DistGrid( dgrid(n), lli=lli(n), lli_z=lli_z(n), global_lli=global_lli(n) ) ! correct name (it defines file to read data) lli(n)%name = trim(region_name(n)) lli_z(n)%name = trim(region_name(n))//'_z' global_lli(n)%name = trim(region_name(n)) end do ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! level definition ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! setup parent level definition: call Init( levi_ec, lme, a_ec, b_ec, status ) ! ecmwf levels IF_NOTOK_RETURN(status=1) ! setup level definition: call Init( levi, levi_ec, echlevs, status ) ! tm half level selection IF_NOTOK_RETURN(status=1) ! determine "old" at, bt of dims module call geomtryv( ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! done ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ status = 0 call goLabel() END SUBROUTINE METEO_INIT_GRIDS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: METEO_DONE_GRIDS ! ! !DESCRIPTION: finalize all grids and levels used for met fields. !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_DONE_GRIDS( status ) ! ! !USES: ! use Grid, only : Done ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Meteo_Done_Grids' integer :: n ! --- begin -------------------------------- call goLabel(rname) ! horizontal (local) and zonal grids do n = 0, nregions_all call Done( lli(n), status ) IF_NOTOK_RETURN(status=1) call Done( lli_z(n), status ) IF_NOTOK_RETURN(status=1) end do ! horizontal (global) grids do n = 1, nregions_all call Done( global_lli(n), status ) IF_NOTOK_RETURN(status=1) end do ! done parent level definition: call Done( levi_ec, status ) IF_NOTOK_RETURN(status=1) ! level definition: call Done( levi, status ) IF_NOTOK_RETURN(status=1) ! done status = 0 call goLabel() END SUBROUTINE METEO_DONE_GRIDS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: METEO_INIT ! ! !DESCRIPTION: Init met fields, i.e. nullify pointers, store shape, and set ! if needed (ie used) according to meteo.rc. !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_INIT( status ) ! ! !USES: ! use GO, only : TrcFile, Init, Done, ReadRc use Binas, only : p_global use TMM, only : Init use dims, only : im, jm, lm, lmax_conv use meteodata, only : Init use global_data, only : rcfile ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Meteo_Init' ! --- local ----------------------------- integer :: region, n integer :: imr, jmr, lmr integer :: halo type(TrcFile) :: rcF integer :: iveg character(len=4) :: sveg integer :: i01, i02, j01, j02 ! --- begin ---------------------------- call goLabel(rname) ! open rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) #ifdef with_tmm_tm5 ! are convection fluxes derived (Tiedtke, sub files) or from IFS (convec files)? call ReadRc( rcF, 'tiedtke', use_tiedtke, status ) IF_NOTOK_RETURN(status=1) #endif ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! meteo database ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! setup interface to TM meteo: call Init( tmmd, rcF, status ) IF_NOTOK_RETURN(status=1) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! define meteo data ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! global mean surface pressure sp_region0 = p_global ! setup meteo fields: not in use, not allocated: do region = 1, nregions_all call Get_DistGrid( dgrid(region), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 ) lmr = lm(region) ! ! ** surface pressure ************************************* ! ! two extra horizontal cells halo = 2 ! end of interval; also reads for sp1 and spm : call Init_MeteoData( sp2_dat(region), 'sp', 'Pa', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','ml','sp'/), region, status ) IF_NOTOK_RETURN(status=1) ! check time interpolation: if ( sp2_dat(region)%tinterp(1:6) /= 'interp' ) then write (gol,'("surface pressure should be interpolated:")'); call goErr write (gol,'(" requested tinterp : ",a)') trim(sp2_dat(region)%tinterp); call goErr call goErr; status=1; return end if ! start of interval (copied from sp2_dat): call Init( sp1_dat(region), 'sp', 'Pa', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! current pressure: call Init( sp_dat(region), 'sp', 'Pa', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! surface pressure at mid of dynamic time interval: call Init( spm_dat(region), 'sp', 'Pa', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! ! ** 3D pressure and mass ************************** ! ! two extra horizontal cells (same as surface pressures) halo = 2 ! pressure at half levels (lm+1): call Init( phlb_dat(region), 'phlb', 'Pa', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! air mass: call Init( m_dat(region), 'm', 'kg', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! ! ** massfluxes ************************************* ! ! ~~ vertical ! no extra cells halo = 0 ! vertical flux (kg/s) call Init_MeteoData( mfw_dat(region), 'mfw', 'kg/s', & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & rcF, (/'* ','ml ','mflux_w'/), region, status ) IF_NOTOK_RETURN(status=1) ! vertical flux (kg/s) : BALANCED ! NOTE: data is copied from mfw, thus use same tinterp ! for correct allocation of data arrays call Init( pw_dat(region), 'pw', 'kg/s', mfw_dat(region)%tinterp, & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! tendency of surface pressure: call Init_MeteoData( tsp_dat(region), 'tsp', 'Pa/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','ml ','mflux_w'/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~ horizontal ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:imr,1:jmr ,1:lmr) in pu_t(0:imr+1,0:jmr+1,0:lmr) ! ^ ^ ^ ^ too large ! ! pv_tmpp --> pv(1:imr,1:jmr+1,1:lmr) in pv_t(0:imr+1,0:jmr+1,0:lmr) ! ^ ^ ^ ^ too large ! ! The extra cells are implemented below as halo cells. ! one extra cell halo = 1 !! east/west flux (kg/s) !call Init( mfu_dat(region), 'mfu', 'kg/s', tinterp_curr, & ! (/0,imr/), (/1,jmr/), halo, (/1,lmr/), & ! sourcekey_curr, write_meteo, status, default=destkey_curr ) !IF_NOTOK_RETURN(status=1) !! south/north flux (kg/s) !call Init( mfv_dat(region), 'mfv', 'kg/s', tinterp_curr, & ! (/1,imr/), (/0,jmr/), halo, (/1,lmr/), & ! sourcekey_curr, write_meteo, status, default=destkey_curr ) !IF_NOTOK_RETURN(status=1) ! east/west flux (kg/s) call Init_MeteoData( mfu_dat(region), 'mfu', 'kg/s', & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & rcF, (/'* ','ml ','mflux_uv'/), region, status ) IF_NOTOK_RETURN(status=1) ! south/north flux (kg/s) call Init_MeteoData( mfv_dat(region), 'mfv', 'kg/s', & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & rcF, (/'* ','ml ','mflux_uv'/), region, status ) IF_NOTOK_RETURN(status=1) !! east/west flux (kg/s) : BALANCED !call Init( pu_dat(region), 'pu', 'kg/s', 'computed', & ! (/0,imr/), (/1,jmr/), halo, (/1,lmr/), & ! 'no-sourcekey', .false., 'no-destkey', status ) !IF_NOTOK_RETURN(status=1) ! !! south/north flux (kg/s) : BALANCED !call Init( pv_dat(region), 'pv', 'kg/s', 'computed', & ! (/1,imr/), (/0,jmr/), halo, (/1,lmr/), & ! 'no-sourcekey', .false., 'no-destkey', status ) !IF_NOTOK_RETURN(status=1) halo = 1 ! east/west flux (kg/s) : BALANCED ! NOTE: data is copied from mfu, thus use same tinterp ! for correct allocation of data arrays call Init( pu_dat(region), 'pu', 'kg/s', mfu_dat(region)%tinterp, & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! south/north flux (kg/s) : BALANCED ! NOTE: data is copied from mfv, thus use same tinterp ! for correct allocation of data arrays call Init( pv_dat(region), 'pv', 'kg/s', mfv_dat(region)%tinterp, & (/i01,i02/), (/j01,j02/), halo, (/0,lmr/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! ! ** temperature ************************************* ! ! no extra cells halo = 0 ! temperature (K) (halo=0) call Init_MeteoData( temper_dat(region), 'T', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','temper'/), region, status ) IF_NOTOK_RETURN(status=1) ! ! ** humidity ************************************* ! ! no extra cells halo = 0 ! humidity (kg/kg) (halo = 0) call Init_MeteoData( humid_dat(region), 'Q', 'kg/kg', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','humid'/), region, status ) IF_NOTOK_RETURN(status=1) ! ! ** computed ************************************* ! halo = 1 ! halo needed for station output in USER_OUTPUT_AEROCOM ! geopotential height(m) (lm+1, halo=0) call Init( gph_dat(region), 'gph', 'm', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! no extra cells halo = 0 ! vertical velocity (Pa/s) (lm+1, halo=0) call Init( omega_dat(region), 'omega', 'Pa/s', 'computed', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr+1/), & 'no-sourcekey', .false., 'no-destkey', status ) IF_NOTOK_RETURN(status=1) ! ! ** clouds ************************************* ! ! no extra cells halo = 0 ! lwc liquid water content (kg/kg) (halo=0) call Init_MeteoData( lwc_dat(region), 'CLWC', 'kg/kg', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','cloud'/), region, status ) IF_NOTOK_RETURN(status=1) ! iwc ice water content (kg/kg) (halo=0) call Init_MeteoData( iwc_dat(region), 'CIWC', 'kg/kg', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','cloud'/), region, status ) IF_NOTOK_RETURN(status=1) ! cc cloud cover (fraction) (halo=0) call Init_MeteoData( cc_dat(region), 'CC', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','cloud'/), region, status ) IF_NOTOK_RETURN(status=1) ! cco cloud cover overhead = above bottom of box (fraction) (halo=0) call Init_MeteoData( cco_dat(region), 'CCO', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','cloud'/), region, status ) IF_NOTOK_RETURN(status=1) ! ccu cloud cover underfeet = below top of box (fraction) (halo=0) call Init_MeteoData( ccu_dat(region), 'CCU', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,lmr/), & rcF, (/'* ','ml ','cloud'/), region, status ) IF_NOTOK_RETURN(status=1) ! ! ** convection ************************************* ! ! no extra cells halo = 0 ! entu entrainement updraft call Init_MeteoData( entu_dat(region), 'eu', 'kg/m2/s', & (/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), & rcF, (/'* ','ml ','convec'/), region, status ) IF_NOTOK_RETURN(status=1) ! entd entrainement downdraft (im,jm,lmax_conv) call Init_MeteoData( entd_dat(region), 'ed', 'kg/m2/s', & (/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), & rcF, (/'* ','ml ','convec'/), region, status ) IF_NOTOK_RETURN(status=1) ! detu detrainement updraft call Init_MeteoData( detu_dat(region), 'du', 'kg/m2/s', & (/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), & rcF, (/'* ','ml ','convec'/), region, status ) IF_NOTOK_RETURN(status=1) ! detd detrainement downdraft call Init_MeteoData( detd_dat(region), 'dd', 'kg/m2/s', & (/i01,i02/), (/j01,j02/), halo, (/1,lmax_conv/), & rcF, (/'* ','ml ','convec'/), region, status ) IF_NOTOK_RETURN(status=1) ! ! *** surface fields ! ! no extra cells halo = 0 ! orography (m*[g]) call Init_MeteoData( oro_dat(region), 'oro', 'm m/s2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.const','sfc.an ','oro '/), region, status ) IF_NOTOK_RETURN(status=1) ! land/sea mask (%) call Init_MeteoData( lsmask_dat(region), 'lsm', '%', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.const','sfc.an ','lsm '/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~~ instantaneous fields ! sea surface temperatue: call Init_MeteoData( sst_dat(region), 'sst', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','sst '/), region, status ) IF_NOTOK_RETURN(status=1) ! 10m u wind (m/s) call Init_MeteoData( u10m_dat(region), 'u10m', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','u10m '/), region, status ) IF_NOTOK_RETURN(status=1) ! 10m v wind (m/s) call Init_MeteoData( v10m_dat(region), 'v10m', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','v10m '/), region, status ) IF_NOTOK_RETURN(status=1) ! 10m wind speed (m/s) call Init_MeteoData( wspd_dat(region), 'wspd', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','wspd '/), region, status ) IF_NOTOK_RETURN(status=1) ! skin reservoir content (m water) ; instant call Init_MeteoData( src_dat(region), 'src', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','src '/), region, status ) IF_NOTOK_RETURN(status=1) ! 2 meter dewpoint temperature (K) ; instant call Init_MeteoData( d2m_dat(region), 'd2m', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','d2m '/), region, status ) IF_NOTOK_RETURN(status=1) ! 2 meter temperature (K) ; instant call Init_MeteoData( t2m_dat(region), 't2m', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','t2m '/), region, status ) IF_NOTOK_RETURN(status=1) ! skin temperature (K) ; instant call Init_MeteoData( skt_dat(region), 'skt', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','skt '/), region, status ) IF_NOTOK_RETURN(status=1) ! boundary layer height (m) ; instant call Init_MeteoData( blh_dat(region), 'blh', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.fc ','blh '/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~~ average field (accumulated) ! surface sensible heat flux (W/m2) ; time aver call Init_MeteoData( sshf_dat(region), 'sshf', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','sshf '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface latent heat flux (W/m2) ; time aver call Init_MeteoData( slhf_dat(region), 'slhf', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','slhf '/), region, status ) IF_NOTOK_RETURN(status=1) ! east-west surface stress (N/m2); time aver call Init_MeteoData( ewss_dat(region), 'ewss', 'N/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ewss '/), region, status ) IF_NOTOK_RETURN(status=1) ! north-south surface stress (N/m2) ; time aver call Init_MeteoData( nsss_dat(region), 'nsss', 'N/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','nsss '/), region, status ) IF_NOTOK_RETURN(status=1) ! convective precipitation (m/s) ; time aver call Init_MeteoData( cp_dat(region), 'cp', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','cp '/), region, status ) IF_NOTOK_RETURN(status=1) ! large scale stratiform precipitation (m/s) ; time aver call Init_MeteoData( lsp_dat(region), 'lsp', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','lsp '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface solar radiation ( W/m2 ) ; time aver call Init_MeteoData( ssr_dat(region), 'ssr', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ssr '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface solar radiation downwards ( W/m2 ) ; time aver call Init_MeteoData( ssrd_dat(region), 'ssrd', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','ssrd '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface thermal radiation ( W/m2 ) ; time aver call Init_MeteoData( str_dat(region), 'str', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','str '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface thermal radiation downwards ( W/m2 ) ; time aver call Init_MeteoData( strd_dat(region), 'strd', 'W/m2', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','strd '/), region, status ) IF_NOTOK_RETURN(status=1) ! snow fall (m water eqv); time aver call Init_MeteoData( sf_dat(region), 'sf', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','sf '/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~~ time averages in grib files tfc+[12,15] etc ! 10m wind gust (m/s) call Init_MeteoData( g10m_dat(region), 'g10m', 'm/s', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.aver','sfc.fc ','g10m '/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~~ in tmpp daily averages ! sea ice: call Init_MeteoData( ci_dat(region), 'ci', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','ci '/), region, status ) IF_NOTOK_RETURN(status=1) ! snow depth (m water eqv); day aver ? call Init_MeteoData( sd_dat(region), 'sd', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','sd '/), region, status ) IF_NOTOK_RETURN(status=1) ! volumetric soil water layer 1 ( m3 water / m3 soil) ; day aver ? call Init_MeteoData( swvl1_dat(region), 'swvl1', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','swvl1 '/), region, status ) IF_NOTOK_RETURN(status=1) ! soil temperature layer 1 (K) call Init_MeteoData( stl1_dat(region), 'stl1', 'K', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.fc ','stl1 '/), region, status ) IF_NOTOK_RETURN(status=1) ! vegetation type (%) ; day aver do iveg = 1, nveg write (sveg,'("tv",i2.2)') iveg call Init_MeteoData( tv_dat(region,iveg), sveg, '%', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status ) IF_NOTOK_RETURN(status=1) end do ! low vegetation cover (0-1) ; day aver call Init_MeteoData( cvl_dat(region), 'cvl', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status ) IF_NOTOK_RETURN(status=1) ! high vegetation cover (0-1) ; day aver call Init_MeteoData( cvh_dat(region), 'cvh', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','veg '/), region, status ) IF_NOTOK_RETURN(status=1) ! albedo ; daily average call Init_MeteoData( albedo_dat(region), 'albedo', '1', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','albedo '/), region, status ) IF_NOTOK_RETURN(status=1) ! surface roughness (ecmwf,ncep) call Init_MeteoData( sr_ecm_dat(region), 'sr', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','sr '/), region, status ) IF_NOTOK_RETURN(status=1) ! ~~~ monthly data ! surface roughness (olsson) ; monthly call Init_MeteoData( sr_ols_dat(region), 'srols', 'm', & (/i01,i02/), (/j01,j02/), halo, (/1,1/), & rcF, (/'* ','sfc ','sfc.inst','sfc.day ','sfc.an ','srols '/), region, status ) IF_NOTOK_RETURN(status=1) end do ! regions ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! done ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! close rcfile: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 call goLabel() END SUBROUTINE METEO_INIT !EOC ! ! Read multiple keys in rcfile to setup meteodata structure. ! The following keys are read: ! ! meteo.tinterp. <-- time interpolation ! tmm.sourcekey.. <-- input file name description ! tmm.output.. <-- write meteo ? ! tmm.destkey.. <-- output file name description ! ! where is first '*' and then set to the grid name, ! and is set to each of the provided keys. ! ! Called for region=1..nregions_all SUBROUTINE INIT_METEODATA( md, name, unit, is, js, halo, ls, & rcF, rcs, region, status ) use GO, only : TRcFile, ReadRc use Dims, only : nregions, nregions_max use MeteoData, only : TMeteoData, Init, Set ! --- in/out ------------------------------------- type(TMeteoData), intent(out) :: md character(len=*), intent(in) :: name, unit integer, intent(in) :: is(2), js(2) integer, intent(in) :: halo integer, intent(in) :: ls(2) type(TRcFile), intent(inout) :: rcF character(len=*), intent(in) :: rcs(:) integer, intent(in) :: region integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/Init_MeteoData' ! --- local ------------------------------------- character(len=10) :: tinterp character(len=256) :: sourcekey logical :: write_meteo character(len=256) :: destkey logical :: used ! --- begin ------------------------------------- ! time interpolation : call ReadRc( rcF, 'meteo.tinterp', rcs, tinterp, status ) IF_NOTOK_RETURN(status=1) ! source filenames: call ReadRc( rcF, 'tmm.sourcekey.*', rcs, sourcekey, status, default='no-sourcekey' ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'tmm.sourcekey.'//trim(lli(region)%name), rcs, sourcekey, status, default=sourcekey ) IF_ERROR_RETURN(status=1) ! write flag: call ReadRc( rcF, 'tmm.output.*', rcs, write_meteo, status, default=.false. ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'tmm.output.'//trim(lli(region)%name), rcs, write_meteo, status, default=write_meteo ) IF_ERROR_RETURN(status=1) ! destination filenames: if ( write_meteo ) then call ReadRc( rcF, 'tmm.destkey.*', rcs, destkey, status, default='no-destkey' ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'tmm.destkey.'//trim(lli(region)%name), rcs, destkey, status, default=destkey ) IF_ERROR_RETURN(status=1) else destkey = 'no-destkey' end if ! define meteo data, ! but should be marked as 'used' to be allocated and filled: call Init( md, name, unit, tinterp, is, js, halo, ls, & sourcekey, write_meteo, destkey, status ) IF_NOTOK_RETURN(status=1) ! read this type of meteo ? ! only regions 1..nregions or the extra fiels above nregions_max ! could be in use: ! [all regions, but do "if test", because nregions may be different from nregions_max] if ( (region <= nregions) .or. (region > nregions_max) ) then call ReadRc( rcF, 'meteo.read.*', rcs, used, status, default=.false. ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'meteo.read.'//trim(lli(region)%name), rcs, used, status, default=used ) IF_ERROR_RETURN(status=1) else used = .false. end if ! in use ? call Set( md, status, used=used ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE INIT_METEODATA ! *** SUBROUTINE METEO_DONE( status ) use TMM, only : Done use Dims, only : nregions_all use meteodata, only : Done ! --- in/out ------------------------------- integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/Meteo_Done' ! --- local ----------------------------- integer :: n integer :: iveg ! --- begin -------------------------------- 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) call Done( wspd_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( g10m_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( src_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( d2m_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( t2m_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( blh_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( sshf_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( slhf_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( ewss_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( nsss_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( cp_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( lsp_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( ssr_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( ssrd_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( str_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( strd_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( skt_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( sd_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( sf_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( swvl1_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( stl1_dat(n), status ) IF_NOTOK_RETURN(status=1) do iveg = 1, nveg call Done( tv_dat(n,iveg), status ) IF_NOTOK_RETURN(status=1) end do call Done( cvl_dat(n), status ) IF_NOTOK_RETURN(status=1) call Done( cvh_dat(n), status ) IF_NOTOK_RETURN(status=1) ! *** end do ! regions ! 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( wspd_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( src_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( d2m_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( t2m_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( skt_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( blh_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( sshf_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( slhf_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( ewss_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( nsss_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( cp_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( lsp_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( ssr_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( ssrd_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( str_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( strd_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( sd_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( sf_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( g10m_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( swvl1_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( stl1_dat(region), status ) IF_NOTOK_RETURN(status=1) do iveg = 1, nveg call Alloc( tv_dat(region,iveg), status ) IF_NOTOK_RETURN(status=1) end do call Alloc( cvl_dat(region), status ) IF_NOTOK_RETURN(status=1) call Alloc( cvh_dat(region), status ) IF_NOTOK_RETURN(status=1) ! *** end do ! regions ! ok status = 0 call goLabel() END SUBROUTINE METEO_ALLOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: METEO_SETUP_MASS ! ! !DESCRIPTION: Set up Mass FLuxes and Surface Pressures !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_SETUP_MASS( tr1, tr2, status, isfirst, check_pressure ) ! ! !USES: ! use go, only : TDate, rTotal, operator(-), wrtgol use go, only : IncrDate, operator(+), Get use grid, only : Match, TllGridInfo, assignment(=), Done use Grid, only : FillMassChange, BalanceMassFluxes, CheckMassBalance use dims, only : nregions, im, jm, lm, parent use dims, only : xcyc use meteodata, only : SetData ! to copy %data and %tr from one MD to another #ifdef with_prism use meteodata, only : TimeInterpolation #endif use restart, only : Restart_Read ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: tr1, tr2 ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status logical, intent(in), optional :: check_pressure logical, intent(in), optional :: isfirst ! ! !REVISION HISTORY: ! ! 12 Mar 2010 - P. Le Sager - Fix when reading restart files. Added ! protex doc. Added comments. ! 9 Jun 2010 - P. Le Sager - Merged with updates for EC-Earth project. ! ! 10 Aug 2010, Arjo Segers ! Reset previous fix since it makes a restart different from a long run. ! Use 'pw_dat' instead of 'mfw_dat' since otherwise the later changed ! while matching a zoom region with its parent, and this would give ! tiny differences during a restart of a zoomed run. ! ! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! ! push of Surf Press is done with sp2 (the only one on which we call ! setup -ie the only one for which %data1 and %data2 matter). Only %data ! of SP and SP1 are updated and used, and not their %data1 and %data2. ! !------------------------------------------------------------------------------ !EOP character(len=*), parameter :: rname = mname//'/Meteo_Setup_Mass' logical :: do_check_pressure, WestBorder, NorthBorder logical :: do_isfirst integer :: n, p integer :: idater(6) real, allocatable :: dm_dt(:,:,:) real :: dt_sec integer :: l, i0, i1, j0, j1, is, js, ie, je real :: tol_rms, tol_diff type(TllGridInfo) :: L_lli real, pointer :: field(:,:), field_parent(:,:) ! work arrays real, allocatable :: islice(:,:), jslice(:,:), bigIslice(:,:), bigJslice(:,:) real, allocatable :: full_pu(:,:,:), full_pv(:,:,:), full_pw(:,:,:), full_dm_dt(:,:,:) #ifdef with_prism integer :: hour1, dhour #endif ! --- begin -------------------------------- call goLabel(rname) ! check pressure ? if ( present(check_pressure) ) then do_check_pressure = check_pressure else do_check_pressure = .false. end if ! EC-EARTH do not check pressure since we would compare pressure from ! restart file at 00 ( as written by previous chunk) with the one received ! at 06 (or whatever the coupling period is) from IFS do_check_pressure = .false. ! initial call ? if ( present(isfirst) ) then do_isfirst = isfirst else do_isfirst = .false. end if ! ! ** MASS FLUXES ************************************* ! do n = 1, nregions_all L_lli = global_lli(n) call Setup_UVW( n, mfu_dat(n), mfv_dat(n), mfw_dat(n), tsp_dat(n),& (/tr1,tr2/), L_lli, 'n', levi, 'w', status) IF_NOTOK_RETURN(status=1) end do ! ! ** SURFACE PRESSURES : SP1, SP ***************************** ! REG: do n = 1, nregions_all if ( .not. sp1_dat(n)%used ) cycle L_lli = global_lli(n) ! Advance 'next' surface pressure (a/k/a sp2%data) to start of ! new interval tr1. If start of a new meteo interval, then data ! is automatically read from file, or recieved from coupler ! with OASIS/prism call Setup( n, sp2_dat(n), (/tr1,tr1/), L_lli, 'n', status ) IF_NOTOK_RETURN(status=1) ! copy SP2 into SP1 (%data and %tr) call SetData( sp1_dat(n), sp2_dat(n), status ) IF_NOTOK_RETURN(status=1) ! GATHER sp1 array (dummy if not root) !----------------- ! ...of parent region, if any: if ( n /= 1 ) then p = parent(n) if (isRoot) then allocate( field_parent(im(p), jm(p)) ) else allocate( field_parent(1,1) ) end if call GATHER( dgrid(p), sp1_dat(p)%data(:,:,1), field_parent, sp1_dat(p)%halo, status ) IF_NOTOK_RETURN(status=1) end if ! ...of current region: if (isRoot) then allocate( field(im(n), jm(n)) ) else allocate( field(1,1) ) end if call GATHER( dgrid(n), sp1_dat(n)%data(:,:,1), field, sp1_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! MATCH surface pressures to ensure mass balance !----------------- if (isRoot) then ! IF global field (i.e first region) : match global region with one-cell ! world value (average global surface pressure), ELSE match with parent if ( n == 1 ) then call Match( 'area-aver', 'n', lli(0), sp_region0, & global_lli(n), field, status ) IF_NOTOK_RETURN(status=1) else call Match( 'area-aver', 'n', global_lli(p), field_parent, & global_lli(n), field, status ) IF_NOTOK_RETURN(status=1) endif end if ! SCATTER sp1 array, and clean !----------------- call SCATTER( dgrid(n), sp1_dat(n)%data(:,:,1), field, sp1_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) deallocate(field) if ( n /= 1 ) then if (associated(field_parent)) deallocate(field_parent) endif ! Set SP !----------------- ! If initial call *OR* beginning of new meteo interval in case of coupled runs, ! then set current surface pressure to just read/advanced sp1. Otherwise, sp ! remains filled with the advected pressure. #ifdef with_prism 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 call Get( tr1, hour=hour1 ) !else partial coupling - not handled here endif ! at begin of dhour interval ? if ( modulo(hour1,dhour) == 0 ) then #else if ( do_isfirst ) then #endif ! write (gol,'(" copy SP1 to SP ...")'); call goPr ! copy sp1 into sp : call SetData( sp_dat(n), sp1_dat(n), status ) IF_NOTOK_RETURN(status=1) ! fill pressure and mass from sp call Pressure_to_Mass( n, status ) IF_NOTOK_RETURN(status=1) #ifndef with_prism ! Eventually replace with fields from restart file, since meteo from hdf meteo ! files is in real(4) while computed pressures and mass are probably in ! real(8). But not for coupled runs, since they receive pressures from ! IFS. Note that this call will do nothing if istart=32. call Restart_Read( status, surface_pressure=.true., pressure=.true., air_mass=.true. ) IF_NOTOK_RETURN(status=1) !AJS>>> don't do this! sp1 contains data interpolated between ! fields received from the archive or the coupled model, ! while sp contains the actual pressure after advection. !! copy sp into sp1 (PLS, 29-3-2010) !call SetData( sp1_dat(n), sp_dat(n), status ) !IF_NOTOK_RETURN(status=1) !<<< #endif end if ! first or new coupling meteo !! fill initial pressure and mass arrays, !! eventually apply cyclic boundaries to mass !call Meteo_SetupMass( n, status ) !IF_NOTOK_RETURN(status=1) ! check 'advected' pressure ? if ( do_check_pressure) then ! compare 'advected' pressure still in sp with just read ! pressure sp1 : diff b/w sp%data and sp1%data call Meteo_CheckPressure( n, status ) IF_NOTOK_RETURN(status=1) end if END DO REG ! regions ! ! ** SURFACE PRESSURES : SP2 ***************************** ! REG2: do n = 1, nregions_all if ( .not. sp2_dat(n)%used ) cycle ! grid and bounds L_lli = global_lli(n) i0 = sp2_dat(n)%is(1) i1 = sp2_dat(n)%is(2) j0 = sp2_dat(n)%js(1) j1 = sp2_dat(n)%js(2) #ifdef with_prism ! sp2 for prism coupler is computed from : sp(t2) = sp(t1) + tsp*(t2-t1) if ( sp2_dat(n)%sourcekey(1:6) == 'prism:' ) then select case ( sp2_dat(n)%tinterp ) case ( 'interp6' ) ; dhour = 6 case ( 'interp3' ) ; dhour = 3 case ( 'interp2' ) ; dhour = 2 case ( 'interp1' ) ; dhour = 1 case default write (gol,'("unsupported time interpolation:")'); call goErr write (gol,'(" md%tinterp : ",a)') trim(sp2_dat(n)%tinterp); call goErr TRACEBACK; status=1; return end select ! current interval [tr1,tr2] at begin of dhour interval ? call Get( tr1, hour=hour1 ) if ( modulo(hour1,dhour) == 0 ) then ! reset sp2_dat%data1 and sp2_dat%data2: ! Read into sp2%data1 : surface pressure received for tr1 ! 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 sp2_dat(n)%filled1 = .false. sp2_dat(n)%filled2 = .false. call Setup( n, sp2_dat(n), (/tr1,tr1/), L_lli, 'n', status ) IF_NOTOK_RETURN(status=1) ! %data2 = %data1 + tsp * dhour*3600.0 !write (gol,'(" compute sp2%data2 from sp2%data1 and sp tendency ...")'); call goPr dt_sec = dhour * 3600.0 ! sec sp2_dat(n)%data2(i0:i1,j0:j1,1) = & sp2_dat(n)%data1(i0:i1,j0:j1,1) + tsp_dat(n)%data(i0:i1,j0:j1,1) * dt_sec sp2_dat(n)%tr2 = tr1 + IncrDate(sec=nint(dt_sec)) endif ! 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) else ! PLS: this one is never used apparently... AJS: it might be used in a partial ! coupling with only some fields exchanged and others read; this was often the ! case during the first coupling experiments, and might be useful for testing ! advance 'next' surface pressure to end of interval: call Setup( n, sp2_dat(n), (/tr2,tr2/), L_lli, 'n', status ) IF_NOTOK_RETURN(status=1) end if ! it is prism sourcekey #else ! advance 'next' surface pressure to end of interval: call Setup( n, sp2_dat(n), (/tr2,tr2/), L_lli, 'n', status ) IF_NOTOK_RETURN(status=1) #endif /* WITH_PRISM */ ! GATHER sp2 array (dummy if not root) !----------------- ! ...of parent region, if any: if ( n /= 1 ) then p = parent(n) if (isRoot) then allocate( field_parent(im(p), jm(p)) ) else allocate( field_parent(1,1) ) end if call GATHER( dgrid(p), sp2_dat(p)%data(:,:,1), field_parent, sp2_dat(p)%halo, status ) IF_NOTOK_RETURN(status=1) end if ! ...of current region: if (isRoot) then allocate( field(im(n), jm(n)) ) else allocate( field(1,1) ) end if call GATHER( dgrid(n), sp2_dat(n)%data(:,:,1), field, sp2_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! MATCH surface pressures to ensure mass balance !----------------- if (isRoot) then ! IF global field (i.e first region) : match global region with one-cell ! world value (average global surface pressure), ELSE match with parent if ( n == 1 ) then call Match( 'area-aver', 'n', lli(0), sp_region0, & global_lli(n), field, status ) IF_NOTOK_RETURN(status=1) else call Match( 'area-aver', 'n', global_lli(p), field_parent, & global_lli(n), field, status ) IF_NOTOK_RETURN(status=1) endif end if ! SCATTER sp2 array, and clean !----------------- call SCATTER( dgrid(n), sp2_dat(n)%data(:,:,1), field, sp2_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) deallocate(field) if ( n /= 1 ) then if (associated(field_parent)) deallocate(field_parent) endif END DO REG2 ! regions #ifndef without_advection ! ! ** MASS BALANCE ***************************** ! ! NOTE: since only the surface pressure gradient is used, ! it is not necessary to use the data1 and data2 arrays do n = 1, nregions_all if ( .not. pu_dat(n)%used ) cycle if ( .not. pv_dat(n)%used ) cycle if ( .not. pw_dat(n)%used ) cycle L_lli = global_lli(n) i0 = sp2_dat(n)%is(1) i1 = sp2_dat(n)%is(2) j0 = sp2_dat(n)%js(1) j1 = sp2_dat(n)%js(2) ! local indices and tile location (is, ie, js, je must be equal to i0, i1, j0, j1 BTW) CALL GET_DISTGRID( dgrid(n), & I_STRT=is, I_STOP=ie, & J_STRT=js, J_STOP=je, & hasWestBorder=WestBorder, hasNorthBorder=NorthBorder) ! length of time step between sp1 and sp2: dt_sec = rTotal( sp2_dat(n)%tr(1) - sp1_dat(n)%tr(1), 'sec' ) ! allocate temporary array: allocate(dm_dt(i0:i1,j0:j1,lm(n))) ! mass change (kg) : call FillMassChange( dm_dt, lli(n), levi, & sp1_dat(n)%data(i0:i1,j0:j1,1), & sp2_dat(n)%data(i0:i1,j0:j1,1), & status ) IF_NOTOK_RETURN(status=1) ! mass tendency (kg/s) : dm_dt = dm_dt / dt_sec ! kg/s ! >>> data1 >>> ! initial guess for balanced fluxes are unbalanced fluxes: pu_dat(n)%data1 = mfu_dat(n)%data1 pu_dat(n)%filled1 = mfu_dat(n)%filled1 pu_dat(n)%tr1 = mfu_dat(n)%tr1 pv_dat(n)%data1 = mfv_dat(n)%data1 pv_dat(n)%filled1 = mfv_dat(n)%filled1 pv_dat(n)%tr1 = mfv_dat(n)%tr1 pw_dat(n)%data1 = mfw_dat(n)%data1 pw_dat(n)%filled1 = mfw_dat(n)%filled1 pw_dat(n)%tr1 = mfw_dat(n)%tr1 !#ifdef with_prism ! EC-Earth 2.4 discussion - Coupling has changed in EC-Earth 3. ! ! 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. ! ! 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. ! !#else ! CHECK INITIAL MASS BALANCE: ! ----------------------------------- ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! tolerance for difference between sp from mass fluxes and sp from tendency: tol_rms = 1.0e-4 ! max rms tol_diff = 1.0e-3 ! max absolute difference CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data1, pu_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data1, pv_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call CheckMassBalance( lli(n), & pu_dat(n)%data1(i0-1:i1, j0:j1 , 1:lm(n) ), & pv_dat(n)%data1( i0:i1, j0:j1+1, 1:lm(n) ), & sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), & sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), & dt_sec, tol_rms, tol_diff, status ) if (status/=0) then write (gol,'("initial mass imbalance too large for region ",i2)') n; call goErr call goErr; status=1; return end if !#endif ! BALANCE HORIZONTAL FLUXES ! ----------------------------------- ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! needs to be done globally... so gather data if (isRoot) then allocate(full_pu( 0:im(n), 1:jm(n), 0:lm(n)) ) ! must have same number of levels as mfu allocate(full_pv( 1:im(n), 1:jm(n)+1, 0:lm(n)) ) allocate(full_pw( 1:im(n), 1:jm(n), 0:lm(n)) ) ! used also as temp arr in comm allocate(full_dm_dt(im(n), jm(n), lm(n)) ) else allocate( full_pu(1,1,1) ) allocate( full_pv(1,1,1) ) allocate( full_pw(1,1,1) ) allocate( full_dm_dt(1,1,1)) end if !for slice scattering allocate(islice(j0:j1,0:lm(n))) allocate(jslice(i0:i1,0:lm(n))) if (isRoot) then allocate(bigIslice(1:jm(n),0:lm(n))) allocate(bigJslice(1:im(n),0:lm(n))) else allocate(bigIslice(1,1)) allocate(bigJslice(1,1)) end if call GATHER( dgrid(n), pu_dat(n)%data1, full_pw, pu_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then full_pu(1:im(n),1:jm(n),:) = full_pw full_pu(0,:,:) = full_pu(im(n),:,:) ! East-West periodicity end if call GATHER( dgrid(n), pv_dat(n)%data1, full_pw, pv_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then full_pv(1:im(n),1:jm(n),:) = full_pw full_pv(:,jm(n)+1,:) = full_pv(:,1,:) ! donut periodicity end if call GATHER( dgrid(n), dm_dt, full_dm_dt, 0, status ) IF_NOTOK_RETURN(status=1) call GATHER( dgrid(n), pw_dat(n)%data1, full_pw, pw_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then call BalanceMassFluxes( global_lli(n), & full_pu(0:im(n),1:jm(n) ,1:lm(n)), & full_pv(1:im(n),1:jm(n)+1,1:lm(n)), & full_pw, full_dm_dt, global_lli(parent(n)), dt_sec, status ) IF_NOTOK_RETURN(status=1) end if call SCATTER( dgrid(n), pw_dat(n)%data1, full_pw, pw_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if(isRoot) full_pw = full_pu(1:im(n),1:jm(n),:) call SCATTER( dgrid(n), pu_dat(n)%data1, full_pw, pu_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! scatter extra column full_pu(0,:,:) - needed only for noncyclic zoom ! region, for others update_halo takes care of it [FIXME: could had a ! test around these 3 lines ] if(isRoot) bigIslice = full_pu(0,1:jm(n),:) CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1) if(WestBorder)pu_dat(n)%data1(0,j0:j1,0:lm(n)) = islice if(isRoot) full_pw = full_pv(1:im(n),1:jm(n),:) call SCATTER( dgrid(n), pv_dat(n)%data1, full_pw, pv_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! Scatter PV(:,jm+1,:) if(isroot) bigJslice=full_pv(1:im(n),jm(n)+1,:) CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n)) if(NorthBorder)pv_dat(n)%data1(i0:i1,jm(n)+1,0:lm(n))=jslice ! CHECK FINAL MASS BALANCE: ! ----------------------------------- ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! tolerance for difference between sp from mass fluxes and sp from tendency: tol_rms = 1.0e-7 ! max rms tol_diff = 1.0e-6 ! max absolute difference CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data1, pu_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data1, pv_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call CheckMassBalance( lli(n), & pu_dat(n)%data1(i0-1:i1, j0:j1 , 1:lm(n) ), & pv_dat(n)%data1( i0:i1, j0:j1+1, 1:lm(n) ), & sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), & sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), & dt_sec, tol_rms, tol_diff, status ) if (status/=0) then write (gol,'("final mass imbalance too large for region ",i2)') n; call goErr call goErr; status=1; return end if !done deallocate(full_pw, full_pu, full_pv, full_dm_dt, bigJslice, bigIslice,& jslice, islice) ! >>> data2 >>> if ( any((/mfu_dat%filled2,mfv_dat%filled2,mfw_dat%filled2/)) ) then 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 ! CHECK INITIAL MASS BALANCE: ! ----------------------------------- ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data2, pu_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data2, pv_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call CheckMassBalance( lli(n), & pu_dat(n)%data2(i0-1:i1, j0:j1 , 1:lm(n) ), & pv_dat(n)%data2( i0:i1, j0:j1+1, 1:lm(n) ), & sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), & sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), & dt_sec, 1.0e-4, 1.0e-3, status ) if (status/=0) then write (gol,'("initial mass imbalance too large for region ",i2)') n; call goErr call goErr; status=1; return end if ! BALANCE HORIZONTAL FLUXES ! ----------------------------------- ! balance horizontal fluxes: ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) if (isRoot) then allocate(full_pu( 0:im(n), 1:jm(n), 0:lm(n)) ) ! must have same number of levels as mfu allocate(full_pv( 1:im(n), 1:jm(n)+1, 0:lm(n)) ) allocate(full_pw( 1:im(n), 1:jm(n), 0:lm(n)) ) ! used also as temp arr in comm allocate(full_dm_dt(im(n), jm(n), lm(n)) ) else allocate( full_pu(1,1,1) ) allocate( full_pv(1,1,1) ) allocate( full_pw(1,1,1) ) allocate( full_dm_dt(1,1,1)) end if !for slice scattering allocate(islice(j0:j1,0:lm(n))) allocate(jslice(i0:i1,0:lm(n))) if (isRoot) then allocate(bigIslice(1:jm(n),0:lm(n))) allocate(bigJslice(1:im(n),0:lm(n))) else allocate(bigIslice(1,1)) allocate(bigJslice(1,1)) end if call GATHER( dgrid(n), pu_dat(n)%data2, full_pw, pu_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then full_pu(1:im(n),1:jm(n),:) = full_pw full_pu(0,:,:) = full_pu(im(n),:,:) ! East-West periodicity end if call GATHER( dgrid(n), pv_dat(n)%data2, full_pw, pv_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then full_pv(1:im(n),1:jm(n),:) = full_pw full_pv(:,jm(n)+1,:) = full_pv(:,1,:) ! donut periodicity end if call GATHER( dgrid(n), dm_dt, full_dm_dt, 0, status ) IF_NOTOK_RETURN(status=1) call GATHER( dgrid(n), pw_dat(n)%data2, full_pw, pw_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if (isRoot) then call BalanceMassFluxes( global_lli(n), & full_pu(0:im(n),1:jm(n) ,1:lm(n)), & full_pv(1:im(n),1:jm(n)+1,1:lm(n)), & full_pw, full_dm_dt, global_lli(parent(n)), dt_sec, status ) IF_NOTOK_RETURN(status=1) end if call SCATTER( dgrid(n), pw_dat(n)%data2, full_pw, pw_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) if(isRoot) full_pw = full_pu(1:im(n),1:jm(n),:) call SCATTER( dgrid(n), pu_dat(n)%data2, full_pw, pu_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! scatter extra column full_pu(0,:,:) - needed only for noncyclic zoom ! regions, for others update_halo takes care of it [FIXME: could had a ! test around these 3 lines ] if(isRoot) bigIslice = full_pu(0,1:jm(n),:) CALL SCATTER_I_BAND( dgrid(n), islice, bigIslice, status, iref=1) if(WestBorder) pu_dat(n)%data2(0,j0:j1,:) = islice if(isRoot) full_pw = full_pv(1:im(n),1:jm(n),:) call SCATTER( dgrid(n), pv_dat(n)%data2, full_pw, pv_dat(n)%halo, status ) IF_NOTOK_RETURN(status=1) ! Scatter PV(:,jm+1,:) if(isroot) bigJslice=full_pv(1:im(n),jm(n)+1,:) CALL SCATTER_J_BAND( dgrid(n), jslice, bigJslice, status, jref=jm(n)) if(NorthBorder)pv_dat(n)%data2(i0:i1,jm(n)+1,:)=jslice ! CHECK FINAL MASS BALANCE: ! ----------------------------------- ! NOTE: strange old indexing: ! pu_tmpp --> pu(0:im(n),1:jm(n) ,1:lm(n)) in pu_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) ! pv_tmpp --> pv(1:im(n),1:jm(n)+1,1:lm(n)) in pv_t(0:im(n)+1,0:jm(n)+1,0:lm(n)) CALL UPDATE_HALO( dgrid(n), pu_dat(n)%data2, pu_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(n), pv_dat(n)%data2, pv_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call CheckMassBalance( lli(n), & pu_dat(n)%data2(i0-1:i1, j0:j1 , 1:lm(n) ), & pv_dat(n)%data2( i0:i1, j0:j1+1, 1:lm(n) ), & sp1_dat(n)%data ( i0:i1, j0:j1 , 1 ), & sp2_dat(n)%data ( i0:i1, j0:j1 , 1 ), & dt_sec, 1.0e-7, 1.0e-6, status ) if (status/=0) then write (gol,'("final mass imbalance too large for region ",i2)') n; call goErr call goErr; status=1; return end if deallocate(full_pw, full_pu, full_pv, full_dm_dt, bigJslice, bigIslice,& jslice, islice) end if ! filled2 ! >>> ! clear deallocate( dm_dt ) end do ! regions #endif /* ADVECTION */ !------------ ! Done !------------ call done(l_lli, status) IF_NOTOK_RETURN(status=1) status = 0 call goLabel() END SUBROUTINE METEO_SETUP_MASS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: METEO_SETUP_OTHER ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_SETUP_OTHER( tr1, tr2, status, isfirst ) ! ! !USES: ! use GO, only : TDate, NewDate, rTotal, wrtgol use GO, only : operator(-), operator(+), operator(/) use GO, only : InterpolFractions use dims, only : nregions, im, jm, lm use dims, only : lmax_conv use dims, only : xcyc use Dims, only : czeta use global_data, only : region_dat #ifndef without_convection use global_data, only : conv_dat #endif use Phys, only : ConvCloudDim ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: tr1, tr2 logical, intent(in), optional :: isfirst ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 9 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Meteo_Setup_Other' logical :: do_isfirst integer :: n, p integer :: i, j, l integer :: lsave, i0, i1, j0, j1 real :: tote, totd, maxe real, pointer :: dxyp(:) type(TDate) :: tmid real :: alfa1, alfa2 integer :: iveg ! --- begin -------------------------------- call goLabel(rname) ! initial call ? if ( present(isfirst) ) then do_isfirst = isfirst else do_isfirst = .false. end if ! ! ** orography ***************************** ! ! read orographies (if necessary): do n = 1, nregions_all #ifdef parallel_cplng call setup( n, oro_dat(n), (/tr1,tr2/), 'n', status ) #else call setup( n, oro_dat(n), (/tr1,tr2/), global_lli(n), 'n', status ) #endif IF_NOTOK_RETURN(status=1) end do ! ! ** spm ************************************** ! do n = 1, nregions ! skip ? if ( .not. spm_dat(n)%used ) cycle ! mid time: tmid = tr1 + ( tr2 - tr1 )/2 ! deterimine weights to sp1 and sp2 : call InterpolFractions( tmid, sp1_dat(n)%tr(1), sp2_dat(n)%tr(1), alfa1, alfa2, status ) IF_NOTOK_RETURN(status=1) call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 ) ! interpolate: spm_dat(n)%data(i0:i1,j0:j1,1) = alfa1 * sp1_dat(n)%data(i0:i1,j0:j1,1) + & alfa2 * sp2_dat(n)%data(i0:i1,j0:j1,1) ! store time: spm_dat(n)%tr = (/tr1,tr2/) end do ! regions ! ! ** omega ************************************** ! do n = 1, nregions_all ! 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 ************************************** ! do n = 1, nregions_all ! ncep meteo requires conversion of virtual temperature using humidity ... if ( (temper_dat(n)%sourcekey(1:4) == 'ncep') .or. (humid_dat(n)%sourcekey(1:4) == 'ncep') ) then ! read temperature and humidity (if necessary): call setup_TQ( n, temper_dat(n), humid_dat(n), (/tr1,tr2/), global_lli(n), levi, status) IF_NOTOK_RETURN(status=1) else ! read temperature (if necessary): #ifdef parallel_cplng call setup( n, temper_dat(n), (/tr1,tr2/), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) ! read humidity (if necessary): call setup( n, humid_dat(n), (/tr1,tr2/), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) #else call setup( n, temper_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) ! read humidity (if necessary): call setup( n, humid_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) #endif end if end do ! regions ! ! ** gph ************************************** ! do n = 1, nregions_all ! re-compute gph from pressure, temperature, and humidity: call compute_gph( n, status ) IF_NOTOK_RETURN(status=1) end do ! ! ** clouds ************************************** ! do n = 1, nregions #ifdef parallel_cplng call setup( n, lwc_dat(n), (/tr1,tr2/), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) call setup( n, iwc_dat(n), (/tr1,tr2/), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) call setup_CloudCovers( n, cc_dat(n), cco_dat(n), ccu_dat(n), (/tr1,tr2/), levi, status) IF_NOTOK_RETURN(status=1) #else call setup( n, lwc_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) call setup( n, iwc_dat(n), (/tr1,tr2/), global_lli(n), 'n', levi, 'n', status) IF_NOTOK_RETURN(status=1) call setup_CloudCovers( n, cc_dat(n), cco_dat(n), ccu_dat(n), (/tr1,tr2/), global_lli(n), levi, status) IF_NOTOK_RETURN(status=1) #endif end do ! ! ** convection ************************************** ! do n = 1, nregions #ifdef parallel_cplng call setup_Convec( n, entu_dat(n), entd_dat(n), detu_dat(n), detd_dat(n), & omega_dat(n), gph_dat(n), (/tr1,tr2/), levi, status ) IF_NOTOK_RETURN(status=1) #else call setup_Convec( n, entu_dat(n), entd_dat(n), detu_dat(n), detd_dat(n), & omega_dat(n), gph_dat(n), (/tr1,tr2/), global_lli(n), levi, status ) IF_NOTOK_RETURN(status=1) #endif end do #ifndef without_convection ! ~~ convective clouds do n = 1, nregions if ( .not. entu_dat(n)%used ) cycle if ( .not. entd_dat(n)%used ) cycle ! update necessary ? if ( any((/entu_dat(n)%changed,entd_dat(n)%changed/)) ) then call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 ) do j = j0, j1 do i = i0, i1 ! compute convective cloud dimensions for this column: call ConvCloudDim( 'u', size(detu_dat(n)%data,3), & detu_dat(n)%data(i,j,:), entd_dat(n)%data(i,j,:),& conv_dat(n)%cloud_base(i,j), & conv_dat(n)%cloud_top (i,j), & conv_dat(n)%cloud_lfs (i,j), & status ) IF_NOTOK_RETURN(status=1) end do end do end if ! changed end do ! regions #endif ! ~~ unit conversion do n = 1, nregions if ( .not. entu_dat(n)%used ) cycle if ( .not. entd_dat(n)%used ) cycle if ( .not. detu_dat(n)%used ) cycle if ( .not. detd_dat(n)%used ) cycle ! update necessary ? if ( any((/ entu_dat(n)%changed, entd_dat(n)%changed, & detu_dat(n)%changed, detd_dat(n)%changed /)) ) then call Get_DistGrid( dgrid(n), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 ) !cmk calculate the rates in kg/gridbox and scale with czeta dxyp => region_dat(n)%dxyp do j = j0, j1 do i = i0, i1 ! kg/m2/s -> kg/gridbox/s * scale_factor entu_dat(n)%data(i,j,:) = entu_dat(n)%data(i,j,:)*dxyp(j)*czeta detu_dat(n)%data(i,j,:) = detu_dat(n)%data(i,j,:)*dxyp(j)*czeta ! ensure netto zero tracer transport by updraught in column ! (add difference between total entrement and detrement ! to level where entrement reaches maximum): tote = sum( entu_dat(n)%data(i,j,:) ) totd = sum( detu_dat(n)%data(i,j,:) ) maxe = entu_dat(n)%data(i,j,1) ! changed: reported by PB feb 2003 lsave = 1 do l = 2, lmax_conv if ( entu_dat(n)%data(i,j,l) > maxe ) then maxe = entu_dat(n)%data(i,j,l) lsave = l end if end do entu_dat(n)%data(i,j,lsave) = entu_dat(n)%data(i,j,lsave) - tote + totd ! kg/m2/s -> kg/gridbox/s * scale_factor entd_dat(n)%data(i,j,:) = entd_dat(n)%data(i,j,:)*dxyp(j)*czeta detd_dat(n)%data(i,j,:) = detd_dat(n)%data(i,j,:)*dxyp(j)*czeta ! ensure netto zero tracer transport by downdraught in column ! (add difference between total entrement and detrement ! to level where entrement reaches maximum): tote = sum( entd_dat(n)%data(i,j,:) ) ! total entrainement totd = sum( detd_dat(n)%data(i,j,:) ) ! total detrainement maxe = 0.0 lsave = lmax_conv do l = 1, lmax_conv if ( entd_dat(n)%data(i,j,l) > maxe ) then maxe = entd_dat(n)%data(i,j,l) lsave = l end if end do entd_dat(n)%data(i,j,lsave) = entd_dat(n)%data(i,j,lsave) - tote + totd end do end do end if ! changed ? end do ! regions ! ! ** surface fields ***************************** ! #ifdef parallel_cplng ! * lsmask call setup( lsmask_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * albedo call setup( albedo_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * sr_ecm call setup( sr_ecm_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) #else ! * lsmask call setup( lsmask_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * albedo call setup( albedo_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * sr_ecm call setup( sr_ecm_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! * sr_ols call setup( sr_ols_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) #ifdef parallel_cplng ! * sea ice call setup( ci_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * sea surface temperature call setup( sst_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * u10m call setup( u10m_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * v10m call setup( v10m_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * windspeed call setup( wspd_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * skin reservoir content call setup( src_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * 2m dewpoint temperature call setup( d2m_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * 2m temperature call setup( t2m_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * slhf call setup( slhf_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * sshf call setup( sshf_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface stress call setup( ewss_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) call setup( nsss_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * convective precipitation call setup( cp_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * large scale stratiform precipitation call setup( lsp_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface solar radiation call setup( ssr_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) call setup( ssrd_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface thermal radiation call setup( str_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) call setup( strd_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * skin temperature call setup( skt_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * boundary layer height #ifndef with_tmm_ecearth call setup( blh_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! * snow fall and depth call setup( sf_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) call setup( sd_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * g10m #ifndef with_tmm_ecearth call setup( g10m_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! * soil water level 1 call setup( swvl1_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * soil temperature level 1 call setup( stl1_dat(1:nregions_all), (/tr1,tr2/), '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(1:nregions_all,iveg), (/tr1,tr2/), '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(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) ! * high vegetation cover call setup( cvh_dat(1:nregions_all), (/tr1,tr2/), 'n', status) IF_NOTOK_RETURN(status=1) #else ! * sea ice call setup( ci_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * sea surface temperature call setup( sst_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * u10m call setup( u10m_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * v10m call setup( v10m_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * windspeed call setup( wspd_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * skin reservoir content call setup( src_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * 2m dewpoint temperature call setup( d2m_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * 2m temperature call setup( t2m_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * slhf call setup( slhf_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * sshf call setup( sshf_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface stress call setup( ewss_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) call setup( nsss_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * convective precipitation call setup( cp_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * large scale stratiform precipitation call setup( lsp_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface solar radiation call setup( ssr_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) call setup( ssrd_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * surface thermal radiation call setup( str_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) call setup( strd_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * skin temperature call setup( skt_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * boundary layer height #ifndef with_tmm_ecearth call setup( blh_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! * snow fall and depth call setup( sf_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) call setup( sd_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * g10m #ifndef with_tmm_ecearth call setup( g10m_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! * soil water level 1 call setup( swvl1_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), '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(1:nregions_all,iveg), (/tr1,tr2/), global_lli(1:nregions_all), '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(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * high vegetation cover call setup( cvh_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) ! * MACC emissions call setup( ch4fire_dat(1:nregions_all), (/tr1,tr2/), global_lli(1:nregions_all), 'n', status) IF_NOTOK_RETURN(status=1) #endif ! ! ** done ******************************************** ! status = 0 call goLabel() END SUBROUTINE METEO_SETUP_OTHER !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SETUPSETUP ! ! !DESCRIPTION: for one met data MD and one time range TR, returns the dates ! at begining and end of the met field interval that ! encompasses TR, and if the data for these dates (%data1 and ! %data2, resp.) must be read or copied. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUPSETUP( md, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) ! ! !USES: ! use GO, only : TDate, NewDate, IncrDate, AnyDate, IsAnyDate, Get, Set, wrtgol use GO, only : rTotal, iTotal use GO, only : operator(+), operator(-), operator(/) use GO, only : operator(==), operator(/=), operator(<), operator(<=) use meteodata, only : TMeteoData use global_data, only : fcmode, tfcday0 ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: tr(2) ! ! !OUTPUT PARAMETERS: ! logical, intent(out) :: data1_read, data1_copy type(TDate), intent(out) :: data1_tref, data1_t1, data1_t2 logical, intent(out) :: data2_read, data2_copy type(TDate), intent(out) :: data2_tref, data2_t1, data2_t2 integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Mar 2010 - P. Le Sager - ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/SetupSetup' integer :: dth, baseh integer :: year, month, day, hour, minu type(TDate) :: tmid type(TDate) :: tc(2) integer :: dth_int type(TDate) :: tprev, tnext real :: dhr ! --- begin ----------------------------- 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_SERIAL ! ! !DESCRIPTION: Fill md%data1 and md%data2 of a 2D met field type (md), with ! data for date tr(1) and tr(2) respectively (and if needed) ! through reading or copying. Also write to disk the met field ! if requested. ! ! Then set md%data according to its type of interpolation (see ! TimeInterpolation in meteodata.F90). ! For constant type, %data => %data1. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_2D_SERIAL( region, md, tr, lli, nuv, status ) ! ! !USES: ! use GO, only : TDate, wrtgol use Grid, only : TllGridInfo use TMM, only : ReadField, Read_SP, Read_SR_OLS, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md ! met field ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number type(TDate), intent(in) :: tr(2) ! dates type(TllGridInfo), intent(in) :: lli ! grid (GLOBAL) character(len=1), intent(in) :: nuv ! staggering ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! return code ! ! !REVISION HISTORY: ! 4 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_2d_serial' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, pointer :: field(:,:) ! work array ! --- begin ----------------------------- call goLabel(rname) ! leave if not in use: if ( .not. md%used ) then call goLabel() status=0; return end if if (okdebug) then write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(md%name),trim(lli%name); call goPr endif ! not changed by default md%changed = .false. !------------------ ! time stuff !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read or copy call SetupSetup( md, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! Read/write primary field ! ------------------------- if ( data1_read ) then ! Switch to global if ( md%ls(1) /= md%ls(2) ) then write (gol,'("SETUP_2D called instead of SETUP_3D, field is 3D:")'); call goErr write (gol, '(" md%ls(1:2) : ",2i3)') md%ls; call goErr status=1; IF_NOTOK_RETURN(status=1) end if ! Need whole region for I/O on root. Dummy else. IF (isRoot) THEN ALLOCATE( field( im(region), jm(region)) ) ELSE ALLOCATE( field(1,1) ) END IF ! Read/write IOroot : IF (isRoot) THEN select case ( md%name ) case ( 'sp', 'sps' ) ! special routine for surface pressure call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) case ( 'srols' ) ! special routine for Olsson surface roughness: call Read_SR_OLS( tmmd, md%sourcekey, & data1_tref, data1_t1, data1_t2, & lli, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) case default ! general field call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, lli, & nuv, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) end select ! write meteofiles if ( md%putout ) then call WriteField( tmmd, md%destkey, & md%tmi1, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli, nuv, FIELD, status ) IF_NOTOK_RETURN(status=1) end if END IF IOroot CALL SCATTER( dgrid(region), md%data1(:,:,1), FIELD, md%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( FIELD ) ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: md%data1 = md%data2 ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. end if ! ------------------------- ! Read/write (or copy or nothing) secondary field ! ------------------------- if ( data2_read ) then ! Need whole region for I/O on root. Dummy else. IF (isRoot) THEN ALLOCATE( field( im(region), jm(region)) ) ELSE ALLOCATE( field(1,1) ) END IF ! Read/write IOroot2: IF (isRoot) THEN select case ( md%name ) case ( 'sp', 'sps' ) ! special routine for surface pressure call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli, FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) case default ! general field call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, lli, & nuv, FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) end select ! write meteo files if ( md%putout ) then call WriteField( tmmd, md%destkey, & md%tmi2, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli, nuv, FIELD, status ) IF_NOTOK_RETURN(status=1) end if END IF IOroot2 CALL SCATTER( dgrid(region), md%data2(:,:,1), FIELD, md%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( FIELD ) ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data from secondary array md%data2 = md%data1 ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 end if ! ------------------------- ! time interpolation ! ------------------------- call TimeInterpolation( md, tr, status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! done ! ------------------------- status = 0 call goLabel() END SUBROUTINE SETUP_2D_SERIAL !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SETUP_2D_PARA ! ! !DESCRIPTION: Same as SETUP_2D_SERIAL, except all processes get field from IFS !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_2D_PARA( region, md, tr, nuv, status ) ! ! !USES: ! use GO, only : TDate, wrtgol use Grid, only : TllGridInfo use TMM, only : ReadField, Read_SP, Read_SR_OLS, WriteField ! use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md ! met field ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number type(TDate), intent(in) :: tr(2) ! dates character(len=1), intent(in) :: nuv ! staggering ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! return code ! ! !REVISION HISTORY: ! 18 Oct 2013 - Ph. Le Sager - v0 ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_2d_para' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 integer :: i1, i2, j1, j2 real, pointer :: field(:,:) ! work array ! --- begin ----------------------------- 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 !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read or copy call SetupSetup( md, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! Read/write primary field ! ------------------------- if ( data1_read ) then ! test if ( md%ls(1) /= md%ls(2) ) then write (gol,'("SETUP_2D called instead of SETUP_3D, field is 3D:")'); call goErr write (gol, '(" md%ls(1:2) : ",2i3)') md%ls; call goErr status=1; IF_NOTOK_RETURN(status=1) end if ! could get those bounds from md% directly call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( field( i1:i2, j1:j2) ) !! bonds are not strictly required, could as well do (i2-i1+1, ..) ! Read/write select case ( md%name ) case ( 'sp', 'sps' ) ! special routine for surface pressure call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli(region), FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) case ( 'srols' ) ! special routine for Olsson surface roughness: call Read_SR_OLS( tmmd, md%sourcekey, & data1_tref, data1_t1, data1_t2, & lli(region), FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) case default ! general field call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, lli(region), & nuv, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) end select md%data1(i1:i2, j1:j2, 1) = field deallocate( field ) ! write meteofiles if ( md%putout ) then write(gol,*)"writing of remapped met field not tested yet.. SKIPPED." ; call goErr TRACEBACK; status=1; return end if ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: md%data1 = md%data2 ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. end if ! ------------------------- ! Read/write (or copy or nothing) secondary field ! ------------------------- if ( data2_read ) then call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( field( i1:i2, j1:j2) ) select case ( md%name ) case ( 'sp', 'sps' ) ! special routine for surface pressure call Read_SP( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli(region), FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) case default ! general field call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, lli(region), & nuv, FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) end select md%data2(i1:i2, j1:j2, 1) = FIELD deallocate( field ) ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data from secondary array md%data2 = md%data1 ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 end if ! ------------------------- ! time interpolation ! ------------------------- call TimeInterpolation( md, tr, status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! done ! ------------------------- status = 0 call goLabel() END SUBROUTINE SETUP_2D_PARA !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SETUP_2D_N_SERIAL ! ! !DESCRIPTION: wrapper around setup_2d to process several regions for the ! same field. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_2D_N_SERIAL( md, tr, lli, nuv, status ) ! ! !USES: ! use GO , only : TDate use Grid , only : TllGridInfo use meteodata, only : TMeteoData ! ! !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: ! 6 Apr 2010 - P. Le Sager - ! ! !REMARKS: ! (1) Attention: we assume that the regions list start from #1. ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_2d_n_serial' integer :: n ! --- begin ----------------------------- ! check ... if ( size(md) /= size(lli) ) then write (gol,'("md and lli arrays should have same size:")' ) ; call goErr write (gol,'(" size md : ",i6)' ) size(md) ; call goErr write (gol,'(" size lli : ",i6)' ) size(lli) ; call goErr TRACEBACK; status=1; return end if ! loop over regions: do n = 1, size(md) if (okdebug) then write (gol,'(" ",a," ",a)') trim(md(n)%name), trim(lli(n)%name); call goPr endif call Setup( n, md(n), tr, lli(n), nuv, status ) IF_NOTOK_RETURN(status=1) end do status = 0 END SUBROUTINE SETUP_2D_N_SERIAL !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SETUP_2D_N_PARA ! ! !DESCRIPTION: wrapper around setup_2d to process several regions for the ! same field. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_2D_N_PARA( md, tr, nuv, status ) ! ! !USES: ! use GO , only : TDate use Grid , only : TllGridInfo use meteodata, only : TMeteoData ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md(:) ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: tr(2) character(len=1), intent(in) :: nuv ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_2d_n_para' integer :: n do n = 1, size(md) if (okdebug) then write (gol,'(" ",a," ",a)') trim(md(n)%name), trim(lli(n)%name); call goPr endif call Setup( n, md(n), tr, nuv, status ) IF_NOTOK_RETURN(status=1) end do status = 0 END SUBROUTINE SETUP_2D_N_PARA !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: SETUP_3D ! ! !DESCRIPTION: same as SETUP_2D, but for 3D fields by accounting for levels !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_3D_SERIAL( region, md, tr, lli, nuv, levi, nw, status ) ! ! !USES: ! use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, ReadField, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md ! met field ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number type(TDate), intent(in) :: tr(2) ! dates type(TllGridInfo), intent(in) :: lli ! grid character(len=1), intent(in) :: nuv ! horiz. staggering type(TLevelInfo), intent(in) :: levi ! levels character(len=1), intent(in) :: nw ! vertical staggering ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! return code ! ! !REVISION HISTORY: ! 4 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_3d_serial' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) real, pointer :: field(:,:,:) ! work array (data) integer :: is(2), js(2) ! work arrays (bounds) ! --- begin ----------------------------- call goLabel(rname) ! leave if not in use: if ( .not. md%used ) then call goLabel() status=0; return end if if (okdebug) then write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(md%name),trim(lli%name); call goPr endif ! not changed by default md%changed = .false. !------------------ ! time stuff !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read or copy call SetupSetup( md, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! Read/write primary field ! ------------------------- if ( data1_read ) then ! Need whole region for I/O on root. Dummy else. Allocate global array for I/O is = (/1,im(region)/) js = (/1,jm(region)/) IF (isRoot) THEN ALLOCATE( FIELD( is(1):is(2), js(1):js(2), md%ls(1):md%ls(2) )) ELSE ALLOCATE( FIELD(1,1,1) ) END IF ! Read/write on root IOroot : IF (isRoot) THEN ! safety check if ( data1_t2 /= data1_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure allocate( tmp_sp( is(1):is(2), js(1):js(2) ) ) ! fill data call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli, nuv, levi, nw, & tmp_sp, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) ! write meteo file if ( md%putout ) then call WriteField( tmmd, md%destkey, & md%tmi1, 'sp', trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli, nuv, levi, nw, & tmp_sp, FIELD, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) END IF IOroot CALL SCATTER( dgrid(region), md%data1, FIELD, md%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( FIELD ) ! data array is filled now md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: md%data1 = md%data2 ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) IF (isRoot) THEN ALLOCATE(field(im(region), jm(region), md%ls(1):md%ls(2))) ELSE ALLOCATE(field(1,1,1)) END IF ! Read/write IOroot2 : IF (isRoot) THEN ! safety check ... if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) ! fill data call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli, nuv, levi, nw, & tmp_sp, FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) ! write meteofiles if ( md%putout ) then call WriteField( tmmd, md%destkey, & md%tmi2, 'sp', trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli, nuv, levi, nw, & tmp_sp, FIELD, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) END IF IOroot2 CALL SCATTER( dgrid(region), md%data2, FIELD, md%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( FIELD ) ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data from secondary array md%data2 = md%data1 ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 end if ! ------------------------- ! time interpolation ! ------------------------- call TimeInterpolation( md, tr, status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! done ! ------------------------- status = 0 call goLabel() END SUBROUTINE SETUP_3D_SERIAL !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: SETUP_3D_PARA ! ! !DESCRIPTION: same as SETUP_3D_SERIAL, except reading is done by every processes. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_3D_PARA( region, md, tr, nuv, levi, nw, status ) ! ! !USES: ! use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, ReadField, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: md ! met field ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number type(TDate), intent(in) :: tr(2) ! dates character(len=1), intent(in) :: nuv ! horiz. staggering type(TLevelInfo), intent(in) :: levi ! levels character(len=1), intent(in) :: nw ! vertical staggering ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! return code ! ! !REVISION HISTORY: ! 18 Oct 2013 - Ph. Le Sager - v0 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_3d_para' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 integer :: i1, i2, j1, j2 real, allocatable :: tmp_sp(:,:) real, pointer :: field(:,:,:) ! work array ! --- begin ----------------------------- 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 !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read or copy call SetupSetup( md, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! Read/write primary field ! ------------------------- if ( data1_read ) then ! could get those bounds from md% directly call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( field( i1:i2, j1:j2, md%ls(1):md%ls(2))) ! safety check if ( data1_t2 /= data1_t1 ) then ! write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr ! call wrtgol( ' data1_t1 : ', data1_t1 ) ; call goErr ! call wrtgol( ' data1_t2 : ', data1_t2 ) ; call goErr ! write (gol,'("please decide what to do with surface pressures ... ")') ; call goErr ! TRACEBACK; status=1; return write (gol,'("WARNING - using instant surface pressure for regridding temporal averaged 3D field ...")'); call goPr end if ! surface pressure allocate( tmp_sp( i1:i2, j1:j2 ) ) ! read data call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data1_tref, data1_t1, data1_t2, & lli(region), nuv, levi, nw, & tmp_sp, FIELD, md%tmi1, status ) IF_NOTOK_RETURN(status=1) md%data1(i1:i2, j1:j2, md%ls(1):md%ls(2)) = field ! write meteo file if ( md%putout ) then write(gol,*)"writing of remapped met field not finished yet.. Sorry." ; call goErr TRACEBACK; status=1; return endif DEALLOCATE( TMP_SP ) DEALLOCATE( FIELD ) ! data array is filled now md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: md%data1 = md%data2 ! data array is filled now: md%filled1 = .true. md%tr1(1) = data1_t1 md%tr1(2) = data1_t2 md%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! could get those bounds from md% directly call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( field( i1:i2, j1:j2, md%ls(1):md%ls(2))) ! safety check ... if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr call wrtgol( ' data2_t1 : ', data2_t1 ) ; call goErr call wrtgol( ' data2_t2 : ', data2_t2 ) ; call goErr write (gol,'("please decide what to do with surface pressures ... ")') ; call goErr TRACEBACK; status=1; return end if ! surface pressure allocate( tmp_sp( i1:i2, j1:j2 ) ) ! read data call ReadField( tmmd, md%sourcekey, trim(md%name), trim(md%unit), & data2_tref, data2_t1, data2_t2, & lli(region), nuv, levi, nw, & tmp_sp, FIELD, md%tmi2, status ) IF_NOTOK_RETURN(status=1) md%data2(i1:i2, j1:j2, md%ls(1):md%ls(2)) = field ! write meteofiles if ( md%putout ) then write(gol,*)"writing of remapped met field not finished yet.. Sorry. SKIPPED." ; call goErr TRACEBACK; status=1; return end if ! clear DEALLOCATE( TMP_SP ) DEALLOCATE( FIELD ) ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data from secondary array md%data2 = md%data1 ! data array is filled now md%filled2 = .true. md%tr2(1) = data2_t1 md%tr2(2) = data2_t2 end if ! ------------------------- ! time interpolation ! ------------------------- call TimeInterpolation( md, tr, status ) IF_NOTOK_RETURN(status=1) ! ------------------------- ! done ! ------------------------- status = 0 call goLabel() END SUBROUTINE SETUP_3D_PARA !EOC ! ************************************************************** ! *** ! *** Specific SETUP routines for MASS FLUXES - Only serial case ! *** since it reads spectral fields from IFS ! *** ! ************************************************************** SUBROUTINE SETUP_UVW( region, md_mfu, md_mfv, md_mfw, md_tsp, tr, lli, nuv, levi, nw, status ) ! Set up MFU and MFV (horizontal fluxes) ! Set up MFW (vertical flux) and TSP (tendency surface pressure) ! Read or copy %data1 and %data2, and get %data through time interpolation use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, TMM_Read_UVW, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! --- in/out ---------------------------------- integer, intent(in) :: region ! region number type(TMeteoData), intent(inout) :: md_mfu type(TMeteoData), intent(inout) :: md_mfv type(TMeteoData), intent(inout) :: md_mfw type(TMeteoData), intent(inout) :: md_tsp type(TDate), intent(in) :: tr(2) ! time range 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_UVW' ! --- local ---------------------------------- logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 logical :: NorthBorder, WestBorder ! tile location real, allocatable :: tmp_spu(:,:) real, allocatable :: tmp_spv(:,:) real, allocatable :: tmp_sp(:,:) ! to read the entire region real, pointer :: wrld_u(:,:,:), wrld_v(:,:,:), wrkarr(:,:,:) real, pointer :: mfw(:,:,:), tsp(:,:) ! work arrays (data) integer, dimension(2) :: is, js, ls integer :: halo, i1, i2, j1, j2 real, allocatable :: bigIslice(:,:), bigJslice(:,:), Islice(:,:), Jslice(:,:) ! --- 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 if (okdebug) then write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(md_mfu%name),trim(lli%name); call goPr write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(md_mfv%name),trim(lli%name); call goPr write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(md_mfw%name),trim(lli%name); call goPr endif ! not changed by default md_mfu%changed = .false. md_mfv%changed = .false. md_mfw%changed = .false. md_tsp%changed = .false. ! local indices and tile location CALL GET_DISTGRID( dgrid(region), & I_STRT=i1, I_STOP=i2, & J_STRT=j1, J_STOP=j2, & hasWestBorder=WestBorder, hasNorthBorder=NorthBorder) !------------------ ! time stuff !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from mfu only) call SetupSetup( md_mfu, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then ! Use fact that mfu and mfv have been allocated with the same bounds and halo ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) ls = md_mfu%ls halo = md_mfu%halo IF (isRoot) THEN ALLOCATE( wrld_u( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) ) ALLOCATE( wrld_v( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) ) ALLOCATE( wrkarr( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) wrld_v = 0. wrld_u = 0. allocate( bigIslice(jm(region),ls(1):ls(2))) allocate( bigJslice(im(region),ls(1):ls(2))) allocate( mfw(is(1):is(2), js(1):js(2), ls(1):ls(2) )) allocate( tsp(is(1):is(2), js(1):js(2)) ) ELSE ALLOCATE(wrld_u(1,1,1), wrld_v(1,1,1), wrkarr(1,1,1)) ALLOCATE( bigIslice(1,1), bigJslice(1,1) ) allocate( mfw(1,1,1), tsp(1,1) ) END IF ALLOCATE( Islice(j1:j2,ls(1):ls(2)) ) ALLOCATE( Jslice(i1:i2,ls(1):ls(2)) ) if (isRoot) then ! only root does IO ! safety check ... if ( data1_t2 /= data1_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_spu(is(1)-1:is(2), js(1):js(2) ) ) allocate( tmp_spv( is(1):is(2), js(1):js(2)+1) ) allocate( tmp_sp ( is(1):is(2), js(1):js(2) ) ) ! 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 TMM_READ_UVW( tmmd, md_mfu%sourcekey, & data1_tref, data1_t1, data1_t2, lli, levi, & tmp_spu, & wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), & md_mfu%tmi1, & tmp_spv, & wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), & md_mfv%tmi1, & tmp_sp, mfw, & tsp, & md_mfw%tmi1, status ) IF_NOTOK_RETURN(status=1) ! write meteofiles if ( md_mfu%putout ) then call WriteField( tmmd, md_mfu%destkey, & md_mfu%tmi1, 'spu', trim(md_mfu%name), trim(md_mfu%unit), & data1_tref, data1_t1, data1_t2, & lli, 'u', levi, 'n', & tmp_spu, wrld_u(is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), & status ) IF_NOTOK_RETURN(status=1) end if if ( md_mfv%putout ) then call WriteField( tmmd, md_mfv%destkey, & md_mfv%tmi1, 'spv', trim(md_mfv%name), trim(md_mfv%unit), & data1_tref, data1_t1, data1_t2, & lli, 'v', levi, 'n', & tmp_spv, wrld_v(is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), & status ) IF_NOTOK_RETURN(status=1) end if if ( md_mfw%putout ) then call WriteField( tmmd, md_mfw%destkey, & md_mfw%tmi1, 'sp', trim(md_mfw%name), trim(md_mfw%unit), & data1_tref, data1_t1, data1_t2, & lli, nuv, levi, nw, & tmp_sp, mfw, status ) IF_NOTOK_RETURN(status=1) end if if ( md_tsp%putout ) then ! use history from mfw ... call WriteField( tmmd, md_tsp%destkey, & md_mfw%tmi1, trim(md_tsp%name), trim(md_tsp%unit), & data1_tref, data1_t1, data1_t2, & lli, nuv, tsp, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_spu ) deallocate( tmp_spv ) deallocate( tmp_sp ) end if ! root ? ! Scatter U if(isRoot) wrkarr = wrld_u(is(1):is(2),js(1):js(2),:) CALL SCATTER( dgrid(region), md_mfu%data1, wrkarr, md_mfu%halo, status) IF_NOTOK_RETURN(status=1) ! manually scatter wrld_u(is(1)-1,:,:). This is needed only with non-cyclic ! zoom regions, since any update_halo will overwrite is(1)-1. [FIXME: could had a ! test around these 3 lines ] !PLS if(isRoot) bigIslice = wrld_u(0,js(1):js(2),:) !PLS CALL SCATTER_I_BAND( dgrid(region), islice, bigIslice, status, iref=1) !PLS if (WestBorder) md_mfu%data1(0,j1:j2,:) = islice ! Scatter V if(isRoot) wrkarr = wrld_v(is(1):is(2),js(1):js(2),:) CALL SCATTER( dgrid(region), md_mfv%data1, wrkarr, md_mfv%halo, status) IF_NOTOK_RETURN(status=1) ! manually SCATTER wrld_v( :, js(2)+1 , :) : NORTH POLE HALO if(isroot) bigJslice=wrld_v(is(1):is(2),jm(region)+1,:) CALL SCATTER_J_BAND( dgrid(region), jslice, bigJslice, status, jref=jm(region)) if (NorthBorder) md_mfv%data1(i1:i2,jm(region)+1,:)=jslice deallocate(wrld_u, wrld_v, wrkarr, bigIslice, bigJslice, Islice, Jslice) ! Scatter W CALL SCATTER( dgrid(region), md_mfw%data1, MFW, md_mfw%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), md_tsp%data1(:,:,1), TSP, md_tsp%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(MFW, TSP) ! data array is filled now md_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. 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_mfu%data1 = md_mfu%data2 md_mfv%data1 = md_mfv%data2 md_mfw%data1 = md_mfw%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. md_mfw%filled1 = .true. md_mfw%tr1(1) = data1_t1 md_mfw%tr1(2) = data1_t2 md_mfw%changed = .true. md_tsp%filled1 = .true. md_tsp%tr1(1) = data1_t1 md_tsp%tr1(2) = data1_t2 md_tsp%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) ls = md_mfu%ls halo = md_mfu%halo IF (isRoot) THEN allocate( wrld_u( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) ) allocate( wrld_v( is(1)-halo:is(2)+halo, js(1)-halo:js(2)+halo, ls(1):ls(2)) ) allocate( wrkarr( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) wrld_v = 0. wrld_u = 0. allocate( bigIslice(jm(region),ls(1):ls(2))) allocate( bigJslice(im(region),ls(1):ls(2))) allocate( mfw(is(1):is(2), js(1):js(2), ls(1):ls(2) )) allocate( tsp(is(1):is(2), js(1):js(2)) ) ELSE ALLOCATE(wrld_u(1,1,1), wrld_v(1,1,1), wrkarr(1,1,1)) ALLOCATE( bigIslice(1,1), bigJslice(1,1) ) allocate( mfw(1,1,1), tsp(1,1) ) END IF ALLOCATE( Islice(j1:j2,ls(1):ls(2)) ) ALLOCATE( Jslice(i1:i2,ls(1):ls(2)) ) if (isRoot) then ! only root does IO if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_spu(is(1)-1:is(2), js(1):js(2) ) ) allocate( tmp_spv( is(1):is(2), js(1):js(2)+1) ) allocate( tmp_sp ( is(1):is(2), js(1):js(2) ) ) ! 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 TMM_READ_UVW( tmmd, md_mfu%sourcekey, & data2_tref, data2_t1, data2_t2, lli, levi, & tmp_spu, & wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), & md_mfu%tmi2, & tmp_spv, & wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), & md_mfv%tmi2, & tmp_sp, MFW, TSP, md_mfw%tmi2, status ) IF_NOTOK_RETURN(status=1) ! write meteofiles if ( md_mfu%putout ) then call WriteField( tmmd, md_mfu%destkey, & md_mfu%tmi2, 'spu', trim(md_mfu%name), trim(md_mfu%unit), & data2_tref, data2_t1, data2_t2, & lli, 'u', levi, 'n', & tmp_spu, wrld_u( is(1)-1:is(2), js(1):js(2), ls(1)+1:ls(2) ), & status ) IF_NOTOK_RETURN(status=1) endif if ( md_mfv%putout ) then call WriteField( tmmd, md_mfv%destkey, & md_mfv%tmi2, 'spv', trim(md_mfv%name), trim(md_mfv%unit), & data2_tref, data2_t1, data2_t2, & lli, 'v', levi, 'n', & tmp_spv, wrld_v( is(1):is(2), js(1):js(2)+1, ls(1)+1:ls(2) ), & status ) IF_NOTOK_RETURN(status=1) end if if ( md_mfw%putout ) then call WriteField( tmmd, md_mfw%destkey, & md_mfw%tmi2, 'sp', trim(md_mfw%name), trim(md_mfw%unit), & data2_tref, data2_t1, data2_t2, & lli, nuv, levi, nw, & tmp_sp, MFW, status ) IF_NOTOK_RETURN(status=1) end if if ( md_tsp%putout ) then ! use history from mfw ... call WriteField( tmmd, md_tsp%destkey, & md_mfw%tmi2, trim(md_tsp%name), trim(md_tsp%unit), & data2_tref, data2_t1, data2_t2, & lli, nuv, TSP, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_spu ) deallocate( tmp_spv ) deallocate( tmp_sp ) end if ! root ! Scatter U if(isRoot) wrkarr = wrld_u(is(1):is(2),js(1):js(2),:) CALL SCATTER( dgrid(region), md_mfu%data2, wrkarr, md_mfu%halo, status) IF_NOTOK_RETURN(status=1) ! important for zoom regions only, since any update_halo will overwrite is(1)-1. [FIXME: could had a ! test around these 3 lines ] !PLS if(isRoot) bigIslice = wrld_u(0,js(1):js(2),:) !PLS CALL SCATTER_I_BAND( dgrid(region), islice, bigIslice, status, iref=1) !PLS if (WestBorder) md_mfu%data2(0,j1:j2,:) = islice ! Scatter V if(isRoot) wrkarr = wrld_v(is(1):is(2),js(1):js(2),:) CALL SCATTER( dgrid(region), md_mfv%data2, wrkarr, md_mfv%halo, status) IF_NOTOK_RETURN(status=1) ! manually SCATTER wrld_v( :, js(2)+1 , :) : NORTH POLE HALO if(isroot) bigJslice=wrld_v(is(1):is(2),jm(region)+1,:) CALL SCATTER_J_BAND( dgrid(region), jslice, bigJslice, status, jref=jm(region)) if (NorthBorder) md_mfv%data2(i1:i2,jm(region)+1,:)=jslice DEALLOCATE(wrld_u, wrld_v, wrkarr, bigIslice, bigJslice, Islice, Jslice) ! Scatter W CALL SCATTER( dgrid(region), md_mfw%data2, MFW, md_mfw%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), md_tsp%data2(:,:,1), TSP, md_tsp%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(MFW, TSP) ! data array is filled now: md_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 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 primary array: md_mfu%data2 = md_mfu%data md_mfv%data2 = md_mfv%data md_mfw%data2 = md_mfw%data1 ! 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 md_mfw%filled2 = .true. md_mfw%tr2(1) = data2_t1 md_mfw%tr2(2) = data2_t2 md_tsp%filled2 = .true. md_tsp%tr2(1) = data2_t1 md_tsp%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( md_mfu, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( md_mfv, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( md_mfw, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( md_tsp, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ status = 0 call goLabel() end subroutine SETUP_UVW ! ************************************************************** ! *** ! *** temperature and humidity ! *** ! ************************************************************** subroutine Setup_TQ( region, md_T, md_Q, tr, lli, levi, status) use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, Read_TQ, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! --- in/out ---------------------------------- integer, intent(in) :: region ! region number type(TMeteoData), intent(inout) :: md_T type(TMeteoData), intent(inout) :: md_Q type(TDate), intent(in) :: tr(2) type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/Setup_TQ' ! --- local ---------------------------------- logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) real, pointer :: T(:,:,:), Q(:,:,:) ! work array integer :: is(2), js(2) ! work arrays (bounds) ! --- begin ----------------------------- 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 !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from T only) call SetupSetup( md_T, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) IF (isRoot) THEN ALLOCATE( T(is(1):is(2), js(1):js(2), md_T%ls(1):md_T%ls(2) )) ALLOCATE( Q(is(1):is(2), js(1):js(2), md_Q%ls(1):md_Q%ls(2) )) ELSE ALLOCATE( T(1,1,1), Q(1,1,1) ) END IF if (isRoot) then ! only root does IO ! safety check ... if ( data1_t2 /= data1_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_sp( is(1):is(2), js(1):js(2) ) ) ! fill data: call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, & data1_tref, data1_t1, data1_t2, lli, levi, & tmp_sp, & T, md_T%tmi1, & Q, md_Q%tmi1, status ) IF_NOTOK_RETURN(status=1) ! write meteofiles ? if ( md_T%putout ) then call WriteField( tmmd, md_T%destkey, & md_T%tmi1, 'sp', trim(md_T%name), trim(md_T%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, 'n', & tmp_sp, T, status ) IF_NOTOK_RETURN(status=1) end if if ( md_Q%putout ) then call WriteField( tmmd, md_Q%destkey, & md_Q%tmi1, 'sp', trim(md_Q%name), trim(md_Q%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, 'n', & tmp_sp, Q, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if ! root ? ! Distribute CALL SCATTER( dgrid(region), md_T%data1, T, md_T%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), md_Q%data1, Q, md_Q%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(T, Q) ! data array is filled now: md_T%filled1 = .true. md_T%tr1(1) = data1_t1 md_T%tr1(2) = data1_t2 md_T%changed = .true. md_Q%filled1 = .true. md_Q%tr1(1) = data1_t1 md_Q%tr1(2) = data1_t2 md_Q%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: md_T%data1 = md_T%data2 md_Q%data1 = md_Q%data2 ! data array is filled now: md_T%filled1 = .true. md_T%tr1(1) = data1_t1 md_T%tr1(2) = data1_t2 md_T%changed = .true. md_Q%filled1 = .true. md_Q%tr1(1) = data1_t1 md_Q%tr1(2) = data1_t2 md_Q%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) IF (isRoot) THEN allocate( T(is(1):is(2), js(1):js(2), md_T%ls(1):md_T%ls(2) )) allocate( Q(is(1):is(2), js(1):js(2), md_Q%ls(1):md_Q%ls(2) )) ELSE allocate( T(1,1,1), Q(1,1,1) ) END IF if (isRoot) then ! only root does IO ! safety check ... if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_sp( is(1):is(2), js(1):js(2)) ) ! fill data: call Read_TQ( tmmd, md_T%sourcekey, md_Q%sourcekey, & data2_tref, data2_t1, data2_t2, lli, levi, & tmp_sp, & T, md_T%tmi2, & Q, md_Q%tmi2, status ) IF_NOTOK_RETURN(status=1) ! write meteofiles ? if ( md_T%putout ) then call WriteField( tmmd, md_T%destkey, & md_T%tmi2, 'sp', trim(md_T%name), trim(md_T%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, 'n', & tmp_sp, T, status ) IF_NOTOK_RETURN(status=1) endif if ( md_Q%putout ) then call WriteField( tmmd, md_Q%destkey, & md_Q%tmi2, 'sp', trim(md_Q%name), trim(md_Q%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, 'n', & tmp_sp, Q, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if ! root CALL SCATTER( dgrid(region), md_T%data2, T, md_T%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), md_Q%data2, Q, md_Q%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(T, Q) ! data array is filled now: md_T%filled2 = .true. md_T%tr2(1) = data2_t1 md_T%tr2(2) = data2_t2 md_Q%filled2 = .true. md_Q%tr2(1) = data2_t1 md_Q%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data from primary array: md_T%data2 = md_T%data1 md_Q%data2 = md_Q%data1 ! data array is filled now: md_T%filled2 = .true. md_T%tr2(1) = data2_t1 md_T%tr2(2) = data2_t2 md_Q%filled2 = .true. md_Q%tr2(1) = data2_t1 md_Q%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( md_T, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( md_Q, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ status = 0 call goLabel() end subroutine Setup_TQ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: METEO_CHECKPRESSURE ! ! !DESCRIPTION: Compute difference b/w sp1_dat (read) and sp_dat (advected), ! and compare to threshold. !\\ !\\ ! !INTERFACE: ! SUBROUTINE METEO_CHECKPRESSURE( n, status ) ! ! !USES: ! use ParTools, only : Par_Reduce use dims, only : idate, newsrun use dims, only : xcyc, im, jm use redgridZoom, only : calc_pdiff #ifdef with_hdf4 use io_hdf, only : io_write2d_32d, DFACC_CREATE #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: n ! region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Meteo_CheckPressure' ! maximum accepted pressure difference: real, parameter :: pdiffmax_treshold = 1.0e2 ! Pa ! --- external ------------------------- integer(4), external :: sfStart, sfEnd ! --- local ----------------------------- real :: pdiffmax, pdiffmax_l integer(4) :: io ! --- begin ------------------------------ call goLabel(rname) ! compare 'advected' pressure with read pressure if ( .not. newsrun ) then ! compute difference between 'advected' pressure sp and read pressure ! sp1, accounting for reduce grid if any call calc_pdiff( n, pdiffmax_l ) ! compute maximum over all pe's call Par_Reduce( pdiffmax_l, 'max', pdiffmax, status, all=.true. ) IF_NOTOK_RETURN(status=1) ! check ... if ( pdiffmax > pdiffmax_treshold ) then write (gol,'("difference between advected and read-in pressure exceeds treshold :")'); call goErr write (gol,'(" max diff. : ",es9.2," [Pa]")') pdiffmax; call goErr write (gol,'(" treshold : ",es9.2," [Pa]")') pdiffmax_treshold; call goErr write (gol,'("pressure arrays saved to local `pressure.hdf`")'); call goErr #ifdef with_hdf4 if (isRoot) then io = sfStart( 'pressure.hdf', DFACC_CREATE ) if ( io > 0 ) then call io_write2d_32d( io, im(n)+4, 'LON', jm(n)+4, 'LAT', sp1_dat(n)%data(:,:,1), 'p' , idate ) call io_write2d_32d( io, im(n)+4, 'LON', jm(n)+4, 'LAT', sp_dat(n)%data(:,:,1), 'pold', idate ) status = sfend(io) else write (gol,'("writing pressures")'); call goErr end if end if ! root #endif TRACEBACK; status=1; return end if ! max diff end if ! no newsrun ! ok status = 0 call goLabel() END SUBROUTINE METEO_CHECKPRESSURE !EOC ! ************************************************************** ! *** ! *** vertical velocity ! *** ! ************************************************************** subroutine Compute_Omega( omega, lli, mfw, status ) use binas , only : grav use grid, only : TllGridInfo, AreaOper use meteodata, only : TMeteoData use tmm, only : SetHistory, AddHistory ! --- in/out ---------------------------------- type(TMeteoData), intent(inout) :: omega ! Pa/s downward type(TllGridInfo), intent(in) :: lli type(TMeteoData), intent(in) :: mfw ! kg/s upward integer, intent(out) :: status ! --- const ----------------------------------- character(len=*), parameter :: rname = mname//'/Compute_Omega' ! --- local ---------------------------------- integer :: l ! --- begin ---------------------------------- ! not in use ? if ( .not. omega%used ) return ! leave if not in use: if ( .not. mfw%used ) then write (gol,'("omega (Pa/s) requires mfw (kg/s)")'); call goErr call goErr; status=1; return end if call goLabel(rname) ! Pa/s = kg/s / m2 * g ! init with mass flux; revert sign from upward to downard, divide by ! gravity accelaration omega%data = - mfw%data * grav ! Pa/s m2 ! loop over levels and divide by cell area (m2) do l = 1, size(omega%data,3) call AreaOper( lli, omega%data(:,:,l), '/', 'm2', status ) IF_NOTOK_RETURN(status=1) end do ! info .. !call SetHistory( omega%tmi, mfw%tmi, status ) !call AddHistory( omega%tmi, 'convert to Pa/s', status ) ! ok status = 0 call goLabel() end subroutine Compute_Omega ! ************************************************************** ! *** ! *** Specific SETUP routine for CONVECTIVE FLUXES ! *** ! ************************************************************** subroutine Setup_Convec_SERIAL( region, entu, entd, detu, detd, omega, gph, & tr, lli, levi, status ) use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, Read_Convec, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! --- in/out ---------------------------------- integer, intent(in) :: region ! region number type(TMeteoData), intent(inout) :: entu, entd, detu, detd type(TMeteoData), intent(in) :: omega, gph type(TDate), intent(in) :: tr(2) type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/SETUP_CONVEC_SERIAL' ! --- local ---------------------------------- logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) ! to differentiate b/w local and global data set real, pointer, dimension(:,:,:) :: L_entu, L_entd, L_detu, L_detd real, pointer :: L_omega(:,:,:), L_gph(:,:,:) integer, dimension(2) :: is, js, ls, auxls integer :: halo ! --- begin ----------------------------- 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 if (okdebug) then write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(entu%name),trim(lli%name); call goPr endif ! 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 ! NOT NEEDED in EC-Earth, since we are using the ec-ll method (read_convec) ! ! omega is required as input: ! if ( .not. omega%used ) then ! write (gol,'("omega should be in use to compute convective stuff")'); call goErr ! call goErr; status=1; return ! end if ! not changed by default entu%changed = .false. entd%changed = .false. detu%changed = .false. detd%changed = .false. !------------------ ! time stuff !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from entu only) call SetupSetup( entu, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then ! Need whole region for I/O on root. Dummy else. is = (/1,im(region)/) js = (/1,jm(region)/) ls = entu%ls auxls = gph%ls IF (isRoot) THEN ! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo=0 allocate( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate(L_gph (im(region), jm(region), auxls(1):auxls(2)) ) allocate(L_omega(im(region), jm(region), auxls(1):auxls(2)) ) ELSE allocate( L_entu(1,1,1), L_entd(1,1,1), L_detu(1,1,1), L_detd(1,1,1)) allocate(L_gph (1,1,1)) allocate(L_omega(1,1,1)) END IF CALL GATHER( dgrid(region), gph%data, L_gph, gph%halo, status) IF_NOTOK_RETURN(status=1) ! NOT NEEDED in EC-Earth, since we are using the ec-ll method (read_convec) ! CALL GATHER( dgrid(region), omega%data, L_omega, omega%halo, status) ! IF_NOTOK_RETURN(status=1) ! Read/write on root IOroot : if (isRoot) then ! safety check ... ! if ( data1_t2 /= data1_t1 ) then !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr !call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr !call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr !call goErr; status=1; return ! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr ! end if ! surface pressure field: allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) ! fill data call Read_Convec( tmmd, entu%sourcekey, & data1_tref, data1_t1, data1_t2, lli, levi, & L_omega, omega%tmi, & L_gph, gph%tmi, & tmp_sp, & L_entu, entu%tmi1, L_entd, entd%tmi1, & L_detu, detu%tmi1, L_detd, detd%tmi1, & status ) IF_NOTOK_RETURN(status=1) ! write meteofiles if ( entu%putout ) then call WriteField( tmmd, entu%destkey, & entu%tmi1, 'sp', trim(entu%name), trim(entu%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, '*', & tmp_sp, L_entu, status ) IF_NOTOK_RETURN(status=1) end if if ( entd%putout ) then call WriteField( tmmd, entd%destkey, & entd%tmi1, 'sp', trim(entd%name), trim(entd%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, '*', & tmp_sp, L_entd, status ) IF_NOTOK_RETURN(status=1) end if if ( detu%putout ) then call WriteField( tmmd, detu%destkey, & detu%tmi1, 'sp', trim(detu%name), trim(detu%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, '*', & tmp_sp, L_detu, status ) IF_NOTOK_RETURN(status=1) end if if ( detd%putout ) then call WriteField( tmmd, detd%destkey, & detd%tmi1, 'sp', trim(detd%name), trim(detd%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, '*', & tmp_sp, L_detd, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if IOroot ! Scatter & clean up CALL SCATTER( dgrid(region), entu%data1, L_entu, entu%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), entd%data1, L_entd, entd%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), detu%data1, L_detu, detu%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), detd%data1, L_detd, detd%halo, status) IF_NOTOK_RETURN(status=1) deallocate(L_entu, L_entd, L_detu, L_detd, L_gph, L_omega) ! data array is filled now: entu%filled1 = .true. entu%tr1(1) = data1_t1 entu%tr1(2) = data1_t2 entu%changed = .true. entd%filled1 = .true. entd%tr1(1) = data1_t1 entd%tr1(2) = data1_t2 entd%changed = .true. detu%filled1 = .true. detu%tr1(1) = data1_t1 detu%tr1(2) = data1_t2 detu%changed = .true. detd%filled1 = .true. detd%tr1(1) = data1_t1 detd%tr1(2) = data1_t2 detd%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: entu%data1 = entu%data2 entd%data1 = entd%data2 detu%data1 = detu%data2 detd%data1 = detd%data2 ! data array is filled now: entu%filled1 = .true. entu%tr1(1) = data1_t1 entu%tr1(2) = data1_t2 entu%changed = .true. entd%filled1 = .true. entd%tr1(1) = data1_t1 entd%tr1(2) = data1_t2 entd%changed = .true. detu%filled1 = .true. detu%tr1(1) = data1_t1 detu%tr1(2) = data1_t2 detu%changed = .true. detd%filled1 = .true. detd%tr1(1) = data1_t1 detd%tr1(2) = data1_t2 detd%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! Need whole region for I/O on root. Dummy else is = (/1,im(1)/) js = (/1,jm(1)/) ls = entu%ls auxls = gph%ls IF (isRoot) THEN ! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo ALLOCATE( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE(L_gph (im(region),jm(region),auxls(1):auxls(2))) ALLOCATE(L_omega (im(region),jm(region),auxls(1):auxls(2))) ELSE ALLOCATE( L_entu(1,1,1), L_entd(1,1,1), L_detu(1,1,1), L_detd(1,1,1)) ALLOCATE( L_gph(1,1,1), L_omega(1,1,1) ) END IF CALL GATHER( dgrid(region), gph%data, L_gph, gph%halo, status) IF_NOTOK_RETURN(status=1) ! NOT NEEDED in EC-Earth, since we are using the ec-ll method (read_convec) ! CALL GATHER( dgrid(region), omega%data, L_omega, omega%halo, status) ! IF_NOTOK_RETURN(status=1) ! Read/write on root IOroot2 : if (isRoot) then ! safety check ... ! if ( data2_t2 /= data2_t1 ) then !write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr !call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr !call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr !write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr !call goErr; status=1; return ! write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr ! end if ! surface pressure field: allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) ! fill data2 call Read_Convec( tmmd, entu%sourcekey, & data2_tref, data2_t1, data2_t2, lli, levi, & L_omega, omega%tmi, & L_gph, gph%tmi, & tmp_sp, & L_entu, entu%tmi2, L_entd, entd%tmi2, & L_detu, detu%tmi2, L_detd, detd%tmi2, & status ) IF_NOTOK_RETURN(status=1) ! write meteofiles ? if ( entu%putout ) then call WriteField( tmmd, entu%destkey, & entu%tmi2, 'sp', trim(entu%name), trim(entu%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, '*', & tmp_sp, L_entu, status ) IF_NOTOK_RETURN(status=1) end if if ( entd%putout ) then call WriteField( tmmd, entd%destkey, & entd%tmi2, 'sp', trim(entd%name), trim(entd%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, '*', & tmp_sp, L_entd, status ) IF_NOTOK_RETURN(status=1) end if if ( detu%putout ) then call WriteField( tmmd, detu%destkey, & detu%tmi2, 'sp', trim(detu%name), trim(detu%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, '*', & tmp_sp, L_detu, status ) IF_NOTOK_RETURN(status=1) end if if ( detd%putout ) then call WriteField( tmmd, detd%destkey, & detd%tmi2, 'sp', trim(detd%name), trim(detd%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, '*', & tmp_sp, L_detd, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if IOroot2 CALL SCATTER( dgrid(region), entu%data2, L_entu, entu%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), entd%data2, L_entd, entd%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), detu%data2, L_detu, detu%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), detd%data2, L_detd, detd%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( L_entu, L_entd, L_detu, L_detd, L_gph, L_omega ) ! data2 array is filled now: entu%filled2 = .true. entu%tr2(1) = data2_t1 entu%tr2(2) = data2_t2 entd%filled2 = .true. entd%tr2(1) = data2_t1 entd%tr2(2) = data2_t2 detu%filled2 = .true. detu%tr2(1) = data2_t1 detu%tr2(2) = data2_t2 detd%filled2 = .true. detd%tr2(1) = data2_t1 detd%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data2 from primary array: entu%data2 = entu%data1 entd%data2 = entd%data1 detu%data2 = detu%data1 detd%data2 = detd%data1 ! data2 array is filled now: entu%filled2 = .true. entu%tr2(1) = data2_t1 entu%tr2(2) = data2_t2 entd%filled2 = .true. entd%tr2(1) = data2_t1 entd%tr2(2) = data2_t2 detu%filled2 = .true. detu%tr2(1) = data2_t1 detu%tr2(2) = data2_t2 detd%filled2 = .true. detd%tr2(1) = data2_t1 detd%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( entu, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( entd, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( detu, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( detd, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ status = 0 call goLabel() END SUBROUTINE SETUP_CONVEC_SERIAL !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: SETUP_CONVEC_PARA ! ! !DESCRIPTION: same as setup_convec_serial_io but with parallel i/o !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_CONVEC_PARA( region, entu, entd, detu, detd, omega, gph, & tr, levi, status ) ! ! !USES: ! use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, Read_Convec, WriteField ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: entu, entd, detu, detd type(TMeteoData), intent(in) :: omega, gph type(TDate), intent(in) :: tr(2) type(TLevelInfo), intent(in) :: levi ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 24 Oct 2013 - Ph. Le Sager - v0 ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/SETUP_CONVEC_PARA' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) ! to differentiate b/w local and global data set real, pointer, dimension(:,:,:) :: L_entu, L_entd, L_detu, L_detd real, pointer :: L_omega(:,:,:), L_gph(:,:,:) integer, dimension(2) :: is, js, ls, auxls integer :: halo integer :: i1, i2, j1, j2 ! --- begin ----------------------------- 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 !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from entu only) call SetupSetup( entu, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! work arrays if (data1_read .or. data2_read) then CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) is = (/i1,i2/) js = (/j1,j2/) ls = entu%ls auxls = gph%ls ! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo=0 allocate( L_entu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_entd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_detu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( L_detd( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) L_gph => gph%data L_omega => omega%data end if !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then !AJS if ( data1_t2 /= data1_t1 ) then !AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr !AJS end if call Read_Convec( tmmd, entu%sourcekey, & data1_tref, data1_t1, data1_t2, lli(region), levi, & L_omega, omega%tmi, & L_gph, gph%tmi, & tmp_sp, & L_entu, entu%tmi1, L_entd, entd%tmi1, & L_detu, detu%tmi1, L_detd, detd%tmi1, & status ) IF_NOTOK_RETURN(status=1) entu%data1(i1:i2,j1:j2,:) = L_entu entd%data1(i1:i2,j1:j2,:) = L_entd detu%data1(i1:i2,j1:j2,:) = L_detu detd%data1(i1:i2,j1:j2,:) = L_detd ! data array is filled now: entu%filled1 = .true. entu%tr1(1) = data1_t1 entu%tr1(2) = data1_t2 entu%changed = .true. entd%filled1 = .true. entd%tr1(1) = data1_t1 entd%tr1(2) = data1_t2 entd%changed = .true. detu%filled1 = .true. detu%tr1(1) = data1_t1 detu%tr1(2) = data1_t2 detu%changed = .true. detd%filled1 = .true. detd%tr1(1) = data1_t1 detd%tr1(2) = data1_t2 detd%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: entu%data1 = entu%data2 entd%data1 = entd%data2 detu%data1 = detu%data2 detd%data1 = detd%data2 ! data array is filled now: entu%filled1 = .true. entu%tr1(1) = data1_t1 entu%tr1(2) = data1_t2 entu%changed = .true. entd%filled1 = .true. entd%tr1(1) = data1_t1 entd%tr1(2) = data1_t2 entd%changed = .true. detu%filled1 = .true. detu%tr1(1) = data1_t1 detu%tr1(2) = data1_t2 detu%changed = .true. detd%filled1 = .true. detd%tr1(1) = data1_t1 detd%tr1(2) = data1_t2 detd%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then !AJS if ( data2_t2 /= data2_t1 ) then !AJS write (gol,'("WARNING - convec for interval, but pressure/gph/etc instant ...")'); call goPr !AJS end if call Read_Convec( tmmd, entu%sourcekey, & data2_tref, data2_t1, data2_t2, lli(region), levi, & L_omega, omega%tmi, & L_gph, gph%tmi, & tmp_sp, & L_entu, entu%tmi2, L_entd, entd%tmi2, & L_detu, detu%tmi2, L_detd, detd%tmi2, & status ) IF_NOTOK_RETURN(status=1) ! write meteofiles ? ! if ( entu%putout ) then ! end if entu%data2(i1:i2,j1:j2,:) = L_entu entd%data2(i1:i2,j1:j2,:) = L_entd detu%data2(i1:i2,j1:j2,:) = L_detu detd%data2(i1:i2,j1:j2,:) = L_detd ! data2 array is filled now: entu%filled2 = .true. entu%tr2(1) = data2_t1 entu%tr2(2) = data2_t2 entd%filled2 = .true. entd%tr2(1) = data2_t1 entd%tr2(2) = data2_t2 detu%filled2 = .true. detu%tr2(1) = data2_t1 detu%tr2(2) = data2_t2 detd%filled2 = .true. detd%tr2(1) = data2_t1 detd%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data2 from primary array: entu%data2 = entu%data1 entd%data2 = entd%data1 detu%data2 = detu%data1 detd%data2 = detd%data1 ! data2 array is filled now: entu%filled2 = .true. entu%tr2(1) = data2_t1 entu%tr2(2) = data2_t2 entd%filled2 = .true. entd%tr2(1) = data2_t1 entd%tr2(2) = data2_t2 detu%filled2 = .true. detu%tr2(1) = data2_t1 detu%tr2(2) = data2_t2 detd%filled2 = .true. detd%tr2(1) = data2_t1 detd%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( entu, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( entd, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( detu, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( detd, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ if (data1_read .or. data2_read) then deallocate(L_entu, L_entd, L_detu, L_detd) deallocate( tmp_sp ) nullify(L_gph, L_omega) end if status = 0 call goLabel() END SUBROUTINE SETUP_CONVEC_PARA !EOC ! ************************************************************** ! *** ! *** Specific SETUP routine for CLOUD COVER ! *** ! ************************************************************** SUBROUTINE SETUP_CLOUDCOVERS_SERIAL( region, cc, cco, ccu, tr, lli, levi, status ) use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, Read_CloudCovers, WriteField use meteodata, only : TMeteoData, TimeInterpolation use dims, only : im, jm ! --- in/out ---------------------------------- integer, intent(in) :: region ! region number type(TMeteoData), intent(inout) :: cc, cco, ccu type(TDate), intent(in) :: tr(2) type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/Setup_CloudCovers_Serial' ! --- local ---------------------------------- logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) ! surface pressure real, pointer, dimension(:,:,:) :: L_cc, L_cco, L_ccu ! work arrays (data) integer :: is(2), js(2), ls(2) ! work arrays (bounds) ! --- begin ----------------------------- 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 if (okdebug) then write (gol,'(" ",a,": ",a," @ ",a)') rname, trim(cc%name),trim(lli%name); call goPr endif ! not changed by default cc%changed = .false. cco%changed = .false. ccu%changed = .false. !------------------ ! time stuff !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from cc only) call SetupSetup( cc, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then ! Allocate global arrays for I/O is = (/1,im(region)/) js = (/1,jm(region)/) ls = cc%ls IF (isRoot) THEN ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ELSE ALLOCATE( L_cc(1,1,1), L_cco(1,1,1), L_ccu(1,1,1) ) END IF ! Read/write on root IOroot : if (isRoot) then ! safety check ... if ( data1_t2 /= data1_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) ! fill data: call Read_CloudCovers( tmmd, cc%sourcekey, & data1_tref, data1_t1, data1_t2, lli, levi, & tmp_sp, & L_cc, cc%tmi1, & L_cco, cco%tmi1, & L_ccu, ccu%tmi1, & status ) IF_NOTOK_RETURN(status=1) ! write meteofiles if ( cc%putout ) then call WriteField( tmmd, cc%destkey, & cc%tmi1, 'sp', trim(cc%name), trim(cc%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, 'n', & tmp_sp, L_cc, status ) IF_NOTOK_RETURN(status=1) end if if ( cco%putout ) then call WriteField( tmmd, cco%destkey, & cco%tmi1, 'sp', trim(cco%name), trim(cco%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, 'n', & tmp_sp, L_cco, status ) IF_NOTOK_RETURN(status=1) end if if ( ccu%putout ) then call WriteField( tmmd, ccu%destkey, & ccu%tmi1, 'sp', trim(ccu%name), trim(ccu%unit), & data1_tref, data1_t1, data1_t2, & lli, 'n', levi, 'n', & tmp_sp, L_ccu, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if IOroot ! Wrap up CALL SCATTER( dgrid(region), cc%data1, L_cc, cc%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), cco%data1, L_cco, cco%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), ccu%data1, L_ccu, ccu%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(L_cc, L_cco, L_ccu) ! data array is filled now: cc%filled1 = .true. cc%tr1(1) = data1_t1 cc%tr1(2) = data1_t2 cc%changed = .true. cco%filled1 = .true. cco%tr1(1) = data1_t1 cco%tr1(2) = data1_t2 cco%changed = .true. ccu%filled1 = .true. ccu%tr1(1) = data1_t1 ccu%tr1(2) = data1_t2 ccu%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: cc%data1 = cc%data2 cco%data1 = cco%data2 ccu%data1 = ccu%data2 ! data array is filled now: cc%filled1 = .true. cc%tr1(1) = data1_t1 cc%tr1(2) = data1_t2 cc%changed = .true. cco%filled1 = .true. cco%tr1(1) = data1_t1 cco%tr1(2) = data1_t2 cco%changed = .true. ccu%filled1 = .true. ccu%tr1(1) = data1_t1 ccu%tr1(2) = data1_t2 ccu%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! Allocate global arrays for I/O is = (/1,im(region)/) js = (/1,jm(region)/) ls = cc%ls IF (isRoot) THEN ! Use the fact that entu, entd, detu, and detd have been allocated with the same bounds and halo ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ELSE ALLOCATE( L_cc(1,1,1), L_cco(1,1,1), L_ccu(1,1,1) ) END IF ! Read/write IOroot2 : IF (isRoot) THEN ! safety check ... if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data2_t1 : ', data2_t1 ); call goErr call wrtgol( ' data2_t2 : ', data2_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr call goErr; status=1; return end if ! surface pressure field: allocate( tmp_sp(is(1):is(2),js(1):js(2)) ) ! fill data2: call Read_CloudCovers( tmmd, cc%sourcekey, data2_tref, & data2_t1, data2_t2, lli, levi, & tmp_sp, & L_cc, cc%tmi2, & L_cco, cco%tmi2, & L_ccu, ccu%tmi2, & status ) IF_NOTOK_RETURN(status=1) ! write meteofiles ? if ( cc%putout ) then call WriteField( tmmd, cc%destkey, & cc%tmi2, 'sp', trim( cc%name), trim( cc%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, 'n', & tmp_sp, L_cc, status ) IF_NOTOK_RETURN(status=1) end if if ( cco%putout ) then call WriteField( tmmd, cco%destkey, & cco%tmi2, 'sp', trim(cco%name), trim(cco%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, 'n', & tmp_sp, L_cco, status ) IF_NOTOK_RETURN(status=1) end if if ( ccu%putout ) then call WriteField( tmmd, ccu%destkey, & ccu%tmi2, 'sp', trim(ccu%name), trim(ccu%unit), & data2_tref, data2_t1, data2_t2, & lli, 'n', levi, 'n', & tmp_sp, L_ccu, status ) IF_NOTOK_RETURN(status=1) end if ! clear deallocate( tmp_sp ) end if IOroot2 ! Wrap up CALL SCATTER( dgrid(region), cc%data2, L_cc, cc%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), cco%data2, L_cco, cco%halo, status) IF_NOTOK_RETURN(status=1) CALL SCATTER( dgrid(region), ccu%data2, L_ccu, ccu%halo, status) IF_NOTOK_RETURN(status=1) DEALLOCATE(L_cc, L_cco, L_ccu) ! data2 array is filled now: cc%filled2 = .true. cc%tr2(1) = data2_t1 cc%tr2(2) = data2_t2 cco%filled2 = .true. cco%tr2(1) = data2_t1 cco%tr2(2) = data2_t2 ccu%filled2 = .true. ccu%tr2(1) = data2_t1 ccu%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data2 from primary array: cc%data2 = cc%data1 cco%data2 = cco%data1 ccu%data2 = ccu%data1 ! data2 array is filled now: cc%filled2 = .true. cc%tr2(1) = data2_t1 cc%tr2(2) = data2_t2 cco%filled2 = .true. cco%tr2(1) = data2_t1 cco%tr2(2) = data2_t2 ccu%filled2 = .true. ccu%tr2(1) = data2_t1 ccu%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( cc, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( cco, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( ccu, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ status = 0 call goLabel() END SUBROUTINE SETUP_CLOUDCOVERS_SERIAL !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: SETUP_CLOUDCOVERS_PARA ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE SETUP_CLOUDCOVERS_PARA( region, cc, cco, ccu, tr, levi, status ) ! ! !USES: ! use GO, only : TDate, wrtgol, operator(/=) use Grid, only : TllGridInfo, TLevelInfo use TMM, only : TMeteoInfo, Read_CloudCovers, WriteField use dims, only : im, jm ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region number ! ! !INPUT/OUTPUT PARAMETERS: ! type(TMeteoData), intent(inout) :: cc, cco, ccu type(TDate), intent(in) :: tr(2) type(TLevelInfo), intent(in) :: levi ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 24 Oct 2013 - Ph. Le Sager - v0 ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/SETUP_CLOUDCOVERS_PARA' logical :: data1_read, data1_copy type(TDate) :: data1_tref, data1_t1, data1_t2 logical :: data2_read, data2_copy type(TDate) :: data2_tref, data2_t1, data2_t2 real, allocatable :: tmp_sp(:,:) ! surface pressure real, pointer, dimension(:,:,:) :: L_cc, L_cco, L_ccu ! work arrays (data) integer :: is(2), js(2), ls(2) ! work arrays (bounds) integer :: i1, i2, j1, j2 ! --- begin ----------------------------- 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 !------------------ ! get time interval of met field and check if data from start and/or end ! of interval must be read (sufficient to setup from cc only) call SetupSetup( cc, tr, & data1_read, data1_copy, data1_tref, data1_t1, data1_t2, & data2_read, data2_copy, data2_tref, data2_t1, data2_t2, & status ) IF_NOTOK_RETURN(status=1) ! work arrays IF (data1_read .OR. data2_read) THEN CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) is = (/i1,i2/) js = (/j1,j2/) ls = cc%ls ALLOCATE( L_cc( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_cco( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( L_ccu( is(1):is(2), js(1):js(2), ls(1):ls(2)) ) ALLOCATE( tmp_sp(is(1):is(2),js(1):js(2)) ) ENDIF !-------------------------- ! read/write primary field !-------------------------- if ( data1_read ) then ! safety check if ( data1_t2 /= data1_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")'); call goErr call wrtgol( ' data1_t1 : ', data1_t1 ); call goErr call wrtgol( ' data1_t2 : ', data1_t2 ); call goErr write (gol,'("please deceide what to do with surface pressures ... ")'); call goErr TRACEBACK; status=1; return end if call Read_CloudCovers( tmmd, cc%sourcekey, & data1_tref, data1_t1, data1_t2, lli(region), levi, & tmp_sp, & L_cc, cc%tmi1, & L_cco, cco%tmi1, & L_ccu, ccu%tmi1, & status ) IF_NOTOK_RETURN(status=1) cc%data1(i1:i2,j1:j2,:) = L_cc cco%data1(i1:i2,j1:j2,:) = L_cco ccu%data1(i1:i2,j1:j2,:) = L_ccu ! data array is filled now: cc%filled1 = .true. cc%tr1(1) = data1_t1 cc%tr1(2) = data1_t2 cc%changed = .true. cco%filled1 = .true. cco%tr1(1) = data1_t1 cco%tr1(2) = data1_t2 cco%changed = .true. ccu%filled1 = .true. ccu%tr1(1) = data1_t1 ccu%tr1(2) = data1_t2 ccu%changed = .true. else if ( data1_copy ) then ! copy data from secondary array: cc%data1 = cc%data2 cco%data1 = cco%data2 ccu%data1 = ccu%data2 ! data array is filled now: cc%filled1 = .true. cc%tr1(1) = data1_t1 cc%tr1(2) = data1_t2 cc%changed = .true. cco%filled1 = .true. cco%tr1(1) = data1_t1 cco%tr1(2) = data1_t2 cco%changed = .true. ccu%filled1 = .true. ccu%tr1(1) = data1_t1 ccu%tr1(2) = data1_t2 ccu%changed = .true. end if !-------------------------- ! read/write secondary field !-------------------------- if ( data2_read ) then ! safety check ... if ( data2_t2 /= data2_t1 ) then write (gol,'("not sure that this routine is correct for time intervals:")') ; call goErr call wrtgol( ' data2_t1 : ', data2_t1 ) ; call goErr call wrtgol( ' data2_t2 : ', data2_t2 ) ; call goErr write (gol,'("please deceide what to do with surface pressures ... ")') ; call goErr call goErr; status=1; return end if call Read_CloudCovers( tmmd, cc%sourcekey, data2_tref, & data2_t1, data2_t2, lli(region), levi, & tmp_sp, & L_cc, cc%tmi2, & L_cco, cco%tmi2, & L_ccu, ccu%tmi2, & status ) IF_NOTOK_RETURN(status=1) cc%data2(i1:i2,j1:j2,:) = L_cc cco%data2(i1:i2,j1:j2,:) = L_cco ccu%data2(i1:i2,j1:j2,:) = L_ccu ! data2 array is filled now: cc%filled2 = .true. cc%tr2(1) = data2_t1 cc%tr2(2) = data2_t2 cco%filled2 = .true. cco%tr2(1) = data2_t1 cco%tr2(2) = data2_t2 ccu%filled2 = .true. ccu%tr2(1) = data2_t1 ccu%tr2(2) = data2_t2 else if ( data2_copy ) then ! copy data2 from primary array: cc%data2 = cc%data1 cco%data2 = cco%data1 ccu%data2 = ccu%data1 ! data2 array is filled now: cc%filled2 = .true. cc%tr2(1) = data2_t1 cc%tr2(2) = data2_t2 cco%filled2 = .true. cco%tr2(1) = data2_t1 cco%tr2(2) = data2_t2 ccu%filled2 = .true. ccu%tr2(1) = data2_t1 ccu%tr2(2) = data2_t2 end if !------------------ ! time interpolation !------------------ call TimeInterpolation( cc, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( cco, tr, status ) IF_NOTOK_RETURN(status=1) call TimeInterpolation( ccu, tr, status ) IF_NOTOK_RETURN(status=1) !------------------ ! done !------------------ if (data1_read .or. data2_read) then deallocate( tmp_sp ) deallocate(L_cc, L_cco, L_ccu) end if status = 0 call goLabel() END SUBROUTINE SETUP_CLOUDCOVERS_PARA !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PRESSURE_TO_MASS ! ! !DESCRIPTION: Get Air Mass: from surface pressure (sp), get pressure at ! box boundaries (so-called half-levels, phlb), and then air ! mass (m_dat). !\\ !\\ ! !INTERFACE: ! SUBROUTINE PRESSURE_TO_MASS( region, status ) ! ! !USES: ! use Binas, only : grav use Grid, only : HPressure !use Grid, only : FillMass use Grid, only : AreaOper use dims, only : im, jm, lm use dims, only : xcyc ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: (old remark: "assume that halo cells in sp have been filled ! correctly..." still valid???) ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Pressure_to_Mass' integer :: l, i0, i1, j0, j1, lmr ! --- begin ---------------------------------- ! Local grid size i0 = sp_dat(region)%is(1) i1 = sp_dat(region)%is(2) j0 = sp_dat(region)%js(1) j1 = sp_dat(region)%js(2) lmr = lm(region) ! Fill pressure boundaries (Pa) if ( phlb_dat(region)%used ) then call HPressure( levi, sp_dat(region)%data(i0:i1, j0:j1, 1), & phlb_dat(region)%data(i0:i1, j0:j1, :), status ) IF_NOTOK_RETURN(status=0) end if ! Fill air mass (kg) if ( m_dat(region)%used ) then !call FillMass( m_dat(region)%data(1:imr,1:jmr,:), lli(region), levi, & ! sp_dat(region)%data(1:imr,1:jmr,1), status ) !IF_NOTOK_RETURN(status=0) ! Pressure difference between top and bottom of layer do l = 1, lmr m_dat(region)%data(:,:,l) = phlb_dat(region)%data(:,:,l) - phlb_dat(region)%data(:,:,l+1) ! Pa end do ! Convert to kg/m2 m_dat(region)%data = m_dat(region)%data / grav ! Pa/g = kg/m2 ! Convert to kg call AreaOper( lli(region), m_dat(region)%data(i0:i1, j0:j1, :), '*', 'm2', status ) ! kg IF_NOTOK_RETURN(status=0) end if ! ok status = 0 END SUBROUTINE PRESSURE_TO_MASS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: MASS_TO_PRESSURE ! ! !DESCRIPTION: get 3D and surface (spm) pressures from 3D Air Mass. !\\ !\\ ! !INTERFACE: ! SUBROUTINE MASS_TO_PRESSURE( region, status ) ! ! !USES: ! use Binas, only : grav use Grid, only : AreaOper use dims, only : lm ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 7 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Mass_to_Pressure' integer :: l, i0, i1, j0, j1, lmr ! --- begin ---------------------------------- ! Local grid size i0 = sp_dat(region)%is(1) i1 = sp_dat(region)%is(2) j0 = sp_dat(region)%js(1) j1 = sp_dat(region)%js(2) lmr = lm(region) ! Fill pressure at half level boundaries: ! o zero in space: phlb_dat(region)%data(:,:,lmr+1) = 0.0 ! kg m/s2 = Pa m2 ! o add for each level pressure gradient: do l = lmr, 1, -1 phlb_dat(region)%data(i0:i1, j0:j1, l) = phlb_dat(region)%data(i0:i1, j0:j1, l+1) & + m_dat(region)%data(i0:i1, j0:j1, l ) * grav ! kg m/s2 = Pa m2 end do ! Divide by grid cell area call AreaOper( lli(region), phlb_dat(region)%data(i0:i1, j0:j1, :), '/', 'm2', status ) ! Pa IF_NOTOK_RETURN(status=0) ! copy surface pressure spm_dat(region)%data(i0:i1, j0:j1, 1) = phlb_dat(region)%data(i0:i1, j0:j1, 1) ! Pa ! ok status = 0 END SUBROUTINE MASS_TO_PRESSURE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COMPUTE_GPH ! ! !DESCRIPTION: compute geopotential height !\\ !\\ ! !INTERFACE: ! SUBROUTINE COMPUTE_GPH( region, status ) ! ! !USES: ! use Dims, only : itau, lm use Dims, only : at, bt use binas, only : grav use datetime, only : tstamp ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/compute_gph' ! --- local ---------------------------------- real,dimension(:,:,:),pointer :: gph, t, q real,dimension(:,:,:),pointer :: ps integer :: i,j,l,i0,i1,j0,j1 real :: tv,pdown,pup ! --- begin ----------------------------- ! leave if not in use: if ( .not. gph_dat(region)%used ) then if (okdebug) then write (gol,'(a," not used on : ",i2)') trim(gph_dat(region)%name),region; call goPr endif status=0; return end if ! other meteo required: if ( (.not. temper_dat(region)%used) .or. (.not. humid_dat(region)%used) & .or. (.not. sp_dat(region)%used) .or. (.not. oro_dat(region)%used)) then write (gol,'("computation of gph requires temper, humid, sp, and oro")'); call goErr TRACEBACK; status=1; return end if ! leave if input did not change: if ( (.not. sp_dat(region)%changed) .and. & (.not. temper_dat(region)%changed) .and. & (.not. humid_dat(region)%changed) ) then if (okdebug) then write (gol,'(a,": not changed for region ",i2)') rname, region; call goPr endif status=0 return end if ! field will be changed ... gph_dat(region)%changed = .true. ! pointers to meteo field ps => sp_dat(region)%data t => temper_dat(region)%data q => humid_dat(region)%data gph => gph_dat(region)%data ! bounds w/o halo (same as: call Get_DistGrid( dgrid(region), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 ) i0 = gph_dat(region)%is(1) i1 = gph_dat(region)%is(2) j0 = gph_dat(region)%js(1) j1 = gph_dat(region)%js(2) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! compute geo potential height ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ gph(i0:i1,j0:j1,1) = oro_dat(region)%data(i0:i1,j0:j1,1)/grav ! oro is stored in g*m do l=1,lm(region)-1 do j=j0,j1 do i=i0,i1 tv = t(i,j,l)*(1. + 0.608*q(i,j,l)) pdown = at(l) + bt(l)*ps(i,j,1) pup = at(l+1) + bt(l+1)*ps(i,j,1) ! rgas in different units! gph(i,j,l+1) = gph(i,j,l) + tv*287.05*alog(pdown/pup)/grav ! note dec 2002 (MK) gph now from 1--->lm+1 end do end do end do !set top of atmosphere at 200 km gph(:,:,lm(region)+1) = 200000.0 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! done ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ nullify( ps ) nullify( t ) nullify( q ) nullify( gph ) status = 0 END SUBROUTINE COMPUTE_GPH !EOC END MODULE METEO