added most_snow, most_snow_param

srcF/sfx_most_snow.f90

+#include "../includeF/sfx_def.fi"
+
+module sfx_most_snow
+    !< @brief MOST surface flux module
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+    use sfx_common
+#endif
+    use sfx_data
+    use sfx_surface
+    use sfx_most_snow_param
+    ! --------------------------------------------------------------------------------
+
+    ! directives list
+    ! --------------------------------------------------------------------------------
+    implicit none
+    private
+    ! --------------------------------------------------------------------------------
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: get_surface_fluxes
+    public :: get_surface_fluxes_vec
+    ! --------------------------------------------------------------------------------
+
+    ! --------------------------------------------------------------------------------
+    type, public :: numericsType
+        integer :: maxiters_charnock = 10      !< maximum (actual) number of iterations in charnock roughness
+        integer :: maxiters_snow = 10          !< maximum (actual) number of iterations in snow roughness
+    end type
+    ! --------------------------------------------------------------------------------
+
+contains
+
+    ! --------------------------------------------------------------------------------
+    subroutine get_surface_fluxes_vec(sfx, sfx2, meteo, numerics, n)
+        !< @brief surface flux calculation for array data
+        !< @details contains C/C++ & CUDA interface
+        ! ----------------------------------------------------------------------------
+        type (sfxDataVecType), intent(inout) :: sfx
+		type (sfxDataAddVecType), intent(inout) :: sfx2
+
+        type (meteoDataVecType), intent(in) :: meteo
+        type (numericsType), intent(in) :: numerics
+        integer, intent(in) :: n
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        type (meteoDataType)  meteo_cell
+        type (sfxDataType) sfx_cell
+		type (sfxDataAddType) sfx2_cell
+        integer i
+        ! ----------------------------------------------------------------------------
+
+        do i = 1, n
+            meteo_cell = meteoDataType(&
+                    h = meteo%h(i), &
+                    U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
+                    z0_m = meteo%z0_m(i))
+
+            call get_surface_fluxes(sfx_cell, sfx2_cell, meteo_cell, numerics)
+
+            call push_sfx_data(sfx, sfx_cell, i)
+            call push_sfx_data_add(sfx2, sfx2_cell, i)
+		end do
+
+    end subroutine get_surface_fluxes_vec
+    ! --------------------------------------------------------------------------------
+
+    ! --------------------------------------------------------------------------------
+    subroutine get_surface_fluxes(sfx, sfx2, meteo, numerics)
+        !< @brief surface flux calculation for single cell
+        !< @details contains C/C++ interface
+        ! ----------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+        use ieee_arithmetic
+#endif
+
+        type (sfxDataType), intent(out) :: sfx
+        type (sfxDataAddType), intent(out) :: sfx2
+        type (meteoDataType), intent(in) :: meteo
+        type (numericsType), intent(in) :: numerics
+        ! ----------------------------------------------------------------------------
+
+        ! --- meteo derived datatype name shadowing
+        ! ----------------------------------------------------------------------------
+        real :: h       !< constant flux layer height [m]
+        real :: U       !< abs(wind speed) at 'h' [m/s]
+        real :: dT      !< difference between potential temperature at 'h' and at surface [K]
+        real :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
+        real :: dQ      !< difference between humidity at 'h' and at surface [g/g]
+        real :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        ! ----------------------------------------------------------------------------
+
+        ! --- local variables
+        ! ----------------------------------------------------------------------------
+        real z0_t               !< thermal roughness [m]
+        real B                  !< = ln(z0_m / z0_t) [n/d]
+        real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
+
+        real u_dyn0             !< dynamic velocity in neutral conditions [m/s]
+        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
+
+        real zeta               !< = z/L [n/d]
+        real Rib                !< bulk Richardson number
+
+        real Udyn, Tdyn, Qdyn   !< dynamic scales
+
+        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]
+
+        ! --- local additional variables
+        ! ----------------------------------------------------------------------------
+		!real phi_m, phi_m
+		real hfx, mfx
+		real S_mean, Lsnow
+	    real z0_s, h_salt
+                
+	
+		integer surface_type    !< surface type = (ocean || land)
+
+
+#ifdef SFX_CHECK_NAN
+        real NaN
+#endif
+        ! ----------------------------------------------------------------------------
+
+#ifdef SFX_CHECK_NAN
+        ! --- checking if arguments are finite
+        if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
+                .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then
+
+            !NaN = ieee_value(0.0, ieee_quiet_nan)   ! setting NaN
+            sfx = sfxDataType(zeta = NaN, Rib = NaN, &
+                    Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
+                    Rib_conv_lim = NaN, &
+                    Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
+            
+            sfx2 = sfxDataAddType(phi_m = NaN, phi_h = NaN, &
+                hfx = NaN, mfx = NaN, Udyn = NaN, S_mean = NaN, &
+                Lsnow = NaN, &
+                z0_s = NaN, h_salt = NaN)
+            
+            
+            return
+        end if
+#endif
+
+        ! --- shadowing names for clarity
+        U = meteo%U
+        Tsemi = meteo%Tsemi
+        dT = meteo%dT
+        dQ = meteo%dQ
+        h = meteo%h
+        z0_m = meteo%z0_m
+
+        ! --- define surface type
+        if (z0_m < 0.0) then
+            surface_type = surface_ocean
+        else if (z0_m == 0.0) then
+            surface_type = surface_snow
+        else
+            surface_type = surface_land
+        end if
+
+        if (surface_type == surface_ocean) then
+            ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
+            call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
+            ! --- define relative height
+            h0_m = h / z0_m
+        endif
+         if (surface_type == surface_snow) then
+            ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
+            call get_snow_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_snow)
+            ! --- define relative height
+            h0_m = h / z0_m
+        endif
+        if (surface_type == surface_land) then
+            ! --- define relative height
+            h0_m = h / z0_m
+            ! --- define dynamic velocity in neutral conditions
+            u_dyn0 = U * kappa / log(h0_m)
+        end if
+
+        ! --- define thermal roughness & B = log(z0_m / z0_h)
+        Re = u_dyn0 * z0_m / nu_air
+        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type,  u_dyn0)
+
+        ! --- define relative height [thermal]
+        h0_t = h / z0_t
+
+        ! --- define Ri-bulk
+        Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
+        
+		
+		! --- get the fluxes
+        ! ----------------------------------------------------------------------------
+        call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+	        Lsnow, S_mean, h_salt, &
+            U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+        ! ----------------------------------------------------------------------------
+
+        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
+
+        ! --- define eddy viscosity & inverse Prandtl number
+        Km = kappa * Cm * U * h / phi_m
+        Pr_t_inv = phi_m / phi_h
+	    ! --- define heat flux and  momentum flux
+		hfx=-rho_air*U*dT*Cm*Ct
+		mfx=-rho_air*Cm*Cm*U*U
+        h_salt=h_s
+
+		! --- setting output
+        sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+                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)
+		 ! --- setting additional output
+		sfx2 = sfxDataAddType(phi_m = phi_m, phi_h = phi_h, &
+                hfx = hfx, mfx = mfx, Udyn = Udyn, S_mean = S_mean, &
+                Lsnow = Lsnow, &
+                z0_s = z0_m, h_salt = h_s)
+				
+    end subroutine get_surface_fluxes
+    ! --------------------------------------------------------------------------------
+
+    !< @brief get dynamic scales
+    ! --------------------------------------------------------------------------------
+    subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+	        Lsnow, S_mean, h_salt, &
+            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(out) :: Lsnow, S_mean             
+        real, intent(out) :: h_salt        
+		
+        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
+        real :: w_snow, sigma_m
+        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
+       
+            if (Udyn>u_thsnow) then
+           
+                call get_S_mean(S_mean, S_salt, h_s, z)
+                call get_sigma_m(sigma_m, rho_air, rho_s)
+                call get_w_snow(w_snow, sigma_m, g, d_s, nu_air)
+                Linv=Linv*((1-S_mean)/(1+sigma_m*S_mean))+(g*w_snow*sigma_m*S_mean/(Udyn**3.0))/(1+sigma_m*S_mean)
+                zeta = z * Linv
+                Lsnow=1/Linv
+				  !write(*,*) S_mean, sigma_m, w_snow
+			      !pause 
+			      !stop
+				
+				endif
+       
+       
+       
+        end do
+
+    end subroutine get_dynamic_scales
+    ! --------------------------------------------------------------------------------
+
+    ! 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 + beta_m * zeta
+            phi_h = 1.0 + beta_h * zeta
+        else
+            phi_m = (1.0 - alpha_m * zeta)**(-0.25)
+            phi_h = (1.0 - alpha_h * zeta)**(-0.5)
+        end if
+
+    end subroutine
+    ! --------------------------------------------------------------------------------
+    subroutine get_S_mean(S_mean, S_salt, h_s, z)
+        !< @brief function for snow consentration
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: S_mean   !< snow consentration
+        real, intent(in) ::  S_salt, h_s, z        
+        ! ----------------------------------------------------------------------------
+        S_mean  = S_salt *  h_s/z
+    end subroutine
+    ! --------------
+    subroutine get_sigma_m(sigma_m, rho_air, rho_s)
+        !< @brief function for 
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: sigma_m   !< s
+        real, intent(in) ::  rho_air, rho_s            
+        ! ----------------------------------------------------------------------------
+        sigma_m  = (rho_s - rho_air)/rho_air
+    end subroutine
+    subroutine get_w_snow(w_snow, sigma_m, g, d_s, nu_air)
+        !< @brief function for smow velosity
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: w_snow   !< 
+        real, intent(in) ::  sigma_m, g, d_s, nu_air            
+        ! ----------------------------------------------------------------------------
+        w_snow  = sigma_m * g * d_s * d_s / (18.0 * nu_air);
+
+    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
+        ! ----------------------------------------------------------------------------
+
+
+        if (zeta >= 0.0) then
+            psi_m = -beta_m * zeta;
+            psi_h = -beta_h * zeta;
+        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
+        ! ----------------------------------------------------------------------------
+
+
+        if (zeta_m >= 0.0) then
+            psi_m = -beta_m * zeta_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
+            psi_h = -beta_h * zeta_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
+
+    end subroutine
+    ! --------------------------------------------------------------------------------
+
+end module sfx_most_snow
\ No newline at end of file

srcF/sfx_most_snow_param.f90

+module sfx_most_snow_param
+    !< @brief MOST BD71 surface flux model parameters
+    !< @details  all in SI units
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+    use sfx_phys_const
+    ! --------------------------------------------------------------------------------
+
+    ! directives list
+    ! --------------------------------------------------------------------------------
+    implicit none
+    ! --------------------------------------------------------------------------------
+
+
+    !< von Karman constant [n/d]
+    real, parameter :: kappa = 0.40
+    !< inverse Prandtl (turbulent) number in neutral conditions [n/d]
+    real, parameter :: Pr_t_0_inv = 1.15
+
+
+    !< stability function coeff. (unstable)
+    real, parameter :: alpha_m = 16.0
+    real, parameter :: alpha_h = 16.0
+
+    !< stability function coeff. (stable)
+    real, parameter :: beta_m = 4.7
+    real, parameter :: beta_h = beta_m
+	
+	!< snow parameters
+    real, parameter :: rho_s =900
+    real, parameter :: d_s=0.0000886
+    real, parameter :: h_s=0.07
+    real, parameter :: u_thsnow=0.25
+    real, parameter :: S_salt=0.0004   
+    real, parameter :: rho_air=1.2
+
+end module sfx_most_snow_param