Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
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
!!! TO DO: consistent psi_m and psi_h needed!!! For now is those of NCAR !!!
!!
MODULE sbcblk_algo_andreas
!!======================================================================
!! *** MODULE sbcblk_algo_andreas ***
!! Computes:
!! * bulk transfer coefficients C_D, C_E and C_H
!! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
!! * the effective bulk wind speed at 10m Ubzu
!! according to Andreas et al. (2015)
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!! Andreas, E.L., Mahrt, L. and Vickers, D. (2015),
!! An improved bulk air–sea surface flux algorithm,
!! including spray‐mediated transfer.
!! Q.J.R. Meteorol. Soc., 141: 642-654. doi:10.1002/qj.2424
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!
!! * bulk transfer coefficients C_D, C_E and C_H
!! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
!! * the effective bulk wind speed at z=zu: Ubzu
!! => all these are used in bulk formulas in sbcblk.F90
!!
!! Using the bulk formulation/param. of Large & Yeager 2008
!!
!! Routine turb_andreas maintained and developed in AeroBulk
!! (https://github.com/brodeau/aerobulk/)
!!
!! ** Author: L. Brodeau, August 2020 / AeroBulk (https://github.com/brodeau/aerobulk)
!!----------------------------------------------------------------------
!! History : 4.x ! 2020-08 (L.Brodeau) Original code
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
IMPLICIT NONE
PRIVATE
!! Important (Brodeau fix):
REAL(wp), PARAMETER :: rRi_max = 0.15_wp ! Bulk Ri above which the algorithm fucks up!
! ! (increasing (>0) Ri means that surface layer increasingly stable and/or wind increasingly weak)
REAL(wp), PARAMETER :: rCs_min = 0.35E-3_wp ! minimum value to tolarate for CE and CH ! Must be larger than "Cx_min" !!!
PUBLIC :: TURB_ANDREAS, psi_m_andreas, psi_h_andreas
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE turb_andreas( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
& Cd, Ch, Ce, t_zu, q_zu, Ubzu, &
& nb_iter, CdN, ChN, CeN )
!!----------------------------------------------------------------------------------
!! *** ROUTINE turb_andreas ***
!!
!! ** Purpose : Computes turbulent transfert coefficients of surface
!! fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
!! Returns the effective bulk wind speed at zu to be used in the bulk formulas
!!
!! INPUT :
!! -------
!! * zt : height for temperature and spec. hum. of air [m]
!! * zu : height for wind speed (usually 10m) [m]
!! * sst : bulk SST [K]
!! * t_zt : potential air temperature at zt [K]
!! * ssq : specific humidity at saturation at SST [kg/kg]
!! * q_zt : specific humidity of air at zt [kg/kg]
!! * U_zu : scalar wind speed at zu [m/s]
!!
!! OUTPUT :
!! --------
!! * Cd : drag coefficient
!! * Ch : sensible heat coefficient
!! * Ce : evaporation coefficient
!! * t_zu : pot. air temperature adjusted at wind height zu [K]
!! * q_zu : specific humidity of air // [kg/kg]
!! * Ubzu : bulk wind speed at zu [m/s]
!!
!!
!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m]
REAL(wp), INTENT(in ) :: zu ! height for U_zu [m]
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin]
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin]
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg]
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg]
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s]
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K]
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg]
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s]
!
INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations
REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN
REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN
REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN
!
INTEGER :: nbit, jit ! iterations...
LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U
!!
REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
REAL(wp), DIMENSION(jpi,jpj) :: z0 ! roughness length (momentum) [m]
REAL(wp), DIMENSION(jpi,jpj) :: UN10 ! Neutral wind speed at zu [m/s]
REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
REAL(wp), DIMENSION(jpi,jpj) :: RiB ! square root of Cd
!!
!!----------------------------------------------------------------------------------
nbit = nb_iter0
IF( PRESENT(nb_iter) ) nbit = nb_iter
l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
Ubzu = MAX( 0.25_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
!! First guess:
UN10 = Ubzu
Cd = 1.1E-3_wp
Ch = 1.1E-3_wp
Ce = 1.1E-3_wp
t_zu = t_zt
q_zu = q_zt
!! First guess of turbulent scales for scalars:
ztmp0 = SQRT(Cd)
t_star = Ch/ztmp0*(t_zu - sst) ! theta*
q_star = Ce/ztmp0*(q_zu - ssq) ! q*
! Bulk Richardson number:
RiB(:,:) = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu )
!! ITERATION BLOCK
DO jit = 1, nbit
WHERE ( RiB < rRi_max )
!! Normal condition case:
u_star = U_STAR_ANDREAS( UN10 )
ELSEWHERE
!! Extremely stable + weak wind !!!
!! => for we force u* to be consistent with minimum value for CD:
!! (otherwize algorithm becomes nonsense...)
u_star = SQRT(Cx_min) * Ubzu ! Cd does not go below Cx_min !
ENDWHERE
!! Stability parameter :
zeta_u = zu*One_on_L( t_zu, q_zu, u_star, t_star, q_star ) ! zu * 1/L
!! Drag coefficient:
ztmp0 = u_star/Ubzu
Cd = MAX( ztmp0*ztmp0 , Cx_min )
!! Roughness length:
z0 = MIN( z0_from_Cd( zu, Cd, ppsi=psi_m_andreas(zeta_u) ) , z0_sea_max )
!! z0t and z0q, based on LKB, just like into COARE 2.5:
ztmp0 = z0 * u_star / visc_air(t_zu) ! Re_r
ztmp1 = z0tq_LKB( 1, ztmp0, z0 ) ! z0t
ztmp2 = z0tq_LKB( 2, ztmp0, z0 ) ! z0q
!! Turbulent scales at zu :
ztmp0 = psi_h_andreas(zeta_u) ! lolo: zeta_u for scalars???
t_star = (t_zu - sst)*vkarmn/(LOG(zu) - LOG(ztmp1) - ztmp0) ! theta* (ztmp1 == z0t in rhs term)
q_star = (q_zu - ssq)*vkarmn/(LOG(zu) - LOG(ztmp2) - ztmp0) ! q* (ztmp2 == z0q in rhs term)
IF( (.NOT. l_zt_equal_zu).AND.( jit > 1 ) ) THEN
!! Re-updating temperature and humidity at zu if zt /= zu:
ztmp0 = zeta_u/zu*zt ! zeta_t
ztmp0 = LOG(zt/zu) + psi_h_andreas(zeta_u) - psi_h_andreas(ztmp0)
t_zu = t_zt - t_star/vkarmn*ztmp0
q_zu = q_zt - q_star/vkarmn*ztmp0
RiB = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu ) !LOLO
ENDIF
!! Update neutral-stability wind at zu:
UN10 = MAX( 0.1_wp , UN10_from_ustar( zu, Ubzu, u_star, psi_m_andreas(zeta_u) ) ) ! UN10
END DO !DO jit = 1, nbit
! Compute transfer coefficients at zu:
ztmp0 = u_star/Ubzu
Cd = MAX( ztmp0*ztmp0 , Cx_min ) ! the earlier use of Cx_min on u* should make use of Cx_min here unnecessary!
ztmp1 = t_zu - sst ; ztmp1 = SIGN( MAX(ABS(ztmp1),1.E-6_wp), ztmp1 ) ! dt_zu
ztmp2 = q_zu - ssq ; ztmp2 = SIGN( MAX(ABS(ztmp2),1.E-9_wp), ztmp2 ) ! dq_zu
Ch = MAX( ztmp0*t_star/ztmp1 , rCs_min )
Ce = MAX( ztmp0*q_star/ztmp2 , rCs_min )
!! Neutral-stability coefficients:
ztmp0 = 1._wp/LOG(zu/z0)
ztmp1 = z0 * u_star / visc_air(t_zu) ! Re_r
IF(PRESENT(CdN)) CdN = vkarmn2*ztmp0*ztmp0
IF(PRESENT(ChN)) ChN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 1, ztmp1, z0 ))
IF(PRESENT(CeN)) CeN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 2, ztmp1, z0 ))
END SUBROUTINE turb_andreas
FUNCTION U_STAR_ANDREAS( pun10 )
!!----------------------------------------------------------------------------------
!! Estimate of the friction velocity as a function of the neutral-stability wind
!! speed at at 10m
!!
!! Origin: Eq.(2.2) of Andreas et al. (2015)
!!
!! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pun10 !: neutral-stability scalar wind speed at 10m (m/s)
REAL(wp), DIMENSION(jpi,jpj) :: u_star_andreas !: friction velocity [m/s]
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: za, zt, zw ! local scalars
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zw = pun10(ji,jj)
za = zw - 8.271_wp
zt = za + SQRT( 0.12_wp*za*za + 0.181_wp )
u_star_andreas(ji,jj) = 0.239_wp + 0.0433_wp * zt
END_2D
END FUNCTION U_STAR_ANDREAS
FUNCTION psi_m_andreas( pzeta )
!!----------------------------------------------------------------------------------
!! Universal profile stability function for momentum
!! TO DO !!!!!!!!!!!!!!!!!!!!!
!! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE
!!
!! pzeta : stability paramenter, z/L where z is altitude measurement
!! and L is M-O length
!!
!! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: psi_m_andreas
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
!
REAL(wp), PARAMETER :: zam = 5._wp ! a_m (just below Eq.(9b)
REAL(wp), PARAMETER :: zbm = zam/6.5_wp ! b_m (just below Eq.(9b)
!
REAL(wp), PARAMETER :: z1o3 = 1._wp/3._wp
REAL(wp), PARAMETER :: zsr3 = SQRT(3._wp)
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zta, zx2, zx, zpsi_unst, zbbm, zpsi_stab, zstab ! local scalars
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and big!)
!
!! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!!
zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5
zx2 = MAX( zx2 , 1._wp )
zx = SQRT(zx2) ! (1 - 16z)^0.25
zpsi_unst = 2._wp*LOG(ABS( (1._wp + zx )*0.5_wp )) &
& + LOG(ABS( (1._wp + zx2)*0.5_wp )) &
& - 2._wp*ATAN(zx) + rpi*0.5_wp
!
!! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(12) Grachev et al 2007]:
zx = ABS(1._wp + zta)**z1o3
zbbm = ABS( (1._wp - zbm)/zbm )**z1o3 ! B_m
!
zpsi_stab = -3.*zam/zbm*(zx - 1._wp) + zam*zbbm/(2.*zbm) * ( &
& 2.*LOG(ABS( ( zx + zbbm )/(1._wp + zbbm ) )) &
& - LOG(ABS( (zx*zx - zx*zbbm + zbbm*zbbm)/(1._wp - zbbm + zbbm*zbbm) )) &
& + 2.*zsr3*( ATAN( (2.*zx - zbbm)/(zsr3*zbbm) ) - ATAN( (2._wp - zbbm)/(zsr3*zbbm) ) ) )
!
!
zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
!
psi_m_andreas(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
!
END_2D
END FUNCTION psi_m_andreas
FUNCTION psi_h_andreas( pzeta )
!!----------------------------------------------------------------------------------
!! Universal profile stability function for temperature and humidity
!!
!! TO DO
!! !! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE
!!
!! pzeta : stability paramenter, z/L where z is altitude measurement
!! and L is M-O length
!!
!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: psi_h_andreas
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
!
REAL(wp), PARAMETER :: zah = 5._wp ! a_h (just below Eq.(9b)
REAL(wp), PARAMETER :: zbh = 5._wp ! b_h (just below Eq.(9b)
REAL(wp), PARAMETER :: zch = 3._wp ! c_h (just below Eq.(9b)
REAL(wp), PARAMETER :: zbbh = SQRT(5._wp) ! B_h (just below Eq.(13)
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zta, zz, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and large!)
!
!! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!!
zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5
zx2 = MAX( zx2 , 1._wp )
zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
!
!! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(13) Grachev et al 2007]:
zz = 2.*zta + zch
zpsi_stab = - 0.5*zbh*LOG(ABS(1._wp + zch*zta + zta*zta)) &
& + (-zah/zbbh + 0.5*zbh*zch/zbbh) &
& *( LOG(ABS((zz - zbbh)/(zz + zbbh))) &
& - LOG(ABS((zch - zbbh)/(zch + zbbh))) )
!
zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
!
psi_h_andreas(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
!
END_2D
END FUNCTION psi_h_andreas
!!======================================================================
END MODULE sbcblk_algo_andreas