source: trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F @ 268

      SUBROUTINE surflayer_interpol(ngrid,nlay,pz0, 
     & pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar)
      IMPLICIT NONE
!=======================================================================
!
!   Subject: interpolation of temperature and velocity norm in the surface layer
!   by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u*, teta*)
!   -------  
!
!   Author: Arnaud Colaitis 05/08/11
!   -------
!
!   Arguments:
!   ----------
!
!   inputs:
!   ------
!     ngrid            size of the horizontal grid
!     pg               gravity (m s -2)
!     pz(ngrid,nlay)   height of layers
!     pu(ngrid,nlay)   u component of the wind
!     pv(ngrid,nlay)   v component of the wind
!     pts(ngrid)       surface temperature
!     ph(ngrid)        potential temperature T*(p/ps)^kappa
!
!
!=======================================================================
!
!-----------------------------------------------------------------------
!   Declarations:
!   -------------

#include "comcstfi.h"

!   Arguments:
!   ----------

      INTEGER, INTENT(IN) :: ngrid,nlay
      REAL, INTENT(IN) :: pz0(ngrid)
      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
      REAL, INTENT(IN) :: wmax(ngrid)
      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
      REAL, INTENT(INOUT) :: z_out                  ! output height (in m above surface)
      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)! interpolated fields at z_out : potential temperature and norm(uv)
      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta*

!   Local:
!   ------

      INTEGER ig

      REAL karman,nu
      DATA karman,nu/.41,0.001/
      SAVE karman,nu

!    Local(2):
!    ---------

      REAL rib(ngrid)  ! Bulk Richardson number
      REAL fm(ngrid) ! stability function for momentum
      REAL fh(ngrid) ! stability function for heat
      REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T

! phim = 1+betam*zeta   or   (1-bm*zeta)**am
! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
      REAL betam, betah, alphah, bm, bh, lambda
! ah and am are assumed to be -0.25 and -0.5 respectively

      REAL cdn(ngrid),chn(ngrid)  ! neutral momentum and heat drag coefficient
      REAL pz0t        ! initial thermal roughness length. (local)
      REAL ric         ! critical richardson number
      REAL reynolds(ngrid)    ! reynolds number for UBL
      REAL prandtl(ngrid)     ! prandtl number for UBL
      REAL pz0tcomp(ngrid)     ! computed z0t
      REAL ite
      REAL residual
      REAL pcdv(ngrid),pcdh(ngrid)
! For output :

      REAL zu2(ngrid)                  ! Large-scale wind at first layer
      REAL L_mo(ngrid)                ! Monin-Obukhov length
!-----------------------------------------------------------------------
!   couche de surface:
!   ------------------
      
      tstar(:)=0.
      ustar(:)=0.
      reynolds(:)=0.

! New formulation (AC) :

! phim = 1+betam*zeta   or   (1-bm*zeta)**am
! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
! am=-0.25, ah=-0.5

      pz0t = 0.     ! for the sake of simplicity
      pz0tcomp(:) = 0.1*pz0(:)
      rib(:)=0.
      pcdv(:)=0.
      pcdh(:)=0.

! this formulation assumes alphah=1., implying betah=betam
! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
      bm=16.            !UBL
      bh=16.            !UBL
      alphah=1.
      betam=5.         !SBL
      betah=5.         !SBL
      lambda=(sqrt(bh/bm))/alphah
      ric=betah/(betam**2)

      DO ig=1,ngrid

         ite=0.
         residual=abs(pz0tcomp(ig)-pz0t)

         do while((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))

         pz0t=pz0tcomp(ig)

         if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then

! Classical Richardson number formulation

!         rib(ig) = (pg/ph(ig,1))*pz(ig,1)*(ph(ig,1)-pts(ig))
!     &           /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))

! Richardson number formulation proposed by D.E. England et al. (1995)

!        IF((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) .lt. 1.) 
!     &            .and. (wmax(ig) .gt. 0.)) THEN
          zu2(ig)=
!     &  (MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
     &  ( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
!     &  pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
!        ELSE
!        zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
!        ENDIF

          rib(ig) = (pg/ph(ig,1))
!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
     &      *sqrt(pz(ig,1)*pz0(ig))
     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
     &      *(ph(ig,1)-pts(ig))/zu2(ig)

!     &  /(MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
!     &  /( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)

         else
         print*,'warning, infinite Richardson at surface'
         print*,pu(ig,1),pv(ig,1)
         rib(ig) = ric          ! traiter ce cas !
         endif

         z1z0=pz(ig,1)/pz0(ig)
         z1z0t=pz(ig,1)/pz0t

         cdn(ig)=karman/log(z1z0)
         cdn(ig)=cdn(ig)*cdn(ig)
         chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) 

! Stable case :
      if (rib(ig) .gt. 0.) then
! From D.E. England et al. (95)
      prandtl(ig)=1.
         if(rib(ig) .lt. ric) then
! Assuming alphah=1. and bh=bm for stable conditions :
            fm(ig)=((ric-rib(ig))/ric)**2
            fh(ig)=((ric-rib(ig))/ric)**2
         else
! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
            fm(ig)=0.
            fh(ig)=0.
         endif
! Unstable case :
      else
! From D.E. England et al. (95)
         fm(ig)=sqrt(1.-lambda*bm*rib(ig))
         fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*
     &                     (1.-lambda*bm*rib(ig))**0.25
         prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/
     &             ((1.-lambda*bh*rib(ig))**0.5)
      endif

       reynolds(ig)=karman*sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2)
     &              *pz0(ig)/(log(z1z0)*nu)
       pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3*
     &              (reynolds(ig)**0.25)*(prandtl(ig)**0.5))

          
      residual = abs(pz0t-pz0tcomp(ig))
      ite = ite+1
!      if(ite .eq. 10) then
!        print*, 'iteration max reached'
!      endif
!      print*, "iteration nnumber, residual",ite,residual

      enddo  ! of while

          pz0t=0.

! Drag computation :

         pcdv(ig)=cdn(ig)*fm(ig)
         pcdh(ig)=chn(ig)*fh(ig)
       
      ENDDO

! Large-scale wind at first layer (without gustiness) and
! u* theta* computation
      DO ig=1,ngrid

         if (rib(ig) .gt. ric) then
           ustar(ig)=0.
           tstar(ig)=0.
         else

!           ustar(ig)=karman*sqrt(fm(ig)*zu2(ig))/(log(pz(ig,1)/pz0(ig)))
!           tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
!     &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))

!simpler definition of u* and teta*, equivalent to the formula above :

            ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))                   
            tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))/sqrt(pcdv(ig))

           if (tstar(ig) .lt. -50) then
             print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
     &               ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
           endif 
         endif
      ENDDO

! Monin Obukhov length :

      DO ig=1,ngrid
         if (rib(ig) .gt. ric) then
            L_mo(ig)=0.
         else
            L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig))  !as defined here, L is positive for SBL, negative for UBL
         endif
      ENDDO

      DO ig=1,ngrid
      IF(z_out .ge. pz(ig,1)) THEN
           z_out=1.
           print*, 'Warning, z_out should be less than the first 
     &          vertical grid-point'
           print*, 'z1 =',pz(ig,1)
           print*, 'z_out=',z_out
           print*, 'z_out has been changed to 1m 
     &                  and computation is resumed'
      ENDIF

      IF(z_out .lt. pz0(ig)) THEN
           z_out=1.
           print*, 'Warning, z_out should be more than the roughness
     &          length'
           print*, 'z0 =',pz0(ig)
           print*, 'z_out=',z_out
           print*, 'z_out has been changed to z0
     &                  and computation is resumed'
      ENDIF
      ENDDO

      print*, 'interpolation of u and teta at z_out=',z_out

      DO ig=1,ngrid
        IF (L_mo(ig) .gt. 0.) THEN
           u_out(ig)=(ustar(ig)/karman)*log(z_out/pz0(ig)) + 
     &            5.*(ustar(ig)/(karman*L_mo(ig)))*(z_out-pz0(ig))
           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
     &                            *log(z_out/pz0tcomp(ig)) +
     &            5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig)))
     &                            *(z_out-pz0tcomp(ig))
        ELSEIF (L_mo(ig) .lt. 0.) THEN

          IF(L_mo(ig) .gt. -1000.) THEN
           
           u_out(ig)=(ustar(ig)/karman)*(
     &  2.*atan((1.-16.*z_out/L_mo(ig))**0.25)
     & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
     & -2.*log(1.+(1.-16.*z_out/L_mo(ig))**0.25)
     & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
     & -   log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
     & +   log(z_out/pz0(ig))
     &                                  )

           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
     &  2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
     & -2.*log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
     & +   log(z_out/pz0tcomp(ig))
     &                                                        )

          ELSE

! We have to treat the case where L is very large and negative,
! corresponding to a very nearly stable atmosphere (but not quite !)
! i.e. teta* <0 and almost zero. We choose the cutoff
! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
! between the two expression for z/L = -1/1000 is -1.54324*10^-9
! (we do that to avoid using r*4 precision, otherwise, we get -inf values)          

           u_out(ig)=(ustar(ig)/karman)*(
     &     (4./L_mo(ig))*(z_out-pz0(ig))
     &   + (20./(L_mo(ig))**2)*(z_out**2-pz0(ig)**2)
     &   + (160./(L_mo(ig))**3)*(z_out**3-pz0(ig)**3)
     &   + log(z_out/pz0(ig))
     &                                   )

           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
     &     (8./L_mo(ig))*(z_out-pz0tcomp(ig))
     &   + (48./(L_mo(ig))**2)*(z_out**2-pz0tcomp(ig)**2)
     &   + (1280./(3.*(L_mo(ig))**3))*(z_out**3-pz0tcomp(ig)**3)
     &   + log(z_out/pz0tcomp(ig))
     &                                                           )

          ENDIF
        ELSE
           u_out(ig)=0.
           Teta_out(ig)=pts(ig)
        ENDIF
      ENDDO

! Usefull diagnostics for the interpolation routine :

         call WRITEDIAGFI(ngrid,'L',
     &              'Monin Obukhov length','m',
     &                         2,L_mo)
         call WRITEDIAGFI(ngrid,'z0T',
     &              'thermal roughness length','m',
     &                         2,pz0tcomp)
         call WRITEDIAGFI(ngrid,'z0',
     &              'roughness length','m',
     &                         2,pz0)


      RETURN
      END