1 | !!! TO DO: consistent psi_m and psi_h needed!!! For now is those of NCAR !!! |
---|
2 | !! |
---|
3 | MODULE sbcblk_algo_andreas |
---|
4 | !!====================================================================== |
---|
5 | !! *** MODULE sbcblk_algo_andreas *** |
---|
6 | !! Computes: |
---|
7 | !! * bulk transfer coefficients C_D, C_E and C_H |
---|
8 | !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed |
---|
9 | !! * the effective bulk wind speed at 10m Ubzu |
---|
10 | !! according to Andreas et al. (2015) |
---|
11 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
12 | !! Andreas, E.L., Mahrt, L. and Vickers, D. (2015), |
---|
13 | !! An improved bulk air–sea surface flux algorithm, |
---|
14 | !! including spray‐mediated transfer. |
---|
15 | !! Q.J.R. Meteorol. Soc., 141: 642-654. doi:10.1002/qj.2424 |
---|
16 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
17 | !! |
---|
18 | !! * bulk transfer coefficients C_D, C_E and C_H |
---|
19 | !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed |
---|
20 | !! * the effective bulk wind speed at z=zu: Ubzu |
---|
21 | !! => all these are used in bulk formulas in sbcblk.F90 |
---|
22 | !! |
---|
23 | !! Using the bulk formulation/param. of Large & Yeager 2008 |
---|
24 | !! |
---|
25 | !! Routine turb_andreas maintained and developed in AeroBulk |
---|
26 | !! (https://github.com/brodeau/aerobulk/) |
---|
27 | !! |
---|
28 | !! ** Author: L. Brodeau, August 2020 / AeroBulk (https://github.com/brodeau/aerobulk) |
---|
29 | !!---------------------------------------------------------------------- |
---|
30 | !! History : 4.x ! 2020-08 (L.Brodeau) Original code |
---|
31 | !!---------------------------------------------------------------------- |
---|
32 | |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | !!---------------------------------------------------------------------- |
---|
35 | USE dom_oce ! ocean space and time domain |
---|
36 | USE phycst ! physical constants |
---|
37 | USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer |
---|
38 | |
---|
39 | IMPLICIT NONE |
---|
40 | PRIVATE |
---|
41 | |
---|
42 | !! Important (Brodeau fix): |
---|
43 | REAL(wp), PARAMETER :: rRi_max = 0.15_wp ! Bulk Ri above which the algorithm fucks up! |
---|
44 | ! ! (increasing (>0) Ri means that surface layer increasingly stable and/or wind increasingly weak) |
---|
45 | REAL(wp), PARAMETER :: rCs_min = 0.35E-3_wp ! minimum value to tolarate for CE and CH ! Must be larger than "Cx_min" !!! |
---|
46 | |
---|
47 | PUBLIC :: TURB_ANDREAS, psi_m_andreas, psi_h_andreas |
---|
48 | |
---|
49 | !! * Substitutions |
---|
50 | # include "do_loop_substitute.h90" |
---|
51 | |
---|
52 | !!---------------------------------------------------------------------- |
---|
53 | CONTAINS |
---|
54 | |
---|
55 | SUBROUTINE turb_andreas( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & |
---|
56 | & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & |
---|
57 | & nb_iter, CdN, ChN, CeN ) |
---|
58 | !!---------------------------------------------------------------------------------- |
---|
59 | !! *** ROUTINE turb_andreas *** |
---|
60 | !! |
---|
61 | !! ** Purpose : Computes turbulent transfert coefficients of surface |
---|
62 | !! fluxes according to Large & Yeager (2004) and Large & Yeager (2008) |
---|
63 | !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu |
---|
64 | !! Returns the effective bulk wind speed at zu to be used in the bulk formulas |
---|
65 | !! |
---|
66 | !! INPUT : |
---|
67 | !! ------- |
---|
68 | !! * zt : height for temperature and spec. hum. of air [m] |
---|
69 | !! * zu : height for wind speed (usually 10m) [m] |
---|
70 | !! * sst : bulk SST [K] |
---|
71 | !! * t_zt : potential air temperature at zt [K] |
---|
72 | !! * ssq : specific humidity at saturation at SST [kg/kg] |
---|
73 | !! * q_zt : specific humidity of air at zt [kg/kg] |
---|
74 | !! * U_zu : scalar wind speed at zu [m/s] |
---|
75 | !! |
---|
76 | !! OUTPUT : |
---|
77 | !! -------- |
---|
78 | !! * Cd : drag coefficient |
---|
79 | !! * Ch : sensible heat coefficient |
---|
80 | !! * Ce : evaporation coefficient |
---|
81 | !! * t_zu : pot. air temperature adjusted at wind height zu [K] |
---|
82 | !! * q_zu : specific humidity of air // [kg/kg] |
---|
83 | !! * Ubzu : bulk wind speed at zu [m/s] |
---|
84 | !! |
---|
85 | !! |
---|
86 | !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
87 | !!---------------------------------------------------------------------------------- |
---|
88 | REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] |
---|
89 | REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] |
---|
90 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin] |
---|
91 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] |
---|
92 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg] |
---|
93 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] |
---|
94 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] |
---|
95 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) |
---|
96 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) |
---|
97 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat) |
---|
98 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] |
---|
99 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] |
---|
100 | REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] |
---|
101 | ! |
---|
102 | INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations |
---|
103 | REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN |
---|
104 | REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN |
---|
105 | REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN |
---|
106 | ! |
---|
107 | INTEGER :: nbit, jit ! iterations... |
---|
108 | LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U |
---|
109 | !! |
---|
110 | REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star |
---|
111 | REAL(wp), DIMENSION(jpi,jpj) :: z0 ! roughness length (momentum) [m] |
---|
112 | REAL(wp), DIMENSION(jpi,jpj) :: UN10 ! Neutral wind speed at zu [m/s] |
---|
113 | REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu |
---|
114 | REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 |
---|
115 | REAL(wp), DIMENSION(jpi,jpj) :: RiB ! square root of Cd |
---|
116 | !! |
---|
117 | !!---------------------------------------------------------------------------------- |
---|
118 | nbit = nb_iter0 |
---|
119 | IF( PRESENT(nb_iter) ) nbit = nb_iter |
---|
120 | |
---|
121 | l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision |
---|
122 | |
---|
123 | 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 |
---|
124 | |
---|
125 | !! First guess: |
---|
126 | UN10 = Ubzu |
---|
127 | Cd = 1.1E-3_wp |
---|
128 | Ch = 1.1E-3_wp |
---|
129 | Ce = 1.1E-3_wp |
---|
130 | t_zu = t_zt |
---|
131 | q_zu = q_zt |
---|
132 | |
---|
133 | !! First guess of turbulent scales for scalars: |
---|
134 | ztmp0 = SQRT(Cd) |
---|
135 | t_star = Ch/ztmp0*(t_zu - sst) ! theta* |
---|
136 | q_star = Ce/ztmp0*(q_zu - ssq) ! q* |
---|
137 | |
---|
138 | ! Bulk Richardson number: |
---|
139 | RiB(:,:) = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu ) |
---|
140 | |
---|
141 | |
---|
142 | !! ITERATION BLOCK |
---|
143 | DO jit = 1, nbit |
---|
144 | |
---|
145 | WHERE ( RiB < rRi_max ) |
---|
146 | !! Normal condition case: |
---|
147 | u_star = U_STAR_ANDREAS( UN10 ) |
---|
148 | ELSEWHERE |
---|
149 | !! Extremely stable + weak wind !!! |
---|
150 | !! => for we force u* to be consistent with minimum value for CD: |
---|
151 | !! (otherwize algorithm becomes nonsense...) |
---|
152 | u_star = SQRT(Cx_min) * Ubzu ! Cd does not go below Cx_min ! |
---|
153 | ENDWHERE |
---|
154 | |
---|
155 | !! Stability parameter : |
---|
156 | zeta_u = zu*One_on_L( t_zu, q_zu, u_star, t_star, q_star ) ! zu * 1/L |
---|
157 | |
---|
158 | !! Drag coefficient: |
---|
159 | ztmp0 = u_star/Ubzu |
---|
160 | |
---|
161 | Cd = MAX( ztmp0*ztmp0 , Cx_min ) |
---|
162 | |
---|
163 | !! Roughness length: |
---|
164 | z0 = MIN( z0_from_Cd( zu, Cd, ppsi=psi_m_andreas(zeta_u) ) , z0_sea_max ) |
---|
165 | |
---|
166 | !! z0t and z0q, based on LKB, just like into COARE 2.5: |
---|
167 | ztmp0 = z0 * u_star / visc_air(t_zu) ! Re_r |
---|
168 | ztmp1 = z0tq_LKB( 1, ztmp0, z0 ) ! z0t |
---|
169 | ztmp2 = z0tq_LKB( 2, ztmp0, z0 ) ! z0q |
---|
170 | |
---|
171 | !! Turbulent scales at zu : |
---|
172 | ztmp0 = psi_h_andreas(zeta_u) ! lolo: zeta_u for scalars??? |
---|
173 | t_star = (t_zu - sst)*vkarmn/(LOG(zu) - LOG(ztmp1) - ztmp0) ! theta* (ztmp1 == z0t in rhs term) |
---|
174 | q_star = (q_zu - ssq)*vkarmn/(LOG(zu) - LOG(ztmp2) - ztmp0) ! q* (ztmp2 == z0q in rhs term) |
---|
175 | |
---|
176 | IF( (.NOT. l_zt_equal_zu).AND.( jit > 1 ) ) THEN |
---|
177 | !! Re-updating temperature and humidity at zu if zt /= zu: |
---|
178 | ztmp0 = zeta_u/zu*zt ! zeta_t |
---|
179 | ztmp0 = LOG(zt/zu) + psi_h_andreas(zeta_u) - psi_h_andreas(ztmp0) |
---|
180 | t_zu = t_zt - t_star/vkarmn*ztmp0 |
---|
181 | q_zu = q_zt - q_star/vkarmn*ztmp0 |
---|
182 | RiB = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu ) !LOLO |
---|
183 | ENDIF |
---|
184 | |
---|
185 | !! Update neutral-stability wind at zu: |
---|
186 | UN10 = MAX( 0.1_wp , UN10_from_ustar( zu, Ubzu, u_star, psi_m_andreas(zeta_u) ) ) ! UN10 |
---|
187 | |
---|
188 | END DO !DO jit = 1, nbit |
---|
189 | |
---|
190 | ! Compute transfer coefficients at zu: |
---|
191 | ztmp0 = u_star/Ubzu |
---|
192 | |
---|
193 | Cd = MAX( ztmp0*ztmp0 , Cx_min ) ! the earlier use of Cx_min on u* should make use of Cx_min here unnecessary! |
---|
194 | |
---|
195 | ztmp1 = t_zu - sst ; ztmp1 = SIGN( MAX(ABS(ztmp1),1.E-6_wp), ztmp1 ) ! dt_zu |
---|
196 | ztmp2 = q_zu - ssq ; ztmp2 = SIGN( MAX(ABS(ztmp2),1.E-9_wp), ztmp2 ) ! dq_zu |
---|
197 | Ch = MAX( ztmp0*t_star/ztmp1 , rCs_min ) |
---|
198 | Ce = MAX( ztmp0*q_star/ztmp2 , rCs_min ) |
---|
199 | |
---|
200 | !! Neutral-stability coefficients: |
---|
201 | ztmp0 = 1._wp/LOG(zu/z0) |
---|
202 | ztmp1 = z0 * u_star / visc_air(t_zu) ! Re_r |
---|
203 | |
---|
204 | IF(PRESENT(CdN)) CdN = vkarmn2*ztmp0*ztmp0 |
---|
205 | IF(PRESENT(ChN)) ChN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 1, ztmp1, z0 )) |
---|
206 | IF(PRESENT(CeN)) CeN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 2, ztmp1, z0 )) |
---|
207 | |
---|
208 | END SUBROUTINE turb_andreas |
---|
209 | |
---|
210 | |
---|
211 | FUNCTION U_STAR_ANDREAS( pun10 ) |
---|
212 | !!---------------------------------------------------------------------------------- |
---|
213 | !! Estimate of the friction velocity as a function of the neutral-stability wind |
---|
214 | !! speed at at 10m |
---|
215 | !! |
---|
216 | !! Origin: Eq.(2.2) of Andreas et al. (2015) |
---|
217 | !! |
---|
218 | !! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
219 | !!---------------------------------------------------------------------------------- |
---|
220 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pun10 !: neutral-stability scalar wind speed at 10m (m/s) |
---|
221 | REAL(wp), DIMENSION(jpi,jpj) :: u_star_andreas !: friction velocity [m/s] |
---|
222 | ! |
---|
223 | INTEGER :: ji, jj ! dummy loop indices |
---|
224 | REAL(wp) :: za, zt, zw ! local scalars |
---|
225 | !!---------------------------------------------------------------------------------- |
---|
226 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
227 | zw = pun10(ji,jj) |
---|
228 | za = zw - 8.271_wp |
---|
229 | zt = za + SQRT( 0.12_wp*za*za + 0.181_wp ) |
---|
230 | u_star_andreas(ji,jj) = 0.239_wp + 0.0433_wp * zt |
---|
231 | END_2D |
---|
232 | END FUNCTION U_STAR_ANDREAS |
---|
233 | |
---|
234 | |
---|
235 | FUNCTION psi_m_andreas( pzeta ) |
---|
236 | !!---------------------------------------------------------------------------------- |
---|
237 | !! Universal profile stability function for momentum |
---|
238 | !! TO DO !!!!!!!!!!!!!!!!!!!!! |
---|
239 | !! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE |
---|
240 | !! |
---|
241 | !! pzeta : stability paramenter, z/L where z is altitude measurement |
---|
242 | !! and L is M-O length |
---|
243 | !! |
---|
244 | !! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
245 | !!---------------------------------------------------------------------------------- |
---|
246 | REAL(wp), DIMENSION(jpi,jpj) :: psi_m_andreas |
---|
247 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
248 | ! |
---|
249 | REAL(wp), PARAMETER :: zam = 5._wp ! a_m (just below Eq.(9b) |
---|
250 | REAL(wp), PARAMETER :: zbm = zam/6.5_wp ! b_m (just below Eq.(9b) |
---|
251 | ! |
---|
252 | REAL(wp), PARAMETER :: z1o3 = 1._wp/3._wp |
---|
253 | REAL(wp), PARAMETER :: zsr3 = SQRT(3._wp) |
---|
254 | ! |
---|
255 | INTEGER :: ji, jj ! dummy loop indices |
---|
256 | REAL(wp) :: zta, zx2, zx, zpsi_unst, zbbm, zpsi_stab, zstab ! local scalars |
---|
257 | !!---------------------------------------------------------------------------------- |
---|
258 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
259 | ! |
---|
260 | zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and big!) |
---|
261 | ! |
---|
262 | !! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!! |
---|
263 | zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5 |
---|
264 | zx2 = MAX( zx2 , 1._wp ) |
---|
265 | zx = SQRT(zx2) ! (1 - 16z)^0.25 |
---|
266 | zpsi_unst = 2._wp*LOG(ABS( (1._wp + zx )*0.5_wp )) & |
---|
267 | & + LOG(ABS( (1._wp + zx2)*0.5_wp )) & |
---|
268 | & - 2._wp*ATAN(zx) + rpi*0.5_wp |
---|
269 | ! |
---|
270 | !! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(12) Grachev et al 2007]: |
---|
271 | zx = ABS(1._wp + zta)**z1o3 |
---|
272 | zbbm = ABS( (1._wp - zbm)/zbm )**z1o3 ! B_m |
---|
273 | ! |
---|
274 | zpsi_stab = -3.*zam/zbm*(zx - 1._wp) + zam*zbbm/(2.*zbm) * ( & |
---|
275 | & 2.*LOG(ABS( ( zx + zbbm )/(1._wp + zbbm ) )) & |
---|
276 | & - LOG(ABS( (zx*zx - zx*zbbm + zbbm*zbbm)/(1._wp - zbbm + zbbm*zbbm) )) & |
---|
277 | & + 2.*zsr3*( ATAN( (2.*zx - zbbm)/(zsr3*zbbm) ) - ATAN( (2._wp - zbbm)/(zsr3*zbbm) ) ) ) |
---|
278 | ! |
---|
279 | ! |
---|
280 | zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 |
---|
281 | ! |
---|
282 | psi_m_andreas(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable |
---|
283 | & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable |
---|
284 | ! |
---|
285 | END_2D |
---|
286 | END FUNCTION psi_m_andreas |
---|
287 | |
---|
288 | |
---|
289 | FUNCTION psi_h_andreas( pzeta ) |
---|
290 | !!---------------------------------------------------------------------------------- |
---|
291 | !! Universal profile stability function for temperature and humidity |
---|
292 | !! |
---|
293 | !! TO DO |
---|
294 | !! !! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE |
---|
295 | !! |
---|
296 | !! pzeta : stability paramenter, z/L where z is altitude measurement |
---|
297 | !! and L is M-O length |
---|
298 | !! |
---|
299 | !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
300 | !!---------------------------------------------------------------------------------- |
---|
301 | REAL(wp), DIMENSION(jpi,jpj) :: psi_h_andreas |
---|
302 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
303 | ! |
---|
304 | REAL(wp), PARAMETER :: zah = 5._wp ! a_h (just below Eq.(9b) |
---|
305 | REAL(wp), PARAMETER :: zbh = 5._wp ! b_h (just below Eq.(9b) |
---|
306 | REAL(wp), PARAMETER :: zch = 3._wp ! c_h (just below Eq.(9b) |
---|
307 | REAL(wp), PARAMETER :: zbbh = SQRT(5._wp) ! B_h (just below Eq.(13) |
---|
308 | ! |
---|
309 | INTEGER :: ji, jj ! dummy loop indices |
---|
310 | REAL(wp) :: zta, zz, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars |
---|
311 | !!---------------------------------------------------------------------------------- |
---|
312 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
313 | ! |
---|
314 | zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and large!) |
---|
315 | ! |
---|
316 | !! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!! |
---|
317 | zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5 |
---|
318 | zx2 = MAX( zx2 , 1._wp ) |
---|
319 | zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) |
---|
320 | ! |
---|
321 | !! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(13) Grachev et al 2007]: |
---|
322 | zz = 2.*zta + zch |
---|
323 | zpsi_stab = - 0.5*zbh*LOG(ABS(1._wp + zch*zta + zta*zta)) & |
---|
324 | & + (-zah/zbbh + 0.5*zbh*zch/zbbh) & |
---|
325 | & *( LOG(ABS((zz - zbbh)/(zz + zbbh))) & |
---|
326 | & - LOG(ABS((zch - zbbh)/(zch + zbbh))) ) |
---|
327 | ! |
---|
328 | zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 |
---|
329 | ! |
---|
330 | psi_h_andreas(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable |
---|
331 | & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable |
---|
332 | ! |
---|
333 | END_2D |
---|
334 | END FUNCTION psi_h_andreas |
---|
335 | |
---|
336 | !!====================================================================== |
---|
337 | END MODULE sbcblk_algo_andreas |
---|