ecearth3/sources/tm5mp/base/diffusion.F90

!#################################################################
!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 minimum value for u* to trap division by zero (MK).
!
! 2012-01 Philippe Le Sager
!   Adapted for Lat-Lon MPI domain decomposition
!
!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
  private

  ! methods
  public  :: Diffusion_Init, Diffusion_Done
  public  :: Calc_Kzz, Read_Diffusion, Write_Diffusion

  ! flag
  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 MeteoData,   only : Set
    use MeteoData,   only : pu_dat,     pv_dat
    use MeteoData,   only : temper_dat, humid_dat, gph_dat
    use MeteoData,   only : oro_dat,    lsmask_dat
    use MeteoData,   only : sr_ecm_dat, sr_ols_dat
    use MeteoData,   only : u10m_dat,   v10m_dat
    use MeteoData,   only : slhf_dat,   sshf_dat
    use MeteoData,   only : ewss_dat,   nsss_dat
    use GO,          only : TrcFile, Init, Done, ReadRc


    ! --- 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( status )

    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 : 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 meteoData,    only : spm_dat, pu_dat, pv_dat, m_dat
    use meteoData,    only : temper_dat, humid_dat, gph_dat
    use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo, dist_arr_stat

    integer, intent(out) :: status
    
    ! --- 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(:,:,:)
    
    integer :: i1, i2, j1, j2

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

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

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


    ! Allocate work arrays, and fill ghost cells
    DO region = 1, nregions

       CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       ALLOCATE(ustar_loc(region)%surf( i1:i2, j1:j2))
       ALLOCATE(   sr_mix(region)%surf( i1:i2, j1:j2))

       CALL UPDATE_HALO( dgrid(region), wind_dat(region)%am, wind_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)
       
       CALL UPDATE_HALO( dgrid(region), wind_dat(region)%bm, wind_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)

       CALL UPDATE_HALO( dgrid(region), wind_dat(region)%cm, wind_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)

       CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)

       CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)
       
    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

    REG: do region = 1,nregions

       ! local grid sizes
       lmr=lm(region)
       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       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   (i1:i2, j1:j2, lmr  ) )
          allocate ( phlb_adj(i1:i2, j1:j2, lmr+1) )
          allocate ( p_adj   (i1:i2, j1:j2       ) )
           
          am   => wind_dat(region)%am
          bm   => wind_dat(region)%bm
          cm   => wind_dat(region)%cm
          dxyp => region_dat(region)%dxyp
          
          ntimes = (nread/(ndyn))*tref(region)
          m_adj  = m_dat(region)%data(i1:i2, j1:j2, 1:lmr)
           
          do i=1,ntimes
             m_adj = m_adj + revert*am(i1-1:i2-1, j1  :j2  , 1:lmr  ) &   !east   
                           - revert*am(i1  :i2  , j1  :j2  , 1:lmr  ) &   !west
                           + revert*bm(i1  :i2  , j1  :j2  , 1:lmr  ) &   !south
                           - revert*bm(i1  :i2  , j1+1:j2+1, 1:lmr  ) &   !north
                           + revert*cm(i1  :i2  , j1  :j2  , 0:lmr-1) &   !lower
                           - revert*cm(i1  :i2  , j1  :j2  , 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 = j1, j2
             do i = i1, i2
                   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 = j1, j2
             do i = i1, i2
                   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(i1:i2, j1:j2, 1:lmr+1) )

          do l = 1, lmr+1
             phlb_t(:,:,l) = at(l) + bt(l)*spm_dat(region)%data(i1:i2, j1:j2, 1)
          end do

          allocate( m_t(i1-2:i2+2, j1-2:j2+2, 1:lmr) )
          m_t = 0.0

          do l = 1, lmr
             do j = j1, j2
                m_t(i1:i2,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, i1, i2, j1, j2)


       if (okdebug) then
!          call dd_field_statistics(region, 'ustar_loc', ustar_loc(region)%surf, i1, j1, status)
          call dist_arr_stat(dgrid(region), 'ustar_loc', ustar_loc(region)%surf, 0, status)
       end if
       
#ifndef without_convection
#ifndef without_diffusion
       ! calculate kzz                                                         
       call bldiff( region, phlb, m, i1, i2, j1, j2 )  
#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 REG

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

    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).
  !
  ! 2012-01 P. Le Sager
  !   adapted for lon-lat MPI domain decomposition
  !------------------------
  
  subroutine dd_calc_ustar(region, i1, i2, j1, j2)

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

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

      integer,intent(in) :: region, i1, i2, j1, j2
      
      ! --- 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

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

      do j=j1, j2
       do i=i1, i2
          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
            !>>> TvN
            ! Over land, the friction velocity ustar 
            ! is calculated from the 10 m wind speed, u10.
            ! In IFS u10 is a diagnostic variable, calculated to be
            ! compatible with "open-terrain" wind speed observations.
            ! Over rough or inhomogeneous terrain (z0M > 0.03 m),
            ! it is calculated using an aerodynamic roughness
            ! length typical for open terrain with grassland (0.03 m)
            ! (see IFS cy31r1 documentation, p. 46).

            ! In the expressions applied in TM5,
            ! the stability profile functions Psi_M have disappeared,
            ! and only the logarithmic terms are kept.
            ! Apart from this, the expression for ustar was incorrect:
            ! 10./sr_land should be 75./sr_land.
            ! This has now been corrected.

            ! Moreover, over islands and near coast lines, 
            ! zeros can occur in sr_ols_dat,
            ! which have to be removed.
            ! However, a lower bound of 1e-2 m = 1 cm
            ! seems too high for smooth surfaces,
            ! like deserts and ice caps.
            ! The minimum positive value in sr_ols_dat
            ! is 2.5e-4 m = 0.025 cm,
            ! which is obtained in parts of Antarctica.
            ! This is likely the result of regridding of
            ! a field with minimum value of 0.1 cm
            ! from 0.5x0.5 to 1x1 degrees.

            ! Please note that with the introduction of CY31R1,
            ! the orographic contribution to the aerodynamic roughness length
            ! in IFS has been replaced by an explicit specification of stress
            ! on model levels due to turbulent orographic form drag,
            ! and the climatological variable previously used
            ! has been replaced by a prognostic variable.
            ! In TM5, the climatological variable is still used,
            ! but it would be better to use the prognostic variable instead.

            ! Note also that the same calculation 
            ! is done in dry_deposition.F90.

            !sr_land=max(1e-2,sr_ols_dat(region)%data(i,j,1))  ! occurs at Islands, etc.
            sr_land=max(1e-3,sr_ols_dat(region)%data(i,j,1))
            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)
            ustar_land=vKarman*wind10m/&
                                    alog(10./sr_help)*alog(75./sr_help)/alog(75./sr_land)
            ! <<< TvN
          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


    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, m, i1, i2, j1, j2 )

    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, pu_dat, pv_dat, gph_dat, temper_dat, humid_dat

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

    integer,intent(in)              :: region, i1, i2, j1, j2
    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(j1:j2)      :: 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

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

    if (okdebug) write(*,*) 'start of bldiff'
    
    allocate(wkvf  ( i1:i2, j1:j2, lm(region))) 
    allocate(zdup2 ( i1:i2, j1:j2, lm(region))) 
    allocate(zrinub( i1:i2, j1:j2, lm(region)))
    allocate(wthm  ( i1:i2, j1:j2, lm(region))) 
    allocate(zthvk ( i1:i2, j1:j2, lm(region)))
    allocate(pf    ( i1:i2, j1:j2, lm(region))) 
    allocate(whm   ( i1:i2, j1:j2, lm(region)))
    allocate(uwind ( i1:i2, j1:j2, lm(region))) 
    allocate(vwind ( i1:i2, j1:j2, lm(region))) 
    allocate(kvh   ( i1:i2, j1:j2, lm(region))) 
    allocate(ustr  ( i1:i2, j1:j2            ))
    allocate(wheat ( i1:i2, j1:j2            )) 
    allocate(wqflx ( i1:i2, j1:j2            )) 
    allocate(wobkl ( i1:i2, j1:j2            )) 
    allocate(wheatv( i1:i2, j1:j2            ))
    allocate(wtseff( i1:i2, j1:j2            ))

    pblh => conv_dat(region)%blh

    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
    dkg => conv_dat(region)%dkg
    

    ! mid-layer pressure levels
    ! potential temperature field
    ! virtual potential temperature field
    !
    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (lm, jsr, jer, isr, ier, region, pf, phlb, t, q, wthm, &
    !PLS $OMP           zthvk, whm, gph) &
    !PLS $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=j1,j2
          do i=i1,i2

             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
    !PLS $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
    !
    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (jm, im, region, wheat, sshf_dat, gph, phlb, slhf_dat, &
    !PLS $OMP           wqflx, ustr, ustar_loc) &
    !PLS $OMP  private (i, j, js, je)
    !      call my_omp_domain (1, jm(region), js, je)
    do j = j1,j2
       do i = i1,i2
          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
    !PLS $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 = j1,j2
       lat(j) = ybeg(region) + 0.5 * yres + yres * (j-1)
    enddo

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

    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (lm, jsr, jer, isr, ier, region, dxx, dyy, pu, pv, m, &
    !PLS $OMP           uwind, vwind) &
    !PLS $OMP  private (i, j, l, jss, jee)
    !      call my_omp_domain (jsr(region), jer(region), jss, jee)
    do l=1,lm(region)
       do j=j1,j2
          do i=i1,i2
             uwind(i,j,l) = dxx(j)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l)
             vwind(i,j,l) = dyy   *(pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l)
          enddo
       enddo
    enddo
    !PLS $OMP END PARALLEL

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

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

    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (lm, jsr, jer, isr, ier, region, gph, zdup2, uwind, &
    !PLS $OMP           vwind, whm, phlb, pf, zthvk, zrinub, wkvf) &
    !PLS $OMP  private (i, j, l, jss, jee, jq, zlambdac, cml2, zdz, arg, &
    !PLS $OMP           zthva, zsstab, zfunst, zfstabh, zkvn)
    !      call my_omp_domain (jsr(region), jer(region), jss, jee)
    do l=1,lm(region)-1
       do j=j1,j2
          do i=i1,i2
             !
             ! 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
    !PLS $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'

    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (lm, jsr, jer, isr, ier, region, wtseff, wthm, q, &
    !PLS $OMP           wheatv, wheat, wqflx, wobkl, ustr, pblh, uwind, vwind, &
    !PLS $OMP           whm ) &
    !PLS $OMP  private (i, j, l, jss, jee, jwork, wrino, vvk, tkv, dtkv, &
    !PLS $OMP           zrino, jq, fmt, wsc, zexcess, thvparc )
    !      call my_omp_domain (jsr(region), jer(region), jss, jee)
    do j=j1,j2
       do i=i1,i2
          ! 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
!PLS $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
    !
    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (lm, jsr, jer, isr, ier, region, pblh, gph, wobkl, &
    !PLS $OMP           wheatv, whm, zrinub, zdup2, kvhmin, wkvf, ustr, &
    !PLS $OMP           wtseff, kvh, zthvk) &
    !PLS $OMP  private (i, j, l, jss, jee, jq, z, zh, zl, zzh, jq1, jq2, jq3, &
    !PLS $OMP           jq4, jq5, jq6, jq7, jqq, cml2, zfstabh, zpblk, zkvh, &
    !PLS $OMP           zslask, fmt, wsc, zwstr, obukov, term1, term2, term3, &
    !PLS $OMP           pr, thvgrad, kvhentr)
    !      call my_omp_domain (jsr(region), jer(region), jss, jee)
    do l=1,lm(region)-1
       do j=j1,j2
          do i=i1,i2

             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
    !PLS $OMP END PARALLEl

    ! calculate vertical diffusion

    !PLS $OMP PARALLEL &
    !PLS $OMP  default (none) &
    !PLS $OMP  shared  (jsr, jer, isr, ier, region, dkg, kvh, m, gph) &
    !PLS $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=j1,j2
          do  i=i1, i2
             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(i1:i2,j1:j2,lmax_conv) = 0.0
    !PLS $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(*,*) 'end of bldiff'

  end subroutine bldiff
#endif /* diffusion */
#endif /* convection */


  !
  ! 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
    !

    !   19 Dec 2012 - P. Le Sager - TM5-MP version, using MDF

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

    use dims,               only : itau, revert, region_name, lm, okdebug
    use global_data,        only : conv_dat
    use datetime,           only : tau2date
    use MDF,                only : MDF_OPEN, MDF_HDF4, MDF_Inq_VarID, MDF_Get_Var, MDF_Close, MDF_READ
    use Partools,           only : isRoot
    use meteodata,          only : global_lli
    use tm5_distgrid,       only : dgrid, scatter

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

    integer, intent(out)        :: status

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

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

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

    character(len=1000)         :: fname
    integer                     :: region, fid, varid, imr, jmr, sz(3)
    integer, dimension(6)       :: idater
    logical                     :: exist
    real,    allocatable        :: global_arr(:,:,:)

    ! --- 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

       ! entire region grid size
       imr = global_lli(region)%nlon
       jmr = global_lli(region)%nlat

       ! for gathering
       sz = shape(conv_dat(region)%dkg)  ! to get 3rd dimension

       if(isRoot)then
          allocate(global_arr(imr,jmr,sz(3)))
       else
          allocate(global_arr(1,1,1))
       end if

       if (isRoot) then
          
          ! 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
             TRACEBACK; return
          end if

          ! Open file
          ! This is a bit too verbose, leading to long output files.
          !       write(*,'(a,": reading ", a)') rname, trim(fname)
          
          call MDF_OPEN(TRIM(fname), MDF_HDF4, MDF_READ, fid, status )       
          IF_NOTOK_RETURN(status=1)

          ! Read dkg
          CALL MDF_Inq_VarID( fid, 'dkg', varid, status )
          IF_ERROR_RETURN(status=1)

          CALL MDF_Get_Var( fid, varid, global_arr, status )
          IF_NOTOK_RETURN(status=1)

       end if

       call scatter( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status)
       IF_NOTOK_RETURN(status=1)
       
       if(isRoot) then
          ! Read blh
          CALL MDF_Inq_VarID( fid, 'blh', varid, status )
          IF_ERROR_RETURN(status=1)

          CALL MDF_Get_Var( fid, varid, global_arr(:,:,1), status )
          IF_NOTOK_RETURN(status=1)
       end if

       call scatter( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status)
       IF_NOTOK_RETURN(status=1)
       
       ! Close file
       if (isRoot) then
          call MDF_Close( fid, status )
          IF_NOTOK_RETURN(status=1)
       end if

       deallocate(global_arr)
       
    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
    !
    !   19 Dec 2012 - P. Le Sager - made a TM5-MP version , that works using MDF
    
    ! --- modules ------------------------------

    use dims,            only : itau, revert, region_name, lm, okdebug
    use global_data,     only : conv_dat
    use datetime,        only : tau2date
    use MDF,             only : MDF_CREATE, MDF_HDF4,  MDF_Def_Dim, MDF_Def_Var, MDF_Close
    use MDF,             only : MDF_DOUBLE, MDF_WRITE, MDF_Put_Var, MDF_EndDef
    use Partools,        only : isRoot
    use meteodata,       only : global_lli
    use tm5_distgrid,    only : dgrid, gather
    ! --- in/out ----------------------------------------------

    integer, intent(out)             :: status

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

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

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

    character(len=1000)      :: fname, targetdir
    integer                  :: region, fid, dkgid, blhid, dimlon, dimlat, dimlev, imr, jmr, sz(3)
    integer, dimension(6)    :: idater
    real,    allocatable     :: global_arr(:,:,:)

    ! --- 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


    do region = 1, nregions

       ! entire region grid size
       imr = global_lli(region)%nlon
       jmr = global_lli(region)%nlat
       
       ! for gathering
       sz = shape(conv_dat(region)%dkg)  ! to get 3rd dimension
       
       if(isRoot)then
          allocate(global_arr(imr,jmr,sz(3)))
       else
          allocate(global_arr(1,1,1))
       end if

       ! Create and define
       if(isRoot) then
          
          ! 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)

          if(okdebug) then
             write(gol,'(a,": creating ", a)') rname, trim(fname); call goPr
          end if

          call MDF_CREATE(TRIM(fname), MDF_HDF4, MDF_WRITE, fid, status )
          IF_NOTOK_RETURN(status=1)
       
          call MDF_Def_Dim( fid, 'lon',    imr,   dimlon, status )
          call MDF_Def_Dim( fid, 'lat',    jmr,   dimlat, status )
          call MDF_Def_Dim( fid, 'levels', sz(3), dimlev, status )
       
          call MDF_Def_Var( fid, 'dkg', MDF_DOUBLE, (/dimlon,dimlat,dimlev/), dkgid, status)
          call MDF_Def_Var( fid, 'blh', MDF_DOUBLE, (/dimlon,dimlat/),        blhid, status)
       
          call MDF_EndDef(fid, status)

       end if

       ! gather and write dkg/blh
       call gather( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) then 
          call MDF_Put_Var( fid, dkgid, global_arr, status)
          IF_NOTOK_RETURN(status=1)
       end if

       call gather( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status)
       IF_NOTOK_RETURN(status=1)

       if(isRoot) then
          call MDF_Put_Var( fid, blhid, global_arr(:,:,1), status)
          IF_NOTOK_RETURN(status=1)

          ! Close file
          call MDF_CLOSE( fid, status )
          IF_NOTOK_RETURN(status=1)

       end if

       deallocate(global_arr)
       
    end do

    status = 0

  end subroutine write_diffusion

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

END MODULE DIFFUSION