Merge branch 'diag_new' of http://tesla.parallel.ru/INMCM60_pbl/INMCM60_sfx into diag_new

diag/pbldia_new_sfx.f90

@@ -2,11 +2,11 @@ module pbldia_new_sfx
 
 
 implicit none 
-real,parameter,private::undef=-999.0,R=287.,appa=0.287
+real,parameter,private::kappa=0.4,R=287.,appa=0.287
 
 contains
 
-! in this subroutine interpolation is performed using integration of universal functions from z1 to z2, 
+! in this module interpolation is performed using integration of universal functions from z1 to z2, 
 ! where z1 is measurement or lowest model level height and z2 is required height
 ! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx)
 !C*    DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND),
@@ -16,8 +16,6 @@ contains
 !C*    ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE)
 !C*    AR2(1)    - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH)
 !C*                HEIGHT OF CONSTANT FLUX LAYER TOP
-!C*    AR2(5)    - ROUGHNESS LENGTH FOR MOMENTUM
-!C*    AR2(6)    - ROUGHNESS LENGTH FOR TEMPERATURE AND HUMIDITY
 !C*    AR2(8)    - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM
 !C*    AR2(9)    - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY
 !C*    ARRAY ARDIN
@@ -29,32 +27,19 @@ contains
 !C*    ARDIN(6)  - = HTEMP
 !C*    OUTPUT DATA (ARRAY ARDOUT):
 !C*    ARDOUT(1)  - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT
-!C*    ARDOUT(2)  - POTENTIAL TEMPERATURE at REQUIRED HEIGHT
-!C*    ARDOUT(3)  - SPECIFIC HUMIDITY at REQUIRED HEIGHT
+!C*    ARDOUT(2)  - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
+!C*    ARDOUT(3)  - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
 !C*    UFWIND     - UNIVERSAL FUNCTION FOR WIND VELOCITY
 !C*    UFTEMP     - UNIVERSAL FUNCTION FOR SCALARS
-!     AR2(1) - NON-DIMENSIONAL CFL HIGHT                            =
-!     AR2(2) - RICHARDSON NUMBER                                    =
-!     AR2(3) - REYNODS NUMBER                                       =
-!     AR2(4) - LN(ZU/ZT)                                            =
-!     AR2(5) - DYNAMICAL ROUGHNESS ZU (M)                           =
-!     AR2(6) - THERMAL   ROUGHNESS ZT (M)                           =
-!     AR2(7) - CRITICAL RICHARDSON NUMBER                           =
-!     AR2(8) - TRANSFER COEFFICIENT FOR MOMENTUM                    =
-!     AR2(9) - TRANSFER COEFFICIENT FR HEAT                         =
-!     AR2(10)- COEFFICIENT OF TURBULENCE (KM) AT CFL HIGHT (M**2/S) =
-!     AR2(11)- ALFT=KT/KM ( KT-COEFFICIENT OF TURBULENCE FOR HEAT)  =
 
 
 subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
 use sfx_sheba, only: &
                 get_psi_sheba => get_psi
 
-      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
-      real,intent(out)::ARDOUT(3) !ARDOUT(4)
-!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
+      real,intent(in):: AR2(11),ARDIN(6)
+      real,intent(out)::ARDOUT(3)
       real,parameter::zetalim = 2. !maximum value of z/L for stable SL
-      real,parameter::kappa = 0.4 
       real psi_m,psi_h,psi_m_hs,psi_h_hs
       real hwind,htemp,ustar,tstar,qstar,&
 & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
@@ -70,8 +55,7 @@ use sfx_sheba, only: &
 
 
      if(zeta.gt.zetalim) zeta=zetalim
-!     if(L.lt.Llim.and.L.gt.0) L=Llim
-      L = HIN/zeta  ! MO length
+      L = HIN/zeta
 
      call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
      call get_psi_sheba(psi_m,psi_h,HWIND/L)
@@ -89,18 +73,16 @@ use sfx_sheba, only: &
       ARDOUT(3) = DQ+ardin(3)
 
 
-      RETURN
-  END subroutine pbldia_new_sheba
+      return
+  end subroutine pbldia_new_sheba
 
 subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
 use sfx_most, only: &
                 get_psi_most => get_psi
 
-      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
-      real,intent(out)::ARDOUT(3) !ARDOUT(4)
-!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
+      real,intent(in):: AR2(11),ARDIN(6)
+      real,intent(out)::ARDOUT(3)
       real,parameter::zetalim = 2. !maximum value of z/L for stable SL
-      real,parameter::kappa = 0.4 
       real psi_m,psi_h,psi_m_hs,psi_h_hs
       real hwind,htemp,ustar,tstar,qstar,&
 & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
@@ -115,8 +97,7 @@ use sfx_most, only: &
       QSTAR = AR2(9) * ARDIN(3)
 
      if(zeta.gt.zetalim) zeta=zetalim
-!     if(L.lt.Llim.and.L.gt.0) L=Llim
-      L = HIN/zeta  ! MO length
+      L = HIN/zeta
 
      call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
      call get_psi_most(psi_m,psi_h,HWIND/L)
@@ -124,10 +105,6 @@ use sfx_most, only: &
      WIND = (USTAR/kappa) * UFWIND
      call get_psi_most(psi_m,psi_h,HTEMP/L)
       UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
-!  if(L>0)then
-!   !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
-!      UFTEMP = ALOG(HTEMP/HIN) 
-!  endif
 
       DTETA = (TSTAR/kappa) * UFTEMP
       DQ = (QSTAR/kappa) * UFTEMP
@@ -137,17 +114,13 @@ use sfx_most, only: &
       ARDOUT(3) = DQ+ardin(3)
 
 
-      RETURN
-  END subroutine pbldia_new_most
+      return
+  end subroutine pbldia_new_most
 
 subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
-!use sfx_esm, only: &
-!                get_psi_esm => get_psi
-      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
-      real,intent(out)::ARDOUT(3) !ARDOUT(4)
-!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
+      real,intent(in):: AR2(11),ARDIN(6)
+      real,intent(out)::ARDOUT(3)
       real,parameter::zetalim = 2. !maximum value of z/L for stable SL
-      real,parameter::kappa = 0.4 
       real psi_m,psi_h,psi_m_hs,psi_h_hs
       real hwind,htemp,ustar,tstar,qstar,&
 & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
@@ -162,8 +135,7 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
       QSTAR = AR2(9) * ARDIN(3)
 
      if(zeta.gt.zetalim) zeta=zetalim
-!     if(L.lt.Llim.and.L.gt.0) L=Llim
-      L = HIN/zeta  ! MO length
+      L = HIN/zeta
 
 
      call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
@@ -171,10 +143,6 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
      WIND = (USTAR/kappa) * UFWIND
      call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
       UFTEMP = psi_h
-!  if(L>0)then
-!   !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
-!      UFTEMP = ALOG(HTEMP/HIN) 
-!  endif
 
       DTETA = (TSTAR/kappa) * UFTEMP
       DQ = (QSTAR/kappa) * UFTEMP
@@ -184,52 +152,53 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
       ARDOUT(3) = DQ+ardin(3)
 
 
-      RETURN
-  END subroutine pbldia_new_esm
+      return
+  end subroutine pbldia_new_esm
 
 
 subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
+
 use sfx_esm_param
-    !< @brief universal functions (momentum) & (heat): neutral case
-    ! ----------------------------------------------------------------------------
-    real, intent(out) :: psi_m, psi_h   !< universal functions
-    real, intent(in) :: z1,z2,L            !
 
-  real an5,D1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2
+real, intent(out) :: psi_m, psi_h
+real, intent(in) :: z1,z2,L
+
+real an5,d1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2
 
       an5 = (Pr_t_inf_inv/Pr_t_0_inv)**4
-      D1=(2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*&
+      d1=(2.0E0*alpha_h-an5*alpha_m-sqrt((an5*alpha_m)**2+4.0E0*an5*alpha_h*&
 &        (alpha_h-alpha_m)))/(2.0E0*alpha_h**2)
 
 
       zeta1 = z1 / L
       zeta2 = z2 / L
 
-      IF(zeta2 .GE. 0.) THEN
+      if(zeta2 .ge. 0.) then
 
       psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
       psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)
 
-      ELSE
+      else
 
-      GMZ2 = SQRT(SQRT(1.E0 - alpha_m * zeta2))
-      GMZ1 = SQRT(SQRT(1.E0 - alpha_m * zeta1))
-      GMST = SQRT(SQRT(1.E0 - alpha_m * D1))
-      GTZ2 = SQRT(1.E0 - alpha_h * zeta2)
-      GTZ1 = SQRT(1.E0 - alpha_h * zeta1)
-      GTST = SQRT(1.E0 - alpha_h * D1)
+      gmz2 = sqrt(sqrt(1.E0 - alpha_m * zeta2))
+      gmz1 = sqrt(sqrt(1.E0 - alpha_m * zeta1))
+      gmst = sqrt(sqrt(1.E0 - alpha_m * D1))
+      gtz2 = sqrt(1.E0 - alpha_h * zeta2)
+      gtz1 = sqrt(1.E0 - alpha_h * zeta1)
+      gtst = sqrt(1.E0 - alpha_h * D1)
 
-      IF(zeta2 .GE. D1) THEN
-      psi_m = FIM(GMZ2) - FIM(GMZ1)
-      psi_h = FIT(GTZ2) - FIT(GTZ1)
-      ELSE
-      psi_m = (3.E0/GMST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIM(GMST) - FIM(GMZ1)
-      psi_h = (3.E0/GTST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIT(GTST) - FIT(GTZ1)
-      ENDIF
-      ENDIF
+      if(zeta2 .ge. d1) then
+      psi_m = fim(gmz2) - fim(gmz1)
+      psi_h = fit(gtz2) - fit(gtz1)
+      else
+      psi_m = (3.E0/gmst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fim(gmst) - fim(gmz1)
+      psi_h = (3.E0/gtst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fit(gtst) - fit(gtz1)
+      endif
+      endif
 
-end subroutine
+end subroutine get_psi_esm1
 
+! functions fim,fit were copied from the original pbldia
 REAL FUNCTION FIM(XX)
   IMPLICIT NONE
   REAL X, XX
@@ -246,6 +215,4 @@ end subroutine
   RETURN
 END function
 
-
-
 end module pbldia_new_sfx
\ No newline at end of file