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

CMakeLists.txt

@@ -88,6 +88,7 @@ set(SOURCES_F
     srcF/sfx_most.f90
     srcF/sfx_most_param.f90
     srcF/sfx_sheba.f90
+    srcF/sfx_sheba_noniterative.f90
     srcF/sfx_sheba_param.f90
     srcF/sfx_sheba_noniterative.f90
     srcF/sfx_sheba_noit_param.f90

srcF/sfx_sheba_noniterative.f90

 #include "../includeF/sfx_def.fi"
 
 module sfx_sheba_noniterative
+<<<<<<< HEAD
     !< @brief main Earth System Model surface flux module
+=======
+    !< @brief SHEBA surface flux module
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
 
     ! modules used
     ! --------------------------------------------------------------------------------
@@ -10,7 +14,12 @@ module sfx_sheba_noniterative
 #endif
     use sfx_data
     use sfx_surface
+<<<<<<< HEAD
     use sfx_sheba_noit_param
+=======
+    use sfx_sheba_param
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
 #if defined(INCLUDE_CXX)
     use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
     use C_FUNC
@@ -27,16 +36,24 @@ module sfx_sheba_noniterative
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+<<<<<<< HEAD
+=======
+    public :: get_psi
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
     type, public :: numericsType
+<<<<<<< HEAD
         integer :: maxiters_convection = 10    !< maximum (actual) number of iterations in convection
+=======
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         integer :: maxiters_charnock = 10      !< maximum (actual) number of iterations in charnock roughness
     end type
     ! --------------------------------------------------------------------------------
 
 #if defined(INCLUDE_CXX)
+<<<<<<< HEAD
     type, BIND(C), public :: sfx_sheba_noit_param_C 
         real(C_FLOAT) :: kappa           
         real(C_FLOAT) :: Pr_t_0_inv
@@ -52,10 +69,27 @@ module sfx_sheba_noniterative
 
     type, BIND(C), public :: sfx_sheba_noit_numericsType_C 
         integer(C_INT) :: maxiters_convection           
+=======
+    type, BIND(C), public :: sfx_sheba_param_C
+        real(C_FLOAT) :: kappa           
+        real(C_FLOAT) :: Pr_t_0_inv
+
+        real(C_FLOAT) :: alpha_m           
+        real(C_FLOAT) :: alpha_h
+        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_numericsType_C 
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         integer(C_INT) :: maxiters_charnock 
     end type
 
     INTERFACE
+<<<<<<< HEAD
         SUBROUTINE c_sheba_noit_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, & 
             name="c_sheba_noit_compute_flux")
             use sfx_data
@@ -70,11 +104,28 @@ module sfx_sheba_noniterative
             type(sfx_sheba_noit_numericsType_C) :: numerics
             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
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
     END INTERFACE
 #endif 
 
 contains
 
+<<<<<<< HEAD
     ! --------------------------------------------------------------------------------
 #if defined(INCLUDE_CXX)
     subroutine set_c_struct_sfx_sheba_noit_param_values(sfx_model_param)
@@ -92,6 +143,25 @@ contains
     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
+#endif
+
+    ! --------------------------------------------------------------------------------
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
     subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
         !< @brief surface flux calculation for array data
         !< @details contains C/C++ & CUDA interface
@@ -112,18 +182,27 @@ contains
         type (meteoDataVecTypeC), target :: meteo_c         !< meteorological data (input)
         type (sfxDataVecTypeC), target :: sfx_c             !< surface flux data (output)
         type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
+<<<<<<< HEAD
         type (sfx_sheba_noit_param_C) :: model_param
         type (sfx_surface_param) :: surface_param
         type (sfx_sheba_noit_numericsType_C) :: numerics_c
         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
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         numerics_c%maxiters_charnock = numerics%maxiters_charnock
 
         phys_constants%Pr_m = Pr_m;
         phys_constants%nu_air = nu_air;
         phys_constants%g = g;
 
+<<<<<<< HEAD
         call set_c_struct_sfx_sheba_noit_param_values(model_param)
         call set_c_struct_sfx_surface_param_values(surface_param)
         call set_meteo_vec_c(meteo, meteo_c)
@@ -136,6 +215,19 @@ contains
         do i = 1, n
 #ifdef SFX_FORCE_DEPRECATED_sheba_CODE
 #else
+=======
+        call set_c_struct_sfx_sheba_param_values(model_param)
+        call set_c_struct_sfx_surface_param_values(surface_param)
+        call set_meteo_vec_c(meteo, meteo_c)
+        call set_sfx_vec_c(sfx, sfx_c)
+
+        meteo_c_ptr = C_LOC(meteo_c)
+        sfx_c_ptr   = C_LOC(sfx_c)
+
+        call c_sheba_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
+#else
+        do i = 1, n
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
@@ -144,10 +236,15 @@ contains
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
             call push_sfx_data(sfx, sfx_cell, i)
+<<<<<<< HEAD
 #endif
         end do
 #endif
 
+=======
+        end do
+#endif
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
     end subroutine get_surface_fluxes_vec
     ! --------------------------------------------------------------------------------
 
@@ -165,6 +262,10 @@ contains
         type (meteoDataType), intent(in) :: meteo
         type (numericsType), intent(in) :: numerics
         ! ----------------------------------------------------------------------------
+<<<<<<< HEAD
+=======
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         ! --- meteo derived datatype name shadowing
         ! ----------------------------------------------------------------------------
         real :: h       !< constant flux layer height [m]
@@ -187,6 +288,7 @@ contains
         real zeta               !< = z/L [n/d]
         real Rib                !< bulk Richardson number
 
+<<<<<<< HEAD
         real zeta_conv_lim      !< z/L critical value for matching free convection limit [n/d]
         real Rib_conv_lim       !< Ri-bulk critical value for matching free convection limit [n/d]
 
@@ -195,12 +297,17 @@ contains
 
         real psi_m, psi_h       !< universal functions (momentum) & (heat) [n/d]
         real psi0_m, psi0_h       !< universal functions (momentum) & (heat) [n/d]
+=======
+        real Udyn, Tdyn, Qdyn   !< dynamic scales
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         real phi_m, phi_h       !< stability functions (momentum) & (heat) [n/d]
 
         real Km                 !< eddy viscosity coeff. at h [m^2/s]
         real Pr_t_inv           !< invese Prandt number [n/d]
 
         real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
+<<<<<<< HEAD
         
         real Udyn, Tdyn
 
@@ -209,6 +316,11 @@ contains
         real fval               !< just a shortcut for partial calculations
 
         real :: C1,A1,A2,lne,lnet,Ribl
+=======
+
+        integer surface_type    !< surface type = (ocean || land)
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
 
 #ifdef SFX_CHECK_NAN
         real NaN
@@ -267,6 +379,7 @@ contains
         ! --- define Ri-bulk
         Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
 
+<<<<<<< HEAD
         ! --- define free convection transition zeta = z/L value
         call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
                 h0_m, h0_t, B)
@@ -346,19 +459,53 @@ contains
             Ct = kappa / psi_h
         end if
         
+=======
+        ! --- get the fluxes
+        ! ----------------------------------------------------------------------------
+        if(Rib > 0)then
+                call get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
+                U, dT, dQ, h, z0_m, z0_t, Rib)
+        else
+                call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+                U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+        end if
+
+        ! ----------------------------------------------------------------------------
+
+        call get_phi(phi_m, phi_h, zeta)
+        ! ----------------------------------------------------------------------------
+
+        ! --- define transfer coeff. (momentum) & (heat)
+        Cm = 0.0
+        if (U > 0.0) then
+            Cm = Udyn / U
+        end if
+        Ct = 0.0
+        if (abs(dT) > 0.0) then
+            Ct = Tdyn / dT
+        end if
+
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
         ! --- define eddy viscosity & inverse Prandtl number
         Km = kappa * Cm * U * h / phi_m
         Pr_t_inv = phi_m / phi_h
 
         ! --- setting output
         sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+<<<<<<< HEAD
             Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
             Rib_conv_lim = Rib_conv_lim, &
             Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+=======
+                Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
+                Rib_conv_lim = 0.0, &
+                Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
 
     end subroutine get_surface_fluxes
     ! --------------------------------------------------------------------------------
 
+<<<<<<< HEAD
     ! convection universal functions shortcuts
     ! --------------------------------------------------------------------------------
     function f_m_conv(zeta)
@@ -581,6 +728,224 @@ contains
 
         ! --- bulk Richardson number
         Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
+=======
+    !< @brief get dynamic scales
+    ! --------------------------------------------------------------------------------
+    subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+            U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: Udyn, Tdyn, Qdyn   !< dynamic scales
+        real, intent(out) :: zeta               !< = z/L
+
+        real, intent(in) :: U                   !< abs(wind speed) at z
+        real, intent(in) :: Tsemi               !< semi-sum of temperature at z and at surface
+        real, intent(in) :: dT, dQ              !< temperature & humidity difference between z and at surface
+        real, intent(in) :: z                   !< constant flux layer height
+        real, intent(in) :: z0_m, z0_t          !< roughness parameters
+        real, intent(in) :: beta                !< buoyancy parameter
+
+        integer, intent(in) :: maxiters         !< maximum number of iterations
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        real, parameter :: gamma = 0.61
+
+        real :: psi_m, psi_h
+        real :: psi0_m, psi0_h
+        real :: Linv
+        integer :: i
+        ! ----------------------------------------------------------------------------
+
+
+        Udyn = kappa * U / log(z / z0_m)
+        Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t)
+        Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t)
+        zeta = 0.0
+
+        ! --- no wind
+        if (Udyn < 1e-5) return
+
+        Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
+        zeta = z * Linv
+
+        ! --- near neutral case
+        if (Linv < 1e-5) return
+
+        do i = 1, maxiters
+
+            call get_psi(psi_m, psi_h, zeta)
+            call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)
+
+            Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
+            Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
+            Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
+
+            if (Udyn < 1e-5) exit
+
+            Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
+            zeta = z * Linv
+        end do
+
+    end subroutine get_dynamic_scales
+
+
+    subroutine get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
+            U, dT, dQ, z, z0_m, z0_t, Rib)
+        ! ----------------------------------------------------------------------------
+        real, parameter  ::  gamma = 2.91, zeta_a = 3.6
+
+        real, intent(out) :: Udyn, Tdyn, Qdyn   !< dynamic scales
+        real, intent(out) :: zeta               !< = z/L
+
+        real, intent(in) :: U                   !< abs(wind speed) at z
+        real, intent(in) :: dT, dQ              !< temperature & humidity difference between z and at surface
+        real, intent(in) :: z                   !< constant flux layer height
+        real, intent(in) :: z0_m, z0_t          !< roughness parameters
+        real, intent(in) :: Rib                 !< bulk Richardson number
+
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        real :: psi_m, psi_h
+        real :: psi0_m, psi0_h
+        real :: C1,A1,A2,lne,lnet,Ribl
+        ! ----------------------------------------------------------------------------
+
+        Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / z) / ((1 - z0_m / z)**2)
+
+        call get_psi(psi_m, psi_h, zeta_a)
+        call get_psi_mh(psi0_m, psi0_h, zeta_a * z0_m / z,  zeta_a * z0_t / z)
+
+        lne = log(z/z0_m)
+        lnet = log(z/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)
+
+        call get_psi(psi_m, psi_h, zeta)
+        call get_psi_mh(psi0_m, psi0_h, zeta * z0_m / z, zeta * z0_t /z)
+
+        Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
+        Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
+        Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
+
+
+    end subroutine get_dynamic_scales_noniterative
+    ! --------------------------------------------------------------------------------
+
+    ! stability functions
+    ! --------------------------------------------------------------------------------
+    subroutine get_phi(phi_m, phi_h, zeta)
+        !< @brief stability functions (momentum) & (heat): neutral case
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: phi_m, phi_h   !< stability functions
+
+        real, intent(in) :: zeta            !< = z/L
+        ! ----------------------------------------------------------------------------
+
+
+        if (zeta >= 0.0) then
+            phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
+            phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
+        else
+            phi_m = (1.0 - alpha_m * zeta)**(-0.25)
+            phi_h = (1.0 - alpha_h * zeta)**(-0.5)
+        end if
+
+    end subroutine
+    ! --------------------------------------------------------------------------------
+
+    ! universal functions
+    ! --------------------------------------------------------------------------------
+    subroutine get_psi(psi_m, psi_h, zeta)
+        !< @brief universal functions (momentum) & (heat): neutral case
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: psi_m, psi_h   !< universal functions
+
+        real, intent(in) :: zeta            !< = z/L
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        real :: x_m, x_h
+        real :: q_m, q_h
+        ! ----------------------------------------------------------------------------
+
+
+        if (zeta >= 0.0) then
+
+            q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
+            q_h = sqrt(c_h * c_h - 4.0)
+
+            x_m = (1.0 + zeta)**(1.0 / 3.0)
+            x_h = zeta
+
+            psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
+                    2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
+                            log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
+                            2.0 * sqrt(3.0) * (&
+                                    atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
+                                            atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
+
+            psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
+                    ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
+                            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)))
+        else
+            x_m = (1.0 - alpha_m * zeta)**(0.25)
+            x_h = (1.0 - alpha_h * zeta)**(0.25)
+
+            psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
+            psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
+        end if
+
+    end subroutine
+
+
+    subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h)
+        !< @brief universal functions (momentum) & (heat): neutral case
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: psi_m, psi_h   !< universal functions
+
+        real, intent(in) :: zeta_m, zeta_h  !< = z/L
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        real :: x_m, x_h
+        real :: q_m, q_h
+        ! ----------------------------------------------------------------------------
+
+
+        if (zeta_m >= 0.0) then
+            q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
+            x_m = (1.0 + zeta_m)**(1.0 / 3.0)
+
+            psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
+                    2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
+                            log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
+                            2.0 * sqrt(3.0) * (&
+                                    atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
+                                            atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
+        else
+            x_m = (1.0 - alpha_m * zeta_m)**(0.25)
+            psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
+        end if
+
+        if (zeta_h >= 0.0) then
+            q_h = sqrt(c_h * c_h - 4.0)
+            x_h = zeta_h
+
+            psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
+                    ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
+                            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)))
+        else
+            x_h = (1.0 - alpha_h * zeta_h)**(0.25)
+            psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
+        end if
+>>>>>>> 9d99a415378a2907d460477f87825d027fae071e
 
     end subroutine
     ! --------------------------------------------------------------------------------