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