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: 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 46176939 : 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 : 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 6371954 : 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 6371954 : 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 46176939 : al2 = log(zref/zTrf)
169 :
170 : !------------------------------------------------------------
171 : ! Initialize
172 : !------------------------------------------------------------
173 :
174 46176939 : cpvir = cp_wv/cp_air-c1 ! defined as cp_wv/cp_air - 1.
175 :
176 46176939 : if (highfreq) then
177 0 : umin = p5 ! minumum allowable wind-ice speed difference of 0.5 m/s
178 : else
179 46176939 : umin = c1 ! minumum allowable wind speed of 1m/s
180 : endif
181 :
182 46176939 : Tref = c0
183 46176939 : Qref = c0
184 46176939 : Uref = c0
185 46176939 : delt = c0
186 46176939 : delq = c0
187 46176939 : shcoef = c0
188 46176939 : 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 46176939 : if (sfctype(1:3)=='ice') then
200 :
201 33852099 : qqq = qqqice ! for qsat
202 33852099 : TTT = TTTice ! for qsat
203 33852099 : Lheat = Lsub ! ice to vapor
204 :
205 33852099 : if (highfreq) then
206 : vmag = max(umin, sqrt( (uatm-uvel)**2 + &
207 0 : (vatm-vvel)**2) )
208 : else
209 33852099 : vmag = max(umin, wind)
210 : endif
211 :
212 33852099 : if (formdrag .and. Cdn_atm > puny) then
213 0 : rdn = sqrt(Cdn_atm)
214 : else
215 33852099 : rdn = vonkar/log(zref/iceruf) ! neutral coefficient
216 33852099 : Cdn_atm = rdn * rdn
217 : endif
218 :
219 12324840 : elseif (sfctype(1:3)=='ocn') then
220 :
221 12324840 : qqq = qqqocn
222 12324840 : TTT = TTTocn
223 12324840 : Lheat = Lvap ! liquid to vapor
224 12324840 : vmag = max(umin, wind)
225 : rdn = sqrt(0.0027_dbl_kind/vmag &
226 12324840 : + .000142_dbl_kind + .0000764_dbl_kind*vmag)
227 :
228 : endif ! sfctype
229 :
230 : !------------------------------------------------------------
231 : ! define some more needed variables
232 : !------------------------------------------------------------
233 :
234 46176939 : TsfK = Tsf + Tffresh ! surface temp (K)
235 46176939 : delt = potT - TsfK ! pot temp diff (K)
236 46176939 : qsat = qqq * exp(-TTT/TsfK) ! saturation humidity (kg/m^3)
237 46176939 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
238 :
239 46176939 : thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K)
240 46176939 : delq = Qa - ssq ! spec hum dif (kg/kg)
241 46176939 : alzm = log(zlvl/zref)
242 46176939 : if (present(zlvs)) then
243 33852099 : alzs = log(zlvs/zref)
244 : else
245 12324840 : alzs = alzm
246 : endif
247 46176939 : 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 46176939 : rhn = rdn
255 46176939 : ren = rdn
256 :
257 : ! ustar,tstar,qstar
258 46176939 : ustar = rdn * vmag
259 46176939 : tstar = rhn * delt
260 46176939 : qstar = ren * delq
261 :
262 : !------------------------------------------------------------
263 : ! iterate to converge on Z/L, ustar, tstar and qstar
264 : !------------------------------------------------------------
265 :
266 46176939 : ustar_prev = c2 * ustar
267 :
268 46176939 : k = 1
269 276153750 : do while (abs(ustar - ustar_prev)/ustar > atmiter_conv .and. k <= natmiter)
270 229976811 : k = k + 1
271 229976811 : 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 229976811 : qstar, Qa)
277 229976811 : if (present(zlvs)) then
278 : hols = compute_stability_parameter(zlvs , thva , &
279 : ustar, tstar, & ! LCOV_EXCL_LINE
280 168580187 : qstar, Qa)
281 : else
282 61396624 : hols = holm
283 : endif
284 :
285 229976811 : call compute_stability_function('momentum', holm, stable, psimh)
286 229976811 : call compute_stability_function('scalar' , hols, stable, psixh)
287 :
288 : ! shift all coeffs to measurement height and stability
289 229976811 : rd = rdn / (c1+rdn/vonkar*(alzm-psimh))
290 229976811 : rh = rhn / (c1+rhn/vonkar*(alzs-psixh))
291 229976811 : re = ren / (c1+ren/vonkar*(alzs-psixh))
292 :
293 : ! update ustar, tstar, qstar using updated, shifted coeffs
294 229976811 : ustar = rd * vmag
295 229976811 : tstar = rh * delt
296 229976811 : qstar = re * delq
297 :
298 : enddo ! end iteration
299 :
300 46176939 : if (calc_strair) then
301 :
302 : ! initialize
303 46176939 : strx = c0
304 46176939 : stry = c0
305 :
306 46176939 : 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 0 : tau = rhoa * rd * rd ! not the stress at zlvl
317 :
318 : ! high frequency momentum coupling following Roberts et al. (2014)
319 0 : strx = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (uatm-uvel)
320 0 : 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 46176939 : tau = rhoa * ustar * rd ! not the stress at zlvl
333 46176939 : strx = tau * uatm
334 46176939 : stry = tau * vatm
335 :
336 : endif
337 :
338 46176939 : 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 46176939 : 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 0 : shcoef = senscoef*cp_air*rhoa*vmag
353 0 : lhcoef = latncoef*Lheat *rhoa*vmag
354 : else ! 'similarity'
355 : !- Monin-Obukhov similarity theory for boundary layer
356 46176939 : shcoef = rhoa * ustar * cp * rh + c1
357 46176939 : lhcoef = rhoa * ustar * Lheat * re
358 : endif
359 :
360 : !------------------------------------------------------------
361 : ! Compute diagnostics: 2m ref T, Q, U
362 : !------------------------------------------------------------
363 :
364 46176939 : hols = hols*zTrf/zlvl
365 46176939 : psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols)
366 : fac = (rh/vonkar) &
367 46176939 : * (alzs + al2 - psixh + psix2)
368 46176939 : Tref = potT - delt*fac
369 46176939 : Tref = Tref - p01*zTrf ! pot temp to temp correction
370 : fac = (re/vonkar) &
371 46176939 : * (alzs + al2 - psixh + psix2)
372 46176939 : Qref = Qa - delq*fac
373 :
374 46176939 : if (highfreq .and. sfctype(1:3)=='ice') then
375 0 : Uref = sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * rd / rdn
376 : else
377 46176939 : Uref = vmag * rd / rdn
378 : endif
379 :
380 46176939 : if (tr_iso .and. sfctype(1:3)=='ice') then
381 0 : Qref_iso(:) = c0
382 0 : do n = 1, n_iso
383 0 : ratio = c0
384 0 : if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2)
385 0 : Qref_iso(n) = Qa_iso(n) - ratio*delq*fac
386 : enddo
387 : endif
388 :
389 46176939 : 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 0 : 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 0 : Lheat ! Lvap or Lsub, depending on surface type
440 :
441 : character(len=*),parameter :: subname='(atmo_boundary_const)'
442 :
443 : !------------------------------------------------------------
444 : ! Initialize
445 : !------------------------------------------------------------
446 :
447 0 : delt = c0
448 0 : delq = c0
449 0 : shcoef = c0
450 0 : lhcoef = c0
451 :
452 0 : if (calc_strair) then
453 :
454 0 : strx = c0
455 0 : stry = c0
456 :
457 : !------------------------------------------------------------
458 : ! momentum flux
459 : !------------------------------------------------------------
460 0 : tau = rhoa * 0.0012_dbl_kind * wind
461 : !AOMIP tau = rhoa * (1.10_dbl_kind + c4*p01*wind) &
462 : !AOMIP * wind * p001
463 0 : strx = tau * uatm
464 0 : stry = tau * vatm
465 :
466 : endif ! calc_strair
467 :
468 : !------------------------------------------------------------
469 : ! define variables that depend on surface type
470 : !------------------------------------------------------------
471 :
472 0 : if (sfctype(1:3)=='ice') then
473 0 : Lheat = Lsub ! ice to vapor
474 0 : elseif (sfctype(1:3)=='ocn') then
475 0 : Lheat = Lvap ! liquid to vapor
476 : endif ! sfctype
477 :
478 : !------------------------------------------------------------
479 : ! potential temperature and specific humidity differences
480 : !------------------------------------------------------------
481 :
482 0 : TsfK = Tsf + Tffresh ! surface temp (K)
483 0 : qsat = qqqocn * exp(-TTTocn/TsfK) ! sat humidity (kg/m^3)
484 0 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
485 :
486 0 : delt= potT - TsfK ! pot temp diff (K)
487 0 : delq= Qa - ssq ! spec hum dif (kg/kg)
488 :
489 : !------------------------------------------------------------
490 : ! coefficients for turbulent flux calculation
491 : !------------------------------------------------------------
492 :
493 0 : shcoef = senscoef*cp_air*rhoa*wind
494 0 : lhcoef = latncoef*Lheat *rhoa*wind
495 :
496 0 : 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 0 : subroutine neutral_drag_coeffs (apnd, hpnd, &
511 : ipnd, & ! LCOV_EXCL_LINE
512 : alvl, vlvl, & ! LCOV_EXCL_LINE
513 : aice, vice, & ! LCOV_EXCL_LINE
514 : vsno, aicen, & ! LCOV_EXCL_LINE
515 : vicen, & ! LCOV_EXCL_LINE
516 : Cdn_ocn, Cdn_ocn_skin, & ! LCOV_EXCL_LINE
517 : Cdn_ocn_floe, Cdn_ocn_keel,& ! LCOV_EXCL_LINE
518 : Cdn_atm, Cdn_atm_skin, & ! LCOV_EXCL_LINE
519 : Cdn_atm_floe, Cdn_atm_pond,& ! LCOV_EXCL_LINE
520 : Cdn_atm_rdg, hfreebd, & ! LCOV_EXCL_LINE
521 : hdraft, hridge, & ! LCOV_EXCL_LINE
522 : distrdg, hkeel, & ! LCOV_EXCL_LINE
523 : dkeel, lfloe, & ! LCOV_EXCL_LINE
524 : dfloe, ncat)
525 :
526 : use icepack_tracers, only: tr_pond
527 :
528 : integer (kind=int_kind), intent(in) :: &
529 : ncat
530 :
531 : real (kind=dbl_kind), dimension (:), intent(in) :: &
532 : apnd ,& ! melt pond fraction of sea ice ! LCOV_EXCL_LINE
533 : hpnd ,& ! mean melt pond depth over sea ice ! LCOV_EXCL_LINE
534 : ipnd ,& ! mean ice pond depth over sea ice in cat n ! LCOV_EXCL_LINE
535 : alvl ,& ! level ice area fraction (of grid cell ?) ! LCOV_EXCL_LINE
536 : vlvl ! level ice mean thickness
537 :
538 : real (kind=dbl_kind), intent(in) :: &
539 : aice , & ! concentration of ice ! LCOV_EXCL_LINE
540 : vice , & ! volume per unit area of ice ! LCOV_EXCL_LINE
541 : vsno ! volume per unit area of snow
542 :
543 : real (kind=dbl_kind), dimension (:), intent(in) :: &
544 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
545 : vicen ! volume per unit area of ice (m)
546 :
547 : real (kind=dbl_kind), intent(out) :: &
548 : hfreebd , & ! freeboard (m) ! LCOV_EXCL_LINE
549 : hdraft , & ! draught of ice + snow column (Stoessel1993) ! LCOV_EXCL_LINE
550 : hridge , & ! ridge height ! LCOV_EXCL_LINE
551 : distrdg , & ! distance between ridges ! LCOV_EXCL_LINE
552 : hkeel , & ! keel depth ! LCOV_EXCL_LINE
553 : dkeel , & ! distance between keels ! LCOV_EXCL_LINE
554 : lfloe , & ! floe length (m) ! LCOV_EXCL_LINE
555 : dfloe , & ! distance between floes ! LCOV_EXCL_LINE
556 : Cdn_ocn , & ! ocean-ice neutral drag coefficient ! LCOV_EXCL_LINE
557 : Cdn_ocn_skin , & ! drag coefficient due to skin drag ! LCOV_EXCL_LINE
558 : Cdn_ocn_floe , & ! drag coefficient due to floe edges ! LCOV_EXCL_LINE
559 : Cdn_ocn_keel , & ! drag coefficient due to keels ! LCOV_EXCL_LINE
560 : Cdn_atm , & ! ice-atmosphere drag coefficient ! LCOV_EXCL_LINE
561 : Cdn_atm_skin , & ! drag coefficient due to skin drag ! LCOV_EXCL_LINE
562 : Cdn_atm_floe , & ! drag coefficient due to floe edges ! LCOV_EXCL_LINE
563 : Cdn_atm_pond , & ! drag coefficient due to ponds ! LCOV_EXCL_LINE
564 : Cdn_atm_rdg ! drag coefficient due to ridges
565 :
566 : real (kind=dbl_kind), parameter :: &
567 : ! [,] = range of values that can be tested
568 : csw = 0.002_dbl_kind ,&! ice-ocn drag coefficient [0.0005,0.005]
569 : csa = 0.0005_dbl_kind,&! ice-air drag coefficient [0.0001,0.001] ! LCOV_EXCL_LINE
570 : mrdg = c20 ,&! screening effect see Lu2011 [5,50] ! LCOV_EXCL_LINE
571 : mrdgo = c10 ,&! screening effect see Lu2011 [5,50] ! LCOV_EXCL_LINE
572 : beta = p5 ,&! power exponent appearing in astar and ! LCOV_EXCL_LINE
573 : ! L=Lmin(A*/(A*-A))**beta [0,1]
574 : Lmin = c8 ,&! min length of floe (m) [5,100]
575 : Lmax = 300._dbl_kind ,&! max length of floe (m) [30,3000] ! LCOV_EXCL_LINE
576 : cfa = p2 ,&! Eq. 12 ratio of local from drag over ! LCOV_EXCL_LINE
577 : ! geometrical parameter [0,1]
578 : cfw = p2 ,&! Eq. 15 ratio of local from drag over
579 : ! geometrical parameter [0,1]
580 : cpa = p2 ,&! Eq. 16 ratio of local form drag over
581 : ! geometrical parameter [0,1]
582 : cra = p2 ,&! Eq. 10 local form drag coefficient [0,1]
583 : crw = p2 ,&! Eq. 11 local form drag coefficient [0,1] ! LCOV_EXCL_LINE
584 : sl = 22._dbl_kind ,&! Sheltering parameter Lupkes2012 [10,30] ! LCOV_EXCL_LINE
585 : lpmin = 2.26_dbl_kind ,&! min pond length (m) see Eq. 17 [1,10] ! LCOV_EXCL_LINE
586 : lpmax = 24.63_dbl_kind ,&! max pond length (m) see Eq. 17 [10,100] ! LCOV_EXCL_LINE
587 : tanar = p4 ,&! 0.25 sail slope = 14 deg [0.4,1] ! LCOV_EXCL_LINE
588 : tanak = p4 ,&! 0.58 keel slope = 30 deg [0.4,1] ! LCOV_EXCL_LINE
589 : phir = 0.8_dbl_kind ,&! porosity of ridges [0.4,1] ! LCOV_EXCL_LINE
590 : phik = 0.8_dbl_kind ,&! porosity of keels [0.4,1] ! LCOV_EXCL_LINE
591 : hkoverhr = c4 ,&! hkeel/hridge ratio [4,8] ! LCOV_EXCL_LINE
592 : dkoverdr = c1 ,&! dkeel/distrdg ratio [1,5] ! LCOV_EXCL_LINE
593 : sHGB = 0.18_dbl_kind ,&! Lupkes2012 Eq. 28, Hanssen1988, ! LCOV_EXCL_LINE
594 : ! Steele1989 suggest instead 0.18
595 : alpha2 = c0 ,&! weight functions for area of
596 : beta2 = p75 ! ridged ice [0,1]
597 :
598 : integer (kind=int_kind) :: &
599 : n ! category index
600 :
601 : real (kind=dbl_kind) :: &
602 : astar, & ! new constant for form drag ! LCOV_EXCL_LINE
603 : ctecaf, & ! constante ! LCOV_EXCL_LINE
604 : ctecwf, & ! constante ! LCOV_EXCL_LINE
605 : sca, & ! wind attenuation function ! LCOV_EXCL_LINE
606 : scw, & ! ocean attenuation function ! LCOV_EXCL_LINE
607 : lp, & ! pond length (m) ! LCOV_EXCL_LINE
608 : ctecar, & ! LCOV_EXCL_LINE
609 : ctecwk, & ! LCOV_EXCL_LINE
610 : ai, aii, & ! ice area and its inverse ! LCOV_EXCL_LINE
611 : ocnrufi, & ! inverse ocean roughness ! LCOV_EXCL_LINE
612 : icerufi, & ! inverse ice roughness ! LCOV_EXCL_LINE
613 0 : tmp1 ! temporary
614 :
615 : real (kind=dbl_kind) :: &
616 : apond , & ! melt pond fraction of grid cell ! LCOV_EXCL_LINE
617 : ardg , & ! ridged ice area fraction of grid cell ! LCOV_EXCL_LINE
618 0 : vrdg ! ridged ice mean thickness
619 :
620 : real (kind=dbl_kind), parameter :: &
621 : ocnruf = 0.000327_dbl_kind ! ocean surface roughness (m)
622 :
623 : real (kind=dbl_kind), parameter :: &
624 : camax = 0.02_dbl_kind , & ! Maximum for atmospheric drag ! LCOV_EXCL_LINE
625 : cwmax = 0.06_dbl_kind ! Maximum for ocean drag
626 :
627 : character(len=*),parameter :: subname='(neutral_drag_coeffs)'
628 :
629 0 : astar = c1/(c1-(Lmin/Lmax)**(c1/beta))
630 :
631 : !-----------------------------------------------------------------
632 : ! Initialize across entire grid
633 : !-----------------------------------------------------------------
634 :
635 0 : ocnrufi = c1/ocnruf ! inverse ocean roughness
636 0 : icerufi = c1/iceruf ! inverse ice roughness
637 0 : hfreebd=c0
638 0 : hdraft =c0
639 0 : hridge =c0
640 0 : distrdg=c0
641 0 : hkeel =c0
642 0 : dkeel =c0
643 0 : lfloe =c0
644 0 : dfloe =c0
645 0 : Cdn_ocn=dragio
646 0 : Cdn_ocn_skin=c0
647 0 : Cdn_ocn_floe=c0
648 0 : Cdn_ocn_keel=c0
649 0 : Cdn_atm = (vonkar/log(zref/iceruf)) * (vonkar/log(zref/iceruf))
650 0 : Cdn_atm_skin=c0
651 0 : Cdn_atm_floe=c0
652 0 : Cdn_atm_pond=c0
653 0 : Cdn_atm_rdg =c0
654 :
655 0 : if (aice > p001) then
656 :
657 0 : Cdn_atm_skin = csa
658 0 : Cdn_ocn_skin = csw
659 :
660 0 : ai = aice
661 0 : aii = c1/ai
662 :
663 : !------------------------------------------------------------
664 : ! Compute average quantities
665 : !------------------------------------------------------------
666 :
667 : ! ponds
668 0 : apond = c0
669 0 : if (tr_pond) then
670 0 : do n = 1,ncat
671 : ! area of pond per unit area of grid cell
672 0 : apond = apond+apnd(n)*aicen(n)
673 : enddo
674 : endif
675 :
676 : ! draft and freeboard (see Eq. 27)
677 0 : hdraft = (rhoi*vice+rhos*vsno)*aii/rhow ! without ponds
678 0 : hfreebd = (vice+vsno)*aii-hdraft
679 :
680 : ! Do not allow draft larger than ice thickness (see Eq. 28)
681 0 : if (hdraft >= vice*aii) then
682 : ! replace excess snow with ice so hi~=hdraft
683 : hfreebd = (hdraft*ai*(c1-rhoi/rhow) + &
684 : (vsno-(vice-hdraft*ai)*rhoi/rhos) * & ! LCOV_EXCL_LINE
685 0 : (c1-rhos/rhow))*aii ! Stoessel1993
686 : endif
687 :
688 : ! floe size parameterization see Eq. 13
689 0 : lfloe = Lmin * (astar / (astar - ai))**beta
690 :
691 : ! distance between floes parameterization see Eq. 14
692 0 : dfloe = lfloe * (c1/sqrt(ai) - c1)
693 :
694 : ! Relate ridge height and distance between ridges to
695 : ! ridged ice area fraction and ridged ice mean thickness
696 : ! Assumes total volume of ridged ice is split into ridges and keels.
697 : ! Then assume total ridges volume over total area of ridges =
698 : ! volume of one average ridge / area of one average ridge
699 : ! Same for keels.
700 :
701 0 : ardg=c0
702 0 : vrdg=c0
703 0 : do n=1,ncat
704 : ! ridged ice area fraction over grid cell
705 0 : ardg=ardg+(c1-alvl(n))*aicen(n)
706 : ! total ridged ice volume per unit grid cell area
707 0 : vrdg=vrdg+(c1-vlvl(n))*vicen(n)
708 : enddo
709 :
710 : ! hridge, hkeel, distrdg and dkeel estimates from CICE for
711 : ! simple triangular geometry
712 0 : if (ardg > p001) then
713 : ! see Eq. 25 and Eq. 26
714 : hridge = vrdg/ardg*c2 &
715 : * (alpha2+beta2*hkoverhr/dkoverdr*tanar/tanak) & ! LCOV_EXCL_LINE
716 0 : / (phir*c1+phik*tanar/tanak*hkoverhr**c2/dkoverdr)
717 : distrdg = c2*hridge*ai/ardg &
718 0 : * (alpha2/tanar+beta2/tanak*hkoverhr/dkoverdr)
719 0 : hkeel = hkoverhr * hridge
720 0 : dkeel = dkoverdr * distrdg
721 :
722 : ! Use the height of ridges relative to the mean freeboard of
723 : ! the pack. Therefore skin drag and ridge drag differ in
724 : ! this code as compared to Tsamados et al. (2014) equations
725 : ! 10 and 18, which reference both to sea level.
726 0 : tmp1 = max(c0,hridge - hfreebd)
727 :
728 : !------------------------------------------------------------
729 : ! Skin drag (atmo)
730 : !------------------------------------------------------------
731 :
732 0 : Cdn_atm_skin = csa*(c1 - mrdg*tmp1/distrdg)
733 0 : Cdn_atm_skin = max(min(Cdn_atm_skin,camax),c0)
734 :
735 : !------------------------------------------------------------
736 : ! Ridge effect (atmo)
737 : !------------------------------------------------------------
738 :
739 0 : if (tmp1 > puny) then
740 0 : sca = c1 - exp(-sHGB*distrdg/tmp1) ! see Eq. 9
741 0 : ctecar = cra*p5
742 : Cdn_atm_rdg = ctecar*tmp1/distrdg*sca* &
743 0 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
744 0 : Cdn_atm_rdg = min(Cdn_atm_rdg,camax)
745 : endif
746 :
747 : ! Use the depth of keels relative to the mean draft of
748 : ! the pack. Therefore skin drag and keel drag differ in
749 : ! this code as compared to Tsamados et al. (2014) equations
750 : ! 11 and 19, which reference both to sea level. In some
751 : ! circumstances, hkeel can be less than hdraft because hkoverhr
752 : ! is constant, and max(c0,...) temporarily addresses this.
753 0 : tmp1 = max(c0,hkeel - hdraft)
754 :
755 : !------------------------------------------------------------
756 : ! Skin drag bottom ice (ocean)
757 : !------------------------------------------------------------
758 :
759 0 : Cdn_ocn_skin = csw * (c1 - mrdgo*tmp1/dkeel)
760 0 : Cdn_ocn_skin = max(min(Cdn_ocn_skin,cwmax), c0)
761 :
762 : !------------------------------------------------------------
763 : ! Keel effect (ocean)
764 : !------------------------------------------------------------
765 :
766 0 : if (tmp1 > puny) then
767 0 : scw = c1 - exp(-sHGB*dkeel/tmp1)
768 0 : ctecwk = crw*p5
769 : Cdn_ocn_keel = ctecwk*tmp1/dkeel*scw* &
770 0 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
771 0 : Cdn_ocn_keel = max(min(Cdn_ocn_keel,cwmax),c0)
772 : endif
773 :
774 : endif ! ardg > 0.001
775 :
776 : !------------------------------------------------------------
777 : ! Floe edge drag effect (atmo)
778 : !------------------------------------------------------------
779 :
780 0 : if (hfreebd > puny) then
781 0 : sca = c1 - exp(-sl*beta*(c1-ai))
782 0 : ctecaf = cfa*p5*(log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2*sca
783 0 : Cdn_atm_floe = ctecaf * hfreebd / lfloe
784 0 : Cdn_atm_floe = max(min(Cdn_atm_floe,camax),c0)
785 : endif
786 :
787 : !------------------------------------------------------------
788 : ! Pond edge effect (atmo)
789 : !------------------------------------------------------------
790 :
791 0 : if (hfreebd > puny) then
792 0 : sca = (apond)**(c1/(zref*beta))
793 0 : lp = lpmin*(1-apond)+lpmax*apond
794 : Cdn_atm_pond = cpa*p5*sca*apond*hfreebd/lp &
795 0 : * (log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2
796 0 : Cdn_atm_pond = min(Cdn_atm_pond,camax)
797 : endif
798 :
799 : !------------------------------------------------------------
800 : ! Floe edge drag effect (ocean)
801 : !------------------------------------------------------------
802 :
803 0 : if (hdraft > puny) then
804 0 : scw = c1 - exp(-sl*beta*(c1-ai))
805 0 : ctecwf = cfw*p5*(log(hdraft*ocnrufi)/log(zref*ocnrufi))**c2*scw
806 0 : Cdn_ocn_floe = ctecwf * hdraft / lfloe
807 0 : Cdn_ocn_floe = max(min(Cdn_ocn_floe,cwmax),c0)
808 : endif
809 :
810 : !------------------------------------------------------------
811 : ! Total drag coefficient (atmo)
812 : !------------------------------------------------------------
813 :
814 0 : Cdn_atm = Cdn_atm_skin + Cdn_atm_floe + Cdn_atm_pond + Cdn_atm_rdg
815 0 : Cdn_atm = min(Cdn_atm,camax)
816 :
817 : !------------------------------------------------------------
818 : ! Total drag coefficient (ocean)
819 : !------------------------------------------------------------
820 :
821 0 : Cdn_ocn = Cdn_ocn_skin + Cdn_ocn_floe + Cdn_ocn_keel
822 0 : Cdn_ocn = min(Cdn_ocn,cwmax)
823 :
824 : endif
825 :
826 0 : end subroutine neutral_drag_coeffs
827 :
828 : !=======================================================================
829 : !autodocument_start icepack_atm_boundary
830 : !
831 :
832 46176939 : subroutine icepack_atm_boundary(sfctype, &
833 : Tsf, potT, & ! LCOV_EXCL_LINE
834 : uatm, vatm, & ! LCOV_EXCL_LINE
835 : wind, zlvl, & ! LCOV_EXCL_LINE
836 : Qa, rhoa, & ! LCOV_EXCL_LINE
837 : strx, stry, & ! LCOV_EXCL_LINE
838 : Tref, Qref, & ! LCOV_EXCL_LINE
839 : delt, delq, & ! LCOV_EXCL_LINE
840 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
841 : Cdn_atm, & ! LCOV_EXCL_LINE
842 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
843 : Qa_iso, Qref_iso, & ! LCOV_EXCL_LINE
844 : uvel, vvel, & ! LCOV_EXCL_LINE
845 : Uref, zlvs)
846 :
847 : character (len=3), intent(in) :: &
848 : sfctype ! ice or ocean
849 :
850 : real (kind=dbl_kind), intent(in) :: &
851 : Tsf , & ! surface temperature of ice or ocean ! LCOV_EXCL_LINE
852 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
853 : uatm , & ! x-direction wind speed (m/s) ! LCOV_EXCL_LINE
854 : vatm , & ! y-direction wind speed (m/s) ! LCOV_EXCL_LINE
855 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
856 : zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) ! LCOV_EXCL_LINE
857 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
858 : rhoa ! air density (kg/m^3)
859 :
860 : real (kind=dbl_kind), intent(inout) :: &
861 : Cdn_atm , & ! neutral drag coefficient ! LCOV_EXCL_LINE
862 : Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff
863 :
864 : real (kind=dbl_kind), intent(inout) :: &
865 : strx , & ! x surface stress (N) ! LCOV_EXCL_LINE
866 : stry ! y surface stress (N)
867 :
868 : real (kind=dbl_kind), intent(inout) :: &
869 : Tref , & ! reference height temperature (K) ! LCOV_EXCL_LINE
870 : Qref , & ! reference height specific humidity (kg/kg) ! LCOV_EXCL_LINE
871 : delt , & ! potential T difference (K) ! LCOV_EXCL_LINE
872 : delq , & ! humidity difference (kg/kg) ! LCOV_EXCL_LINE
873 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
874 : lhcoef ! transfer coefficient for latent heat
875 :
876 : real (kind=dbl_kind), intent(in), dimension(:), optional :: &
877 : Qa_iso ! specific isotopic humidity (kg/kg)
878 :
879 : real (kind=dbl_kind), intent(inout), dimension(:), optional :: &
880 : Qref_iso ! reference specific isotopic humidity (kg/kg)
881 :
882 : real (kind=dbl_kind), intent(in), optional :: &
883 : uvel , & ! x-direction ice speed (m/s) ! LCOV_EXCL_LINE
884 : vvel , & ! y-direction ice speed (m/s) ! LCOV_EXCL_LINE
885 : zlvs ! atm level height for scalars (if different than zlvl) (m)
886 :
887 : real (kind=dbl_kind), intent(out), optional :: &
888 : Uref ! reference height wind speed (m/s)
889 :
890 : !autodocument_end
891 :
892 : ! local variables
893 :
894 : real (kind=dbl_kind) :: &
895 6371954 : l_uvel, l_vvel, l_Uref
896 :
897 : logical (kind=log_kind), save :: &
898 : first_call_ice = .true. ! first call flag
899 :
900 : character(len=*),parameter :: subname='(icepack_atm_boundary)'
901 :
902 : !------------------------------------------------------------
903 : ! Check optional arguments
904 : ! Need separate first_call flags for 'ice' and 'ocn' sfctype
905 : !------------------------------------------------------------
906 :
907 46176939 : if (sfctype == 'ice') then
908 33852099 : if (icepack_chkoptargflag(first_call_ice)) then
909 44 : if (tr_iso) then
910 0 : if (.not.(present(Qa_iso).and.present(Qref_iso))) then
911 0 : call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T')
912 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
913 0 : return
914 : endif
915 : endif
916 : endif
917 : endif
918 :
919 46176939 : l_uvel = c0
920 46176939 : l_vvel = c0
921 46176939 : l_Uref = c0
922 46176939 : if (present(uvel)) then
923 33852099 : l_uvel = uvel
924 : endif
925 46176939 : if (present(vvel)) then
926 33852099 : l_vvel = vvel
927 : endif
928 :
929 46176939 : Cdn_atm_ratio_n = c1
930 :
931 46176939 : if (trim(atmbndy) == 'constant') then
932 : call atmo_boundary_const (sfctype, calc_strair, &
933 : uatm, vatm, & ! LCOV_EXCL_LINE
934 : wind, rhoa, & ! LCOV_EXCL_LINE
935 : strx, stry, & ! LCOV_EXCL_LINE
936 : Tsf, potT, & ! LCOV_EXCL_LINE
937 : Qa, & ! LCOV_EXCL_LINE
938 : delt, delq, & ! LCOV_EXCL_LINE
939 0 : lhcoef, shcoef )
940 0 : if (icepack_warnings_aborted(subname)) return
941 : else
942 : call atmo_boundary_layer (sfctype, &
943 : calc_strair, formdrag, & ! LCOV_EXCL_LINE
944 : Tsf, potT, & ! LCOV_EXCL_LINE
945 : uatm, vatm, & ! LCOV_EXCL_LINE
946 : wind, zlvl, & ! LCOV_EXCL_LINE
947 : Qa, rhoa, & ! LCOV_EXCL_LINE
948 : strx, stry, & ! LCOV_EXCL_LINE
949 : Tref, Qref, & ! LCOV_EXCL_LINE
950 : delt, delq, & ! LCOV_EXCL_LINE
951 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
952 : Cdn_atm, & ! LCOV_EXCL_LINE
953 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
954 : Qa_iso=Qa_iso, & ! LCOV_EXCL_LINE
955 : Qref_iso=Qref_iso, & ! LCOV_EXCL_LINE
956 : uvel=l_uvel, vvel=l_vvel,& ! LCOV_EXCL_LINE
957 46176939 : Uref=l_Uref, zlvs=zlvs )
958 46176939 : if (icepack_warnings_aborted(subname)) return
959 : endif ! atmbndy
960 :
961 46176939 : if (present(Uref)) then
962 33852099 : Uref = l_Uref
963 : endif
964 :
965 46176939 : if (sfctype == 'ice') then
966 33852099 : first_call_ice = .false.
967 : endif
968 :
969 : end subroutine icepack_atm_boundary
970 :
971 : !=======================================================================
972 :
973 398556998 : function compute_stability_parameter(zlvl , thva , &
974 : ustar, tstar, & ! LCOV_EXCL_LINE
975 : qstar, Qa) & ! LCOV_EXCL_LINE
976 46413088 : result(hol)
977 :
978 : real (kind=dbl_kind), intent(in) :: &
979 : zlvl , & ! atm level height (m) ! LCOV_EXCL_LINE
980 : thva , & ! virtual temperature (K) ! LCOV_EXCL_LINE
981 : ustar , & ! turbulent scale for momentum ! LCOV_EXCL_LINE
982 : tstar , & ! turbulent scale for temperature ! LCOV_EXCL_LINE
983 : qstar , & ! turbulent scale for humidity ! LCOV_EXCL_LINE
984 : Qa ! specific humidity (kg/kg)
985 :
986 : real (kind=dbl_kind) :: &
987 : hol ! H (at zlvl) over L
988 :
989 : character(len=*),parameter :: subname='(compute_stability_parameter)'
990 :
991 : hol = vonkar * gravit * zlvl &
992 : * (tstar/thva & ! LCOV_EXCL_LINE
993 : + qstar/(c1/zvir+Qa)) & ! LCOV_EXCL_LINE
994 398556998 : / ustar**2
995 398556998 : hol = sign( min(abs(hol),c10), hol)
996 :
997 398556998 : end function compute_stability_parameter
998 :
999 : !=======================================================================
1000 :
1001 459953622 : subroutine compute_stability_function(qty, hol, stable, psi)
1002 :
1003 : character (len=*), intent(in) :: &
1004 : qty ! 'momentum' or 'scalar'
1005 :
1006 : real (kind=dbl_kind), intent(in) :: &
1007 : hol ! H over L
1008 :
1009 : real (kind=dbl_kind), intent(out) :: &
1010 : psi , & ! stability function at hol ! LCOV_EXCL_LINE
1011 : stable ! unit step function at hol
1012 :
1013 : ! local variables
1014 :
1015 : real (kind=dbl_kind) :: &
1016 : psi_stable , & ! stable stability funcion at hol ! LCOV_EXCL_LINE
1017 62709046 : psi_unstable ! unstable stability funcion at hol
1018 :
1019 : character(len=*),parameter :: subname='(compute_stability_function)'
1020 :
1021 459953622 : stable = p5 + sign(p5 , hol)
1022 :
1023 : psi_stable = -(0.7_dbl_kind*hol &
1024 : + 0.75_dbl_kind*(hol-14.3_dbl_kind) & ! LCOV_EXCL_LINE
1025 459953622 : * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
1026 :
1027 459953622 : if(trim(qty) == 'momentum') then
1028 229976811 : psi_unstable = psi_momentum_unstable(hol)
1029 229976811 : elseif(trim(qty) == 'scalar') then
1030 229976811 : psi_unstable = psi_scalar_unstable(hol)
1031 : else
1032 0 : call icepack_warnings_add(subname//' incorrect qty: ' // qty)
1033 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1034 : endif
1035 :
1036 459953622 : psi = psi_stable*stable + (c1 - stable)*psi_unstable
1037 :
1038 459953622 : end subroutine compute_stability_function
1039 :
1040 : !------------------------------------------------------------
1041 : ! Define functions
1042 : !------------------------------------------------------------
1043 :
1044 : !=======================================================================
1045 :
1046 229976811 : real(kind=dbl_kind) function psi_momentum_unstable(hol)
1047 :
1048 : real(kind=dbl_kind), intent(in) :: hol
1049 :
1050 31354523 : real(kind=dbl_kind) :: xd
1051 :
1052 229976811 : xd = capital_X(hol)
1053 :
1054 : psi_momentum_unstable = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
1055 229976811 : - c2*atan(xd) + pih
1056 :
1057 229976811 : end function psi_momentum_unstable
1058 :
1059 : !=======================================================================
1060 :
1061 276153750 : real(kind=dbl_kind) function psi_scalar_unstable(hol)
1062 :
1063 : real(kind=dbl_kind), intent(in) :: hol
1064 :
1065 37726477 : real(kind=dbl_kind) :: xd
1066 :
1067 276153750 : xd = capital_X(hol)
1068 :
1069 276153750 : psi_scalar_unstable = c2 * log((c1 + xd*xd)/c2)
1070 :
1071 276153750 : end function psi_scalar_unstable
1072 :
1073 : !=======================================================================
1074 :
1075 506130561 : real(kind=dbl_kind) function capital_X(hol)
1076 :
1077 : real(kind=dbl_kind), intent(in) :: hol
1078 :
1079 506130561 : capital_X = sqrt(max(sqrt(abs(c1 - c16*hol)) , c1))
1080 :
1081 506130561 : end function capital_X
1082 :
1083 : !=======================================================================
1084 :
1085 : end module icepack_atmo
1086 :
1087 : !=======================================================================
|
