Newer
Older
MODULE drag3
USE param
type, public :: meteoDataType
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)
end type
type, public :: meteoDataVecType
real, allocatable :: h(:) ! constant flux layer height [m]
real, allocatable :: U(:) ! abs(wind speed) at 'h' [m/s]
real, allocatable :: dT(:) ! difference between potential temperature at 'h' and at surface [K]
real, allocatable :: Tsemi(:) ! semi-sum of potential temperature at 'h' and at surface [K]
real, allocatable :: dQ(:) ! difference between humidity at 'h' and at surface [g/g]
real, allocatable :: z0_m(:) ! surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
type, public :: sfxDataType
real :: zeta ! = z/L [n/d]
real :: Rib ! bulk Richardson number [n/d]
real :: Re ! Reynolds number [n/d]
real :: lnzuzt
real :: z0_m ! aerodynamic roughness [m]
real :: z0_t ! thermal roughness [m]
real :: Rib_conv_lim ! Ri-bulk convection critical value [n/d]
real :: cm
real :: ch
real :: ct
real :: ckt
type, public :: numericsType
integer :: maxiters_convection = 10 ! maximum (actual) number of iterations in convection
integer :: maxiters_charnock = 10 ! maximum (actual) number of iterations in charnock roughness
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
real, dimension (numst) :: masout_zl
real, dimension (numst) :: masout_ri
real, dimension (numst) :: masout_re
real, dimension (numst) :: masout_lnzuzt
real, dimension (numst) :: masout_zu
real, dimension (numst) :: masout_ztout
real, dimension (numst) :: masout_rith
real, dimension (numst) :: masout_cm
real, dimension (numst) :: masout_ch
real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
type (data_lutyp) lu1
integer i,numst
do i = 1, numst
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 surf_flux(meteo_cell, sfx_cell, numerics, lu1)
masout_zl(i)=sfx_cell%zeta
masout_ri(i)=sfx_cell%Rib
masout_re(i)=sfx_cell%Re
masout_lnzuzt(i)=sfx_cell%lnzuzt
masout_zu(i)=sfx_cell%z0_m
masout_ztout(i)=sfx_cell%z0_t
masout_rith(i)=sfx_cell%Rib_conv_lim
masout_cm(i)=sfx_cell%cm
masout_ch(i)=sfx_cell%ch
masout_ct(i)=sfx_cell%ct
masout_ckt(i)=sfx_cell%ckt
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
END SUBROUTINE surf_flux_vec
function f_m_conv(zeta)
real f_m_conv
real zeta
f_m_conv = (1.0 - alpha_m * zeta)**0.25
end function f_m_conv
function f_h_conv(zeta)
real f_h_conv
real zeta
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
! stability functions (neutral)
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
! ----------------------------------------------------------------------------
zeta = 0.0
psi_m = log(h0_m)
psi_h = log(h0_t) / Pr_t_0_inv
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
! --------------------------------------------------------------------------------
! stability functions (stable)
! --------------------------------------------------------------------------------
subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: Rib_coeff
real :: psi0_m, psi0_h
real :: phi
real :: c
! ----------------------------------------------------------------------------
psi0_m = alog(h0_m)
psi0_h = B / psi0_m
Rib_coeff = beta_m * Rib
c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
(2.0 * beta_m * (1.0 - Rib_coeff))
phi = beta_m * zeta
psi_m = psi0_m + phi
psi_h = (psi0_m + B) / Pr_t_0_inv + phi
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection-semistrong)
! --------------------------------------------------------------------------------
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, x1, y0, y1
integer :: i
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t)
if (abs(B) < 1.0e-10) psi_h = psi_m
zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
y1 = (1.0 - alpha_m * zeta)**0.25e0
x1 = sqrt(1.0 - alpha_h_fix * zeta)
y0 = (1.0 - alpha_m * zeta0_m)**0.25e0
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = max(y0, 1.000001e0)
x0 = max(x0, 1.000001e0)
psi_m = log((y1 - 1.0e0)*(y0 + 1.0e0)/((y1 + 1.0e0)*(y0 - 1.0e0))) + 2.0e0*(atan(y1) - atan(y0))
psi_h = log((x1 - 1.0e0)*(x0 + 1.0e0)/((x1 + 1.0e0)*(x0 - 1.0e0))) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection)
! --------------------------------------------------------------------------------
subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, x10, y10, &
h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: zeta_conv_lim
real, intent(in) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, y0
real :: p1, p0
real :: a1, b1, c, f
integer :: i
! ----------------------------------------------------------------------------
p1 = 2.0 * atan(y10) + log((y10 - 1.0) / (y10 + 1.0))
p0 = log((x10 - 1.0) / (x10 + 1.0))
zeta = zeta_conv_lim
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
a1 = (zeta_conv_lim / zeta)**(1.0 / 3.0)
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = (1.0 - alpha_m * zeta0_m)**0.25
c = log((x0 + 1.0)/(x0 - 1.0))
b1 = -2.0*atan(y0) + log((y0 + 1.0)/(y0 - 1.0))
f = 3.0 * (1.0 - a1)
psi_m = f / y10 + p1 + b1
psi_h = (f / x10 + p0 + c) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! define convection limit
! --------------------------------------------------------------------------------
subroutine get_convection_lim(zeta_lim, Rib_lim, x10, y10, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim, Rib_lim
real, intent(out) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: an4, an5
real :: psi_m, psi_h
! ----------------------------------------------------------------------------
an5=(Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (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)
y10 = f_m_conv(zeta_lim)
x10 = f_h_conv(zeta_lim)
! ......definition of r-prim......
an4 = zeta_lim / h0_m
an5 = zeta_lim / h0_t
if (abs(B) < 1.0e-10) an5=an4
an5 = sqrt(1.0 - alpha_h_fix * an5)
an4 = (1.0 - alpha_m * an4)**0.25
psi_m = 2.0 * (atan(y10) - atan(an4)) + alog((y10 - 1.0) * (an4 + 1.0)/((y10 + 1.0) * (an4 - 1.0)))
psi_h = alog((x10 - 1.0) * (an5 + 1.0)/((x10 + 1.0) * (an5 - 1.0))) / Pr_t_0_inv
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! --------------------------------------------------------------------------------
! charnock roughness
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_m, u_dyn0
real, intent(in) :: U, h
integer, intent(in) :: maxiters
real :: U10m
real :: a1, c1, y1, c1min, f
integer :: i, j
! ----------------------------------------------------------------------------
U10m = U
a1 = 0.0e0
y1 = 25.0e0
c1min = log(h_charnock) / kappa
do i = 1, maxiters
f = c1_charnock - 2.0 * log(U10m)
do j = 1, maxiters
c1 = (f + 2.0 * log(y1)) / kappa
! looks like the check should use U10 instead of U
! but note that a1 should be set = 0 in cycle beforehand
if(U <= 8.0e0) a1 = log(1.0 + c2_charnock * ((y1 / U10m)**3)) / kappa
c1 = max(c1 - a1, c1min)
y1 = c1
end do
z0_m = h_charnock * exp(-c1 * kappa)
z0_m = max(z0_m,0.000015e0)
U10m = U*alog(h_charnock/z0_m)/(alog(h/z0_m))
end do
! --- define dynamic velocity in neutral conditions
u_dyn0 = U10m / c1
end subroutine
! --------------------------------------------------------------------------------
SUBROUTINE surf_flux(meteo, sfx, numerics, lu)
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
type (data_lutyp) lu
real h ! surface flux layer height [m]
real z0_m, z0_t ! aerodynamic & 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 Re ! roughness Reynolds number = u_dyn0 * z0 / nu [n/d]
real u_dyn0 ! dynamic velocity in neutral conditions [m/s]
real zeta ! = z/L [n/d]
real zeta_conv_lim ! z/L critical value for matching free convection limit [n/d]
real Rib ! bulk Richardson number
real Rib_conv_lim ! Ri-bulk critical value for matching free convection limit [n/d]
real psi_m, psi_h ! universal functions [momentum] & [thermal]
real U ! wind velocity at h [m/s]
real Tsemi !
integer lu_indx
! output ... !
real an4, o, am
real c4, c0, an0
! ...
! local unnamed !
real x10, y10
real al
real q4, t4
! .....
! just local variables
real f1, a1
! ....
integer is_ocean
Tsemi = meteo%Tsemi
U = meteo%U
t4 = meteo%dT
q4 = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
if(z0_m < 0.0) then
!> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
is_ocean = 1
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
else
! ......parameters from viscosity sublayer......
is_ocean = 0
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
Re = u_dyn0 * z0_m / nu_air
if(Re <= Re_rough_min) then
B = B1_rough * alog(B3_rough * Re) + B2_rough
else
! B4 takes into account Re value at z' ~ O(10) z0
B = B4_rough * (Re**B2_rough)
! --- apply max restriction based on surface type
if (is_ocean == 1) then
B = min(B, B_max_ocean)
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
! --- define relative height [thermal]
h0_t = h / z0_t
! ......humidity stratification and ri-number......
al = g / Tsemi
Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
h0_m, h0_t, B)
if (Rib > 0.0e0) then
! ......stable stratification......
!write (*,*) 'stable'
Rib = min(Rib, Rib_max)
call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
f1 = beta_m * zeta
am = 1.0 + f1
o = 1.0/Pr_t_0_inv + f1
else if (Rib < Rib_conv_lim) then
! ......strong instability.....
!write (*,*) 'instability'
call get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)
a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
am = a1 / y10
o = a1 / (Pr_t_0_inv * x10)
else if (Rib > -0.001e0) then
! ......nearly neutral......
write (*,*) 'neutral'
call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
am = 1.0e0
o = 1.0e0 / Pr_t_0_inv
else
! ......week and semistrong instability......
!write (*,*) 'semistrong'
call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
am = (1.0 - alpha_m * zeta)**(-0.25)
o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
end if
! ......computation of cu,co,k(h),alft
c4=kappa/psi_m
c0=kappa/psi_h
an4=kappa*c4*U*h/am
! ......exit......
sfx%zeta = zeta
sfx%Rib = Rib
sfx%Re = Re
sfx%lnzuzt=B
sfx%z0_m = z0_m
sfx%z0_t = z0_t
sfx%Rib_conv_lim = Rib_conv_lim
sfx%cm=c4
sfx%ch=c0
sfx%ct=an4
sfx%ckt=an0
END SUBROUTINE surf_flux
END MODULE drag3