MODULE screenp_mod

! This module contains some procedures for calculation of the first
! guess of temperature, specific humidity and wind at a reference level
! coefficients for turbulent diffusion at surface

  IMPLICIT NONE

CONTAINS

!****************************************************************************************

!r original routine

      SUBROUTINE screenp(klon, knon, nsrf, &
                     speed, tair, qair, &
                     ts, qsurf, rugos, lmon, &
                     ustar, testar, qstar, zref, &
                     delu, delte, delq)
      IMPLICIT NONE
!-------------------------------------------------------------------------

! Objet : calcul "predicteur" des anomalies du vent, de la temperature 
!         potentielle et de l'humidite relative au niveau de reference zref et 
!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
!         a partir des relations de Dyer-Businger.

! Reference : Hess, Colman et McAvaney (1995)

! I. Musat, 01.07.2002
!-------------------------------------------------------------------------

! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
! knon----input-I- nombre de points pour un type de surface
! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
! speed---input-R- module du vent au 1er niveau du modele
! tair----input-R- temperature de l'air au 1er niveau du modele
! qair----input-R- humidite relative au 1er niveau du modele
! ts------input-R- temperature de l'air a la surface
! qsurf---input-R- humidite relative a la surface
! rugos---input-R- rugosite
! lmon----input-R- longueur de Monin-Obukov
! ustar---input-R- facteur d'echelle pour le vent
! testar--input-R- facteur d'echelle pour la temperature potentielle
! qstar---input-R- facteur d'echelle pour l'humidite relative
! zref----input-R- altitude de reference

! delu----input-R- anomalie du vent par rapport au 1er niveau
! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
! delq----input-R- anomalie de l'humidite relative par rapport a la surface

      INTEGER, INTENT(IN) :: klon, knon, nsrf
      REAL, DIMENSION(klon), INTENT(IN) :: speed, tair, qair
      REAL, DIMENSION(klon), INTENT(IN) :: ts, qsurf, rugos
      DOUBLE PRECISION, DIMENSION(klon), INTENT(IN) :: lmon
      REAL, DIMENSION(klon), INTENT(IN) :: ustar, testar, qstar
      REAL, INTENT(IN) :: zref

      REAL, DIMENSION(klon), INTENT(OUT) :: delu, delte, delq

!-------------------------------------------------------------------------
! Variables locales et constantes :
      REAL, PARAMETER :: RKAR=0.40
      INTEGER :: i
      REAL :: xtmp, xtmp0
!-------------------------------------------------------------------------
      DO i = 1, knon

        IF (lmon(i)>=0.) THEN

! STABLE CASE

          IF (speed(i)>1.5.AND.lmon(i)<=1.0                        &
                        .AND. rugos(i)<=1.0) THEN
            delu(i) = (ustar(i)/RKAR)* &
                      (log(zref/(rugos(i))+1.) + &
                      min(5.d0, 5.0 *(zref - rugos(i))/lmon(i)))
            delte(i) = (testar(i)/RKAR)* &
                       (log(zref/(rugos(i))+1.) + &
                       min(5.d0, 5.0 * (zref - rugos(i))/lmon(i)))
            delq(i) = (qstar(i)/RKAR)* &
                      (log(zref/(rugos(i))+1.) + &
                      min(5.d0, 5.0 * (zref - rugos(i))/lmon(i)))
          ELSE
            delu(i)  = 0.1 * speed(i)
            delte(i) = 0.1 * (tair(i) - ts(i) )
            delq(i)  = 0.1 * (max(qair(i),0.0) - max(qsurf(i),0.0))
          ENDIF
        ELSE  

! UNSTABLE CASE

          IF (speed(i)>5.0.AND.abs(lmon(i))<=50.0) THEN
            xtmp = (1. - 16. * (zref/lmon(i)))**(1./4.)
            xtmp0 = (1. - 16. * (rugos(i)/lmon(i)))**(1./4.)
            delu(i) = (ustar(i)/RKAR)* &
                      (log(zref/(rugos(i))+1.) & 
                      - 2.*log(0.5*(1. + xtmp)) &
                      + 2.*log(0.5*(1. + xtmp0)) &
                      - log(0.5*(1. + xtmp*xtmp)) &
                      + log(0.5*(1. + xtmp0*xtmp0)) &
                      + 2.*atan(xtmp) - 2.*atan(xtmp0))
            delte(i) = (testar(i)/RKAR)* &
                       (log(zref/(rugos(i))+1.) &
                       - 2.0 * log(0.5*(1. + xtmp*xtmp)) & 
                       + 2.0 * log(0.5*(1. + xtmp0*xtmp0)))
            delq(i)  = (qstar(i)/RKAR)* &
                       (log(zref/(rugos(i))+1.) &
                       - 2.0 * log(0.5*(1. + xtmp*xtmp)) & 
                       + 2.0 * log(0.5*(1. + xtmp0*xtmp0)))
          ELSE
            delu(i)  = 0.5 * speed(i)
            delte(i) = 0.5 * (tair(i) - ts(i) )
            delq(i)  = 0.5 * (max(qair(i),0.0) - max(qsurf(i),0.0))
          ENDIF
        ENDIF

      ENDDO

      END SUBROUTINE screenp

      SUBROUTINE screenpn(klon, knon, nsrf, &
                     speed, tair, qair, &
                     ts, qsurf, rugos, zri1, &
                     zref, &
                     delu, delte, delq)
      IMPLICIT NONE
!-------------------------------------------------------------------------

! Objet : calcul "predicteur" des anomalies du vent, de la temperature 
!         potentielle et de l'humidite relative au niveau de reference zref et 
!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
!         a partir des relations de Dyer-Businger.

! Reference : Hess, Colman et McAvaney (1995)

! I. Musat, 01.07.2002
!-------------------------------------------------------------------------

! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
! knon----input-I- nombre de points pour un type de surface
! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
! speed---input-R- module du vent au 1er niveau du modele
! tair----input-R- temperature de l'air au 1er niveau du modele
! qair----input-R- humidite relative au 1er niveau du modele
! ts------input-R- temperature de l'air a la surface
! qsurf---input-R- humidite relative a la surface
! rugos---input-R- rugosite
! zref----input-R- altitude de reference

! delu----input-R- anomalie du vent par rapport au 1er niveau
! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
! delq----input-R- anomalie de l'humidite relative par rapport a la surface

      INTEGER, INTENT(IN) :: klon, knon, nsrf
      REAL, DIMENSION(klon), INTENT(IN) :: speed, tair, qair
      REAL, DIMENSION(klon), INTENT(IN) :: ts, qsurf, rugos
      REAL, DIMENSION(klon), INTENT(IN) :: zri1
      REAL, INTENT(IN) :: zref

      REAL, DIMENSION(klon), INTENT(OUT) :: delu, delte, delq

!-------------------------------------------------------------------------
! Variables locales et constantes :
      REAL, PARAMETER :: RKAR=0.40
      INTEGER :: i
      REAL :: xtmp, xtmp0
!-------------------------------------------------------------------------
      DO i = 1, knon

       IF (zri1(i)>=0.) THEN
          delu(i)  = 0.1 * speed(i)
          delte(i) = 0.1 * (tair(i) - ts(i) )
          delq(i)  = 0.1 * (max(qair(i),0.0) - max(qsurf(i),0.0))
       ELSE
          delu(i)  = 0.5 * speed(i)
          delte(i) = 0.5 * (tair(i) - ts(i) )
          delq(i)  = 0.5 * (max(qair(i),0.0) - max(qsurf(i),0.0))
       ENDIF

      ENDDO

      END SUBROUTINE screenpn
END MODULE screenp_mod