1 | MODULE sbcblk_algo_ice_an05 |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE sbcblk_algo_ice_an05 *** |
---|
4 | !! Computes turbulent components of surface fluxes over sea-ice |
---|
5 | !! |
---|
6 | !! Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results. |
---|
7 | !! Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7 |
---|
8 | !! |
---|
9 | !! * bulk transfer coefficients C_D, C_E and C_H |
---|
10 | !! * air temp. and spec. hum. adjusted from zt (usually 2m) to zu (usually 10m) if needed |
---|
11 | !! * the "effective" bulk wind speed at zu: Ub (including gustiness contribution in unstable conditions) |
---|
12 | !! => all these are used in bulk formulas in sbcblk.F90 |
---|
13 | !! |
---|
14 | !! Routine turb_ice_an05 maintained and developed in AeroBulk |
---|
15 | !! (https://github.com/brodeau/aerobulk/) |
---|
16 | !! |
---|
17 | !! Author: Laurent Brodeau, Summer 2020 |
---|
18 | !! |
---|
19 | !!---------------------------------------------------------------------- |
---|
20 | USE par_kind, ONLY: wp |
---|
21 | USE par_oce, ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls, ntsi, ntsj, ntei, ntej |
---|
22 | USE lib_mpp, ONLY: ctl_stop ! distribued memory computing library |
---|
23 | USE phycst ! physical constants |
---|
24 | USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer |
---|
25 | |
---|
26 | IMPLICIT NONE |
---|
27 | PRIVATE |
---|
28 | |
---|
29 | PUBLIC :: turb_ice_an05 |
---|
30 | |
---|
31 | INTEGER , PARAMETER :: nbit = 8 ! number of itterations |
---|
32 | |
---|
33 | !! * Substitutions |
---|
34 | # include "do_loop_substitute.h90" |
---|
35 | !!---------------------------------------------------------------------- |
---|
36 | CONTAINS |
---|
37 | |
---|
38 | SUBROUTINE turb_ice_an05( zt, zu, Ts_i, t_zt, qs_i, q_zt, U_zu, & |
---|
39 | & Cd_i, Ch_i, Ce_i, t_zu_i, q_zu_i, & |
---|
40 | & CdN, ChN, CeN, xz0, xu_star, xL, xUN10 ) |
---|
41 | !!---------------------------------------------------------------------- |
---|
42 | !! *** ROUTINE turb_ice_an05 *** |
---|
43 | !! |
---|
44 | !! ** Purpose : Computes turbulent transfert coefficients of surface |
---|
45 | !! fluxes according to: |
---|
46 | !! Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results. |
---|
47 | !! Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7 |
---|
48 | !! |
---|
49 | !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu |
---|
50 | !! Returns the effective bulk wind speed at zu to be used in the bulk formulas |
---|
51 | !! |
---|
52 | !! INPUT : |
---|
53 | !! ------- |
---|
54 | !! * zt : height for temperature and spec. hum. of air [m] |
---|
55 | !! * zu : height for wind speed (usually 10m) [m] |
---|
56 | !! * Ts_i : surface temperature of sea-ice [K] |
---|
57 | !! * t_zt : potential air temperature at zt [K] |
---|
58 | !! * qs_i : saturation specific humidity at temp. Ts_i over ice [kg/kg] |
---|
59 | !! * q_zt : specific humidity of air at zt [kg/kg] |
---|
60 | !! * U_zu : scalar wind speed at zu [m/s] |
---|
61 | !! |
---|
62 | !! OUTPUT : |
---|
63 | !! -------- |
---|
64 | !! * Cd_i : drag coefficient over sea-ice |
---|
65 | !! * Ch_i : sensible heat coefficient over sea-ice |
---|
66 | !! * Ce_i : sublimation coefficient over sea-ice |
---|
67 | !! * t_zu_i : pot. air temp. adjusted at zu over sea-ice [K] |
---|
68 | !! * q_zu_i : spec. hum. of air adjusted at zu over sea-ice [kg/kg] |
---|
69 | !! |
---|
70 | !! OPTIONAL OUTPUT: |
---|
71 | !! ---------------- |
---|
72 | !! * CdN : neutral-stability drag coefficient |
---|
73 | !! * ChN : neutral-stability sensible heat coefficient |
---|
74 | !! * CeN : neutral-stability evaporation coefficient |
---|
75 | !! * xz0 : return the aerodynamic roughness length (integration constant for wind stress) [m] |
---|
76 | !! * xu_star : return u* the friction velocity [m/s] |
---|
77 | !! * xL : return the Obukhov length [m] |
---|
78 | !! * xUN10 : neutral wind speed at 10m [m/s] |
---|
79 | !! |
---|
80 | !! ** Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
81 | !!---------------------------------------------------------------------------------- |
---|
82 | REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] |
---|
83 | REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] |
---|
84 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: Ts_i ! ice surface temperature [Kelvin] |
---|
85 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] |
---|
86 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: qs_i ! sat. spec. hum. at ice/air interface [kg/kg] |
---|
87 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! spec. air humidity at zt [kg/kg] |
---|
88 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] |
---|
89 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Cd_i ! drag coefficient over sea-ice |
---|
90 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ch_i ! transfert coefficient for heat over ice |
---|
91 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ce_i ! transfert coefficient for sublimation over ice |
---|
92 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: t_zu_i ! pot. air temp. adjusted at zu [K] |
---|
93 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: q_zu_i ! spec. humidity adjusted at zu [kg/kg] |
---|
94 | !!---------------------------------------------------------------------------------- |
---|
95 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CdN |
---|
96 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: ChN |
---|
97 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CeN |
---|
98 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xz0 ! Aerodynamic roughness length [m] |
---|
99 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xu_star ! u*, friction velocity |
---|
100 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xL ! zeta (zu/L) |
---|
101 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xUN10 ! Neutral wind at zu |
---|
102 | !!---------------------------------------------------------------------------------- |
---|
103 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Ubzu |
---|
104 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp0, ztmp1, ztmp2 ! temporary stuff |
---|
105 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z0, dt_zu, dq_zu |
---|
106 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: u_star, t_star, q_star |
---|
107 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: znu_a !: Nu_air = kinematic viscosity of air |
---|
108 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_u, zeta_t ! stability parameter at height zu |
---|
109 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z0tq |
---|
110 | !! |
---|
111 | INTEGER :: jit |
---|
112 | LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U |
---|
113 | !! |
---|
114 | LOGICAL :: lreturn_cdn=.FALSE., lreturn_chn=.FALSE., lreturn_cen=.FALSE. |
---|
115 | LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. |
---|
116 | !! |
---|
117 | CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ice_an05@sbcblk_algo_ice_an05.f90' |
---|
118 | !!---------------------------------------------------------------------------------- |
---|
119 | ALLOCATE ( Ubzu(jpi,jpj), u_star(jpi,jpj), t_star(jpi,jpj), q_star(jpi,jpj), & |
---|
120 | & zeta_u(jpi,jpj), dt_zu(jpi,jpj), dq_zu(jpi,jpj), & |
---|
121 | & znu_a(jpi,jpj), ztmp1(jpi,jpj), ztmp2(jpi,jpj), & |
---|
122 | & z0(jpi,jpj), z0tq(jpi,jpj,2), ztmp0(jpi,jpj) ) |
---|
123 | |
---|
124 | lreturn_cdn = PRESENT(CdN) |
---|
125 | lreturn_chn = PRESENT(ChN) |
---|
126 | lreturn_cen = PRESENT(CeN) |
---|
127 | lreturn_z0 = PRESENT(xz0) |
---|
128 | lreturn_ustar = PRESENT(xu_star) |
---|
129 | lreturn_L = PRESENT(xL) |
---|
130 | lreturn_UN10 = PRESENT(xUN10) |
---|
131 | |
---|
132 | l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) |
---|
133 | IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) |
---|
134 | |
---|
135 | !! Scalar wind speed cannot be below 0.2 m/s |
---|
136 | Ubzu = MAX( U_zu, wspd_thrshld_ice ) |
---|
137 | |
---|
138 | !! First guess of temperature and humidity at height zu: |
---|
139 | t_zu_i = MAX( t_zt , 100._wp ) ! who knows what's given on masked-continental regions... |
---|
140 | q_zu_i = MAX( q_zt , 0.1e-6_wp ) ! " |
---|
141 | |
---|
142 | !! Air-Ice differences (and we don't want it to be 0!) |
---|
143 | dt_zu = t_zu_i - Ts_i ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) |
---|
144 | dq_zu = q_zu_i - qs_i ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) |
---|
145 | |
---|
146 | znu_a = visc_air(t_zu_i) ! Air viscosity (m^2/s) at zt given from temperature in (K) |
---|
147 | |
---|
148 | !! Very crude first guesses of z0: |
---|
149 | z0 = 8.0E-4_wp |
---|
150 | |
---|
151 | !! Crude first guess of turbulent scales |
---|
152 | u_star = 0.035_wp*Ubzu*LOG( 10._wp/z0 )/LOG( zu/z0 ) |
---|
153 | z0 = rough_leng_m( u_star , znu_a ) |
---|
154 | |
---|
155 | DO jit = 1, 2 |
---|
156 | u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)) , 1.E-9 ) |
---|
157 | z0 = rough_leng_m( u_star , znu_a ) |
---|
158 | END DO |
---|
159 | |
---|
160 | z0tq = rough_leng_tq( z0, u_star , znu_a ) |
---|
161 | t_star = dt_zu*vkarmn/(LOG(zu/z0tq(:,:,1))) |
---|
162 | q_star = dq_zu*vkarmn/(LOG(zu/z0tq(:,:,2))) |
---|
163 | |
---|
164 | |
---|
165 | !! ITERATION BLOCK |
---|
166 | DO jit = 1, nbit |
---|
167 | |
---|
168 | !!Inverse of Obukov length (1/L) : |
---|
169 | ztmp0 = One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star) ! 1/L == 1/[Obukhov length] |
---|
170 | ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) |
---|
171 | |
---|
172 | !! Stability parameters "zeta" : |
---|
173 | zeta_u = zu*ztmp0 |
---|
174 | zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u ) |
---|
175 | IF( .NOT. l_zt_equal_zu ) THEN |
---|
176 | zeta_t = zt*ztmp0 |
---|
177 | zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) |
---|
178 | END IF |
---|
179 | |
---|
180 | !! Roughness lengthes z0, z0t, & z0q : |
---|
181 | z0 = rough_leng_m ( u_star , znu_a ) |
---|
182 | z0tq = rough_leng_tq( z0, u_star , znu_a ) |
---|
183 | |
---|
184 | !! Turbulent scales at zu : |
---|
185 | ztmp0 = psi_h_ice(zeta_u) |
---|
186 | t_star = dt_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,1)) - ztmp0) |
---|
187 | q_star = dq_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,2)) - ztmp0) |
---|
188 | u_star = MAX( Ubzu*vkarmn/(LOG(zu) - LOG(z0(:,:)) - psi_m_ice(zeta_u)) , 1.E-9 ) |
---|
189 | |
---|
190 | IF( .NOT. l_zt_equal_zu ) THEN |
---|
191 | !! Re-updating temperature and humidity at zu if zt /= zu : |
---|
192 | ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_ice(zeta_t) |
---|
193 | t_zu_i = t_zt - t_star/vkarmn*ztmp1 |
---|
194 | q_zu_i = q_zt - q_star/vkarmn*ztmp1 |
---|
195 | dt_zu = t_zu_i - Ts_i ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) |
---|
196 | dq_zu = q_zu_i - qs_i ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) |
---|
197 | END IF |
---|
198 | |
---|
199 | END DO !DO jit = 1, nbit |
---|
200 | |
---|
201 | ! compute transfer coefficients at zu : |
---|
202 | ztmp0 = u_star/Ubzu |
---|
203 | Cd_i = ztmp0*ztmp0 |
---|
204 | Ch_i = ztmp0*t_star/dt_zu |
---|
205 | Ce_i = ztmp0*q_star/dq_zu |
---|
206 | |
---|
207 | IF( lreturn_cdn .OR. lreturn_chn .OR. lreturn_cen ) ztmp0 = 1._wp/LOG( zu/z0(:,:) ) |
---|
208 | IF( lreturn_cdn ) CdN = vkarmn2*ztmp0*ztmp0 |
---|
209 | IF( lreturn_chn ) ChN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,1)) |
---|
210 | IF( lreturn_cen ) CeN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,2)) |
---|
211 | |
---|
212 | IF( lreturn_z0 ) xz0 = z0 |
---|
213 | IF( lreturn_ustar ) xu_star = u_star |
---|
214 | IF( lreturn_L ) xL = 1./One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star) |
---|
215 | IF( lreturn_UN10 ) xUN10 = u_star/vkarmn*LOG(10./z0) |
---|
216 | |
---|
217 | DEALLOCATE ( Ubzu, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu, z0, z0tq, znu_a, ztmp0, ztmp1, ztmp2 ) |
---|
218 | IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) |
---|
219 | |
---|
220 | END SUBROUTINE turb_ice_an05 |
---|
221 | |
---|
222 | |
---|
223 | |
---|
224 | FUNCTION rough_leng_m( pus , pnua ) |
---|
225 | !!---------------------------------------------------------------------------------- |
---|
226 | !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 19) |
---|
227 | !! |
---|
228 | !! Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
229 | !!---------------------------------------------------------------------------------- |
---|
230 | REAL(wp), DIMENSION(jpi,jpj) :: rough_leng_m ! roughness length over sea-ice [m] |
---|
231 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! u* = friction velocity [m/s] |
---|
232 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua ! kinematic viscosity of air [m^2/s] |
---|
233 | !! |
---|
234 | INTEGER :: ji, jj ! dummy loop indices |
---|
235 | REAL(wp) :: zus, zz |
---|
236 | !!---------------------------------------------------------------------------------- |
---|
237 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
238 | zus = MAX( pus(ji,jj) , 1.E-9_wp ) |
---|
239 | |
---|
240 | zz = (zus - 0.18_wp) / 0.1_wp |
---|
241 | |
---|
242 | rough_leng_m(ji,jj) = 0.135*pnua(ji,jj)/zus + 0.035*zus*zus/grav*( 5.*EXP(-zz*zz) + 1._wp ) ! Eq.(19) Andreas et al., 2005 |
---|
243 | END_2D |
---|
244 | !! |
---|
245 | END FUNCTION rough_leng_m |
---|
246 | |
---|
247 | FUNCTION rough_leng_tq( pz0, pus , pnua ) |
---|
248 | !!---------------------------------------------------------------------------------- |
---|
249 | !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 22) |
---|
250 | !! => which still relies on Andreas 1987 ! |
---|
251 | !! |
---|
252 | !! Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
253 | !!---------------------------------------------------------------------------------- |
---|
254 | REAL(wp), DIMENSION(jpi,jpj,2) :: rough_leng_tq ! temp.,hum. roughness lengthes over sea-ice [m] |
---|
255 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 ! roughness length [m] |
---|
256 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! u* = friction velocity [m/s] |
---|
257 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua ! kinematic viscosity of air [m^2/s] |
---|
258 | !! |
---|
259 | INTEGER :: ji, jj ! dummy loop indices |
---|
260 | REAL(wp) :: zz0, zus, zre, zsmoot, ztrans, zrough |
---|
261 | REAL(wp) :: zb0, zb1, zb2, zlog, zlog2, zlog_z0s_on_z0 |
---|
262 | !!---------------------------------------------------------------------------------- |
---|
263 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
264 | zz0 = pz0(ji,jj) |
---|
265 | zus = MAX( pus(ji,jj) , 1.E-9_wp ) |
---|
266 | zre = MAX( zus*zz0/pnua(ji,jj) , 0._wp ) ! Roughness Reynolds number |
---|
267 | |
---|
268 | !! *** TABLE 1 of Andreas et al. 2005 *** |
---|
269 | zsmoot = 0._wp ; ztrans = 0._wp ; zrough = 0._wp |
---|
270 | IF ( zre <= 0.135_wp ) THEN ! Smooth flow condition (R* <= 0.135): |
---|
271 | zsmoot = 1._wp |
---|
272 | ELSEIF( zre < 2.5_wp ) THEN ! Transition (0.135 < R* < 2.5) |
---|
273 | ztrans = 1._wp |
---|
274 | ELSE ! Rough ( R* > 2.5) |
---|
275 | zrough = 1._wp |
---|
276 | ENDIF |
---|
277 | |
---|
278 | zlog = LOG(zre) |
---|
279 | zlog2 = zlog*zlog |
---|
280 | |
---|
281 | !! z0t: |
---|
282 | zb0 = zsmoot*1.25_wp + ztrans*0.149_wp + zrough*0.317_wp |
---|
283 | zb1 = - ztrans*0.550_wp - zrough*0.565_wp |
---|
284 | zb2 = - zrough*0.183_wp |
---|
285 | zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2 |
---|
286 | rough_leng_tq(ji,jj,1) = zz0 * EXP( zlog_z0s_on_z0 ) |
---|
287 | |
---|
288 | !! z0q: |
---|
289 | zb0 = zsmoot*1.61_wp + ztrans*0.351_wp + zrough*0.396_wp |
---|
290 | zb1 = - ztrans*0.628_wp - zrough*0.512_wp |
---|
291 | zb2 = - zrough*0.180_wp |
---|
292 | zlog = LOG(zre) |
---|
293 | zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2 |
---|
294 | rough_leng_tq(ji,jj,2) = zz0 * EXP( zlog_z0s_on_z0 ) |
---|
295 | |
---|
296 | END_2D |
---|
297 | !! |
---|
298 | END FUNCTION rough_leng_tq |
---|
299 | |
---|
300 | |
---|
301 | |
---|
302 | FUNCTION psi_m_ice( pzeta ) |
---|
303 | !!---------------------------------------------------------------------------------- |
---|
304 | !! ** Purpose: compute the universal profile stability function for momentum |
---|
305 | !! |
---|
306 | !! |
---|
307 | !! Andreas et al 2005 == Jordan et al. 1999 |
---|
308 | !! |
---|
309 | !! Psi: |
---|
310 | !! Unstable => Paulson 1970 |
---|
311 | !! Stable => Holtslag & De Bruin 1988 |
---|
312 | !! |
---|
313 | !! pzeta : stability paramenter, z/L where z is altitude |
---|
314 | !! measurement and L is M-O length |
---|
315 | !! |
---|
316 | !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
317 | !!---------------------------------------------------------------------------------- |
---|
318 | REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ice |
---|
319 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
320 | ! |
---|
321 | INTEGER :: ji, jj ! dummy loop indices |
---|
322 | REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab |
---|
323 | !!---------------------------------------------------------------------------------- |
---|
324 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! |
---|
325 | zta = pzeta(ji,jj) |
---|
326 | ! |
---|
327 | ! Unstable stratification: |
---|
328 | zx = ABS(1._wp - 16._wp*zta)**.25 ! (16 here, not 15!) |
---|
329 | |
---|
330 | zpsi_u = LOG( (1._wp + zx*zx)/2. ) & ! Eq.(30) Jordan et al. 1999 |
---|
331 | & + 2.*LOG( (1._wp + zx )/2. ) & |
---|
332 | & - 2.*ATAN( zx ) + 0.5*rpi |
---|
333 | |
---|
334 | ! Stable stratification: |
---|
335 | zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp ) ! Eq.(33) Jordan et al. 1999 |
---|
336 | |
---|
337 | !! Combine: |
---|
338 | zstab = 0.5 + SIGN(0.5_wp, zta) |
---|
339 | psi_m_ice(ji,jj) = (1._wp - zstab) * zpsi_u & ! Unstable (zta < 0) |
---|
340 | & + zstab * zpsi_s ! Stable (zta > 0) |
---|
341 | ! |
---|
342 | END_2D |
---|
343 | END FUNCTION psi_m_ice |
---|
344 | |
---|
345 | |
---|
346 | FUNCTION psi_h_ice( pzeta ) |
---|
347 | !!---------------------------------------------------------------------------------- |
---|
348 | !! ** Purpose: compute the universal profile stability function for |
---|
349 | !! temperature and humidity |
---|
350 | !! |
---|
351 | !! |
---|
352 | !! Andreas et al 2005 == Jordan et al. 1999 |
---|
353 | !! |
---|
354 | !! Psi: |
---|
355 | !! Unstable => Paulson 1970 |
---|
356 | !! Stable => Holtslag & De Bruin 1988 |
---|
357 | !! |
---|
358 | !! pzeta : stability paramenter, z/L where z is altitude |
---|
359 | !! measurement and L is M-O length |
---|
360 | !! |
---|
361 | !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
362 | !!---------------------------------------------------------------------------------- |
---|
363 | REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ice |
---|
364 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
365 | ! |
---|
366 | INTEGER :: ji, jj ! dummy loop indices |
---|
367 | REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab |
---|
368 | !!---------------------------------------------------------------------------------- |
---|
369 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! |
---|
370 | zta = pzeta(ji,jj) |
---|
371 | ! |
---|
372 | ! Unstable stratification: |
---|
373 | zx = ABS(1._wp - 16._wp*zta)**.25 ! (16 here, not 15!) |
---|
374 | |
---|
375 | zpsi_u = 2.*LOG( (1._wp + zx*zx)/2. ) ! Eq.(31) Jordan et al. 1999 |
---|
376 | |
---|
377 | ! Stable stratification (identical to Psi_m!): |
---|
378 | zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp ) ! Eq.(33) Jordan et al. 1999 |
---|
379 | |
---|
380 | !! Combine: |
---|
381 | zstab = 0.5 + SIGN(0.5_wp, zta) |
---|
382 | psi_h_ice(ji,jj) = (1._wp - zstab) * zpsi_u & ! Unstable (zta < 0) |
---|
383 | & + zstab * zpsi_s ! Stable (zta > 0) |
---|
384 | ! |
---|
385 | END_2D |
---|
386 | END FUNCTION psi_h_ice |
---|
387 | |
---|
388 | !!====================================================================== |
---|
389 | END MODULE sbcblk_algo_ice_an05 |
---|