some changes to pbldia_new and sfx_sheba_noniterative

diag/compile_simple.sh

@@ -2,6 +2,7 @@
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_phys_const.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_esm_param.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba_param.f90
+gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba_noit_param.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_most_param.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_log_param.f90
 
@@ -10,6 +11,7 @@ gfortran -c -cpp -Wuninitialized ../srcF/sfx_surface.f90
 
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_esm.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba.f90
+gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba_noniterative.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_most.f90
 gfortran -c -cpp -Wuninitialized ../srcF/sfx_log.f90
 
@@ -17,5 +19,5 @@ gfortran -c -cpp -Wuninitialized pbldia_new_sfx.f90
 
 gfortran -c -cpp -Wuninitialized diag_sfc_simple.f90
 
-gfortran -o a.exe sfx_phys_const.o sfx_sheba_param.o sfx_esm_param.o sfx_most_param.o sfx_log_param.o sfx_data.o sfx_surface.o sfx_sheba.o sfx_esm.o sfx_most.o sfx_log.o pbldia_new_sfx.o diag_sfc_simple.o
+gfortran -o a.exe sfx_phys_const.o sfx_sheba_param.o sfx_sheba_noit_param.o sfx_esm_param.o sfx_most_param.o sfx_log_param.o sfx_data.o sfx_surface.o sfx_sheba.o sfx_sheba_noniterative.o sfx_esm.o sfx_most.o sfx_log.o pbldia_new_sfx.o diag_sfc_simple.o
 rm *.o *.mod

diag/diag_sfc_simple.f90

@@ -8,6 +8,9 @@ use sfx_data
 use sfx_sheba, only: &
         get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
         numericsType_sheba => numericsType
+use sfx_sheba_noniterative, only: &
+        get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
+        numericsType_sheba_noit => numericsType
 use sfx_esm, only: &
         get_surface_fluxes_vec_esm => get_surface_fluxes_vec, &
         numericsType_esm => numericsType
@@ -26,6 +29,9 @@ use sfx_phys_const
 type(meteoDataVecType) :: meteo         !< meteorological data (input)
 type(meteoDataType) :: meteo_cell
 type(sfxDataVecType) :: sfx             !< surface fluxes (output)
+#ifdef SHEBA_NOIT
+type(numericsType_sheba_noit) :: numerics_sheba_noit      !< surface flux module (SHEBA) numerics parameters
+#endif
 #ifdef SHEBA
 type(numericsType_sheba) :: numerics_sheba      !< surface flux module (SHEBA) numerics parameters
 #endif
@@ -54,14 +60,14 @@ integer j,k
 
 ! user-defined input parameters for wind, temperature and humidity interpolation in the surface layer
 !mandatory:
-u = 5.0
+u = 10.
 Tin = 288.0
 Ts  = 283.0
 Qin = 0.005
 Qs  = 0.008
 z0m = 0.03  !momentum roughness length
 pot_temp = .true. !if inout temperatures are potential
-h = 25.      ! input height
+h = 50.      ! input height
 zw = 10     ! height to which wind interpolation is performed
 zt = 2     ! height to wich temperature and humidity interpolation is performed
 
@@ -133,6 +139,9 @@ meteo%z0_m = z0m
 #ifdef SHEBA
 call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, 1)
 #endif
+#ifdef SHEBA_NOIT
+call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, 1)
+#endif
 #ifdef ESM
 call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, 1)
 #endif
@@ -152,13 +161,14 @@ ar2(9) = sfx%Ct(1)
 else
 
 ar2(1) = h / L
-ar2(5) = z0m
-ar2(6) = z0h
+!ar2(5) = z0m
+!ar2(6) = z0h
 ar2(8) = cm
 ar2(9) = ct
 
 endif
 
+print *, sfx%zeta,-1*1005.*sfx%Cm(1)*sfx%Ct(1)*u*dt*1.3
 
 ardin(1) = u
 ardin(2) = dt
@@ -175,6 +185,9 @@ ardin(6) = zt
 #ifdef SHEBA
 call pbldia_new_sheba(ar2(:),ardin(:),ardout(:))
 #endif
+#ifdef SHEBA_NOIT
+call pbldia_new_sheba_noit(ar2(:),ardin(:),ardout(:))
+#endif
 #ifdef MOST
 call pbldia_new_most(ar2(:),ardin(:),ardout(:))
 #endif

diag/pbldia_new_sfx.f90

@@ -84,13 +84,78 @@ use sfx_sheba, only: &
       return
   end subroutine pbldia_new_sheba
 
+subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT)
+use sfx_sheba_noniterative, only: &
+                get_psi_stable_sheba => get_psi_stable
+
+      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 psi_m,psi_h,psi_m_hs,psi_h_hs
+      real hwind,htemp,ustar,tstar,qstar,&
+& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
+
+      HWIND = ARDIN(5)
+      HTEMP = ARDIN(6)
+      HIN = ardin(4)
+      zeta = AR2(1)
+
+      USTAR = AR2(8) * ARDIN(1)
+      TSTAR = AR2(9) * ARDIN(2)
+      QSTAR = AR2(9) * ARDIN(3)
+
+     if(zeta.ne.0)then
+
+     if(zeta.gt.zetalim) zeta=zetalim
+      L = HIN/zeta
+
+     if(zeta.gt.0)then ! stable stratification
+
+
+     call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
+     call get_psi_stable_sheba(psi_m,psi_h,HWIND/L,HWIND/L)
+     UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+     call get_psi_stable_sheba(psi_m,psi_h,HTEMP/L,HTEMP/L)
+     UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
+
+     else    ! unstable stratification (functions from ESM)
+
+     call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
+     UFWIND = psi_m
+     call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
+     UFTEMP = psi_h
+
+     endif
+
+     else
+
+! neutral stratification (z/L=0)
+
+      UFWIND = ALOG(HWIND/HIN)
+      UFTEMP = ALOG(HTEMP/HIN) 
+
+     endif
+
+      WIND = (USTAR/kappa) * UFWIND
+      DTETA = (TSTAR/kappa) * UFTEMP
+      DQ = (QSTAR/kappa) * UFTEMP
+
+
+      ARDOUT(1) = ARDIN(1)+WIND
+      ARDOUT(2) = DTETA+ardin(2)
+      ARDOUT(3) = DQ+ardin(3)
+
+
+      return
+  end subroutine pbldia_new_sheba_noit
+
 subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
 use sfx_most, only: &
                 get_psi_most => get_psi
 
       real,intent(in):: AR2(11),ARDIN(6)
       real,intent(out)::ARDOUT(3)
-      real,parameter::zetalim = 1. !maximum value of z/L for stable SL
+      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
       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
@@ -137,7 +202,7 @@ use sfx_most, only: &
 subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
       real,intent(in):: AR2(11),ARDIN(6)
       real,intent(out)::ARDOUT(3)
-      real,parameter::zetalim = 1. !maximum value of z/L for stable SL
+      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
       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

srcF/sfx_sheba_noniterative.f90

@@ -27,6 +27,7 @@ module sfx_sheba_noniterative
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+    public :: get_psi_stable
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -48,6 +49,11 @@ module sfx_sheba_noniterative
         real(C_FLOAT) :: Rib_max
         real(C_FLOAT) :: gamma
         real(C_FLOAT) :: zeta_a
+        real(C_FLOAT) ::  a_m
+        real(C_FLOAT) ::  b_m
+        real(C_FLOAT) ::  a_h
+        real(C_FLOAT) ::  b_h
+        real(C_FLOAT) ::  c_h
     end type
 
     type, BIND(C), public :: sfx_sheba_noit_numericsType_C 
@@ -70,21 +76,6 @@ module sfx_sheba_noniterative
             type(sfx_phys_constants) :: constants
         END SUBROUTINE c_sheba_noit_compute_flux
 
-        SUBROUTINE c_sheba_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, & 
-            name="c_sheba_compute_flux")
-            use sfx_data
-            use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
-            Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C
-            implicit none
-            INTEGER(C_INT) :: grid_size
-            type(C_PTR), value :: sfx
-            type(C_PTR), value :: meteo
-            type(sfx_sheba_param_C) :: model_param
-            type(sfx_surface_param) :: surface_param
-            type(sfx_sheba_numericsType_C) :: numerics
-            type(sfx_phys_constants) :: constants
-        END SUBROUTINE c_sheba_compute_flux
-
     END INTERFACE
 #endif 
 
@@ -104,23 +95,12 @@ contains
         sfx_model_param%Rib_max = Rib_max
         sfx_model_param%gamma = gamma
         sfx_model_param%zeta_a = zeta_a
-    end subroutine set_c_struct_sfx_sheba_noit_param_values
-#endif
-
-#if defined(INCLUDE_CXX)
-    subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
-        type (sfx_sheba_param_C), intent(inout) :: sfx_model_param
-        sfx_model_param%kappa = kappa
-        sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
-
-        sfx_model_param%alpha_m = alpha_m
-        sfx_model_param%alpha_h = alpha_h
         sfx_model_param%a_m = a_m
         sfx_model_param%b_m = b_m
         sfx_model_param%a_h = a_h
         sfx_model_param%b_h = b_h
         sfx_model_param%c_h = c_h
-    end subroutine set_c_struct_sfx_sheba_param_values
+    end subroutine set_c_struct_sfx_sheba_noit_param_values
 #endif
 
     ! --------------------------------------------------------------------------------
@@ -150,12 +130,6 @@ contains
         type (sfx_phys_constants) :: phys_constants
 
         numerics_c%maxiters_convection = numerics%maxiters_convection
-
-        type (sfx_sheba_param_C) :: model_param
-        type (sfx_surface_param) :: surface_param
-        type (sfx_sheba_numericsType_C) :: numerics_c
-        type (sfx_phys_constants) :: phys_constants
-
         numerics_c%maxiters_charnock = numerics%maxiters_charnock
 
         phys_constants%Pr_m = Pr_m;
@@ -245,8 +219,6 @@ contains
 
         real fval               !< just a shortcut for partial calculations
 
-        real :: C1,A1,A2,lne,lnet,Ribl
-
 
 #ifdef SFX_CHECK_NAN
         real NaN
@@ -316,21 +288,9 @@ contains
 
             !   --- restrict bulk Ri value
             !   *: note that this value is written to output
-            Rib = min(Rib, Rib_max)
-
-            Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
-
-            call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
-            call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h,  zeta_a * z0_t / h)
-
-            lne = log(h/z0_m)
-            lnet = log(h/z0_t)
-            C1 = (lne**2)/lnet
-            A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
-&           / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
-            A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
+!            Rib = min(Rib, Rib_max)
 
-            zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
+            call get_zeta(zeta, Rib, h, z0_m, z0_t)
 
             call get_psi_stable(psi_m, psi_h, zeta, zeta)
             call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
@@ -397,6 +357,29 @@ contains
     end subroutine get_surface_fluxes
     ! --------------------------------------------------------------------------------
 
+    subroutine get_zeta(zeta, Rib, h, z0_m, z0_t)
+    real,intent(out) :: zeta
+    real,intent(in) :: Rib, h, z0_m, z0_t
+
+    real :: Ribl, C1, A1, A2, lne, lnet
+    real :: psi_m, psi_h, psi0_m, psi0_h
+
+        Ribl = (Rib * Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
+
+        call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
+        call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h,  zeta_a * z0_t / h)
+
+        lne = log(h/z0_m)
+        lnet = log(h/z0_t)
+        C1 = (lne**2)/lnet
+        A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
+&           / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
+        A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
+
+        zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
+
+    end subroutine get_zeta
+
     ! convection universal functions shortcuts
     ! --------------------------------------------------------------------------------
     function f_m_conv(zeta)
@@ -469,7 +452,7 @@ contains
                             log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
                                     log((c_h - q_h) / (c_h + q_h)))
 
-    end subroutine
+    end subroutine get_psi_stable
 
 
     subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)