    subroutine surface_layer (kstep, nlon, nlat, h_atm, p_atm, ps, u_atm, v_atm, t_atm, q_atm, tvirt_atm, roscd, &
 rd, cpd, g, ep, rdrcp, tskin, qskin, alp0, fmask, fice, &
 rgm, rgq, rgmd, fsnow, ustar, tstar, &
 kturb_surf_m, kturb_surf_h, kturb_surf_q, hflux, qflux, &
 nstep_sl_filter, myid, kturb_surf_h_mem, kturb_surf_q_mem)

! Computes parameters and variables of the turbulent exchanging in the surface layer
! Computes heat and specific humidity fluxes at the ground
! Fluxes are positive upward

#ifdef chem
  use mod_bol2chem
#endif

! Input variables:

integer :: kstep, nlon, nlat, nstep_sl_filter, myid
real :: alp0
real, dimension(nlon,nlat) :: h_atm, p_atm, ps, u_atm, v_atm, t_atm, q_atm, tvirt_atm, roscd, tskin, qskin, &
 fmask, fice, rgm, rgq, rgmd, fsnow

! Ouput variables:

real, dimension(nlon,nlat) :: ustar, tstar, &
 kturb_surf_m, kturb_surf_h, kturb_surf_q, hflux, qflux
real, dimension(nlon,nlat,60) :: kturb_surf_h_mem, kturb_surf_q_mem

! Internal variables:

    integer, parameter :: n_std_lev_atm = 2
    real(4), parameter :: zak=.4, zgam=16., zaholt=1., zbholt=2./3., zcholt=5., zdholt=0.35
    real(4), dimension(n_std_lev_atm) :: psim_sl, psih_sl, psiq_sl
    real(4) :: ztetav_bottom, zlev_bottom
    real(4), dimension(nlon,nlat) :: kturb_surf_h_filter, kturb_surf_q_filter
    real(4), dimension(n_std_lev_atm) :: std_lev_atm = (/2., 10./)

!  Businger functions

    psium(zz1,zz2) = log( (1.+zz1)**2*(1.+zz1**2)/((1.+zz2)**2*(1.+zz2**2)) ) -2.*(atan(zz1)-atan(zz2))
    psiuh(zz1,zz2) = 2.*log((1.+zz1**2)/(1.+zz2**2))

!  Holtslag functions

    psism(zz1,zz2)=-zaholt*zz1-zbholt*(zz1-zcholt/zdholt)*exp(-zdholt*zz1)  &
                   +zaholt*zz2+zbholt*(zz2-zcholt/zdholt)*exp(-zdholt*zz2)
    psish(zz1,zz2)=-(1.+2./3.*zaholt*zz1)**1.5-zbholt*(zz1-zcholt/zdholt)*exp(-zdholt*zz1) &
                   +(1.+2./3.*zaholt*zz2)**1.5+zbholt*(zz2-zcholt/zdholt)*exp(-zdholt*zz2)

!--------------------------------------------
!  loop in latitude to the end of the routine
!--------------------------------------------
 OPEN (51,FILE="surf_rich_num_01.dat",STATUS="UNKNOWN")

    jstart = 1
    jend = nlat

#ifdef globo
    if (ip_s.eq.ip_null) jstart = 1
    if (ip_n.eq.ip_null) jend = nlat
#endif

    do jlat = jstart, jend
    jlatp1 = min (jlat+1, nlat) ! necessaria per globo

!  conversion factor from t to theta
!  definition of virtual theta
!  specific humidity at saturation
!  half-level height

!------------------------------------------------------------------------------------------
!  turbulent fluxes at the ground (positive upward) - surface layer - Monin Obukhov theory
!------------------------------------------------------------------------------------------

    do jlon = 1, nlon

    zzpp = p_atm(jlon,jlat)
    zconv = exp(rdrcp*(alp0-log(zzpp)))
    ztetav_bottom = tvirt_atm(jlon,jlat)*zconv

    zza   = h_atm(jlon,jlat) + rgm(jlon,jlat) !  sigma=1 is located at z=rgm
    zua   = u_atm(jlon,jlat)
    zva   = v_atm(jlon,jlat)
    zmod2 = zua**2 + zva**2 + .07
    zmod  = sqrt(zmod2)

!  virtual potential temperature computed with skin temperature

    zconvg = exp(rdrcp*(alp0-log(ps(jlon,jlat))))*(1.+ep*qskin(jlon,jlat))
    ztevg = tskin(jlon,jlat)*zconvg

!  bulk Richardson number

    ztebar =.5*(ztevg+ztetav_bottom)
    beta   = g/ztebar
    zri    = zza*beta*(ztetav_bottom-ztevg)/zmod2
   WRITE (51,'(100E16.8)') float(kstep)*60./3600., zri

    if(fmask(jlon,jlat).ge..5.and.fice(jlon,jlat).lt.0.1) then

!  computation of Charnock roughness

      zchar = 5.e-4
      zcoch1 = .0185*zak**2/g*zmod2
      do jiter = 1, 5
      zcoch2 = zza/zchar
      zchar = zcoch1/log(zcoch2)**2
      enddo
      zrgm = zchar

! roughness lengths over sea interpolated logarithmically between 0.<Ri<0.25 with Large & Pond values

       if(zri.ge.0.25) then
       zrgt = 2.2e-9
       zrgq = zrgt
       elseif(zri.lt.0.) then
       zrgt = 5.e-5
       zrgq = 9.e-5
       else
       zrgt = 3.80227e-6*exp(-zri*29.819387)
       zrgq = zrgt
       endif

    else                     ! soil or sea ice > 0.1
      zrgm = rgm(jlon,jlat)
      zrgt = rgq(jlon,jlat)
      zrgm = zrgm*(1.-0.5*fsnow(jlon,jlat)*0.4/(0.4+zrgm))  ! reduction over snow
      zrgt = zrgt*(1.-0.5*fsnow(jlon,jlat)*0.4/(0.4+zrgt))  ! reduction over snow
      zrgq = zrgt
    endif

    rgmd(jlon,jlat) = zrgm

#ifdef chem
      bolchem_rgm(jlon,jlat)=zrgm
      bolchem_rgq(jlon,jlat)=zrgq
#endif

    zalzam = log(zza/zrgm)
    zalzat = log(zza/zrgt)
    zalzaq = log(zza/zrgq)

    if(zri.gt.0.) then         ! Holtslag functions

      zal = .2 + 4.*zri
      zpsim = psism(zal, zal*zrgm/zza)
      zpsih = psish(zal, zal*zrgt/zza)
      do jiter = 1, 10
      zalm = zal
      zal = zri*(zalzam-zpsim)**2/(zalzat-zpsih)   ! za/L from eq (6.47) of Businger
      error = abs(zal-zalm)/zal
      if (error.lt.1.e-2) go to 1
      zpsim = psism(zal, zal*zrgm/zza)
      zpsih = psish(zal, zal*zrgt/zza)
      enddo
 1    continue
      zpsiq = zpsih  ! because rgt=rgq in the stable case
      ustar(jlon,jlat) = zak*zmod/(zalzam-zpsim)
      tstar(jlon,jlat) = zak/(zalzat-zpsih)*(ztetav_bottom-ztevg)
      zqstar           = zak/(zalzaq-zpsiq)*(q_atm(jlon,jlat)-qskin(jlon,jlat))
      zcdm = ustar(jlon,jlat)**2/zmod
      zcdt = ustar(jlon,jlat)*zak/(zalzat-zpsih)
      zcdq = zcdt
           zpsim10 = psism(zal*(10.+zrgm)/zza, zal*zrgm/zza)  ! to compute 10 m wind
           zpsih2  = psish(zal*(2. +zrgt)/zza, zal*zrgt/zza)  ! to compute 2 m temp.
           zpsiq2  = zpsih2    ! because rgt=rgq in the stable case

    else                       ! Businger functions

      zpsim = 0.
      zpsih = 0.
      do jiter = 1, 4
      zal = zri*(zalzam-zpsim)**2/(zalzat-zpsih)   ! za/L from eq (6.47) of Businger
      zz1 = (1.-zgam*zal)**.25
      zx2 = (1.-zgam*zal*zrgm/zza)**.25
      zy2 = (1.-zgam*zal*zrgt/zza)**.25
      zpsim = psium(zz1,zx2)
      zpsih = psiuh(zz1,zy2)
      enddo
      zz2 = (1.-zgam*zal*zrgq/zza)**.25
      zpsiq = psiuh(zz1,zz2)
      ustar(jlon,jlat) = zak*zmod/(zalzam-zpsim)
      tstar(jlon,jlat) = zak/(zalzat-zpsih)*(ztetav_bottom-ztevg)
      zqstar           = zak/(zalzaq-zpsiq)*(q_atm(jlon,jlat)-qskin(jlon,jlat))
      zcdm = ustar(jlon,jlat)**2/zmod
      zcdt = ustar(jlon,jlat)*zak/(zalzat-zpsih)
      zcdq = ustar(jlon,jlat)*zak/(zalzaq-zpsiq)
           zz1 = (1.-zgam*zal*(10.+zrgm)/zza)**.25   ! for 10 m wind
           zpsim10 = psium(zz1, zx2)                 ! for 10 m wind
           zz1 = (1.-zgam*zal*(2. +zrgt)/zza)**.25   ! for 2 m t
           zpsih2  = psiuh(zz1, zy2)                 ! for 2 m t
           zz1 = (1.-zgam*zal*(2. +zrgq)/zza)**.25   ! for 2 m q
           zpsiq2  = psiuh(zz1, zz2)                 ! for 2 m q

    endif

!--------------------------------------------------------------------------------------------------
! Output variables
!--------------------------------------------------------------------------------------------------

    zlev_bottom = h_atm(jlon,jlat)

    kturb_surf_m(jlon,jlat) = abs( -zcdm*zlev_bottom )
    kturb_surf_h(jlon,jlat) = abs( (-zcdt/zconvg)*zlev_bottom )
    kturb_surf_q(jlon,jlat) = abs( -zcdq*zlev_bottom)

    if (kturb_surf_h(jlon,jlat) > 1.) then
!      zzz=max(-0.08*kturb_surf_h(jlon,jlat)+1.08, 0.03) 
      zzz=max(-0.10*kturb_surf_h(jlon,jlat)+1.10, 0.02) 
      kturb_surf_h(jlon,jlat) = kturb_surf_h(jlon,jlat)*zzz ! poisk
    else
      kturb_surf_h(jlon,jlat) = kturb_surf_h(jlon,jlat)*2.
    endif
!    kturb_surf_q(jlon,jlat) = min(kturb_surf_q(jlon,jlat)*0.8, 7.) ! poisk
    kturb_surf_q(jlon,jlat) = min(kturb_surf_q(jlon,jlat)*0.5, 7.) ! poisk

!--------------------------------------------------------------------------------------------------

!  coefficients for comput. of surface fluxes of heat and q

    zros = ps(jlon,jlat)/(rd*tskin(jlon,jlat)*(1.+ep*qskin(jlon,jlat)))
    roscd(jlon,jlat)=  .25*zros*zcdq + .75*roscd(jlon,jlat)

!--------------------------------------------------------------------------------------------------

#ifdef chem
  invra(jlon,jlat) = zcdt
  rhos(jlon,jlat) = zros
#endif

!--------------------------------------------------------------------------------------------------

 enddo ! end of loop in longitude
 enddo ! end of loop in latitude

!--------------------------------------------------------------------------------------------------

! Time filter for turbulent coefficient (heat and water vapour) in the surface layer

 if (nstep_sl_filter > 1) then

   if (kstep > nstep_sl_filter) then
     do n = 1,nstep_sl_filter-1
       kturb_surf_h_mem(:,:,n) = kturb_surf_h_mem(:,:,n+1) 
       kturb_surf_q_mem(:,:,n) = kturb_surf_q_mem(:,:,n+1) 
     enddo
     kturb_surf_h_mem(:,:,nstep_sl_filter) = kturb_surf_h(:,:)
     kturb_surf_q_mem(:,:,nstep_sl_filter) = kturb_surf_q(:,:)
     kturb_surf_h_filter(:,:) = 0.
     kturb_surf_q_filter(:,:) = 0.
     do n = 1,nstep_sl_filter
       kturb_surf_h_filter(:,:) = kturb_surf_h_filter(:,:) + kturb_surf_h_mem(:,:,n)
       kturb_surf_q_filter(:,:) = kturb_surf_q_filter(:,:) + kturb_surf_q_mem(:,:,n)
     enddo
     kturb_surf_h_filter(:,:) = kturb_surf_h_filter(:,:) / float(nstep_sl_filter)
     kturb_surf_q_filter(:,:) = kturb_surf_q_filter(:,:) / float(nstep_sl_filter)
     kturb_surf_h(:,:) = kturb_surf_h_filter(:,:)
     kturb_surf_q(:,:) = kturb_surf_q_filter(:,:)
  !   if (kstep == nstep_sl_filter+1.or.mod( (kstep*int(dtstep)), 3600) == 0) then
  !     do jlat = 2, nlat-1
  !     do jlon = 2, nlon-1
  !       zmin = minval( kturb_surf_h_mem(jlon,jlat,1:nstep_sl_filter) )
  !       zmax = maxval( kturb_surf_h_mem(jlon,jlat,1:nstep_sl_filter) )
  !       if (zmax > 1..and.zmin < 1.e-2) then
  !         kturb_surf_h(jlon,jlat) = kturb_surf_h_filter(jlon,jlat)
  !         kturb_surf_q(jlon,jlat) = kturb_surf_q_filter(jlon,jlat)
  !       endif
  !     enddo 
  !     enddo 
  !   endif
   else
     kturb_surf_h_mem(:,:,kstep) = kturb_surf_h(:,:)
     kturb_surf_q_mem(:,:,kstep) = kturb_surf_q(:,:)
   endif

 endif

!  fluxes at new time step (qflux is redefined in subr. soil but only over land)
!  (note that to effectively suppress h and/or q fluxes to the atmosphere, it is necessary
!  to set zct(:) and/or zcq(:) to zero at the beginning of vdiff - but hflux and qflux
!  still affect heat and moisture fluxes into soil and sea)

 do jlat = 1, nlat
 do jlon = 1, nlon
   zzpp = p_atm(jlon,jlat)
   zconv = exp(rdrcp*(alp0-log(zzpp)))
   ztetav_bottom = tvirt_atm(jlon,jlat)*zconv
   zconvg = exp(rdrcp*(alp0-log(ps(jlon,jlat))))*(1.+ep*qskin(jlon,jlat))
   ztevg = tskin(jlon,jlat)*zconvg
   zlev_bottom = h_atm(jlon,jlat)
   zros = ps(jlon,jlat)/(rd*tskin(jlon,jlat)*(1.+ep*qskin(jlon,jlat)))
   hflux(jlon,jlat) = kturb_surf_h(jlon,jlat)*(zros*cpd)/zlev_bottom*(ztevg-ztetav_bottom)
   qflux(jlon,jlat) = kturb_surf_q(jlon,jlat)*zros/zlev_bottom*(qskin(jlon,jlat)-q_atm(jlon,jlat))
 enddo
 enddo

    return
    end subroutine surface_layer
