ecearth3/sources/tm5mp/proj/ecearth/meteo.F90

!### 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.<param>              <-- time interpolation
  !   tmm.sourcekey.<grid>.<param>       <-- input file name description
  !   tmm.output.<grid>.<param>          <-- write meteo ?
  !   tmm.destkey.<grid>.<param>         <-- output file name description
  !
  ! where <grid> is first '*' and then set to the grid name,
  ! and <param> 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