Line data Source code
1 : !=======================================================================
2 :
3 : ! Atmospheric boundary interface (stability based flux calculations)
4 :
5 : ! author: Elizabeth C. Hunke, LANL
6 : !
7 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
8 : ! 2004: Block structure added by William Lipscomb
9 : ! 2006: Converted to free source form (F90) by Elizabeth Hunke
10 : ! 2013: Form drag routine added (neutral_drag_coeffs) by David Schroeder
11 : ! 2014: Adjusted form drag and added high frequency coupling by Andrew Roberts
12 :
13 : module icepack_atmo
14 :
15 : use icepack_kinds
16 : use icepack_parameters, only: c0, c1, c2, c4, c5, c8, c10
17 : use icepack_parameters, only: c16, c20, p001, p01, p2, p4, p5, p75, puny
18 : use icepack_parameters, only: senscoef, latncoef
19 : use icepack_parameters, only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
20 : use icepack_parameters, only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
21 : use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
22 : use icepack_parameters, only: atmbndy, calc_strair, formdrag
23 : use icepack_parameters, only: icepack_chkoptargflag
24 : use icepack_tracers, only: ncat, n_iso
25 : use icepack_tracers, only: tr_iso
26 : use icepack_warnings, only: warnstr, icepack_warnings_add
27 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
28 :
29 : implicit none
30 :
31 : private
32 : public :: neutral_drag_coeffs, &
33 : icepack_atm_boundary
34 :
35 : !=======================================================================
36 :
37 : contains
38 :
39 : !=======================================================================
40 :
41 : ! Compute coefficients for atm/ice fluxes, stress, and reference
42 : ! temperature and humidity. NOTE:
43 : ! (1) All fluxes are positive downward,
44 : ! (2) Here, tstar = (WT)/U*, and qstar = (WQ)/U*,
45 : ! (3a) wind speeds should all be above a minimum speed (eg. 1.0 m/s).
46 : !
47 : ! ASSUME:
48 : ! The saturation humidity of air at T(K): qsat(T) (kg/m**3)
49 : !
50 : ! Code originally based on CSM1
51 :
52 3558446561 : subroutine atmo_boundary_layer (sfctype, &
53 : calc_strair, formdrag, & ! LCOV_EXCL_LINE
54 : Tsf, potT, & ! LCOV_EXCL_LINE
55 : uatm, vatm, & ! LCOV_EXCL_LINE
56 : wind, zlvl, & ! LCOV_EXCL_LINE
57 : Qa, rhoa, & ! LCOV_EXCL_LINE
58 : strx, stry, & ! LCOV_EXCL_LINE
59 : Tref, Qref, & ! LCOV_EXCL_LINE
60 : delt, delq, & ! LCOV_EXCL_LINE
61 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
62 : Cdn_atm, & ! LCOV_EXCL_LINE
63 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
64 3558446561 : Qa_iso, Qref_iso, & ! LCOV_EXCL_LINE
65 : uvel, vvel, & ! LCOV_EXCL_LINE
66 : Uref, zlvs )
67 :
68 : use icepack_parameters, only: highfreq, natmiter, atmiter_conv
69 :
70 : character (len=3), intent(in) :: &
71 : sfctype ! ice or ocean
72 :
73 : logical (kind=log_kind), intent(in) :: &
74 : calc_strair, & ! if true, calculate wind stress components ! LCOV_EXCL_LINE
75 : formdrag ! if true, calculate form drag
76 :
77 : real (kind=dbl_kind), intent(in) :: &
78 : Tsf , & ! surface temperature of ice or ocean ! LCOV_EXCL_LINE
79 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
80 : uatm , & ! x-direction wind speed (m/s) ! LCOV_EXCL_LINE
81 : vatm , & ! y-direction wind speed (m/s) ! LCOV_EXCL_LINE
82 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
83 : zlvl , & ! atm level height (m) ! LCOV_EXCL_LINE
84 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
85 : rhoa ! air density (kg/m^3)
86 :
87 : real (kind=dbl_kind), intent(inout) :: &
88 : Cdn_atm ! neutral drag coefficient
89 :
90 : real (kind=dbl_kind), intent(inout) :: &
91 : Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff
92 :
93 : real (kind=dbl_kind), intent(inout) :: &
94 : strx , & ! x surface stress (N) ! LCOV_EXCL_LINE
95 : stry ! y surface stress (N)
96 :
97 : real (kind=dbl_kind), intent(inout) :: &
98 : Tref , & ! reference height temperature (K) ! LCOV_EXCL_LINE
99 : Qref , & ! reference height specific humidity (kg/kg) ! LCOV_EXCL_LINE
100 : delt , & ! potential T difference (K) ! LCOV_EXCL_LINE
101 : delq , & ! humidity difference (kg/kg) ! LCOV_EXCL_LINE
102 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
103 : lhcoef ! transfer coefficient for latent heat
104 :
105 : real (kind=dbl_kind), intent(in), dimension(:), optional :: &
106 : Qa_iso ! specific isotopic humidity (kg/kg)
107 :
108 : real (kind=dbl_kind), intent(inout), dimension(:), optional :: &
109 : Qref_iso ! reference specific isotopic humidity (kg/kg)
110 :
111 : real (kind=dbl_kind), intent(in) :: &
112 : uvel , & ! x-direction ice speed (m/s) ! LCOV_EXCL_LINE
113 : vvel ! y-direction ice speed (m/s)
114 :
115 : real (kind=dbl_kind), intent(out) :: &
116 : Uref ! reference height wind speed (m/s)
117 :
118 : real (kind=dbl_kind), intent(in), optional :: &
119 : zlvs ! atm level height (scalar quantities) (m)
120 :
121 : ! local variables
122 :
123 : integer (kind=int_kind) :: &
124 : k,n ! iteration index
125 :
126 : real (kind=dbl_kind) :: &
127 : TsfK , & ! surface temperature in Kelvin (K) ! LCOV_EXCL_LINE
128 : psimh , & ! stability function at zlvl (momentum) ! LCOV_EXCL_LINE
129 : tau , & ! stress at zlvl ! LCOV_EXCL_LINE
130 : fac , & ! interpolation factor ! LCOV_EXCL_LINE
131 : al2 , & ! ln(z10 /zTrf) ! LCOV_EXCL_LINE
132 : psix2 , & ! stability function at zTrf (heat and water) ! LCOV_EXCL_LINE
133 : ssq , & ! sat surface humidity (kg/kg) ! LCOV_EXCL_LINE
134 : qqq , & ! for qsat, dqsfcdt ! LCOV_EXCL_LINE
135 : TTT , & ! for qsat, dqsfcdt ! LCOV_EXCL_LINE
136 : qsat , & ! the saturation humidity of air (kg/m^3) ! LCOV_EXCL_LINE
137 : Lheat , & ! Lvap or Lsub, depending on surface type ! LCOV_EXCL_LINE
138 : umin ! minimum wind speed (m/s)
139 :
140 : real (kind=dbl_kind) :: &
141 : ustar , & ! ustar (m/s) ! LCOV_EXCL_LINE
142 : ustar_prev , & ! ustar_prev (m/s) ! LCOV_EXCL_LINE
143 : tstar , & ! tstar ! LCOV_EXCL_LINE
144 : qstar , & ! qstar ! LCOV_EXCL_LINE
145 : ratio , & ! ratio ! LCOV_EXCL_LINE
146 : rdn , & ! sqrt of neutral exchange coefficient (momentum) ! LCOV_EXCL_LINE
147 : rhn , & ! sqrt of neutral exchange coefficient (heat) ! LCOV_EXCL_LINE
148 : ren , & ! sqrt of neutral exchange coefficient (water) ! LCOV_EXCL_LINE
149 : rd , & ! sqrt of exchange coefficient (momentum) ! LCOV_EXCL_LINE
150 : re , & ! sqrt of exchange coefficient (water) ! LCOV_EXCL_LINE
151 : rh , & ! sqrt of exchange coefficient (heat) ! LCOV_EXCL_LINE
152 : vmag , & ! surface wind magnitude (m/s) ! LCOV_EXCL_LINE
153 : alzm , & ! ln(zlvl /z10) ! LCOV_EXCL_LINE
154 : alzs , & ! ln(zlvs /z10) (if zlvs present) ! LCOV_EXCL_LINE
155 : thva , & ! virtual temperature (K) ! LCOV_EXCL_LINE
156 : cp , & ! specific heat of moist air ! LCOV_EXCL_LINE
157 : holm , & ! H (at zlvl ) over L ! LCOV_EXCL_LINE
158 : hols , & ! H (at zlvs ) over L (if zlvs present) ! LCOV_EXCL_LINE
159 : stable, & ! stability factor ! LCOV_EXCL_LINE
160 : cpvir , & ! defined as cp_wv/cp_air - 1. ! LCOV_EXCL_LINE
161 : psixh ! stability function at zlvl (at zlvs if present) (heat and water)
162 :
163 : real (kind=dbl_kind), parameter :: &
164 : zTrf = c2 ! reference height for air temp (m)
165 :
166 : character(len=*),parameter :: subname='(atmo_boundary_layer)'
167 :
168 3558446561 : al2 = log(zref/zTrf)
169 :
170 : !------------------------------------------------------------
171 : ! Initialize
172 : !------------------------------------------------------------
173 :
174 3558446561 : cpvir = cp_wv/cp_air-c1 ! defined as cp_wv/cp_air - 1.
175 :
176 3558446561 : if (highfreq) then
177 91761810 : umin = p5 ! minumum allowable wind-ice speed difference of 0.5 m/s
178 : else
179 3466684751 : umin = c1 ! minumum allowable wind speed of 1m/s
180 : endif
181 :
182 3558446561 : Tref = c0
183 3558446561 : Qref = c0
184 3558446561 : Uref = c0
185 3558446561 : delt = c0
186 3558446561 : delq = c0
187 3558446561 : shcoef = c0
188 3558446561 : lhcoef = c0
189 :
190 : !------------------------------------------------------------
191 : ! Compute turbulent flux coefficients, wind stress, and
192 : ! reference temperature and humidity.
193 : !------------------------------------------------------------
194 :
195 : !------------------------------------------------------------
196 : ! define variables that depend on surface type
197 : !------------------------------------------------------------
198 :
199 3558446561 : if (sfctype(1:3)=='ice') then
200 :
201 1596104129 : qqq = qqqice ! for qsat
202 1596104129 : TTT = TTTice ! for qsat
203 1596104129 : Lheat = Lsub ! ice to vapor
204 :
205 1596104129 : if (highfreq) then
206 : vmag = max(umin, sqrt( (uatm-uvel)**2 + &
207 70577610 : (vatm-vvel)**2) )
208 : else
209 1525526519 : vmag = max(umin, wind)
210 : endif
211 :
212 1596104129 : if (formdrag .and. Cdn_atm > puny) then
213 16672872 : rdn = sqrt(Cdn_atm)
214 : else
215 1579431257 : rdn = vonkar/log(zref/iceruf) ! neutral coefficient
216 1579431257 : Cdn_atm = rdn * rdn
217 : endif
218 :
219 1962342432 : elseif (sfctype(1:3)=='ocn') then
220 :
221 1962342432 : qqq = qqqocn
222 1962342432 : TTT = TTTocn
223 1962342432 : Lheat = Lvap ! liquid to vapor
224 1962342432 : vmag = max(umin, wind)
225 : rdn = sqrt(0.0027_dbl_kind/vmag &
226 1962342432 : + .000142_dbl_kind + .0000764_dbl_kind*vmag)
227 :
228 : endif ! sfctype
229 :
230 : !------------------------------------------------------------
231 : ! define some more needed variables
232 : !------------------------------------------------------------
233 :
234 3558446561 : TsfK = Tsf + Tffresh ! surface temp (K)
235 3558446561 : delt = potT - TsfK ! pot temp diff (K)
236 3558446561 : qsat = qqq * exp(-TTT/TsfK) ! saturation humidity (kg/m^3)
237 3558446561 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
238 :
239 3558446561 : thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K)
240 3558446561 : delq = Qa - ssq ! spec hum dif (kg/kg)
241 3558446561 : alzm = log(zlvl/zref)
242 3558446561 : if (present(zlvs)) then
243 1594049606 : alzs = log(zlvs/zref)
244 : else
245 1964396955 : alzs = alzm
246 : endif
247 3558446561 : cp = cp_air*(c1 + cpvir*ssq)
248 :
249 : !------------------------------------------------------------
250 : ! first estimate of Z/L and ustar, tstar and qstar
251 : !------------------------------------------------------------
252 :
253 : ! neutral coefficients, z/L = 0.0
254 3558446561 : rhn = rdn
255 3558446561 : ren = rdn
256 :
257 : ! ustar,tstar,qstar
258 3558446561 : ustar = rdn * vmag
259 3558446561 : tstar = rhn * delt
260 3558446561 : qstar = ren * delq
261 :
262 : !------------------------------------------------------------
263 : ! iterate to converge on Z/L, ustar, tstar and qstar
264 : !------------------------------------------------------------
265 :
266 3558446561 : ustar_prev = c2 * ustar
267 :
268 3558446561 : k = 1
269 21223544357 : do while (abs(ustar - ustar_prev)/ustar > atmiter_conv .and. k <= natmiter)
270 17665097796 : k = k + 1
271 17665097796 : ustar_prev = ustar
272 :
273 : ! compute stability & evaluate all stability functions
274 : holm = compute_stability_parameter(zlvl , thva , & ! LCOV_EXCL_LINE
275 : ustar, tstar, & ! LCOV_EXCL_LINE
276 17665097796 : qstar, Qa)
277 17665097796 : if (present(zlvs)) then
278 : hols = compute_stability_parameter(zlvs , thva , &
279 : ustar, tstar, & ! LCOV_EXCL_LINE
280 7843326626 : qstar, Qa)
281 : else
282 9821771170 : hols = holm
283 : endif
284 :
285 17665097796 : call compute_stability_function('momentum', holm, stable, psimh)
286 17665097796 : call compute_stability_function('scalar' , hols, stable, psixh)
287 :
288 : ! shift all coeffs to measurement height and stability
289 17665097796 : rd = rdn / (c1+rdn/vonkar*(alzm-psimh))
290 17665097796 : rh = rhn / (c1+rhn/vonkar*(alzs-psixh))
291 17665097796 : re = ren / (c1+ren/vonkar*(alzs-psixh))
292 :
293 : ! update ustar, tstar, qstar using updated, shifted coeffs
294 17665097796 : ustar = rd * vmag
295 17665097796 : tstar = rh * delt
296 17665097796 : qstar = re * delq
297 :
298 : enddo ! end iteration
299 :
300 3558446561 : if (calc_strair) then
301 :
302 : ! initialize
303 3558446561 : strx = c0
304 3558446561 : stry = c0
305 :
306 3558446561 : if (highfreq .and. sfctype(1:3)=='ice') then
307 :
308 : !------------------------------------------------------------
309 : ! momentum flux for high frequency coupling (RASM/CESM)
310 : !------------------------------------------------------------
311 : ! tau = rhoa * rd * rd
312 : ! strx = tau * |Uatm-U| * (uatm-u)
313 : ! stry = tau * |Uatm-U| * (vatm-v)
314 : !------------------------------------------------------------
315 :
316 70577610 : tau = rhoa * rd * rd ! not the stress at zlvl
317 :
318 : ! high frequency momentum coupling following Roberts et al. (2014)
319 70577610 : strx = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (uatm-uvel)
320 70577610 : stry = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (vatm-vvel)
321 :
322 : else
323 :
324 : !------------------------------------------------------------
325 : ! momentum flux
326 : !------------------------------------------------------------
327 : ! tau = rhoa * ustar * ustar
328 : ! strx = tau * uatm / vmag
329 : ! stry = tau * vatm / vmag
330 : !------------------------------------------------------------
331 :
332 3487868951 : tau = rhoa * ustar * rd ! not the stress at zlvl
333 3487868951 : strx = tau * uatm
334 3487868951 : stry = tau * vatm
335 :
336 : endif
337 :
338 3558446561 : Cdn_atm_ratio_n = rd * rd / rdn / rdn
339 :
340 : endif ! calc_strair
341 :
342 : !------------------------------------------------------------
343 : ! coefficients for turbulent flux calculation
344 : !------------------------------------------------------------
345 : ! add windless coefficient for sensible heat flux
346 : ! as in Jordan et al (JGR, 1999)
347 : !------------------------------------------------------------
348 :
349 3558446561 : if (trim(atmbndy) == 'mixed') then
350 : !- Use constant coefficients for sensible and latent heat fluxes
351 : ! similar to atmo_boundary_const but using vmag instead of wind
352 6369800 : shcoef = senscoef*cp_air*rhoa*vmag
353 6369800 : lhcoef = latncoef*Lheat *rhoa*vmag
354 : else ! 'similarity'
355 : !- Monin-Obukhov similarity theory for boundary layer
356 3552076761 : shcoef = rhoa * ustar * cp * rh + c1
357 3552076761 : lhcoef = rhoa * ustar * Lheat * re
358 : endif
359 :
360 : !------------------------------------------------------------
361 : ! Compute diagnostics: 2m ref T, Q, U
362 : !------------------------------------------------------------
363 :
364 3558446561 : hols = hols*zTrf/zlvl
365 3558446561 : psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols)
366 : fac = (rh/vonkar) &
367 3558446561 : * (alzs + al2 - psixh + psix2)
368 3558446561 : Tref = potT - delt*fac
369 3558446561 : Tref = Tref - p01*zTrf ! pot temp to temp correction
370 : fac = (re/vonkar) &
371 3558446561 : * (alzs + al2 - psixh + psix2)
372 3558446561 : Qref = Qa - delq*fac
373 :
374 3558446561 : if (highfreq .and. sfctype(1:3)=='ice') then
375 70577610 : Uref = sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * rd / rdn
376 : else
377 3487868951 : Uref = vmag * rd / rdn
378 : endif
379 :
380 3558446561 : if (tr_iso .and. sfctype(1:3)=='ice') then
381 63097468 : Qref_iso(:) = c0
382 63097468 : do n = 1, n_iso
383 47323101 : ratio = c0
384 47323101 : if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2)
385 63097468 : Qref_iso(n) = Qa_iso(n) - ratio*delq*fac
386 : enddo
387 : endif
388 :
389 3558446561 : end subroutine atmo_boundary_layer
390 :
391 : !=======================================================================
392 :
393 : ! Compute coefficients for atm/ice fluxes, stress
394 : ! NOTE: \\
395 : ! (1) all fluxes are positive downward, \\
396 : ! (2) reference temperature and humidity are NOT computed
397 :
398 20039172 : subroutine atmo_boundary_const (sfctype, calc_strair, &
399 : uatm, vatm, & ! LCOV_EXCL_LINE
400 : wind, rhoa, & ! LCOV_EXCL_LINE
401 : strx, stry, & ! LCOV_EXCL_LINE
402 : Tsf, potT, & ! LCOV_EXCL_LINE
403 : Qa, & ! LCOV_EXCL_LINE
404 : delt, delq, & ! LCOV_EXCL_LINE
405 : lhcoef, shcoef )
406 :
407 : character (len=3), intent(in) :: &
408 : sfctype ! ice or ocean
409 :
410 : logical (kind=log_kind), intent(in) :: &
411 : calc_strair ! if true, calculate wind stress components
412 :
413 : real (kind=dbl_kind), intent(in) :: &
414 : Tsf , & ! surface temperature of ice or ocean ! LCOV_EXCL_LINE
415 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
416 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
417 : uatm , & ! x-direction wind speed (m/s) ! LCOV_EXCL_LINE
418 : vatm , & ! y-direction wind speed (m/s) ! LCOV_EXCL_LINE
419 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
420 : rhoa ! air density (kg/m^3)
421 :
422 : real (kind=dbl_kind), intent(inout):: &
423 : strx , & ! x surface stress (N) ! LCOV_EXCL_LINE
424 : stry ! y surface stress (N)
425 :
426 : real (kind=dbl_kind), intent(out):: &
427 : delt , & ! potential T difference (K) ! LCOV_EXCL_LINE
428 : delq , & ! humidity difference (kg/kg) ! LCOV_EXCL_LINE
429 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
430 : lhcoef ! transfer coefficient for latent heat
431 :
432 : ! local variables
433 :
434 : real (kind=dbl_kind) :: &
435 : TsfK, & ! surface temperature in Kelvin (K) ! LCOV_EXCL_LINE
436 : qsat, & ! the saturation humidity of air (kg/m^3) ! LCOV_EXCL_LINE
437 : ssq , & ! sat surface humidity (kg/kg) ! LCOV_EXCL_LINE
438 : tau, & ! stress at zlvl ! LCOV_EXCL_LINE
439 : Lheat ! Lvap or Lsub, depending on surface type
440 :
441 : character(len=*),parameter :: subname='(atmo_boundary_const)'
442 :
443 : !------------------------------------------------------------
444 : ! Initialize
445 : !------------------------------------------------------------
446 :
447 20039172 : delt = c0
448 20039172 : delq = c0
449 20039172 : shcoef = c0
450 20039172 : lhcoef = c0
451 :
452 20039172 : if (calc_strair) then
453 :
454 6367836 : strx = c0
455 6367836 : stry = c0
456 :
457 : !------------------------------------------------------------
458 : ! momentum flux
459 : !------------------------------------------------------------
460 6367836 : tau = rhoa * 0.0012_dbl_kind * wind
461 : !AOMIP tau = rhoa * (1.10_dbl_kind + c4*p01*wind) &
462 : !AOMIP * wind * p001
463 6367836 : strx = tau * uatm
464 6367836 : stry = tau * vatm
465 :
466 : endif ! calc_strair
467 :
468 : !------------------------------------------------------------
469 : ! define variables that depend on surface type
470 : !------------------------------------------------------------
471 :
472 20039172 : if (sfctype(1:3)=='ice') then
473 3082116 : Lheat = Lsub ! ice to vapor
474 16957056 : elseif (sfctype(1:3)=='ocn') then
475 16957056 : Lheat = Lvap ! liquid to vapor
476 : endif ! sfctype
477 :
478 : !------------------------------------------------------------
479 : ! potential temperature and specific humidity differences
480 : !------------------------------------------------------------
481 :
482 20039172 : TsfK = Tsf + Tffresh ! surface temp (K)
483 20039172 : qsat = qqqocn * exp(-TTTocn/TsfK) ! sat humidity (kg/m^3)
484 20039172 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
485 :
486 20039172 : delt= potT - TsfK ! pot temp diff (K)
487 20039172 : delq= Qa - ssq ! spec hum dif (kg/kg)
488 :
489 : !------------------------------------------------------------
490 : ! coefficients for turbulent flux calculation
491 : !------------------------------------------------------------
492 :
493 20039172 : shcoef = senscoef*cp_air*rhoa*wind
494 20039172 : lhcoef = latncoef*Lheat *rhoa*wind
495 :
496 20039172 : end subroutine atmo_boundary_const
497 :
498 : !=======================================================================
499 :
500 : ! Neutral drag coefficients for ocean and atmosphere also compute the
501 : ! intermediate necessary variables ridge height, distance, floe size
502 : ! based upon Tsamados et al. (2014), JPO, DOI: 10.1175/JPO-D-13-0215.1.
503 : ! Places where the code varies from the paper are commented.
504 : !
505 : ! authors: Michel Tsamados, CPOM
506 : ! David Schroeder, CPOM
507 : !
508 : ! changes: Andrew Roberts, NPS (RASM/CESM coupling and documentation)
509 :
510 14602944 : subroutine neutral_drag_coeffs (apondn, &
511 14602944 : alvl, vlvl, & ! LCOV_EXCL_LINE
512 : aice, vice, & ! LCOV_EXCL_LINE
513 14602944 : vsno, aicen, & ! LCOV_EXCL_LINE
514 14602944 : vicen, & ! LCOV_EXCL_LINE
515 : Cdn_ocn, Cdn_ocn_skin, & ! LCOV_EXCL_LINE
516 : Cdn_ocn_floe, Cdn_ocn_keel,& ! LCOV_EXCL_LINE
517 : Cdn_atm, Cdn_atm_skin, & ! LCOV_EXCL_LINE
518 : Cdn_atm_floe, Cdn_atm_pond,& ! LCOV_EXCL_LINE
519 : Cdn_atm_rdg, hfreebd, & ! LCOV_EXCL_LINE
520 : hdraft, hridge, & ! LCOV_EXCL_LINE
521 : distrdg, hkeel, & ! LCOV_EXCL_LINE
522 : dkeel, lfloe, & ! LCOV_EXCL_LINE
523 : dfloe)
524 :
525 : use icepack_tracers, only: tr_pond
526 :
527 : real (kind=dbl_kind), dimension (:), intent(in) :: &
528 : apondn ,& ! melt pond fraction of sea ice category ! LCOV_EXCL_LINE
529 : alvl ,& ! level ice area fraction (of grid cell ?) ! LCOV_EXCL_LINE
530 : vlvl ! level ice mean thickness
531 :
532 : real (kind=dbl_kind), intent(in) :: &
533 : aice , & ! concentration of ice ! LCOV_EXCL_LINE
534 : vice , & ! volume per unit area of ice ! LCOV_EXCL_LINE
535 : vsno ! volume per unit area of snow
536 :
537 : real (kind=dbl_kind), dimension (:), intent(in) :: &
538 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
539 : vicen ! volume per unit area of ice (m)
540 :
541 : real (kind=dbl_kind), intent(out) :: &
542 : hfreebd , & ! freeboard (m) ! LCOV_EXCL_LINE
543 : hdraft , & ! draught of ice + snow column (Stoessel1993) ! LCOV_EXCL_LINE
544 : hridge , & ! ridge height ! LCOV_EXCL_LINE
545 : distrdg , & ! distance between ridges ! LCOV_EXCL_LINE
546 : hkeel , & ! keel depth ! LCOV_EXCL_LINE
547 : dkeel , & ! distance between keels ! LCOV_EXCL_LINE
548 : lfloe , & ! floe length (m) ! LCOV_EXCL_LINE
549 : dfloe , & ! distance between floes ! LCOV_EXCL_LINE
550 : Cdn_ocn , & ! ocean-ice neutral drag coefficient ! LCOV_EXCL_LINE
551 : Cdn_ocn_skin , & ! drag coefficient due to skin drag ! LCOV_EXCL_LINE
552 : Cdn_ocn_floe , & ! drag coefficient due to floe edges ! LCOV_EXCL_LINE
553 : Cdn_ocn_keel , & ! drag coefficient due to keels ! LCOV_EXCL_LINE
554 : Cdn_atm , & ! ice-atmosphere drag coefficient ! LCOV_EXCL_LINE
555 : Cdn_atm_skin , & ! drag coefficient due to skin drag ! LCOV_EXCL_LINE
556 : Cdn_atm_floe , & ! drag coefficient due to floe edges ! LCOV_EXCL_LINE
557 : Cdn_atm_pond , & ! drag coefficient due to ponds ! LCOV_EXCL_LINE
558 : Cdn_atm_rdg ! drag coefficient due to ridges
559 :
560 : real (kind=dbl_kind), parameter :: &
561 : ! [,] = range of values that can be tested
562 : csw = 0.002_dbl_kind ,&! ice-ocn drag coefficient [0.0005,0.005]
563 : csa = 0.0005_dbl_kind,&! ice-air drag coefficient [0.0001,0.001] ! LCOV_EXCL_LINE
564 : mrdg = c20 ,&! screening effect see Lu2011 [5,50] ! LCOV_EXCL_LINE
565 : mrdgo = c10 ,&! screening effect see Lu2011 [5,50] ! LCOV_EXCL_LINE
566 : beta = p5 ,&! power exponent appearing in astar and ! LCOV_EXCL_LINE
567 : ! L=Lmin(A*/(A*-A))**beta [0,1]
568 : Lmin = c8 ,&! min length of floe (m) [5,100]
569 : Lmax = 300._dbl_kind ,&! max length of floe (m) [30,3000] ! LCOV_EXCL_LINE
570 : cfa = p2 ,&! Eq. 12 ratio of local from drag over ! LCOV_EXCL_LINE
571 : ! geometrical parameter [0,1]
572 : cfw = p2 ,&! Eq. 15 ratio of local from drag over
573 : ! geometrical parameter [0,1]
574 : cpa = p2 ,&! Eq. 16 ratio of local form drag over
575 : ! geometrical parameter [0,1]
576 : cra = p2 ,&! Eq. 10 local form drag coefficient [0,1]
577 : crw = p2 ,&! Eq. 11 local form drag coefficient [0,1] ! LCOV_EXCL_LINE
578 : sl = 22._dbl_kind ,&! Sheltering parameter Lupkes2012 [10,30] ! LCOV_EXCL_LINE
579 : lpmin = 2.26_dbl_kind ,&! min pond length (m) see Eq. 17 [1,10] ! LCOV_EXCL_LINE
580 : lpmax = 24.63_dbl_kind ,&! max pond length (m) see Eq. 17 [10,100] ! LCOV_EXCL_LINE
581 : tanar = p4 ,&! 0.25 sail slope = 14 deg [0.4,1] ! LCOV_EXCL_LINE
582 : tanak = p4 ,&! 0.58 keel slope = 30 deg [0.4,1] ! LCOV_EXCL_LINE
583 : phir = 0.8_dbl_kind ,&! porosity of ridges [0.4,1] ! LCOV_EXCL_LINE
584 : phik = 0.8_dbl_kind ,&! porosity of keels [0.4,1] ! LCOV_EXCL_LINE
585 : hkoverhr = c4 ,&! hkeel/hridge ratio [4,8] ! LCOV_EXCL_LINE
586 : dkoverdr = c1 ,&! dkeel/distrdg ratio [1,5] ! LCOV_EXCL_LINE
587 : sHGB = 0.18_dbl_kind ,&! Lupkes2012 Eq. 28, Hanssen1988, ! LCOV_EXCL_LINE
588 : ! Steele1989 suggest instead 0.18
589 : alpha2 = c0 ,&! weight functions for area of
590 : beta2 = p75 ! ridged ice [0,1]
591 :
592 : integer (kind=int_kind) :: &
593 : n ! category index
594 :
595 : real (kind=dbl_kind) :: &
596 : astar, & ! new constant for form drag ! LCOV_EXCL_LINE
597 : ctecaf, & ! constante ! LCOV_EXCL_LINE
598 : ctecwf, & ! constante ! LCOV_EXCL_LINE
599 : sca, & ! wind attenuation function ! LCOV_EXCL_LINE
600 : scw, & ! ocean attenuation function ! LCOV_EXCL_LINE
601 : lp, & ! pond length (m) ! LCOV_EXCL_LINE
602 : ctecar, & ! LCOV_EXCL_LINE
603 : ctecwk, & ! LCOV_EXCL_LINE
604 : ai, aii, & ! ice area and its inverse ! LCOV_EXCL_LINE
605 : ocnrufi, & ! inverse ocean roughness ! LCOV_EXCL_LINE
606 : icerufi, & ! inverse ice roughness ! LCOV_EXCL_LINE
607 : tmp1 ! temporary
608 :
609 : real (kind=dbl_kind) :: &
610 : apond , & ! melt pond fraction of grid cell ! LCOV_EXCL_LINE
611 : ardg , & ! ridged ice area fraction of grid cell ! LCOV_EXCL_LINE
612 : vrdg ! ridged ice mean thickness
613 :
614 : real (kind=dbl_kind), parameter :: &
615 : ocnruf = 0.000327_dbl_kind ! ocean surface roughness (m)
616 :
617 : real (kind=dbl_kind), parameter :: &
618 : camax = 0.02_dbl_kind , & ! Maximum for atmospheric drag ! LCOV_EXCL_LINE
619 : cwmax = 0.06_dbl_kind ! Maximum for ocean drag
620 :
621 : character(len=*),parameter :: subname='(neutral_drag_coeffs)'
622 :
623 14602944 : astar = c1/(c1-(Lmin/Lmax)**(c1/beta))
624 :
625 : !-----------------------------------------------------------------
626 : ! Initialize across entire grid
627 : !-----------------------------------------------------------------
628 :
629 14602944 : ocnrufi = c1/ocnruf ! inverse ocean roughness
630 14602944 : icerufi = c1/iceruf ! inverse ice roughness
631 14602944 : hfreebd=c0
632 14602944 : hdraft =c0
633 14602944 : hridge =c0
634 14602944 : distrdg=c0
635 14602944 : hkeel =c0
636 14602944 : dkeel =c0
637 14602944 : lfloe =c0
638 14602944 : dfloe =c0
639 14602944 : Cdn_ocn=dragio
640 14602944 : Cdn_ocn_skin=c0
641 14602944 : Cdn_ocn_floe=c0
642 14602944 : Cdn_ocn_keel=c0
643 14602944 : Cdn_atm = (vonkar/log(zref/iceruf)) * (vonkar/log(zref/iceruf))
644 14602944 : Cdn_atm_skin=c0
645 14602944 : Cdn_atm_floe=c0
646 14602944 : Cdn_atm_pond=c0
647 14602944 : Cdn_atm_rdg =c0
648 :
649 14602944 : if (aice > p001) then
650 :
651 3324194 : Cdn_atm_skin = csa
652 3324194 : Cdn_ocn_skin = csw
653 :
654 3324194 : ai = aice
655 3324194 : aii = c1/ai
656 :
657 : !------------------------------------------------------------
658 : ! Compute average quantities
659 : !------------------------------------------------------------
660 :
661 : ! ponds
662 3324194 : apond = c0
663 3324194 : if (tr_pond) then
664 19945164 : do n = 1,ncat
665 : ! area of pond per unit area of grid cell
666 19945164 : apond = apond+apondn(n)*aicen(n)
667 : enddo
668 : endif
669 :
670 : ! draft and freeboard (see Eq. 27)
671 3324194 : hdraft = (rhoi*vice+rhos*vsno)*aii/rhow ! without ponds
672 3324194 : hfreebd = (vice+vsno)*aii-hdraft
673 :
674 : ! Do not allow draft larger than ice thickness (see Eq. 28)
675 3324194 : if (hdraft >= vice*aii) then
676 : ! replace excess snow with ice so hi~=hdraft
677 : hfreebd = (hdraft*ai*(c1-rhoi/rhow) + &
678 : (vsno-(vice-hdraft*ai)*rhoi/rhos) * & ! LCOV_EXCL_LINE
679 336 : (c1-rhos/rhow))*aii ! Stoessel1993
680 : endif
681 :
682 : ! floe size parameterization see Eq. 13
683 3324194 : lfloe = Lmin * (astar / (astar - ai))**beta
684 :
685 : ! distance between floes parameterization see Eq. 14
686 3324194 : dfloe = lfloe * (c1/sqrt(ai) - c1)
687 :
688 : ! Relate ridge height and distance between ridges to
689 : ! ridged ice area fraction and ridged ice mean thickness
690 : ! Assumes total volume of ridged ice is split into ridges and keels.
691 : ! Then assume total ridges volume over total area of ridges =
692 : ! volume of one average ridge / area of one average ridge
693 : ! Same for keels.
694 :
695 3324194 : ardg=c0
696 3324194 : vrdg=c0
697 19945164 : do n=1,ncat
698 : ! ridged ice area fraction over grid cell
699 16620970 : ardg=ardg+(c1-alvl(n))*aicen(n)
700 : ! total ridged ice volume per unit grid cell area
701 19945164 : vrdg=vrdg+(c1-vlvl(n))*vicen(n)
702 : enddo
703 :
704 : ! hridge, hkeel, distrdg and dkeel estimates from CICE for
705 : ! simple triangular geometry
706 3324194 : if (ardg > p001) then
707 : ! see Eq. 25 and Eq. 26
708 : hridge = vrdg/ardg*c2 &
709 : * (alpha2+beta2*hkoverhr/dkoverdr*tanar/tanak) & ! LCOV_EXCL_LINE
710 966669 : / (phir*c1+phik*tanar/tanak*hkoverhr**c2/dkoverdr)
711 : distrdg = c2*hridge*ai/ardg &
712 966669 : * (alpha2/tanar+beta2/tanak*hkoverhr/dkoverdr)
713 966669 : hkeel = hkoverhr * hridge
714 966669 : dkeel = dkoverdr * distrdg
715 :
716 : ! Use the height of ridges relative to the mean freeboard of
717 : ! the pack. Therefore skin drag and ridge drag differ in
718 : ! this code as compared to Tsamados et al. (2014) equations
719 : ! 10 and 18, which reference both to sea level.
720 966669 : tmp1 = max(c0,hridge - hfreebd)
721 :
722 : !------------------------------------------------------------
723 : ! Skin drag (atmo)
724 : !------------------------------------------------------------
725 :
726 966669 : Cdn_atm_skin = csa*(c1 - mrdg*tmp1/distrdg)
727 966669 : Cdn_atm_skin = max(min(Cdn_atm_skin,camax),c0)
728 :
729 : !------------------------------------------------------------
730 : ! Ridge effect (atmo)
731 : !------------------------------------------------------------
732 :
733 966669 : if (tmp1 > puny) then
734 966669 : sca = c1 - exp(-sHGB*distrdg/tmp1) ! see Eq. 9
735 966669 : ctecar = cra*p5
736 : Cdn_atm_rdg = ctecar*tmp1/distrdg*sca* &
737 966669 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
738 966669 : Cdn_atm_rdg = min(Cdn_atm_rdg,camax)
739 : endif
740 :
741 : ! Use the depth of keels relative to the mean draft of
742 : ! the pack. Therefore skin drag and keel drag differ in
743 : ! this code as compared to Tsamados et al. (2014) equations
744 : ! 11 and 19, which reference both to sea level. In some
745 : ! circumstances, hkeel can be less than hdraft because hkoverhr
746 : ! is constant, and max(c0,...) temporarily addresses this.
747 966669 : tmp1 = max(c0,hkeel - hdraft)
748 :
749 : !------------------------------------------------------------
750 : ! Skin drag bottom ice (ocean)
751 : !------------------------------------------------------------
752 :
753 966669 : Cdn_ocn_skin = csw * (c1 - mrdgo*tmp1/dkeel)
754 966669 : Cdn_ocn_skin = max(min(Cdn_ocn_skin,cwmax), c0)
755 :
756 : !------------------------------------------------------------
757 : ! Keel effect (ocean)
758 : !------------------------------------------------------------
759 :
760 966669 : if (tmp1 > puny) then
761 966669 : scw = c1 - exp(-sHGB*dkeel/tmp1)
762 966669 : ctecwk = crw*p5
763 : Cdn_ocn_keel = ctecwk*tmp1/dkeel*scw* &
764 966669 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
765 966669 : Cdn_ocn_keel = max(min(Cdn_ocn_keel,cwmax),c0)
766 : endif
767 :
768 : endif ! ardg > 0.001
769 :
770 : !------------------------------------------------------------
771 : ! Floe edge drag effect (atmo)
772 : !------------------------------------------------------------
773 :
774 3324194 : if (hfreebd > puny) then
775 3324194 : sca = c1 - exp(-sl*beta*(c1-ai))
776 3324194 : ctecaf = cfa*p5*(log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2*sca
777 3324194 : Cdn_atm_floe = ctecaf * hfreebd / lfloe
778 3324194 : Cdn_atm_floe = max(min(Cdn_atm_floe,camax),c0)
779 : endif
780 :
781 : !------------------------------------------------------------
782 : ! Pond edge effect (atmo)
783 : !------------------------------------------------------------
784 :
785 3324194 : if (hfreebd > puny) then
786 3324194 : sca = (apond)**(c1/(zref*beta))
787 3324194 : lp = lpmin*(1-apond)+lpmax*apond
788 : Cdn_atm_pond = cpa*p5*sca*apond*hfreebd/lp &
789 3324194 : * (log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2
790 3324194 : Cdn_atm_pond = min(Cdn_atm_pond,camax)
791 : endif
792 :
793 : !------------------------------------------------------------
794 : ! Floe edge drag effect (ocean)
795 : !------------------------------------------------------------
796 :
797 3324194 : if (hdraft > puny) then
798 3324194 : scw = c1 - exp(-sl*beta*(c1-ai))
799 3324194 : ctecwf = cfw*p5*(log(hdraft*ocnrufi)/log(zref*ocnrufi))**c2*scw
800 3324194 : Cdn_ocn_floe = ctecwf * hdraft / lfloe
801 3324194 : Cdn_ocn_floe = max(min(Cdn_ocn_floe,cwmax),c0)
802 : endif
803 :
804 : !------------------------------------------------------------
805 : ! Total drag coefficient (atmo)
806 : !------------------------------------------------------------
807 :
808 3324194 : Cdn_atm = Cdn_atm_skin + Cdn_atm_floe + Cdn_atm_pond + Cdn_atm_rdg
809 3324194 : Cdn_atm = min(Cdn_atm,camax)
810 :
811 : !------------------------------------------------------------
812 : ! Total drag coefficient (ocean)
813 : !------------------------------------------------------------
814 :
815 3324194 : Cdn_ocn = Cdn_ocn_skin + Cdn_ocn_floe + Cdn_ocn_keel
816 3324194 : Cdn_ocn = min(Cdn_ocn,cwmax)
817 :
818 : endif
819 :
820 14602944 : end subroutine neutral_drag_coeffs
821 :
822 : !=======================================================================
823 : !autodocument_start icepack_atm_boundary
824 : !
825 :
826 3578485733 : subroutine icepack_atm_boundary(sfctype, &
827 : Tsf, potT, & ! LCOV_EXCL_LINE
828 : uatm, vatm, & ! LCOV_EXCL_LINE
829 : wind, zlvl, & ! LCOV_EXCL_LINE
830 : Qa, rhoa, & ! LCOV_EXCL_LINE
831 : strx, stry, & ! LCOV_EXCL_LINE
832 : Tref, Qref, & ! LCOV_EXCL_LINE
833 : delt, delq, & ! LCOV_EXCL_LINE
834 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
835 : Cdn_atm, & ! LCOV_EXCL_LINE
836 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
837 3578485733 : Qa_iso, Qref_iso, & ! LCOV_EXCL_LINE
838 : uvel, vvel, & ! LCOV_EXCL_LINE
839 : Uref, zlvs)
840 :
841 : character (len=3), intent(in) :: &
842 : sfctype ! ice or ocean
843 :
844 : real (kind=dbl_kind), intent(in) :: &
845 : Tsf , & ! surface temperature of ice or ocean ! LCOV_EXCL_LINE
846 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
847 : uatm , & ! x-direction wind speed (m/s) ! LCOV_EXCL_LINE
848 : vatm , & ! y-direction wind speed (m/s) ! LCOV_EXCL_LINE
849 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
850 : zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) ! LCOV_EXCL_LINE
851 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
852 : rhoa ! air density (kg/m^3)
853 :
854 : real (kind=dbl_kind), intent(inout) :: &
855 : Cdn_atm , & ! neutral drag coefficient ! LCOV_EXCL_LINE
856 : Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff
857 :
858 : real (kind=dbl_kind), intent(inout) :: &
859 : strx , & ! x surface stress (N) ! LCOV_EXCL_LINE
860 : stry ! y surface stress (N)
861 :
862 : real (kind=dbl_kind), intent(inout) :: &
863 : Tref , & ! reference height temperature (K) ! LCOV_EXCL_LINE
864 : Qref , & ! reference height specific humidity (kg/kg) ! LCOV_EXCL_LINE
865 : delt , & ! potential T difference (K) ! LCOV_EXCL_LINE
866 : delq , & ! humidity difference (kg/kg) ! LCOV_EXCL_LINE
867 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
868 : lhcoef ! transfer coefficient for latent heat
869 :
870 : real (kind=dbl_kind), intent(in), dimension(:), optional :: &
871 : Qa_iso ! specific isotopic humidity (kg/kg)
872 :
873 : real (kind=dbl_kind), intent(inout), dimension(:), optional :: &
874 : Qref_iso ! reference specific isotopic humidity (kg/kg)
875 :
876 : real (kind=dbl_kind), intent(in), optional :: &
877 : uvel , & ! x-direction ice speed (m/s) ! LCOV_EXCL_LINE
878 : vvel , & ! y-direction ice speed (m/s) ! LCOV_EXCL_LINE
879 : zlvs ! atm level height for scalars (if different than zlvl) (m)
880 :
881 : real (kind=dbl_kind), intent(out), optional :: &
882 : Uref ! reference height wind speed (m/s)
883 :
884 : !autodocument_end
885 :
886 : ! local variables
887 :
888 : real (kind=dbl_kind) :: &
889 : l_uvel, l_vvel, l_Uref
890 :
891 : logical (kind=log_kind), save :: &
892 : first_call_ice = .true. ! first call flag
893 :
894 : character(len=*),parameter :: subname='(icepack_atm_boundary)'
895 :
896 : !------------------------------------------------------------
897 : ! Check optional arguments
898 : ! Need separate first_call flags for 'ice' and 'ocn' sfctype
899 : !------------------------------------------------------------
900 :
901 3578485733 : if (sfctype == 'ice') then
902 1599186245 : if (icepack_chkoptargflag(first_call_ice)) then
903 3968 : if (tr_iso) then
904 138 : if (.not.(present(Qa_iso).and.present(Qref_iso))) then
905 0 : call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T')
906 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
907 0 : return
908 : endif
909 : endif
910 : endif
911 : endif
912 :
913 3578485733 : l_uvel = c0
914 3578485733 : l_vvel = c0
915 3578485733 : l_Uref = c0
916 3578485733 : if (present(uvel)) then
917 1599186245 : l_uvel = uvel
918 : endif
919 3578485733 : if (present(vvel)) then
920 1599186245 : l_vvel = vvel
921 : endif
922 :
923 3578485733 : Cdn_atm_ratio_n = c1
924 :
925 3578485733 : if (trim(atmbndy) == 'constant') then
926 : call atmo_boundary_const (sfctype, calc_strair, &
927 : uatm, vatm, & ! LCOV_EXCL_LINE
928 : wind, rhoa, & ! LCOV_EXCL_LINE
929 : strx, stry, & ! LCOV_EXCL_LINE
930 : Tsf, potT, & ! LCOV_EXCL_LINE
931 : Qa, & ! LCOV_EXCL_LINE
932 : delt, delq, & ! LCOV_EXCL_LINE
933 20039172 : lhcoef, shcoef )
934 20039172 : if (icepack_warnings_aborted(subname)) return
935 : else
936 : call atmo_boundary_layer (sfctype, &
937 : calc_strair, formdrag, & ! LCOV_EXCL_LINE
938 : Tsf, potT, & ! LCOV_EXCL_LINE
939 : uatm, vatm, & ! LCOV_EXCL_LINE
940 : wind, zlvl, & ! LCOV_EXCL_LINE
941 : Qa, rhoa, & ! LCOV_EXCL_LINE
942 : strx, stry, & ! LCOV_EXCL_LINE
943 : Tref, Qref, & ! LCOV_EXCL_LINE
944 : delt, delq, & ! LCOV_EXCL_LINE
945 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
946 : Cdn_atm, & ! LCOV_EXCL_LINE
947 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
948 : Qa_iso=Qa_iso, & ! LCOV_EXCL_LINE
949 : Qref_iso=Qref_iso, & ! LCOV_EXCL_LINE
950 : uvel=l_uvel, vvel=l_vvel,& ! LCOV_EXCL_LINE
951 3558446561 : Uref=l_Uref, zlvs=zlvs )
952 3558446561 : if (icepack_warnings_aborted(subname)) return
953 : endif ! atmbndy
954 :
955 3578485733 : if (present(Uref)) then
956 1599186245 : Uref = l_Uref
957 : endif
958 :
959 3578485733 : if (sfctype == 'ice') then
960 1599186245 : first_call_ice = .false.
961 : endif
962 :
963 : end subroutine icepack_atm_boundary
964 :
965 : !=======================================================================
966 :
967 25508424422 : function compute_stability_parameter(zlvl , thva , &
968 : ustar, tstar, & ! LCOV_EXCL_LINE
969 : qstar, Qa) & ! LCOV_EXCL_LINE
970 : result(hol)
971 :
972 : real (kind=dbl_kind), intent(in) :: &
973 : zlvl , & ! atm level height (m) ! LCOV_EXCL_LINE
974 : thva , & ! virtual temperature (K) ! LCOV_EXCL_LINE
975 : ustar , & ! turbulent scale for momentum ! LCOV_EXCL_LINE
976 : tstar , & ! turbulent scale for temperature ! LCOV_EXCL_LINE
977 : qstar , & ! turbulent scale for humidity ! LCOV_EXCL_LINE
978 : Qa ! specific humidity (kg/kg)
979 :
980 : real (kind=dbl_kind) :: &
981 : hol ! H (at zlvl) over L
982 :
983 : character(len=*),parameter :: subname='(compute_stability_parameter)'
984 :
985 : hol = vonkar * gravit * zlvl &
986 : * (tstar/thva & ! LCOV_EXCL_LINE
987 : + qstar/(c1/zvir+Qa)) & ! LCOV_EXCL_LINE
988 25508424422 : / ustar**2
989 25508424422 : hol = sign( min(abs(hol),c10), hol)
990 :
991 25508424422 : end function compute_stability_parameter
992 :
993 : !=======================================================================
994 :
995 35330195592 : subroutine compute_stability_function(qty, hol, stable, psi)
996 :
997 : character (len=*), intent(in) :: &
998 : qty ! 'momentum' or 'scalar'
999 :
1000 : real (kind=dbl_kind), intent(in) :: &
1001 : hol ! H over L
1002 :
1003 : real (kind=dbl_kind), intent(out) :: &
1004 : psi , & ! stability function at hol ! LCOV_EXCL_LINE
1005 : stable ! unit step function at hol
1006 :
1007 : ! local variables
1008 :
1009 : real (kind=dbl_kind) :: &
1010 : psi_stable , & ! stable stability funcion at hol ! LCOV_EXCL_LINE
1011 : psi_unstable ! unstable stability funcion at hol
1012 :
1013 : character(len=*),parameter :: subname='(compute_stability_function)'
1014 :
1015 35330195592 : stable = p5 + sign(p5 , hol)
1016 :
1017 : psi_stable = -(0.7_dbl_kind*hol &
1018 : + 0.75_dbl_kind*(hol-14.3_dbl_kind) & ! LCOV_EXCL_LINE
1019 35330195592 : * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
1020 :
1021 35330195592 : if(trim(qty) == 'momentum') then
1022 17665097796 : psi_unstable = psi_momentum_unstable(hol)
1023 17665097796 : elseif(trim(qty) == 'scalar') then
1024 17665097796 : psi_unstable = psi_scalar_unstable(hol)
1025 : else
1026 0 : call icepack_warnings_add(subname//' incorrect qty: ' // qty)
1027 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1028 : endif
1029 :
1030 35330195592 : psi = psi_stable*stable + (c1 - stable)*psi_unstable
1031 :
1032 35330195592 : end subroutine compute_stability_function
1033 :
1034 : !------------------------------------------------------------
1035 : ! Define functions
1036 : !------------------------------------------------------------
1037 :
1038 : !=======================================================================
1039 :
1040 17665097796 : real(kind=dbl_kind) function psi_momentum_unstable(hol)
1041 :
1042 : real(kind=dbl_kind), intent(in) :: hol
1043 :
1044 : real(kind=dbl_kind) :: xd
1045 :
1046 17665097796 : xd = capital_X(hol)
1047 :
1048 : psi_momentum_unstable = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
1049 17665097796 : - c2*atan(xd) + pih
1050 :
1051 17665097796 : end function psi_momentum_unstable
1052 :
1053 : !=======================================================================
1054 :
1055 21223544357 : real(kind=dbl_kind) function psi_scalar_unstable(hol)
1056 :
1057 : real(kind=dbl_kind), intent(in) :: hol
1058 :
1059 : real(kind=dbl_kind) :: xd
1060 :
1061 21223544357 : xd = capital_X(hol)
1062 :
1063 21223544357 : psi_scalar_unstable = c2 * log((c1 + xd*xd)/c2)
1064 :
1065 21223544357 : end function psi_scalar_unstable
1066 :
1067 : !=======================================================================
1068 :
1069 38888642153 : real(kind=dbl_kind) function capital_X(hol)
1070 :
1071 : real(kind=dbl_kind), intent(in) :: hol
1072 :
1073 38888642153 : capital_X = sqrt(max(sqrt(abs(c1 - c16*hol)) , c1))
1074 :
1075 38888642153 : end function capital_X
1076 :
1077 : !=======================================================================
1078 :
1079 : end module icepack_atmo
1080 :
1081 : !=======================================================================
|
