From d2bc97ad8fe3b31b9ee626403c8249af851558b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=D0=92=D0=B8=D0=BA=D1=82=D0=BE=D1=80=D0=B8=D1=8F=20=D0=A1?=
 =?UTF-8?q?=D1=83=D1=8F=D0=B7=D0=BE=D0=B2=D0=B0?=
 <viktoriasuazova@MacBook-Pro-Viktoria.local>
Date: Fri, 13 Dec 2024 16:50:06 +0300
Subject: [PATCH] added coast1,coast2

---
 srcF/sfx_config.f90          |   6 +-
 srcF/sfx_surface.f90         |  37 ++++++++--
 srcF/sfx_z0m_all_surface.f90 | 138 ++++++++++++++++++++++++++++++++++-
 3 files changed, 170 insertions(+), 11 deletions(-)

diff --git a/srcF/sfx_config.f90 b/srcF/sfx_config.f90
index b2a2620..68df2c5 100644
--- a/srcF/sfx_config.f90
+++ b/srcF/sfx_config.f90
@@ -167,8 +167,8 @@ contains
         dataset%filename = get_dataset_filename(id)
         dataset%nmax = 0
 
-        dataset%surface = surface_ice
-        dataset%surface_type = surface_ice
+        dataset%surface = surface_lake
+        dataset%surface_type = surface_lake
         dataset%z0_h = -1.0
         dataset%lai = 1.0
         dataset%depth = 10.0
@@ -194,7 +194,7 @@ contains
         else if (id == dataset_toga) then
             ! *: check & fix values
             dataset%h = 15.0
-            dataset%surface = surface_ocean
+            dataset%surface = surface_lake
             dataset%z0_m = -1.0
         end if        
        write (*,*) dataset%surface, 'dataset%surface'
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index 8f4ebe6..ed09eb1 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -53,7 +53,9 @@ module sfx_surface
     integer, public, parameter :: z0m_map_id = 3      
     integer, public, parameter :: z0m_user = 4 
     integer, public, parameter :: z0m_and = 5 
-
+    integer, public, parameter :: z0m_cst1 = 6
+    integer, public, parameter :: z0m_cst2 = 7
+    integer, public, parameter :: z0m_cst3 = 8 
 
     integer, public, parameter :: z0t_kl_water = 0 
     integer, public, parameter :: z0t_kl_land = 1   
@@ -75,6 +77,9 @@ module sfx_surface
     character(len = 16), parameter :: z0m_tag_map = 'map'
     character(len = 16), parameter :: z0m_tag_user = 'user'
     character(len = 16), parameter :: z0m_tag_and = 'andreas'
+    character(len = 16), parameter :: z0m_tag_cst1 = 'coast1'
+    character(len = 16), parameter :: z0m_tag_cst2 = 'coast2'
+    character(len = 16), parameter :: z0m_tag_cst3 = 'coast3'
 
 
     character(len = 16), parameter :: z0t_tag_kl_water = 'kl_water'
@@ -93,7 +98,7 @@ module sfx_surface
 
     integer, public, parameter :: ocean_z0m_id = z0m_ch     !< ocean surface
     integer, public, parameter :: land_z0m_id = z0m_map_id        !< land surface
-    integer, public, parameter :: lake_z0m_id = z0m_fe        !< lake surface
+    integer, public, parameter :: lake_z0m_id = z0m_map_id        !< lake surface
     integer, public, parameter :: snow_z0m_id = z0m_map_id        !< snow covered surface    
     integer, public, parameter :: forest_z0m_id = z0m_map_id     !< forest csurface  
     integer, public, parameter :: usersf_z0m_id = z0m_ch       !< user surface  
@@ -101,12 +106,12 @@ module sfx_surface
 
 
     integer, public, parameter :: ocean_z0t_id = z0t_kl_water     !< ocean surface
-    integer, public, parameter :: land_z0t_id = z0t_mix        !< land surface
-    integer, public, parameter :: lake_z0t_id = z0t_cz        !< lake surface
-    integer, public, parameter :: snow_z0t_id = z0t_zi         !< snow covered surface    
-    integer, public, parameter :: forest_z0t_id = z0t_ot      !< forest csurface  
+    integer, public, parameter :: land_z0t_id = z0t_kl_land        !< land surface
+    integer, public, parameter :: lake_z0t_id = z0t_kl_land        !< lake surface
+    integer, public, parameter :: snow_z0t_id = z0t_kl_land         !< snow covered surface    
+    integer, public, parameter :: forest_z0t_id = z0t_kl_land      !< forest csurface  
     integer, public, parameter :: usersf_z0t_id = z0t_kl_water      !< user surface  
-    integer, public, parameter :: ice_z0t_id = z0t_zi      !< user surface  
+    integer, public, parameter :: ice_z0t_id = z0t_kl_land     !< user surface  
     
     ! --------------------------------------------------------------------------------
     real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
@@ -210,6 +215,12 @@ contains
             z0m_id = z0m_map_id
         else if (trim(tag_z0m) == trim(z0m_tag_user)) then
             z0m_id = z0m_user
+        else if (trim(tag_z0m) == trim(z0m_tag_cst1)) then
+            z0m_id = z0m_cst1
+        else if (trim(tag_z0m) == trim(z0m_tag_cst2)) then
+            z0m_id = z0m_cst2
+        else if (trim(tag_z0m) == trim(z0m_tag_cst3)) then
+            z0m_id = z0m_cst3
         end if
 
     end function
@@ -230,6 +241,12 @@ contains
             tag_z0m = z0m_tag_map
         else if (z0m_id == z0m_user) then
             tag_z0m = z0m_tag_user
+        else if (z0m_id == z0m_cst1) then
+            tag_z0m = z0m_tag_cst1
+        else if (z0m_id == z0m_cst2) then
+            tag_z0m = z0m_tag_cst2
+        else if (z0m_id == z0m_cst3) then
+            tag_z0m = z0m_tag_cst3
         end if 
      end function
 
@@ -393,6 +410,12 @@ contains
         call get_dynamic_roughness_and(z0_m, u_dyn0, U, h, z0m_map)
    else if (z0m_id == z0m_user) then
        write(*, *) 'z0m_user'
+    else if (z0m_id == z0m_cst1) then
+        call get_dynamic_roughness_coast1(z0_m, u_dyn0, U, h, maxiters)
+    else if (z0m_id == z0m_cst2) then
+        call get_dynamic_roughness_coast2(z0_m, u_dyn0, U, h, maxiters)
+    else if (z0m_id == z0m_cst3) then
+        call get_dynamic_roughness_coast3(z0_m, u_dyn0, U, h, maxiters)
 end if 
      
 end subroutine
diff --git a/srcF/sfx_z0m_all_surface.f90 b/srcF/sfx_z0m_all_surface.f90
index 4d6f2ed..1f8b69c 100644
--- a/srcF/sfx_z0m_all_surface.f90
+++ b/srcF/sfx_z0m_all_surface.f90
@@ -8,6 +8,9 @@ module sfx_z0m_all_surface
     public :: get_dynamic_roughness_ow
     public :: get_dynamic_roughness_fetch
     public :: get_dynamic_roughness_and
+    public :: get_dynamic_roughness_coast1
+    public :: get_dynamic_roughness_coast2
+    public :: get_dynamic_roughness_coast3
     ! --------------------------------------------------------------------------------
 
 
@@ -252,6 +255,139 @@ subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
             end subroutine
          ! --------------------------------------------------------------------------------
 
-
+            subroutine get_dynamic_roughness_coast1(z0_m, u_dyn0, U, h, maxiters)
+                ! ----------------------------------------------------------------------------
+                real, intent(out) :: z0_m           !< aerodynamic roughness [m]
+                real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
+        
+                real, intent(in) :: h               !< constant flux layer height [m]
+                real, intent(in) :: U               !< abs(wind speed) [m/s]
+                integer, intent(in) :: maxiters     !< maximum number of iterations
+                ! ----------------------------------------------------------------------------
+        ! Taylor&Yelland formulation
+                ! --- local variables
+                real :: Uc, c, c_min, a, b                  ! wind speed at h_charnock [m/s] & Cd^(-1)
+        
+               real :: pi=3.14159265
+               real :: hs, Tp, Lp
+               integer :: i, j
+                ! ----------------------------------------------------------------------------
+        
+                Uc = U
+                a = 0
+                b = 0
+                c_min = log(h_charnock) / kappa
+        
+                do i = 1, maxiters
+                  hs = 0.0248*(Uc**2.) !hs is the significant wave height
+                  Tp = 0.729*max(Uc,0.1) !Tp dominant wave period
+                  Lp = g*Tp**2/(2*pi) !Lp is the wavelength of the dominant wave
+                    do j = 1, maxiters
+                        z0_m = 1200.*hs*(hs/Lp)**4.5
+                        c = log(h_charnock / z0_m) / kappa
+                        if (Uc <= 8.0e0) then
+                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
+                        c = max(c - a, c_min)
+                        b = c
+                        z0_m = h_charnock * exp(-kappa*c)
+                        end if
+                    end do
+                z0_m = max( z0_m, 1.27e-7)  !These max/mins were suggested by
+                z0_m = min( z0_m, 2.85e-3)  !Davis et al. (2008)
+                Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
+                end do
+                u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
+        !       c = Uc/u_dyn0
+        !       write (*,*) 'out1', u_dyn0
+        
+            end subroutine
+        ! Hwang
+            subroutine get_dynamic_roughness_coast2(z0_m, u_dyn0, U, h, maxiters)
+                ! ----------------------------------------------------------------------------
+                real, intent(out) :: z0_m           !< aerodynamic roughness [m]
+                real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
+        
+                real, intent(in) :: h               !< constant flux layer height [m]
+                real, intent(in) :: U               !< abs(wind speed) [m/s]
+                integer, intent(in) :: maxiters     !< maximum number of iterations
+                ! ----------------------------------------------------------------------------
+                ! --- local variables
+                real :: Uc, c, c1, a, b, c_min                  ! wind speed at h_charnock [m/s] & Cd^(-1)
+               integer :: i, j
+                ! ----------------------------------------------------------------------------
+                Uc = U
+                a = 0
+                b = 0
+                c_min = log(h_charnock) / kappa
+        
+                do i = 1, maxiters
+                    c1 = 0.016 * Uc**2
+                    do j = 1, maxiters
+                        c = 1.0 / (0.01*sqrt(8.058 + 0.967 * Uc - c1))
+                        if (Uc <= 8.0e0) then
+                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
+                        c = max(c - a, c_min)
+                        b = c
+                        end if
+                    end do
+                    z0_m = h_charnock * exp(-kappa*c)
+                    z0_m = max(z0_m, 0.000015e0)
+                    Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
+                end do
+                ! --- define dynamic velocity in neutral conditions
+                u_dyn0 = Uc / c
+        !       write (*,*) 'out1', u_dyn0
+            ! --------------------------------------------------------------------------------
+            end subroutine
+        ! Zhao et al 2015 10^3*z0=15.6*u_dyn0**2/g + 10**(-2)
+            subroutine get_dynamic_roughness_coast3(z0_m, u_dyn0, U, h, maxiters)
+                   ! ----------------------------------------------------------------------------
+                real, intent(out) :: z0_m           !< aerodynamic roughness [m]
+                real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
+        
+                real, intent(in) :: h               !< constant flux layer height [m]
+                real, intent(in) :: U               !< abs(wind speed) [m/s]
+                integer, intent(in) :: maxiters     !< maximum number of iterations
+                ! ----------------------------------------------------------------------------
+                ! --- local variables
+                real :: Uc, u_dyn, z0_m0                  ! wind speed at h_charnock [m/s]
+        
+                real :: c, c_min, a, b
+                real :: f1
+        
+                integer :: i, j
+                ! ----------------------------------------------------------------------------
+                Uc = U
+                u_dyn = U / 28.0
+                a = 0
+                b = 0
+                z0_m0=0.082
+                c_min = log(h_charnock) / kappa
+                    if (Uc >8.00 .AND. Uc <= 35.0e0) then
+                    do i = 1, maxiters
+                        f1 = 15.6*u_dyn**2/g
+                        z0_m0 = 0.001 * (f1 + 0.01)
+                        c = log(h_charnock / z0_m0) / kappa
+                    end do
+                    end if
+                    if (Uc <= 8.0e0) then
+                    do j = 1, maxiters
+                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
+                        c = max(c - a, c_min)
+                        b = c
+                        z0_m0 = h_charnock * exp(-c * kappa)
+                    end do
+                    else
+                        z0_m0 = 2.82
+                    end if
+                    z0_m = max(z0_m0, 0.000015e0)
+                    Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
+        
+                ! --- define dynamic velocity in neutral conditions
+        !        u_dyn0 = Uc / c
+                u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
+        
+            end subroutine
+        
    
     end module  sfx_z0m_all_surface
\ No newline at end of file
-- 
GitLab