!#################################################################
!BOP
!
! !MODULE:  diffusion
!
! !DESCRIPTION:
!
! Diffusion routines.
!
!
! !REVISION HISTORY: -----------------------------
!
! ????-?? Bram Bregman.
!   Adjusted for TM5 by
!
! 2005-04 ??
!   Adjoint version
!
! 2010-08 Arjo Segers
!   Debugged computation of u* over sea: 
!   abs. surface stress should be devided by air density.
!   Set minum value for u* to trap division by zero (MK).
!
!EOP
!#################################################################
!
! global input fields:
! phlb(:,:,lm+1)          ! half level pressure [Pa]
! gph(:,:,lm+1)           ! height of half-levels [m]    !CMK note gph starts NOW at 1!
! t                       ! temperature [K]
! q                       ! specific humidity [kg/kg]
! uwind                   ! mean wind W-E [m/s]
! vwind                   ! mean wind S-N [m/s]
! sshf                    ! surface sensible heat flux [W / m^2]
! slhf                    ! surface latent heat flux [W / m^2]
! ustar (ustar_loc)          ! velocity scale
!
! global output fields:
! kvh                     ! eddy diffusivity for heat at top [m^2/s]
! pblh                    ! boundary layer height[m]
! dkg                                     ! vertical diffusion coefficient [ks s-1]
!
!### 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"
!
!#################################################################

module diffusion

  use GO, only : gol, goPr, goErr
  
  use dims        , only : nregions
  use global_types, only : emis_data

  implicit none
  
  ! --- in/out ------------------------------------

  private

  public  ::  Diffusion_Init, Diffusion_Done
  public  ::  Calc_Kzz, Read_Diffusion, Write_Diffusion
  
  public  :: diffusion_write
  
  ! --- const --------------------------------------

  character(len=*), parameter  ::  mname = 'diffusion'

  logical                      ::  diffusion_write = .false.
  character(len=2000)          ::  diffusion_dir = ''
  
  ! --- local ------------------------------------

  type(emis_data), dimension(nregions), target  ::  ustar_loc, sr_mix


contains


  ! ================================================================
  
  
  subroutine Diffusion_Init( status )
  
    use global_data, only : rcfile
    use Meteo      , only : Set
    use Meteo      , only : pu_dat, pv_dat
    use Meteo      , only : temper_dat, humid_dat, gph_dat
    use Meteo      , only : oro_dat, lsmask_dat
    use Meteo      , only : sr_ecm_dat, sr_ols_dat
    use Meteo      , only : u10m_dat, v10m_dat
    use Meteo      , only : slhf_dat, sshf_dat
    use Meteo      , only : ewss_dat, nsss_dat
    use GO         , only : TrcFile, Init, Done, ReadRc
    use toolbox    , only : escape_tm                                                                                                                                                                   



    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Diffusion_Init'
    
    ! --- local -------------------------------------
    
    integer           ::  region
    type(TrcFile)               :: rcF        
    
    ! --- begin --------------------------------
    
    ! loop over all (zoom) regions:
    do region = 1, nregions
      ! meteo used by diffusion
      call Set(     pu_dat(region), status, used=.true. )
      call Set(     pv_dat(region), status, used=.true. )
      call Set( temper_dat(region), status, used=.true. )
      call Set(  humid_dat(region), status, used=.true. )
      call Set(    gph_dat(region), status, used=.true. )
      call Set(    oro_dat(region), status, used=.true. )
      call Set( lsmask_dat(region), status, used=.true. )
      call Set( sr_ecm_dat(region), status, used=.true. )
      call Set( sr_ols_dat(region), status, used=.true. )
      call Set(   u10m_dat(region), status, used=.true. )
      call Set(   v10m_dat(region), status, used=.true. )
      call Set(   sshf_dat(region), status, used=.true. )
      call Set(   slhf_dat(region), status, used=.true. )
      call Set(   ewss_dat(region), status, used=.true. )
      call Set(   nsss_dat(region), status, used=.true. )
    end do
        
    call Init( rcF, rcfile , status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc(rcF, 'diffusion.write', diffusion_write, status, default = .false.)
    IF_ERROR_RETURN(status=1)

    if (diffusion_write) then
       call ReadRc(rcF, 'diffusion.dir', diffusion_dir, status)
       IF_NOTOK_RETURN(status=1)
    endif
    call Done(rcF, status)
    IF_NOTOK_RETURN(status=1)
                                                  

    ! ok
    status = 0
    
  end subroutine Diffusion_Init
  
  
  ! ***
  
  
  subroutine Diffusion_Done( status )
  
    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Diffusion_Done'

    ! --- begin --------------------------------
    
    ! nothing to be done

    ! ok
    status = 0
    
  end subroutine Diffusion_Done
  
  
  
  ! ================================================================
  

  !------------------------------------
  !
  ! Purpose:
  ! -------
  ! this subroutine reads and prepares all fields needed for the calculation of kzz
  !
  ! External
  ! --------
  ! dd_get_3_hourly_surface_fields
  ! dd_calc_ustar_raero_rb
  ! dd_coarsen_fields
  !
  ! Reference
  ! ---------
  ! Ganzeveld and Lelieveld (1996) and references therein.
  ! Adjusted for TM5, Bram Bregman, August 2003
  !
  !

  subroutine calc_kzz

    use binas,       only: grav
    use dims,        only: im, jm, lmax_conv, idate, lm, okdebug
    use dims,        only: nread, ndyn, tref, at, bt
    use dims,        only: isr,jsr,ier,jer
    use dims,        only: revert
    use global_data, only: mass_dat, wind_dat, region_dat
#ifndef without_convection
#ifndef without_diffusion
    use global_data, only: conv_dat
#endif
#endif
    use meteo      , only: spm_dat
    use meteo      , only: pu_dat, pv_dat
    use meteo      , only: temper_dat, humid_dat, gph_dat
    use meteodata  , only : m_dat
#ifdef MPI
    use mpi_const  , only: myid, root_t, ierr, com_trac, my_real, lmloc
#endif

    ! --- in/out ----------------------------

    ! --- local ------------------------------

    character(len=10)                            :: c_time
    real,dimension(:,:,:),pointer                :: phlb,gph,t,q,m,pu,pv
    real,dimension(:,:,:),pointer                :: am,bm,cm ! fluxes in kg/timestep(region)
#ifndef without_convection
#ifndef without_diffusion
    real,dimension(:,:,:),pointer                :: dkg     !vert diff (ks/s)
#endif
#endif
    real,dimension(:),pointer                    :: dxyp
    integer,parameter                            :: avg_field=1,nlon360=360,nlat180=180
    integer                                      :: i,region, nsend, ntimes
    real, dimension(:,:,:),allocatable,target    :: m_adj, phlb_adj
    real, dimension(:,:),allocatable             :: p_adj
    integer                                      :: imr, jmr, lmr, j,l

    real, allocatable, target                    :: phlb_t(:,:,:)
    real, allocatable, target                    :: m_t(:,:,:)
    
    ! --- begin ----------------------------------

    if (okdebug) write(*,*) 'start of calc_kzz'

#ifdef MPI
    if(myid == root_t) then    ! reading /calculations by root!
#endif

      do region = 1, nregions
        allocate(ustar_loc(region)%surf(im(region),jm(region)))
        allocate(sr_mix(region)%surf(im(region),jm(region)))
      enddo                                                                                                                           

      ! -- for output of time header on debug fields

      write(c_time,'(i4,3(i2))') idate(1:4)   !CMK corrected
      do i=1,10
        if (c_time(i:i)==' ') c_time(i:i)='0'
      enddo

      !
      ! surface data sets may have the following characteristics
      ! 1. instanteneous. The time written in the file is valid from t-1.5 to t+1.5
      ! 2. Accumulated. The time written in the file is valid from t to t+3
      ! 3. Daily averaged. The time on the file is always 00 hours, and valid until t+24
      !
      ! ***   It is essential to understand this timing when reading and applying these sets.
      !
      !
      ! fd2mk it has to be decided to 'save' fields like vgrat_low and daily average surface fields
      ! for later use, or to read it every 3 hours time step

      do region = 1,nregions

        ! local grid sizes
        imr = im(region) ; jmr = jm(region) ; lmr = lm(region)
        
        if (revert == -1) then 
           !mkadj
           !mass and pressure are at t=t+dt in adjoint
           !these should be brought back to t to get
           !exactly the same diffusion as in forward run
           allocate ( m_adj(1:imr,1:jmr,1:lmr) )
           allocate ( phlb_adj(imr,jmr,lmr+1))
           allocate ( p_adj(imr,jmr))
           am => wind_dat(region)%am_t
           bm => wind_dat(region)%bm_t
           cm => wind_dat(region)%cm_t
           dxyp => region_dat(region)%dxyp
           ntimes = (nread/(ndyn))*tref(region)
           m_adj = m_dat(region)%data(1:imr,1:jmr,1:lmr)
           do i=1,ntimes
              m_adj = m_adj + revert*am(0:imr-1,1:jmr  ,1:lmr  ) &   !east   
                            - revert*am(1:imr  ,1:jmr  ,1:lmr  ) &   !west
                            + revert*bm(1:imr  ,1:jmr  ,1:lmr  ) &   !south
                            - revert*bm(1:imr  ,2:jmr+1,1:lmr  ) &   !north
                            + revert*cm(1:imr  ,1:jmr  ,0:lmr-1) &   !lower
                            - revert*cm(1:imr  ,1:jmr  ,1:lmr  )     !upper
           enddo
           p_adj(:,:) = 0.0
           do l=1,lmr
              ! note that on the EDGES masses are coarse>
              ! diffusion iis applied only on core zoom>
              do j=jsr(region), jer(region)  
                 do i=isr(region), ier(region)
                    p_adj(i,j) = p_adj(i,j) + m_adj(i,j,l)*grav/dxyp(j)
                 end do
              end do
           end do
           do l=1,lmr+1
              do j=jsr(region), jer(region)  
                 do i=isr(region), ier(region)
                    phlb_adj(i,j,l) = at(l)+bt(l)*p_adj(i,j)
                 end do
              end do
           end do
           deallocate(p_adj)
           nullify(am)
           nullify(bm)
           nullify(cm)
           nullify(dxyp)
           
           phlb => phlb_adj
           m => m_adj
        else
           ! Bug fix 2006-09-21 (AJS)
           ! The two half level pressure arrays phbl_t and phlb_k
           ! are sometimes not updated together when running in parallel
           ! with zoom regions. To avoid problems, pressure and mass
           ! are computed here given the surface pressure in the middle
           ! of the current time interval, stored in spm_dat.
           !--
           !phlb => mass_dat(region)%phlb_t
           !m => mass_dat(region)%m_t
           !--
           allocate( phlb_t(1:imr,1:jmr,1:lmr+1) )
           do l = 1, lmr+1
             phlb_t(:,:,l) = at(l) + bt(l)*spm_dat(region)%data(1:imr,1:jmr,1)
           end do
           allocate( m_t(-1:imr+2,-1:jmr+2,1:lmr) )
           m_t = 0.0
           do l = 1, lmr
             do j = 1, jmr
               m_t(1:imr,j,l) = ( phlb_t(:,j,l) - phlb_t(:,j,l+1) )*region_dat(region)%dxyp(j)/grav
             end do
           end do
           phlb => phlb_t
           m    => m_t
           !--
        endif
        pu  =>     pu_dat(region)%data
        pv  =>     pv_dat(region)%data
        gph =>    gph_dat(region)%data
        t   => temper_dat(region)%data
        q   =>  humid_dat(region)%data
#ifndef without_convection
#ifndef without_diffusion
        dkg => conv_dat(region)%dkg
#endif
#endif

        ! ustar_loc
        call dd_calc_ustar(region)

#ifndef without_convection
#ifndef without_diffusion
        ! calculate kzz                                                         
        call bldiff( region, phlb, gph, t, q, m, pu, pv, dkg )  
#endif
#endif                                                                                   

        nullify(phlb)
        nullify(m)
        !--
        deallocate( phlb_t )
        deallocate( m_t )
        !--
        nullify(pu)
        nullify(pv)
#ifndef without_convection
#ifndef without_diffusion
        nullify(dkg)
#endif
#endif
        nullify(t)
        nullify(q)
        nullify(gph)
        if(revert == -1) then
           deallocate(m_adj)
           deallocate(phlb_adj)
        endif

      end do

#ifdef MPI
    end if   !myid == root_t
#endif

#ifdef MPI
#ifndef without_convection
#ifndef without_diffusion
    do region = 1, nregions
      nsend = im(region)*jm(region)*lmax_conv 
      call mpi_bcast(conv_dat(region)%dkg ,nsend, my_real, root_t, com_trac, ierr)
    end do
#endif
#endif

    if ( myid == root_t ) then
#endif
      if (okdebug) print *, 'vertical diffusion kzz (dkg) broadcasted' 

      do region = 1,nregions
        deallocate(sr_mix(region)%surf)
        deallocate(ustar_loc(region)%surf)
      end do

#ifdef MPI
    endif
#endif

    if (okdebug) write(*,*) 'end of calc_kzz'

  end subroutine calc_kzz


    
    !------------------------------------------------------------------------------
    !                    TM5                                                      !
    !------------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:  dd_calc_ustar
    !
    ! !DESCRIPTION: 
    !
    !   Calculate friction velocity (ustar)
    ! 
    !   Uses the log normal wind profile information stored by ECMWF in 10 meter wind
    !   to estimate a local ustar over land
    !   uses Charnock equation over sea to estimate ustar
    !   aerodynamic resistance is calculated from heat fluxes and ustar using a constant reference height
    !   sub laminar resistance from ustar
    !
    !   Reference:
    !    o M.Z. Jacobson, "Fundamentals of Atmospheric Modeling", sections 8.3-8.4
    !
    ! !REVISION HISTORY:
    !
    ! 2003 ??; Ge Verver (personal communication, 2003)
    !
    ! 2010-08 Arjo Segers
    !   Debugged computation of u* over sea: 
    !   abs. surface stress should be devided by air density.
    !   Set minum value for u* to trap division by zero (MK).
    !
    !------------------------

    subroutine dd_calc_ustar(region)

      use binas, only : grav, vKarman
      use dims, only : im, jm
      use dims, only : okdebug
      use meteo, only : lsmask_dat, sr_ecm_dat, sr_ols_dat
      use meteo, only : u10m_dat, v10m_dat
      use meteo, only : slhf_dat, sshf_dat
      use meteo, only : ewss_dat, nsss_dat

      !--- in/out ------------------------------------

      integer,intent(in)              :: region
      
      ! --- const -------------------------------------

      ! Garret relation:
      real,parameter  :: alfa_garret    = 0.11
      ! air density at seal level and 293 K :
      real, parameter :: rho_air = 1.2  ! kg/m3

      ! Charnock relation:
      real,parameter  :: alfa_charnock  = 0.018
      ! kinematic viscosity of air at about 300 K :
      real,parameter  :: nu_air         = 1.5e-5    ! m2/s
      
      real, parameter :: Href=30.     ! constant reference height for calculations of raero
                                      ! some constants specific to calculation of ustar and raero
      real,parameter  :: rhoCp          = 1231.0
      real,parameter  :: rhoLv          = 3013000.0
      
      real, parameter :: rz0 =2.0

      !--- local -------------------------------------

      integer             ::  i,j
      real                ::  tau_z_abs
      real                ::  ustar_sea,ustar_land,xland,sr_sea,sr_land,sr_help,wind10m
      real, allocatable   ::  field(:,:)

      ! --- begin ----------------------------------

      allocate(field(im(region),jm(region)))

      do j=1,jm(region)
        do i=1,im(region)
          xland=lsmask_dat(region)%data(i,j,1)/100.    ! 0-1

          ! SEA:
          !
          ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling",
          !   section "8.3 Friction velocity" :
          ! The friction velocity :
          !   u* = sqrt( |tau_z|/rho_a )   [m/s]
          ! where:
          !   tau_z = (ewss,nsss)       surface stress [N/m2=kg/m/s2]
          !   rho_air                   air density [kg/m3] ; about 1.2  at sea level and 293 K
          !
          tau_z_abs = sqrt(ewss_dat(region)%data(i,j,1)**2+&
                           nsss_dat(region)%data(i,j,1)**2  )   ! N/m2
          ustar_sea = sqrt(tau_z_abs/rho_air)      ! m/s
          !
          ! limitation to avoid division by zero:
          ustar_sea = max( 0.001, ustar_sea )    ! m/s ;  minium of 1 mm/s
          !
          ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling",
          !   section "8.4 Surface roughness lengths" :
          ! "For smooth surfaces, such as over a smooth ocean with low wind speeds (Garret):
          !   z0m = 0.11 nu_a / u*
          ! where nu_a is the kinematic viscosity of air.
          ! Over a rough ocean with high wind speeds:
          !   z0m = alpha_c (u*)**2 / grav
          ! which is the 'Charnock relation',
          ! where alpha_c ~ 0.016 is the 'Charnock constant'."
          !
          ! Here the sum of these two parameterizations is used,
          ! where the first part is dominant for u*<0.1 and the second for u*>0.1 ;
          ! minimum value for u* prevents division by zero:
          !
          sr_sea = alfa_garret * nu_air / ustar_sea + &
                   alfa_charnock * ustar_sea**2 / grav

          !
          ! LAND
          ! calculate the 'local' surface roughness for momentum that is consistent with 10 m wind
          ! and the Olsson data base, we assume that the windspeed at 75 m is independent of
          ! surface roughness
          !
          wind10m = sqrt( u10m_dat(region)%data(i,j,1)**2 + v10m_dat(region)%data(i,j,1)**2 ) 

          if ( xland > 0.0 ) then
            sr_land=max(1e-2,sr_ols_dat(region)%data(i,j,1))  ! occurs at Islands, etc.
            sr_help=min(sr_ecm_dat(region)%data(i,j,1),0.03)
            !fd ustar_land=vKarman*wind10m(i,j)/alog(10./sr_help) !ustar consistent with 'clipped' large scale roughness
            ustar_land=vKarman*wind10m/&
                                    alog(10./sr_help)*alog(75./sr_help)/alog(10./sr_land)
          else
            sr_land = 0.0
            ustar_land = 0.0
          endif

          ! interpolate between land and sea for mixed pixels:
          ustar_loc(region)%surf(i,j) = xland*ustar_land+(1-xland)*ustar_sea
          sr_mix   (region)%surf(i,j) = xland*   sr_land+(1-xland)*   sr_sea

        end do !i
      end do !j
      field = ustar_loc(region)%surf

      if (okdebug) call dd_field_statistics('ustar_loc',field,im(region),jm(region))
      !if (okdebug) call dumpfield(0,'ustar_loc'//c_time//'.hdf',field,'ustar_loc')

      deallocate(field)

      if (okdebug) write(*,*) 'end calc_ustar'
 
    end subroutine dd_calc_ustar
    
    
    ! ***
    

#ifndef without_convection
#ifndef without_diffusion
    !------------------------------------------------------------------
    ! Procedure:
    ! (1) calc_Kv: Diffusion coefficients full atmosphere from shear
    ! (2) pblhght: Determine the boundary layer height
    ! (3) difcoef: Replace diffusion coefficients in the boundary layer
    !------------------------------------------------------------------
    ! subroutine adjusted for TM5, Bram Bregman, August 2003.
    !-------------------------IO---------------------------------------

    subroutine bldiff( region, phlb, gph, t, q, m, pu, pv, dkg )

      use binas,   only: ae, cp_air, Rgas, grav, Lvap
      use binas,   only: xmair
      use binas,   only: vkarman
      use dims,    only: im, jm, lm, lmax_conv
      use dims,    only: okdebug
      use dims,    only: gtor, dx, dy, ybeg, xref, yref
      use dims,    only: isr, ier, jsr, jer
      use toolbox, only: dumpfield
      use global_data, only : conv_dat
      use meteodata  , only : slhf_dat, sshf_dat
      use omp_partools,only : my_omp_domain

      ! --- in/out ----------------------------------------------
      
      integer,intent(in)              :: region
      real,dimension(:,:,:),pointer   :: phlb, gph, t, q, m, pu, pv, dkg

      ! --- const ----------------------------------
      
      ! minimum value diffusion coefficient
      real, parameter          :: ckmin = 1.e-15

      real, parameter          :: sffrac=  0.1
      real, parameter          :: onet   = 1./3.
      real, parameter          :: rair   = Rgas*1.e3/xmair
      real, parameter          :: ccpq   = 1869.46/cp_air -1.
      real, parameter          :: epsilo = rair/461.51
      real, parameter          :: apzero = 100000. ! reference pressure (usually 1000 hPa)
      real, parameter          :: ccon   = sffrac*vkarman
      real, parameter          :: betam=15.,betah=15.
      real, parameter          :: binm   = betam*sffrac
      real, parameter          :: binh   = betah*sffrac
      real, parameter          :: fr_excess = 1.
      real, parameter          :: fak = fr_excess * 8.5
      real, parameter          :: fakn = fr_excess * 7.2

      real, parameter          :: zkappa = rair / cp_air
      real, parameter          :: zcrdq = 1.0/epsilo - 1.0

      real, parameter          :: tiny     = 1.e-9   ! bound wind**2 to avoid dividing by 0
      real, parameter          :: zacb     = 100.    !factor in ustr-shearproduction term
      real, parameter          :: ricr     = 0.3     ! critical richardson number

      real, parameter          :: ssfrac=0.1

      ! --- local ------------------------------------

      real, dimension(:,:,:), allocatable:: wkvf,&     ! outer layer Kv (at the top of each layer)
                                            zdup2,&    ! vertical shear squared
                                            zrinub,&   ! Richardson number
                                            wthm,&     ! potential temperature
                                            zthvk,&    ! virtual potential temperature
                                            pf,&       ! full level pressure
                                            whm, &     ! height of layer centers
                                            uwind,vwind,& ! wind velocities [m s-1]
                                            kvh                     ! 

      real, dimension(:,:), allocatable  :: ustr,&     ! velocity scale     
                                            wheat,&    ! surface sensible heat flux [K m/s]
                                            wqflx,&    ! surface latent heat flux [kg m/ (kg s)]
                                            wobkl,&    ! Monin Obukhov length
                                            wheatv,&   ! surface virtual heat flux [K m/s]
                                            wtseff   ! theta_v at lowest level [K]

      integer             :: i, j, l
      real                :: intfac   !cmk BUG was integer!

      ! variables for the calculation of free atmosphere vertical diffusion coefficients wkvf
      real                        :: cml2,zlambdac,arg,zdz,zthva,zfunst, zfstabh,zsstab,zkvn,dyy
      real,dimension(jm(region))  :: dxx,lat
      integer                     :: jq! jqif function

      ! variabales for the calculation of the boundary layer height
      real                        :: wrino      ! Richardson number
      real                        :: thvparc    ! parcel virtual potential temperature
      real                        :: fmt,wsc,tkv,zrino,vvk,dtkv,zexcess
      integer                     :: jwork 
      real,dimension(:,:),pointer   :: pblh         ! planetary boundary layer height [m]

      ! variables to calculate difusion coefficient in the boundary layer
      !real                        :: wsc,cml2
      real                        :: z,zwstr,zkvh,kvhmin,zm,zp,zslask,zsl1, &
                                     zh,zzh,zl,pr,zpblk,zfstabh1,zfstabh2,&
                                     term1,term2,term3,obukov
      integer                     :: jq1,jq2,jq3,jq4,jq5,jck,jqq, jq6, jq7

      real                        :: yres
      real                        :: xres

      real                        :: thvgrad, kvhentr
      ! Openmp parameters
      integer                     :: js, je, jss, jee
      
      ! --- begin ------------------------------------------

      if (okdebug) write(*,*) 'start of bldiff'

      allocate(wkvf(im(region),jm(region),lm(region))) 
      allocate(zdup2(im(region),jm(region),lm(region))) 
      allocate(zrinub(im(region),jm(region),lm(region)))
      allocate(wthm(im(region),jm(region),lm(region))) 
      allocate(zthvk(im(region),jm(region),lm(region)))
      allocate(pf(im(region),jm(region),lm(region))) 
      allocate(whm(im(region),jm(region),lm(region)))
      allocate(uwind(im(region),jm(region),lm(region))) 
      allocate(vwind(im(region),jm(region),lm(region))) 
      allocate(kvh(im(region),jm(region),lm(region))) 
      allocate(ustr(im(region),jm(region)))
      allocate(wheat(im(region),jm(region))) 
      allocate(wqflx(im(region),jm(region))) 
      allocate(wobkl(im(region),jm(region))) 
      allocate(wheatv(im(region),jm(region)))
      allocate(wtseff(im(region),jm(region)))
      
      pblh => conv_dat(region)%blh

      ! mid-layer pressure levels
      ! potential temperature field
      ! virtual potential temperature field
      !
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, pf, phlb, t, q, wthm, &
!$OMP           zthvk, whm, gph) &
!$OMP  private (i, j, l, jss, jee, intfac)
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do l=1,lm(region)
         do j=jss, jee
            do i=isr(region), ier(region)
               pf(i,j,l)    = (phlb(i,j,l) + phlb(i,j,l+1) )*0.5
               wthm(i,j,l)  = t(i,j,l) * ( apzero/pf(i,j,l) )**zkappa
               zthvk(i,j,l) = wthm(i,j,l) * (1.0 + zcrdq*q(i,j,l))

               ! mid-layer heights (note: gph array is for levels 1:lm+1) 

               intfac = (phlb(i,j,l)-pf(i,j,l)) / (phlb(i,j,l)-phlb(i,j,l+1))
               ! cmk: intfac is now real!
               whm(i,j,l) = (gph(i,j,l)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1)  -gph(i,j,1)) * intfac
               ! cmk bug corrected:
               ! cmk old: whm(i,j,l) = (gph(i,j,l+1)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1) -gph(i,j,1)) * intfac
            end do
         end do
      end do
!$OMP END PARALLEL

      ! convert units of sensible and latent heat fluxes
      ! sshf: W/m2 / (J/kg/K) / (kg/m3) =  J/s/m2 /J kg K /kg m3 =  K m/s
      ! slhf: W/m2 / (J/kg) / (kg/m3) = J/s/m2 /J kg /kg m3 = m/s
      ! ustr is local velocity scale,  NOT from ECMWF surface stresses
      !
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (region, wheat, sshf_dat, gph, phlb, slhf_dat, &
!$OMP           wqflx, ustr, ustar_loc) &
!$OMP  private (i, j, js, je)
      call my_omp_domain (1, jm(region), js, je)
      do j = js, je
         do i = 1, im(region)
            wheat(i,j) = -sshf_dat(region)%data(i,j,1) * (gph(i,j,2) - &
                         gph(i,j,1)) * grav / (phlb(i,j,1) -           &
                         phlb(i,j,2)) / cp_air
            wqflx(i,j) = -slhf_dat(region)%data(i,j,1) * (gph(i,j,2) - &
                         gph(i,j,1)) * grav / (phlb(i,j,1) -           &
                         phlb(i,j,2)) / lvap
            ustr(i,j)  = max (ustar_loc(region)%surf(i,j), 0.01)
         end do
      end do
!$OMP END PARALLEL

      !
      !-----------------------------------------------------------------
      ! compute the free atmosphere vertical diffusion coefficients wkvf
      ! height-dependent asymptotic mixing length implemented
      ! Holtslag and Boville, 1993, J. Climate, 6, 1825-1842.
      !------------------------------------------------------------------

      ! first calculate wind velocities [m s-1] from the mass fluxes

      yres = dy/yref(region)
      xres = dx/xref(region)

      do j = 1,jm(region)
          lat(j) = ybeg(region) + 0.5 * yres + yres * (j-1)
      enddo

      dxx(:) = ae * xres * gtor * cos(lat(:)*gtor)
      dyy = ae * yres * gtor

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, dxx, dyy, pu, pv, m, &
!$OMP           uwind, vwind) &
!$OMP  private (i, j, l, jss, jee)
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do l=1,lm(region)
        do j=jss,jee
          do i=isr(region), ier(region)
             uwind(i,j,l) = dxx(j)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l)
             ! cmk Bug repared, 06/2005 yres --> dyy
             vwind(i,j,l) = dyy*  (pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l)
          enddo
        enddo
      enddo
!$OMP END PARALLEL

      !    
      ! Initialize gradient Richardson number and free tropospheric k values
      !

      if (okdebug) write(*,*) '  Richardson number'

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, gph, zdup2, uwind, &
!$OMP           vwind, whm, phlb, pf, zthvk, zrinub, wkvf) &
!$OMP  private (i, j, l, jss, jee, jq, zlambdac, cml2, zdz, arg, &
!$OMP           zthva, zsstab, zfunst, zfstabh, zkvn)
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do l=1,lm(region)-1
        do j=jss,jee
          do i=isr(region), ier(region)
            !
            ! if (delta gph < 1000) then jq=1
            jq = jqif( 1000.0, (gph(i,j,l+1)-gph(i,j,1)) )
            zlambdac=jq*300.+(1-jq)*(30. + 270. * exp (1.-(gph(i,j,l+1)-gph(i,j,1))/1000.))
            cml2=(1./( 1./zlambdac + 1./(vkarman*(gph(i,j,l+1)-gph(i,j,1))) ))**2.
            !
            !  vertical shear squared,
            !  min value of (delta v)**2 prevents zero shear
            !
            zdup2(i,j,l) = (uwind(i,j,l+1)-uwind(i,j,l))**2 + &
                                               (vwind(i,j,l+1)-vwind(i,j,l))**2
            zdup2(i,j,l) = max(zdup2(i,j,l),1.e-10)
            zdz   = whm(i,j,l+1) - whm(i,j,l)
            zdup2(i,j,l) = zdup2(i,j,l)/(zdz**2)
            !
            ! static stability (use virtual potential temperature)
            ! dv now calculated at interface height (lineair in pressure)
            !
            arg=(log(phlb(i,j,l+1))-log(pf(i,j,l)))/(log(pf(i,j,l+1))-log(pf(i,j,l)))
            zthva = zthvk(i,j,l) + arg * (zthvk(i,j,l+1)-zthvk(i,j,l))
            zsstab = grav/zthva*(zthvk(i,j,l+1) - zthvk(i,j,l))/zdz
            !
            ! gradient Richardson number
            !
            zrinub(i,j,l)  = zsstab/zdup2(i,j,l)
            !
            ! stability functions
            !
            zfunst = max(1. - 18.*zrinub(i,j,l),0.)
            zfunst = sqrt(zfunst)
            !
            zfstabh = 1.0 / (1.0 + 10.0*zrinub(i,j,l)*(1.0+8.0*zrinub(i,j,l)))
            !
            ! neutral diffusion coefficient
            !
            zkvn = cml2 * sqrt(zdup2(i,j,l))
            !
            ! correction with stability functions
            !
            jq = jqif( zrinub(i,j,l),0.0 )  ! thus: if (zrinub>0) then jq=1
            wkvf(i,j,l)=jq*(max(ckmin,zkvn*zfstabh))+(1-jq)*(max(ckmin,zkvn*zfunst))

          end do
        end do
      end do
!$OMP END PARALLEL
      !
      ! Determine the boundary layer height
      !--------------------------------------------------------------------------------
      ! based on: Vogelezang and Holtslag, 1996, Bound. Layer Meteorol., 81, 245-269.
      !--------------------------------------------------------------------------------

      if (okdebug) write(*,*) '  boundary layer height'
      
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, wtseff, wthm, q, &
!$OMP           wheatv, wheat, wqflx, wobkl, ustr, pblh, uwind, vwind, &
!$OMP           whm ) &
!$OMP  private (i, j, l, jss, jee, jwork, wrino, vvk, tkv, dtkv, &
!$OMP           zrino, jq, fmt, wsc, zexcess, thvparc )
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do j=jss,jee
        do i=isr(region),ier(region)
          ! compute bottom level virtual temperature and heat flux and Monin-Obukhov Length
          wtseff(i,j)  = wthm(i,j,1)*(1. + zcrdq*q(i,j,1))
          wheatv(i,j)  = wheat(i,j) + zcrdq*wthm(i,j,1)*wqflx(i,j)
          wobkl(i,j)   = -wtseff(i,j)*ustr(i,j)**3 / (grav*vkarman*(wheatv(i,j)+sign(1.e-10,wheatv(i,j))) )
          ! initialize
          pblh(i,j)   = 0.0
          jwork  = 1
          wrino  = 0.0

          do l=1,lm(region)

            vvk   = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + &
                                     zacb*ustr(i,j)*ustr(i,j)
            vvk   = max(vvk,tiny)

            tkv   = wthm(i,j,l)*(1. + zcrdq*q(i,j,l))
            dtkv  = tkv-wtseff(i,j)
            !
            ! Bulk Richardson number
            !
            zrino = grav*dtkv * (whm(i,j,l)-whm(i,j,1)) / (wtseff(i,j)*vvk)
            zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero
            zrino = jwork*zrino

            jq    = jqif(ricr,zrino)  ! thus: if (zrino < ricr) then jq=1

            !
            ! calculate pblh from linear interpolation in Ri(bulk)
            !
            if ( l == 1 ) then
              pblh(i,j) = jq*pblh(i,j) + (1-jq) *  &
                    (                (ricr-wrino)/(zrino-wrino)*(whm(i,j,l)             ) )
            else
              pblh(i,j) = jq*pblh(i,j) + (1-jq) *  &
                    ( whm(i,j,l-1) + (ricr-wrino)/(zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1)) )
            end if
            !
            ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0 
            ! and pblh does not change anymore
            !
            jwork = jq*jwork
            !
            ! if pblh is found already avoid dividing by zero next level by a faked value
            ! of wrino (jwork=0 already)
            !
            wrino = zrino+(1-jwork)*0.1

          end do

          pblh(i,j)  = (1-jwork)*pblh(i,j) + whm(i,j,lm(region))*jwork

          !
          ! second calculation of pblh including excess surface temperature using
          ! a convective velocity scale wsc (wm, not equal to w*!!!), using result
          ! of the first calculation of pblh
          !
          jq         = jqif(wheatv(i,j),0.) ! thus: if (wheatv(i,j) > 0.) then jq=1
          jwork = jq
          fmt = ( jq*(1.0 - binm*pblh(i,j)/wobkl(i,j)) + (1.0-jq) )**onet
          wsc     =  ustr(i,j)*fmt
          zexcess =  wheatv(i,j)*fak/wsc
          !
          ! improve pblh estimate under convective conditions using
          ! convective temperature excess (zexcess)
          !
          thvparc = wtseff(i,j)+jq*zexcess
          !
          jwork  = 1
          wrino  = 0.0
          !
          do l=2,lm(region)
            !
            ! Bulk Richardson number:
            !  if stable corrected with extra shear term
            !  if unstable now also corrected for temperature excess
            !
            vvk   = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + &
                                      zacb*ustr(i,j)*ustr(i,j)
            vvk   = max(vvk,tiny)
            tkv   = wthm(i,j,l)*(1. + zcrdq*q(i,j,l))
            zrino = grav*(tkv-thvparc) * (whm(i,j,l)-whm(i,j,1))/(vvk*wtseff(i,j))
            zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero
            zrino = zrino*jwork

            jq    = jqif(ricr,zrino)  ! thus: if (zrino < ricr) then jq=0
            !
            ! calculate pblh from linear interpolation in Ri(bulk)
            !
            pblh(i,j)  = jq*pblh(i,j)+(1-jq)* (whm(i,j,l-1)+(ricr-wrino)/ &
                                                (zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1)))
            !
            ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0 
            ! and pblh does not change anymore
            !
            jwork = jq*jwork
            !
            ! if pblh is found already avoid dividing by zero next level by a faked value
            ! of wrino (jwork=0 already)
            !
            wrino = zrino+(1-jwork)*0.1

          end do

          pblh(i,j) = (1-jwork)*pblh(i,j) + jwork*whm(i,j,lm(region))
          pblh(i,j) = max(pblh(i,j), 100.0) ! set minimum value for pblh

        end do
      end do
!$OMP END PARALLEL


      !
      ! Diffusion coefficients in the surface layer
      !
      !--------------------------------------------------------------------------------
      ! the revised LTG scheme (Beljaars and Viterbo, 1998)
      ! Modification - April 1, 1999: prevent zpblk to become too small in very stable conditions
      !--------------------------------------------------------------------------------

      if (okdebug) write(*,*) '  diffusion coeff'

      ! initialize output array with minimum value
      kvhmin = 0.1
      !
      ! Loop to calculate the diffusivity
      !
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, pblh, gph, wobkl, &
!$OMP           wheatv, whm, zrinub, zdup2, kvhmin, wkvf, ustr, &
!$OMP           wtseff, kvh, zthvk) &
!$OMP  private (i, j, l, jss, jee, jq, z, zh, zl, zzh, jq1, jq2, jq3, &
!$OMP           jq4, jq5, jq6, jq7, jqq, cml2, zfstabh, zpblk, zkvh, &
!$OMP           zslask, fmt, wsc, zwstr, obukov, term1, term2, term3, &
!$OMP           pr, thvgrad, kvhentr)
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do l=1,lm(region)-1
        do j=jss,jee
          do i=isr(region), ier(region)

            jq  = jqif(pblh(i,j),(gph(i,j,l+1)-gph(i,j,1))) ! if top of level hh < pblh jq=1 inside BL
            z   = gph(i,j,l+1)-gph(i,j,1)
            zh  = jq*z/pblh(i,j) ! only defined in BL (jq=1)

            ! zl must not be zero
            zl  = jq*z/wobkl(i,j)+(1-jq) ! z/L outside BL = 1
            zzh = jq*(1.-zh)**2           ! (1-z/h)^2 only in BL

            jq1 = jqif(wheatv(i,j),0.)    ! jq1 is 1 in unstable case
            jq2 = jq*(1-jq1)              ! jq2 is 1 in stable BL
            jq3 = jq*jq1                  ! jq3 is 1 in unstable BL

            jq6 = jqif(pblh(i,j), whm(i,j,l))
            jq7 = jqif(whm(i,j,l+1), pblh(i,j))
            jq6 = jq1*jq6*jq7  ! jq6 is 1 at at the kvh-level that is located between the
                               ! mid-levels that surround pblh

            jq1 = jqif(sffrac,zh)         ! jq1 is 1 in surface layer
            jq4 = jq3*jq1                 ! jq4 is 1 in unstable surface layer
            jq5 = jq3*(1-jq1)             ! jq5 is 1 outside unstable surface layer
            !
            ! calculate coefficients for momentum, the values for heat are obtained by dividing by Prantl number
            !
            ! stable and neutral:
            !
            cml2    = (1./( 1./(vkarman*z) + 1./450. ))**2.
            jqq     = jqif(zrinub(i,j,l),0.)
            !
            ! stability function for momentum
            !
            zfstabh  = jqq/(1.+10.*zrinub(i,j,l)*sqrt(1.+jqq*zrinub(i,j,l))) + (1-jqq)

            zpblk   = cml2 * sqrt(zdup2(i,j,l))*zfstabh
            zkvh    = jq2*max(zpblk,kvhmin)+(1-jq2)*wkvf(i,j,l) ! take free troposp. values outside BL
            !
            ! unstable case in surface layer
            !
            zslask = 1.-betah*zl
            zslask = jq4*zslask+(1-jq4)
            zpblk  = (1-jq4)*zkvh+jq4*(ustr(i,j)*vkarman*z*zzh*zslask**(1./3.))
            !
            ! unstable case above surface layer
            !
            ! evaluate stability function for momentum at height z = 0.1*pblh (top of surface layer)
            !
            zslask = 1.-ssfrac*betah*pblh(i,j)/wobkl(i,j)
            zslask = jq5*zslask+(1-jq5)
            fmt    = zslask**(1./3.)
            !
            ! use surface layer definition for w_m
            wsc    = ustr(i,j)*fmt
            zpblk  = (1-jq5)*zpblk + jq5*wsc*vkarman*z*zzh
            !
            ! Determine Prandt Number
            !
            ! Pr-number is assumed to be constant throughout the CBL, i.e. in surface and mixed layer
            ! NOTE: this is different from the original formulation
            ! zwstr not used if wheatv<0.
            !
            zwstr  = (abs(wheatv(i,j))*grav*pblh(i,j)/wtseff(i,j))**(1./3.)         ! bh

            obukov = jq3*wobkl(i,j)-(1-jq3)
            term1 = (1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.)
            term2 = sqrt(1.-ssfrac*betah*pblh(i,j)/obukov)
            term3 = ccon*fakn*zwstr/(ustr(i,j)*(1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.))

            pr = jq3*(term1/term2+term3)+(1-jq3)
            !
            ! NOTE that kvh applies to the top of the level
            !
            kvh(i,j,l) = jq3*max(zpblk/pr,kvhmin)+(1-jq3)*zkvh

            ! if in the entrainment zone, override with prescribed entrainment
            thvgrad = (zthvk(i,j,l+1)-zthvk(i,j,l))/(whm(i,j,l+1)-whm(i,j,l))
            kvhentr = 0.2*wheatv(i,j)/thvgrad
            kvh(i,j,l) = jq6*kvhentr + (1-jq6)*kvh(i,j,l)

          end do
        end do
      end do
!$OMP END PARALLEl

      ! calculate vertical diffusion

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, dkg, kvh, m, gph) &
!$OMP  private (i, j, l, jss, jee)
      call my_omp_domain (jsr(region), jer(region), jss, jee)
      do  l=1,lmax_conv-1
        do  j=jss,jee
          do  i=isr(region), ier(region)
            dkg(i,j,l)=max(0.,kvh(i,j,l))*2.*(m(i,j,l+1)+m(i,j,l))/(gph(i,j,l+2)-gph(i,j,l))**2
            ! CMK dec 2002 ----> gph changed to 1--->lm+1 boundaries
          enddo
        enddo
      enddo
      dkg(:,jss:jee,lmax_conv) = 0.0
!$OMP END PARALLEL

      deallocate(wkvf)
      deallocate(zdup2)
      deallocate(zrinub)
      deallocate(wthm)
      deallocate(zthvk)
      deallocate(pf)
      deallocate(whm)
      deallocate(ustr)
      deallocate(wheat)
      deallocate(wqflx)
      deallocate(wobkl)
      deallocate(wheatv)
      deallocate(wtseff)
      deallocate(uwind)
      deallocate(vwind)
      deallocate(kvh)

      nullify(pblh)

      if (okdebug) write(*,*) 'start of bldiff'

    end subroutine bldiff
#endif
#endif
    
!  end subroutine calc_kzz


  ! ***


  !-----------------------------------------------
  !
  ! Purpose:
  ! -------
  ! this subroutine calculate the min,mean,max of a real field
  !
  ! External
  ! --------
  ! none
  !
  ! Reference
  ! ---------
  ! none
  !------------------------------------

  subroutine dd_field_statistics(name,field,ix,jx)

    ! --- in/out -------------------------------
    
    character(len=*),intent(in)        :: name
    integer,intent(in)                 :: ix,jx
    real,dimension(ix,jx),intent(in)   :: field
    
    ! --- local --------------------------------
    
    integer :: i,j,ntel_non_zero
    real :: maxf,minf,meanf,mean_non_zero
    
    ! --- begin -------------------------------
    
    maxf=-1e20
    minf=1e12
    meanf=0.
    mean_non_zero=0.
    ntel_non_zero=0
    do i=1,ix
    do j=1,jx
      meanf=meanf+field(i,j)
      maxf=max(maxf,field(i,j))
      minf=min(minf,field(i,j))
      if (field(i,j).ne.0) then
        ntel_non_zero=ntel_non_zero+1
        mean_non_zero=mean_non_zero+field(i,j)
      endif
    enddo
    enddo
    meanf=meanf/ix/jx
    if (ntel_non_zero.gt.0) mean_non_zero= mean_non_zero/ntel_non_zero
    write(6,'(a10,4(a3,1x,1pe10.3))') name,'max',maxf,'min',minf,'mean',meanf,'mn0',mean_non_zero

  end subroutine dd_field_statistics


  ! ****


  !------------------------------------
  !
  ! Purpose:
  ! -------
  ! this subroutine calculate the min,mean,max of a field of int8
  !
  ! External
  ! --------
  ! none
  !
  ! Reference
  ! ---------
  ! none
  !------------------------------------

  subroutine dd_field_statistics_i8(name,field,ix,jx)

    ! --- in/out ----------------------------

    character(len=*),intent(in)                   :: name
    integer,intent(in)                            :: ix,jx
    integer(kind=1),dimension(ix,jx),intent(in)   :: field

    ! --- local ----------------------------

    integer :: i,j,ntel_non_zero
    real :: maxf,minf,meanf,mean_non_zero
    
    ! --- begin ---------------------------
    
    maxf=-1e20
    minf=1e12
    meanf=0.
    mean_non_zero=0.
    ntel_non_zero=0
    do i=1,ix
      do j=1,jx
        meanf=meanf+field(i,j)
        maxf=max(maxf,real(field(i,j)))
        minf=min(minf,real(field(i,j)))
        if (field(i,j).ne.0) then
          ntel_non_zero=ntel_non_zero+1
          mean_non_zero=mean_non_zero+field(i,j)
        endif
      enddo
    enddo
    meanf=meanf/ix/jx
    if (ntel_non_zero.gt.0) mean_non_zero= mean_non_zero/ntel_non_zero
    write(6,'(a10,4(a3,1x,1pe10.3))') name,'max',maxf,'min',minf,'mean',meanf,'mn0',mean_non_zero

  end subroutine dd_field_statistics_i8
  
  
  ! ***


  !
  ! define vectorizable 'if' function
  ! if xxxz>yyyz jqif=1 else jqif=0
  !

  integer function jqif(xxxz,yyyz)

    ! --- in/out ---------------------------
    
    real, intent(in)    ::  xxxz, yyyz

    ! --- begin ----------------------------
    
    !
    !jqif = nint( 0.5 - sign( 0.5, -dim(xxxz,yyyz)) )
    !
    jqif = nint( 0.5 + sign( 0.5, xxxz-yyyz ) )

  end function jqif

  !=====================================================================================================
  !=====================================================================================================   

  subroutine read_diffusion( status )

    !
    ! Read vertical diffusion from file and store in conv_dat%dkg
    !
    ! Output:
    !  o status      0 = ok
    !               -1 = file not available for at least one region
    !                1 = read error
    !

    ! --- modules ------------------------------

    use dims,                    only : itau, revert, region_name, lm
    use global_data,             only : conv_dat
    use datetime,                only : tau2date
    use file_hdf,                only : Init, ReadData, Done, THdfFile, TSds
    use toolbox,                 only : escape_tm

    ! --- in/out ----------------------------------------------

    integer, intent(out)             :: status

    ! --- const -----------------------------------------------

    character(len=*), parameter      :: rname = mname//', read_diffusion'

    ! --- local -----------------------------------------------

    character(len=1000)              :: fname
    integer                          :: region
    integer, dimension(6)            :: idater
    logical                          :: exist
    type(THdfFile)                   :: hdf
    type(TSds)                       :: sds

    ! --- begin -----------------------------------------------

    if ( revert == 1 ) then
       call tau2date( itau, idater )
    else
       ! read field from 3 hours before
       call tau2date( itau - 3600*3, idater )
    end if

    do region = 1, nregions

       ! Create file name
       write(fname,'(2a,i4.4,a,i2.2,a,i2.2,3a,i2,a,i4,5i2.2,a)') trim(diffusion_dir), & 
            '/',idater(1),'/',idater(2),'/dkg/',idater(3),'/dkg_', &
            trim(region_name(region)), '_z', lm(region), '_', idater, '.hdf'

       inquire( file=fname, exist=exist )

       if ( .not. exist ) then
          write(*,'(a,": dkg file not found for region",i2)') rname, region
          status = -1
          return
       end if

       ! Open file
       ! This is a bit too verbose, leading to long output files.
       !       write(*,'(a,": reading ", a)') rname, trim(fname)

       call Init( hdf, fname, 'read', status )
       IF_NOTOK_RETURN(status=1)

       ! Read dkg
       call Init( sds, hdf, 'dkg', status )
       IF_NOTOK_RETURN(status=1)

       call ReadData( sds, conv_dat(region)%dkg ,status)

       call Done( sds, status )
       IF_NOTOK_RETURN(status=1)

       ! Read blh
       call Init( sds, hdf, 'blh', status )
       IF_NOTOK_RETURN(status=1)

       call ReadData( sds, conv_dat(region)%blh , status)

       call Done( sds, status )
       IF_NOTOK_RETURN(status=1)

       ! Close file
       call Done( hdf, status )
       IF_NOTOK_RETURN(status=1)

    end do

    status = 0

  end subroutine read_diffusion

  !=====================================================================================================
  !=====================================================================================================   

  subroutine write_diffusion( status )

    !
    ! Write vertical diffusion to files (for all regions)
    !
    ! Output:
    !  o status      0 = ok, else not ok
    !

    ! --- modules ------------------------------

    use dims,                    only : itau, revert, region_name, lm
    use global_data,             only : conv_dat
    use datetime,                only : tau2date
    use file_hdf,                only : Init, WriteData, Done, THdfFile, TSds
    use toolbox                , only : escape_tm
    use Partools,                only : myid, root

    ! --- in/out ----------------------------------------------

    integer, intent(out)             :: status

    ! --- const -----------------------------------------------

    character(len=*), parameter      :: rname = mname//', write_diffusion'

    ! --- local -----------------------------------------------

    character(len=1000)               :: fname,targetdir
    integer                          :: region
    integer, dimension(6)            :: idater
    type(THdfFile)                   :: hdf
    type(TSds)                       :: sds

    ! --- begin -----------------------------------------------

    if ( revert == 1 ) then
       call tau2date( itau, idater )
    else
       ! read field from 3 hours before
       call tau2date( itau - 3600*3, idater )
    end if

    status=0

    if (myid /= root) return !WP! only root PE writes fields

    do region = 1, nregions

       ! Create file name
       write(fname,'(2a,i4.4,a,i2.2,a,i2.2,3a,i2,a,i4,5i2.2,a)') trim(diffusion_dir), & 
            '/',idater(1),'/',idater(2),'/dkg/',idater(3),'/dkg_', &
            trim(region_name(region)), '_z', lm(region), '_', idater, '.hdf'

       ! Create directory if necessary
       write(targetdir,'(2a,i4.4,a,i2.2,a,i2.2)') trim(diffusion_dir), & 
            '/',idater(1),'/',idater(2),'/dkg/',idater(3)
       call system('mkdir -p ' // targetdir)

       ! Create file
       write(*,'(a,": creating ", a)') rname, trim(fname)
       call Init( hdf, fname, 'create', status )
       IF_NOTOK_RETURN(status=1)

       ! Write dkg
       call Init( sds, hdf, 'dkg', shape(conv_dat(region)%dkg), 'real(8)', status )
       IF_NOTOK_RETURN(status=1)

       call WriteData( sds, conv_dat(region)%dkg ,status)

       call Done( sds, status )
       IF_NOTOK_RETURN(status=1)


       ! Write blh
       call Init( sds, hdf, 'blh', shape(conv_dat(region)%blh), 'real(8)', status )
       IF_NOTOK_RETURN(status=1)

       call WriteData( sds, conv_dat(region)%blh,status )

       call Done( sds, status )
       IF_NOTOK_RETURN(status=1)

       ! Close file
       call Done( hdf, status )
       IF_NOTOK_RETURN(status=1)

    end do

    status = 0

  end subroutine write_diffusion

  !=====================================================================================================
  !=====================================================================================================   

  
end module diffusion

