flxsurf3bx.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!-------------------------------------- LICENCE BEGIN ------------------------------------
!Environment Canada - Atmospheric Science and Technology License/Disclaimer,
!                     version 3; Last Modified: May 7, 2008.
!This is free but copyrighted software; you can use/redistribute/modify it under the terms
!of the Environment Canada - Atmospheric Science and Technology License/Disclaimer
!version 3 or (at your option) any later version that should be found at:
!http://collaboration.cmc.ec.gc.ca/science/rpn.comm/license.html
!
!This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
!without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!See the above mentioned License/Disclaimer for more details.
!You should have received a copy of the License/Disclaimer along with this software;
!if not, you can write to: EC-RPN COMM Group, 2121 TransCanada, suite 500, Dorval (Quebec),
!CANADA, H9P 1J3; or send e-mail to service.rpn@ec.gc.ca
!-------------------------------------- LICENCE END --------------------------------------
!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
!!!S/P  FLXSURF3
!
      SUBROUTINE flxsurf3bx(CMU, CTU, RIB, FTEMP, FVAP, ILMO,           &
     &                   ue, fcor, ta , qa , zu, zt, va,                &
     &                   tg , qg , h , z0 , z0t,                        &
     &                   lzz0, lzz0t, fm, fh, n )
!
      USE yomhook   ,ONLY : lhook,   dr_hook
      USE parkind1  ,ONLY : jprb
!RJ: added modi, after freeform conversion
      USE modi_vslog
!
      IMPLICIT NONE
      INTEGER :: n
      REAL :: cmu(n),ctu(n),rib(n),fcor(n),ilmo(n)
      REAL :: ftemp(n),fvap(n),ta(n),qa(n),zu(n),zt(n),va(n)
      REAL :: tg(n),qg(n),h(n),z0(n),ue(n)
      REAL :: z0t(n),lzz0(n),lzz0t(n)
      REAL :: fm(n),fh(n)
!
!Author
!          Y.Delage (Jul 1990)
!Revision
! 001      G. Pellerin (Jun 94) New function for unstable case
! 002      G. Pellerin (Jui 94) New formulation for stable case
! 003      B. Bilodeau (Nov 95) Replace VK by KARMAN
! 004      M. Desgagne (Dec 95) Add safety code in function ff
!                               and ensures that RIB is non zero
! 005      R. Sarrazin (Jan 96) Correction for H
! 006      C. Girard (Nov 95) - Diffuse T instead of Tv
! 007      G. Pellerin (Feb 96) Revised calculation for H (stable)
! 008      G. Pellerin (Feb 96) Remove corrective terms to CTU
! 009      Y. Delage and B. Bilodeau (Jul 97) - Cleanup
! 010      Y. Delage (Feb 98) - Addition of HMIN
! 011      D. Talbot and Y. Delage (Jan 02) -
!             Correct bug of zero divide by dg in loop 35
! 012      Y. Delage (Oct 03) - Set top of surface layer at ZU +Z0
!                   - Output UE instead of UE**2 and rename subroutine
!                   - Change iteration scheme for stable case
!                   - Introduce log-linear profile for near-neutral stable cases
!                   - set VAMIN inside flxsurf and initialise ILMO and H
!                   - Put stability functions into local functions via stabfunc.h
! 013      Y. Delage (Sep 04) - Input of wind and temperature/humidity
!                                at different levels
! 014      R. McTaggart-Cowan and B. Bilodeau (May 2006) -
!             Clean up stabfunc.h
! 015      L. Spacek (Dec 07) - Correction of the log-linear profile
!                               Double precision for rib calculations
!
!Object
!          to calculate surface layer transfer coefficients and fluxes
!
!Arguments
!
!          - Output -
! CMU      transfer coefficient of momentum times UE
! CTU      transfer coefficient of temperature times UE
! RIB      bulk Richardson number
! FTEMP    temperature flux
! FVAP     vapor flux
! ILMO     (1/length of Monin-Obukov)
! UE       friction velocity
! H        height of the boundary layer
! FM       momentum stability function
! FH       heat stability function
! LZZ0     log ((zu+z0)/z0)
! LZZ0T    log ((zt+z0)/z0t)
!
!          - Input -
! FCOR     Coriolis factor
! ZU       height of wind input (measured from model base at topo height + Z0)
! ZT       height of temperature and humidity input
! TA       potential temperature at ZT
! QA       specific humidity at ZT
! VA       wind speed at ZU
! TG       surface temperature
! QG       specific humidity at the surface
! Z0       roughness length for momentum      flux calculations
! Z0T      roughness length for heat/moisture flux calculations
! N        horizontal dimension
!
!
!RJ #include "surfcon.h"
!RJ       LOGICAL :: INIT
!     PHYSICAL CONSTANTS
       REAL,PARAMETER :: cpd      =.100546e+4        ! J K-1 kg-1    ! specific heat of dry air
       REAL,PARAMETER :: cpv      =.186946e+4        ! J K-1 kg-1    ! specific heat of water vapour
       REAL,PARAMETER :: rgasd    =.28705e+3         ! J K-1 kg-1    ! gas constant for dry air
       REAL,PARAMETER :: rgasv    =.46151e+3         ! J K-1 kg-1    ! gas constant for water vapour
       REAL,PARAMETER :: trpl     =.27316e+3         ! K             ! triple point of water
       REAL,PARAMETER :: tcdk     =.27315e+3         !               ! conversion from kelvin to celsius
       REAL,PARAMETER :: rauw     =.1e+4             !               ! density of liquid H2O
       REAL,PARAMETER :: eps1     =.62194800221014   !               ! RGASD/RGASV
       REAL,PARAMETER :: eps2     =.3780199778986    !               ! 1 - EPS1
       REAL,PARAMETER :: delta    =.6077686814144    !               ! 1/EPS1 - 1
       REAL,PARAMETER :: cappa    =.28549121795      !               ! RGASD/CPD
       REAL,PARAMETER :: tgl      =.27316e+3         ! K             ! ice temperature in the atmosphere
       REAL,PARAMETER :: consol   =.1367e+4          ! W m-2         ! solar constant
       REAL,PARAMETER :: grav     =.980616e+1        ! M s-2         ! gravitational acceleration
       REAL,PARAMETER :: rayt     =.637122e+7        ! M             ! mean radius of the earth
       REAL,PARAMETER :: stefan   =.566948e-7        ! J m-2 s-1 K-4 ! Stefan-Boltzmann constant
       REAL,PARAMETER :: pi       =.314159265359e+1  !               ! PI constant = ACOS(-1)
       REAL,PARAMETER :: omega    =.7292e-4          ! s-1           ! angular speed of rotation of the earth
       REAL,PARAMETER :: knams    =.514791           !               ! conversion from knots to m/s
       REAL,PARAMETER :: stlo     =.6628486583943e-3 ! K s2 m-2      ! Schuman-Newell Lapse Rate
       REAL,PARAMETER :: karman   =.35               !               ! Von Karman constant
       REAL,PARAMETER :: ric      =.2                !               ! Critical Richardson number
       REAL,PARAMETER :: chlc     =.2501e+7          ! J kg-1        ! latent heat of condensation
       REAL,PARAMETER :: chlf     =.334e+6           ! J kg-1        ! latent heat of fusion
       REAL,PARAMETER :: t1s      =.27316e+3         ! K             ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: t2s      =.25816e+3         ! K             ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: aw       =.3135012829948e+4 !               ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: bw       =.2367075766316e+1 !               ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: ai       =.2864887713087e+4 !               ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: bi       =.166093131502     !               ! constant used to calculate L/Cp in fcn HTVOCP
       REAL,PARAMETER :: slp      =.6666666666667e-1 !               ! constant used to calculate L/Cp in fcn HTVOCP

!RJ #include "consphy.h"
!     INITIALIZES THE CONSTANTS FOR THE COMMONS OF THE FLXSURF3 ROUTINE FROM
!     CANADIAN METEOROLOGICAL CENTER
      REAL,PARAMETER :: as    = 12.
      REAL,PARAMETER :: asx   = 5.
      REAL,PARAMETER :: ci    = 40.
      REAL,PARAMETER :: bs    = 1.0
      REAL,PARAMETER :: beta  = 1.0
      REAL,PARAMETER :: factn = 1.2
      REAL,PARAMETER :: hmin  = 30.
      REAL,PARAMETER :: angmax= 0.85
      REAL,PARAMETER :: rac3  = sqrt(3.)
!
!*
!
      INTEGER,PARAMETER :: jdbl=8
!
      INTEGER :: j
      INTEGER :: it
      INTEGER,PARAMETER :: itmax = 3
      REAL,PARAMETER :: hmax = 1500.0
      REAL,PARAMETER :: cormin = 0.7e-4
      REAL,PARAMETER :: epsln = 1.0e-05
      REAL,PARAMETER :: vamin = 0.1
      REAL :: cm,ct,zp
      REAL :: f,g,dg
      REAL :: hi,he,hs,unsl
      REAL(KIND=JDBL) :: dthv,tva,tvs
      REAL :: hl,u
      REAL :: cs,xx,xx0,yy,yy0
      REAL :: zb,dd,ilmox
      REAL :: df,zz,betsasx
      REAL :: aa,bb,cc
      REAL(KIND=JPRB) :: zhook_handle
!
      df(zz)=(1-zz*hi)*sqrt(1+4*as*beta*unsl*zz/(1-zz*hi))
      cs=as*2.5
      betsasx=1./asx
!
      IF (lhook) CALL dr_hook('FLXSURF3BX',0,zhook_handle)
!
      DO j=1,n
        lzz0(j)=1+zu(j)/z0(j)
        lzz0t(j)=(zt(j)+z0(j))/z0t(j)
      ENDDO
!
      call vslog(lzz0t,lzz0t,n)
      call vslog(lzz0 ,lzz0 ,n)
!
      DO j=1,n
!
!  CALCULATE THE RICHARDSON NUMBER
        zp=zu(j)**2/(zt(j)+z0(j)-z0t(j))
        u=max(vamin,va(j))
        tva=(1.0_jdbl+delta*qa(j))*ta(j)
        tvs=(1.0_jdbl+delta*qg(j))*tg(j)
        dthv=tva-tvs
        rib(j)=grav/(tvs+0.5_jdbl*dthv)*zp*dthv/(u*u)
        if (rib(j)>=0.0_jdbl) rib(j) = max(rib(j), epsln)
        if (rib(j)<0.0_jdbl) rib(j) = min(rib(j),-epsln)
!
!  FIRST APPROXIMATION TO ILMO
        IF(rib(j)>0.0_jdbl)  THEN
           fm(j)=lzz0(j)+cs*rib(j)/max(2*z0(j),1.0_jdbl)
           fh(j)=beta*(lzz0t(j)+cs*rib(j))/                             &
     &          max(sqrt(z0(j)*z0t(j)),1.0_jdbl)
           ilmo(j)=rib(j)*fm(j)*fm(j)/(zp*fh(j))
           f=max(abs(fcor(j)),cormin)
           h(j)=bs*sqrt(karman*u/(ilmo(j)*f*fm(j)))
        ELSE
           fm(j)=lzz0(j)-min(0.7_jdbl+log(1-rib(j)),lzz0(j)-1)
           fh(j)=beta*(lzz0t(j)-min(0.7_jdbl+log(1-rib(j)),lzz0t(j)-1))
           ilmo(j)=rib(j)*fm(j)*fm(j)/(zp*fh(j))
        ENDIF
      ENDDO
!
! - - - - - - - - -  BEGINNING OF ITERATION LOOP - - - - - - - - - - -
      DO 35 it=1,itmax
      DO 35 j=1,n
        u=max(vamin,va(j))
        zp=zu(j)**2/(zt(j)+z0(j)-z0t(j))
        IF(rib(j)>0.0_jdbl)  THEN
!----------------------------------------------------------------------
!  STABLE CASE
        ilmo(j)=max(epsln,ilmo(j))
        hl=(zu(j)+10*z0(j))*factn
        f=max(abs(fcor(j)),cormin)
        hs=bs*sqrt(karman*u/(ilmo(j)*f*fm(j)))
        h(j)=max(hmin,hs,hl,factn/(4*as*beta*ilmo(j)))
        hi=1/h(j)
!CDIR IEXPAND
        fm(j)=lzz0(j)+psi(zu(j)+z0(j),hi,ilmo(j))-psi(z0(j),hi,ilmo(j))
!CDIR IEXPAND
        fh(j)=beta*(lzz0t(j)+psi(zt(j)+z0(j),hi,ilmo(j))-psi(z0t(j),hi, &
     &       ilmo(j)))
        unsl=ilmo(j)
        dg=-zp*fh(j)/(fm(j)*fm(j))*(1+beta*(df(zt(j)+z0(j))-df(z0t(j)))/&
     &     (2*fh(j))-(df(zu(j)+z0(j))-df(z0(j)))/fm(j))
!----------------------------------------------------------------------
!  UNSTABLE CASE
      ELSE
        ilmo(j)=min(0.,ilmo(j))
!CDIR IEXPAND
        fm(j)=fmi(zu(j)+z0(j),z0(j),lzz0(j),ilmo(j),xx,xx0)
!CDIR IEXPAND
        fh(j)=fhi(zt(j)+z0(j),z0t(j),lzz0t(j),ilmo(j),yy,yy0)
         dg=-zp*fh(j)/(fm(j)*fm(j))*(1+beta/fh(j)*(1/yy-1/yy0)-2/fm(j)* &
     &               (1/xx-1/xx0))
      ENDIF
!----------------------------------------------------------------------
      IF(it<itmax) THEN
             g=rib(j)-fh(j)/(fm(j)*fm(j))*zp*ilmo(j)
             ilmo(j)=ilmo(j)-g/dg
      ENDIF
   35 CONTINUE
! - - - - - -  - - - END OF ITERATION LOOP - - - - - - - - - - - - - -
!
      DO 80 j=1,n
      u=max(vamin,va(j))
      if(asx<as) then
!----------------------------------------------------------------------
!  CALCULATE ILMO AND STABILITY FUNCTIONS FROM LOG-LINEAR PROFILE
!     (SOLUTION OF A QUADRATIC EQATION)
!
        zb=zu(j)/(zt(j)+z0(j)-z0t(j))
!  DISCRIMINANT
        dd=(beta*lzz0t(j)*zb)**2-4*rib(j)*asx*lzz0(j)*                  &
     &       (beta*lzz0t(j)*zb-lzz0(j))
        if(rib(j)>0.0_jdbl.and.rib(j)<betsasx.and.dd>=0.) then
!  COEFFICIENTS
           aa=asx*asx*rib(j)-asx
           bb=-beta*lzz0t(j)*zb+2*rib(j)*asx*lzz0(j)
           cc=rib(j)*lzz0(j)**2
!  SOLUTION
           if(bb>=0)then
              ilmox=(-bb-sqrt(dd))/(2*zu(j)*aa)
           else
              ilmox=2*cc/(zu(j)*(-bb+sqrt(dd)))
           endif
           if(ilmox<ilmo(j)) then
              ilmo(j)=ilmox
              fm(j)=lzz0(j)+asx*zu(j)*ilmox
              fh(j)=beta*lzz0t(j)+asx*(zt(j)+z0(j)-z0t(j))*ilmox
           endif
        endif
      endif
!----------------------------------------------------------------------
        cm=karman/fm(j)
        ct=karman/fh(j)
        ue(j)=u*cm
        cmu(j)=cm*ue(j)
        ctu(j)=ct*ue(j)
        if (rib(j)>0.0_jdbl) then
!          stable case
              h(j)=min(h(j),hmax)
        else
!          unstable case
              f=max(abs(fcor(j)),cormin)
              he=max(hmin,0.3_jdbl*ue(j)/f)
              h(j)=min(he,hmax)
        endif
        ftemp(j)=-ctu(j)*(ta(j)-tg(j))
        fvap(j)=-ctu(j)*(qa(j)-qg(j))
   80 CONTINUE
      IF (lhook) CALL dr_hook('FLXSURF3BX',1,zhook_handle)
      CONTAINS
!RJ: inlining directly
!RJ #include "stabfunc2.h"
!
!   Internal function FMI
!   Stability function for momentum in the unstable regime (ilmo<0)
!   Reference: Delage Y. and Girard C. BLM 58 (19-31) Eq. 19
!
      FUNCTION fmi(Z2,Z02,LZZ02,ILMO2,X,X0)
      IMPLICIT NONE
!
      REAL              :: fmi
      REAL, INTENT(IN ) :: z2,z02,lzz02,ilmo2
      REAL, INTENT(OUT) :: x,x0
!
      x =(1-ci*z2 *beta*ilmo2)**(0.16666666)
      x0=(1-ci*z02*beta*ilmo2)**(0.16666666)
      fmi=lzz02+log((x0+1)**2*sqrt(x0**2-x0+1)*(x0**2+x0+1)**1.5        &
     &               /((x+1)**2*sqrt(x**2-x+1)*(x**2+x+1)**1.5))        &
     &              +rac3*atan(rac3*((x**2-1)*x0-(x0**2-1)*x)/          &
     &              ((x0**2-1)*(x**2-1)+3*x*x0))
!
      RETURN
      END FUNCTION fmi
!
!   Internal function FHI
!   Stability function for heat and moisture in the unstable regime (ilmo<0)
!   Reference: Delage Y. and Girard C. BLM 58 (19-31) Eq. 17
!
      FUNCTION fhi(Z2,Z0T2,LZZ0T2,ILMO2,Y,Y0)
      IMPLICIT NONE
!
      REAL              :: fhi
      REAL, INTENT(IN ) :: z2,z0t2,lzz0t2,ilmo2
      REAL, INTENT(OUT) :: y,y0
!
      y =(1-ci*z2  *beta*ilmo2)**(0.33333333)
      y0=(1-ci*z0t2*beta*ilmo2)**(0.33333333)
      fhi=beta*(lzz0t2+1.5*log((y0**2+y0+1)/(y**2+y+1))+rac3*           &
     &        atan(rac3*2*(y-y0)/((2*y0+1)*(2*y+1)+3)))
!
      RETURN
      END FUNCTION fhi
!
!   Internal function psi
!   Stability function for momentum in the stable regime (unsl>0)
!   Reference :  Y. Delage, BLM, 82 (p23-48) (Eqs.33-37)
!
      FUNCTION psi(Z2,HI2,ILMO2)
      IMPLICIT NONE
!
      REAL :: psi
      REAL :: a,b,c,d
      REAL, INTENT(IN ) :: ilmo2,z2,hi2
!
      d = 4*as*beta*ilmo2
      c = d*hi2 - hi2**2
      b = d - 2*hi2
      a = sqrt(1 + b*z2 - c*z2**2)
      psi = 0.5 * (a-z2*hi2-log(1+b*z2*0.5+a)-                          &
     &            b/(2*sqrt(c))*asin((b-2*c*z2)/d))
!
      RETURN
      END FUNCTION psi
!
      END SUBROUTINE flxsurf3bx