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 10161340 : subroutine atmo_boundary_layer (sfctype, &
53 : calc_strair, formdrag, &
54 : Tsf, potT, &
55 : uatm, vatm, &
56 : wind, zlvl, &
57 : Qa, rhoa, &
58 : strx, stry, &
59 : Tref, Qref, &
60 : delt, delq, &
61 : lhcoef, shcoef, &
62 : Cdn_atm, &
63 : Cdn_atm_ratio_n, &
64 10161340 : Qa_iso, Qref_iso, &
65 : uvel, vvel, &
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
75 : formdrag ! if true, calculate form drag
76 :
77 : real (kind=dbl_kind), intent(in) :: &
78 : Tsf , & ! surface temperature of ice or ocean
79 : potT , & ! air potential temperature (K)
80 : uatm , & ! x-direction wind speed (m/s)
81 : vatm , & ! y-direction wind speed (m/s)
82 : wind , & ! wind speed (m/s)
83 : zlvl , & ! atm level height (m)
84 : Qa , & ! specific humidity (kg/kg)
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)
95 : stry ! y surface stress (N)
96 :
97 : real (kind=dbl_kind), intent(inout) :: &
98 : Tref , & ! reference height temperature (K)
99 : Qref , & ! reference height specific humidity (kg/kg)
100 : delt , & ! potential T difference (K)
101 : delq , & ! humidity difference (kg/kg)
102 : shcoef , & ! transfer coefficient for sensible heat
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)
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)
128 : psimh , & ! stability function at zlvl (momentum)
129 : tau , & ! stress at zlvl
130 : fac , & ! interpolation factor
131 : al2 , & ! ln(z10 /zTrf)
132 : psix2 , & ! stability function at zTrf (heat and water)
133 : ssq , & ! sat surface humidity (kg/kg)
134 : qqq , & ! for qsat, dqsfcdt
135 : TTT , & ! for qsat, dqsfcdt
136 : qsat , & ! the saturation humidity of air (kg/m^3)
137 : Lheat , & ! Lvap or Lsub, depending on surface type
138 : umin ! minimum wind speed (m/s)
139 :
140 : real (kind=dbl_kind) :: &
141 : ustar , & ! ustar (m/s)
142 : ustar_prev , & ! ustar_prev (m/s)
143 : tstar , & ! tstar
144 : qstar , & ! qstar
145 : ratio , & ! ratio
146 : rdn , & ! sqrt of neutral exchange coefficient (momentum)
147 : rhn , & ! sqrt of neutral exchange coefficient (heat)
148 : ren , & ! sqrt of neutral exchange coefficient (water)
149 : rd , & ! sqrt of exchange coefficient (momentum)
150 : re , & ! sqrt of exchange coefficient (water)
151 : rh , & ! sqrt of exchange coefficient (heat)
152 : vmag , & ! surface wind magnitude (m/s)
153 : alzm , & ! ln(zlvl /z10)
154 : alzs , & ! ln(zlvs /z10) (if zlvs present)
155 : thva , & ! virtual temperature (K)
156 : cp , & ! specific heat of moist air
157 : holm , & ! H (at zlvl ) over L
158 : hols , & ! H (at zlvs ) over L (if zlvs present)
159 : stable, & ! stability factor
160 : cpvir , & ! defined as cp_wv/cp_air - 1.
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 10161340 : al2 = log(zref/zTrf)
169 :
170 : !------------------------------------------------------------
171 : ! Initialize
172 : !------------------------------------------------------------
173 :
174 10161340 : cpvir = cp_wv/cp_air-c1 ! defined as cp_wv/cp_air - 1.
175 :
176 10161340 : if (highfreq) then
177 390544 : umin = p5 ! minumum allowable wind-ice speed difference of 0.5 m/s
178 : else
179 9770796 : umin = c1 ! minumum allowable wind speed of 1m/s
180 : endif
181 :
182 10161340 : Tref = c0
183 10161340 : Qref = c0
184 10161340 : Uref = c0
185 10161340 : delt = c0
186 10161340 : delq = c0
187 10161340 : shcoef = c0
188 10161340 : 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 10161340 : if (sfctype(1:3)=='ice') then
200 :
201 8287900 : qqq = qqqice ! for qsat
202 8287900 : TTT = TTTice ! for qsat
203 8287900 : Lheat = Lsub ! ice to vapor
204 :
205 8287900 : if (highfreq) then
206 : vmag = max(umin, sqrt( (uatm-uvel)**2 + &
207 318148 : (vatm-vvel)**2) )
208 : else
209 7969752 : vmag = max(umin, wind)
210 : endif
211 :
212 8287900 : if (formdrag .and. Cdn_atm > puny) then
213 216628 : rdn = sqrt(Cdn_atm)
214 : else
215 8071272 : rdn = vonkar/log(zref/iceruf) ! neutral coefficient
216 8071272 : Cdn_atm = rdn * rdn
217 : endif
218 :
219 1873440 : elseif (sfctype(1:3)=='ocn') then
220 :
221 1873440 : qqq = qqqocn
222 1873440 : TTT = TTTocn
223 1873440 : Lheat = Lvap ! liquid to vapor
224 1873440 : vmag = max(umin, wind)
225 : rdn = sqrt(0.0027_dbl_kind/vmag &
226 1873440 : + .000142_dbl_kind + .0000764_dbl_kind*vmag)
227 :
228 : endif ! sfctype
229 :
230 : !------------------------------------------------------------
231 : ! define some more needed variables
232 : !------------------------------------------------------------
233 :
234 10161340 : TsfK = Tsf + Tffresh ! surface temp (K)
235 10161340 : delt = potT - TsfK ! pot temp diff (K)
236 10161340 : qsat = qqq * exp(-TTT/TsfK) ! saturation humidity (kg/m^3)
237 10161340 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
238 :
239 10161340 : thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K)
240 10161340 : delq = Qa - ssq ! spec hum dif (kg/kg)
241 10161340 : alzm = log(zlvl/zref)
242 10161340 : if (present(zlvs)) then
243 0 : alzs = log(zlvs/zref)
244 : else
245 10161340 : alzs = alzm
246 : endif
247 10161340 : 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 10161340 : rhn = rdn
255 10161340 : ren = rdn
256 :
257 : ! ustar,tstar,qstar
258 10161340 : ustar = rdn * vmag
259 10161340 : tstar = rhn * delt
260 10161340 : qstar = ren * delq
261 :
262 : !------------------------------------------------------------
263 : ! iterate to converge on Z/L, ustar, tstar and qstar
264 : !------------------------------------------------------------
265 :
266 10161340 : ustar_prev = c2 * ustar
267 :
268 10161340 : k = 1
269 60786995 : do while (abs(ustar - ustar_prev)/ustar > atmiter_conv .and. k <= natmiter)
270 50625655 : k = k + 1
271 50625655 : ustar_prev = ustar
272 :
273 : ! compute stability & evaluate all stability functions
274 : holm = compute_stability_parameter(zlvl , thva , &
275 : ustar, tstar, &
276 50625655 : qstar, Qa)
277 50625655 : if (present(zlvs)) then
278 : hols = compute_stability_parameter(zlvs , thva , &
279 : ustar, tstar, &
280 0 : qstar, Qa)
281 : else
282 50625655 : hols = holm
283 : endif
284 :
285 50625655 : call compute_stability_function('momentum', holm, stable, psimh)
286 50625655 : call compute_stability_function('scalar' , hols, stable, psixh)
287 :
288 : ! shift all coeffs to measurement height and stability
289 50625655 : rd = rdn / (c1+rdn/vonkar*(alzm-psimh))
290 50625655 : rh = rhn / (c1+rhn/vonkar*(alzs-psixh))
291 50625655 : re = ren / (c1+ren/vonkar*(alzs-psixh))
292 :
293 : ! update ustar, tstar, qstar using updated, shifted coeffs
294 50625655 : ustar = rd * vmag
295 50625655 : tstar = rh * delt
296 50625655 : qstar = re * delq
297 :
298 : enddo ! end iteration
299 :
300 10161340 : if (calc_strair) then
301 :
302 : ! initialize
303 10161340 : strx = c0
304 10161340 : stry = c0
305 :
306 10161340 : 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 318148 : tau = rhoa * rd * rd ! not the stress at zlvl
317 :
318 : ! high frequency momentum coupling following Roberts et al. (2014)
319 318148 : strx = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (uatm-uvel)
320 318148 : 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 9843192 : tau = rhoa * ustar * rd ! not the stress at zlvl
333 9843192 : strx = tau * uatm
334 9843192 : stry = tau * vatm
335 :
336 : endif
337 :
338 10161340 : 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 10161340 : 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 10161340 : shcoef = rhoa * ustar * cp * rh + c1
357 10161340 : lhcoef = rhoa * ustar * Lheat * re
358 : endif
359 :
360 : !------------------------------------------------------------
361 : ! Compute diagnostics: 2m ref T, Q, U
362 : !------------------------------------------------------------
363 :
364 10161340 : hols = hols*zTrf/zlvl
365 10161340 : psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols)
366 : fac = (rh/vonkar) &
367 10161340 : * (alzs + al2 - psixh + psix2)
368 10161340 : Tref = potT - delt*fac
369 10161340 : Tref = Tref - p01*zTrf ! pot temp to temp correction
370 : fac = (re/vonkar) &
371 10161340 : * (alzs + al2 - psixh + psix2)
372 10161340 : Qref = Qa - delq*fac
373 :
374 10161340 : if (highfreq .and. sfctype(1:3)=='ice') then
375 318148 : Uref = sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * rd / rdn
376 : else
377 9843192 : Uref = vmag * rd / rdn
378 : endif
379 :
380 10161340 : if (tr_iso .and. sfctype(1:3)=='ice') then
381 901808 : Qref_iso(:) = c0
382 901808 : do n = 1, n_iso
383 676356 : ratio = c0
384 676356 : if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2)
385 901808 : Qref_iso(n) = Qa_iso(n) - ratio*delq*fac
386 : enddo
387 : endif
388 :
389 10161340 : 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 588714 : subroutine atmo_boundary_const (sfctype, calc_strair, &
399 : uatm, vatm, &
400 : wind, rhoa, &
401 : strx, stry, &
402 : Tsf, potT, &
403 : Qa, &
404 : delt, delq, &
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
415 : potT , & ! air potential temperature (K)
416 : Qa , & ! specific humidity (kg/kg)
417 : uatm , & ! x-direction wind speed (m/s)
418 : vatm , & ! y-direction wind speed (m/s)
419 : wind , & ! wind speed (m/s)
420 : rhoa ! air density (kg/m^3)
421 :
422 : real (kind=dbl_kind), intent(inout):: &
423 : strx , & ! x surface stress (N)
424 : stry ! y surface stress (N)
425 :
426 : real (kind=dbl_kind), intent(out):: &
427 : delt , & ! potential T difference (K)
428 : delq , & ! humidity difference (kg/kg)
429 : shcoef , & ! transfer coefficient for sensible heat
430 : lhcoef ! transfer coefficient for latent heat
431 :
432 : ! local variables
433 :
434 : real (kind=dbl_kind) :: &
435 : TsfK, & ! surface temperature in Kelvin (K)
436 : qsat, & ! the saturation humidity of air (kg/m^3)
437 : ssq , & ! sat surface humidity (kg/kg)
438 : tau, & ! stress at zlvl
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 588714 : delt = c0
448 588714 : delq = c0
449 588714 : shcoef = c0
450 588714 : lhcoef = c0
451 :
452 588714 : if (calc_strair) then
453 :
454 588714 : strx = c0
455 588714 : stry = c0
456 :
457 : !------------------------------------------------------------
458 : ! momentum flux
459 : !------------------------------------------------------------
460 588714 : tau = rhoa * 0.0012_dbl_kind * wind
461 : !AOMIP tau = rhoa * (1.10_dbl_kind + c4*p01*wind) &
462 : !AOMIP * wind * p001
463 588714 : strx = tau * uatm
464 588714 : stry = tau * vatm
465 :
466 : endif ! calc_strair
467 :
468 : !------------------------------------------------------------
469 : ! define variables that depend on surface type
470 : !------------------------------------------------------------
471 :
472 588714 : if (sfctype(1:3)=='ice') then
473 490038 : Lheat = Lsub ! ice to vapor
474 98676 : elseif (sfctype(1:3)=='ocn') then
475 98676 : Lheat = Lvap ! liquid to vapor
476 : endif ! sfctype
477 :
478 : !------------------------------------------------------------
479 : ! potential temperature and specific humidity differences
480 : !------------------------------------------------------------
481 :
482 588714 : TsfK = Tsf + Tffresh ! surface temp (K)
483 588714 : qsat = qqqocn * exp(-TTTocn/TsfK) ! sat humidity (kg/m^3)
484 588714 : ssq = qsat / rhoa ! sat surf hum (kg/kg)
485 :
486 588714 : delt= potT - TsfK ! pot temp diff (K)
487 588714 : delq= Qa - ssq ! spec hum dif (kg/kg)
488 :
489 : !------------------------------------------------------------
490 : ! coefficients for turbulent flux calculation
491 : !------------------------------------------------------------
492 :
493 588714 : shcoef = senscoef*cp_air*rhoa*wind
494 588714 : lhcoef = latncoef*Lheat *rhoa*wind
495 :
496 588714 : 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 96528 : subroutine neutral_drag_coeffs (apondn, &
511 96528 : alvl, vlvl, &
512 : aice, vice, &
513 96528 : vsno, aicen, &
514 96528 : vicen, &
515 : Cdn_ocn, Cdn_ocn_skin, &
516 : Cdn_ocn_floe, Cdn_ocn_keel,&
517 : Cdn_atm, Cdn_atm_skin, &
518 : Cdn_atm_floe, Cdn_atm_pond,&
519 : Cdn_atm_rdg, hfreebd, &
520 : hdraft, hridge, &
521 : distrdg, hkeel, &
522 : dkeel, lfloe, &
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
529 : alvl ,& ! level ice area fraction (of grid cell ?)
530 : vlvl ! level ice mean thickness
531 :
532 : real (kind=dbl_kind), intent(in) :: &
533 : aice , & ! concentration of ice
534 : vice , & ! volume per unit area of ice
535 : vsno ! volume per unit area of snow
536 :
537 : real (kind=dbl_kind), dimension (:), intent(in) :: &
538 : aicen , & ! concentration of ice
539 : vicen ! volume per unit area of ice (m)
540 :
541 : real (kind=dbl_kind), intent(out) :: &
542 : hfreebd , & ! freeboard (m)
543 : hdraft , & ! draught of ice + snow column (Stoessel1993)
544 : hridge , & ! ridge height
545 : distrdg , & ! distance between ridges
546 : hkeel , & ! keel depth
547 : dkeel , & ! distance between keels
548 : lfloe , & ! floe length (m)
549 : dfloe , & ! distance between floes
550 : Cdn_ocn , & ! ocean-ice neutral drag coefficient
551 : Cdn_ocn_skin , & ! drag coefficient due to skin drag
552 : Cdn_ocn_floe , & ! drag coefficient due to floe edges
553 : Cdn_ocn_keel , & ! drag coefficient due to keels
554 : Cdn_atm , & ! ice-atmosphere drag coefficient
555 : Cdn_atm_skin , & ! drag coefficient due to skin drag
556 : Cdn_atm_floe , & ! drag coefficient due to floe edges
557 : Cdn_atm_pond , & ! drag coefficient due to ponds
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]
564 : mrdg = c20 ,&! screening effect see Lu2011 [5,50]
565 : mrdgo = c10 ,&! screening effect see Lu2011 [5,50]
566 : beta = p5 ,&! power exponent appearing in astar and
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]
570 : cfa = p2 ,&! Eq. 12 ratio of local from drag over
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]
578 : sl = 22._dbl_kind ,&! Sheltering parameter Lupkes2012 [10,30]
579 : lpmin = 2.26_dbl_kind ,&! min pond length (m) see Eq. 17 [1,10]
580 : lpmax = 24.63_dbl_kind ,&! max pond length (m) see Eq. 17 [10,100]
581 : tanar = p4 ,&! 0.25 sail slope = 14 deg [0.4,1]
582 : tanak = p4 ,&! 0.58 keel slope = 30 deg [0.4,1]
583 : phir = 0.8_dbl_kind ,&! porosity of ridges [0.4,1]
584 : phik = 0.8_dbl_kind ,&! porosity of keels [0.4,1]
585 : hkoverhr = c4 ,&! hkeel/hridge ratio [4,8]
586 : dkoverdr = c1 ,&! dkeel/distrdg ratio [1,5]
587 : sHGB = 0.18_dbl_kind ,&! Lupkes2012 Eq. 28, Hanssen1988,
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
597 : ctecaf, & ! constante
598 : ctecwf, & ! constante
599 : sca, & ! wind attenuation function
600 : scw, & ! ocean attenuation function
601 : lp, & ! pond length (m)
602 : ctecar, &
603 : ctecwk, &
604 : ai, aii, & ! ice area and its inverse
605 : ocnrufi, & ! inverse ocean roughness
606 : icerufi, & ! inverse ice roughness
607 : tmp1 ! temporary
608 :
609 : real (kind=dbl_kind) :: &
610 : apond , & ! melt pond fraction of grid cell
611 : ardg , & ! ridged ice area fraction of grid cell
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
619 : cwmax = 0.06_dbl_kind ! Maximum for ocean drag
620 :
621 : character(len=*),parameter :: subname='(neutral_drag_coeffs)'
622 :
623 96528 : astar = c1/(c1-(Lmin/Lmax)**(c1/beta))
624 :
625 : !-----------------------------------------------------------------
626 : ! Initialize across entire grid
627 : !-----------------------------------------------------------------
628 :
629 96528 : ocnrufi = c1/ocnruf ! inverse ocean roughness
630 96528 : icerufi = c1/iceruf ! inverse ice roughness
631 96528 : hfreebd=c0
632 96528 : hdraft =c0
633 96528 : hridge =c0
634 96528 : distrdg=c0
635 96528 : hkeel =c0
636 96528 : dkeel =c0
637 96528 : lfloe =c0
638 96528 : dfloe =c0
639 96528 : Cdn_ocn=dragio
640 96528 : Cdn_ocn_skin=c0
641 96528 : Cdn_ocn_floe=c0
642 96528 : Cdn_ocn_keel=c0
643 96528 : Cdn_atm = (vonkar/log(zref/iceruf)) * (vonkar/log(zref/iceruf))
644 96528 : Cdn_atm_skin=c0
645 96528 : Cdn_atm_floe=c0
646 96528 : Cdn_atm_pond=c0
647 96528 : Cdn_atm_rdg =c0
648 :
649 96528 : if (aice > p001) then
650 :
651 72392 : Cdn_atm_skin = csa
652 72392 : Cdn_ocn_skin = csw
653 :
654 72392 : ai = aice
655 72392 : aii = c1/ai
656 :
657 : !------------------------------------------------------------
658 : ! Compute average quantities
659 : !------------------------------------------------------------
660 :
661 : ! ponds
662 72392 : apond = c0
663 72392 : if (tr_pond) then
664 0 : do n = 1,ncat
665 : ! area of pond per unit area of grid cell
666 0 : apond = apond+apondn(n)*aicen(n)
667 : enddo
668 : endif
669 :
670 : ! draft and freeboard (see Eq. 27)
671 72392 : hdraft = (rhoi*vice+rhos*vsno)*aii/rhow ! without ponds
672 72392 : hfreebd = (vice+vsno)*aii-hdraft
673 :
674 : ! Do not allow draft larger than ice thickness (see Eq. 28)
675 72392 : 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) * &
679 0 : (c1-rhos/rhow))*aii ! Stoessel1993
680 : endif
681 :
682 : ! floe size parameterization see Eq. 13
683 72392 : lfloe = Lmin * (astar / (astar - ai))**beta
684 :
685 : ! distance between floes parameterization see Eq. 14
686 72392 : 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 72392 : ardg=c0
696 72392 : vrdg=c0
697 434352 : do n=1,ncat
698 : ! ridged ice area fraction over grid cell
699 361960 : ardg=ardg+(c1-alvl(n))*aicen(n)
700 : ! total ridged ice volume per unit grid cell area
701 434352 : 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 72392 : if (ardg > p001) then
707 : ! see Eq. 25 and Eq. 26
708 : hridge = vrdg/ardg*c2 &
709 : * (alpha2+beta2*hkoverhr/dkoverdr*tanar/tanak) &
710 0 : / (phir*c1+phik*tanar/tanak*hkoverhr**c2/dkoverdr)
711 : distrdg = c2*hridge*ai/ardg &
712 0 : * (alpha2/tanar+beta2/tanak*hkoverhr/dkoverdr)
713 0 : hkeel = hkoverhr * hridge
714 0 : 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 0 : tmp1 = max(c0,hridge - hfreebd)
721 :
722 : !------------------------------------------------------------
723 : ! Skin drag (atmo)
724 : !------------------------------------------------------------
725 :
726 0 : Cdn_atm_skin = csa*(c1 - mrdg*tmp1/distrdg)
727 0 : Cdn_atm_skin = max(min(Cdn_atm_skin,camax),c0)
728 :
729 : !------------------------------------------------------------
730 : ! Ridge effect (atmo)
731 : !------------------------------------------------------------
732 :
733 0 : if (tmp1 > puny) then
734 0 : sca = c1 - exp(-sHGB*distrdg/tmp1) ! see Eq. 9
735 0 : ctecar = cra*p5
736 : Cdn_atm_rdg = ctecar*tmp1/distrdg*sca* &
737 0 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
738 0 : 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 0 : tmp1 = max(c0,hkeel - hdraft)
748 :
749 : !------------------------------------------------------------
750 : ! Skin drag bottom ice (ocean)
751 : !------------------------------------------------------------
752 :
753 0 : Cdn_ocn_skin = csw * (c1 - mrdgo*tmp1/dkeel)
754 0 : Cdn_ocn_skin = max(min(Cdn_ocn_skin,cwmax), c0)
755 :
756 : !------------------------------------------------------------
757 : ! Keel effect (ocean)
758 : !------------------------------------------------------------
759 :
760 0 : if (tmp1 > puny) then
761 0 : scw = c1 - exp(-sHGB*dkeel/tmp1)
762 0 : ctecwk = crw*p5
763 : Cdn_ocn_keel = ctecwk*tmp1/dkeel*scw* &
764 0 : (log(tmp1*icerufi)/log(zref*icerufi))**c2
765 0 : 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 72392 : if (hfreebd > puny) then
775 72392 : sca = c1 - exp(-sl*beta*(c1-ai))
776 72392 : ctecaf = cfa*p5*(log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2*sca
777 72392 : Cdn_atm_floe = ctecaf * hfreebd / lfloe
778 72392 : Cdn_atm_floe = max(min(Cdn_atm_floe,camax),c0)
779 : endif
780 :
781 : !------------------------------------------------------------
782 : ! Pond edge effect (atmo)
783 : !------------------------------------------------------------
784 :
785 72392 : if (hfreebd > puny) then
786 72392 : sca = (apond)**(c1/(zref*beta))
787 72392 : lp = lpmin*(1-apond)+lpmax*apond
788 : Cdn_atm_pond = cpa*p5*sca*apond*hfreebd/lp &
789 72392 : * (log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2
790 72392 : Cdn_atm_pond = min(Cdn_atm_pond,camax)
791 : endif
792 :
793 : !------------------------------------------------------------
794 : ! Floe edge drag effect (ocean)
795 : !------------------------------------------------------------
796 :
797 72392 : if (hdraft > puny) then
798 72392 : scw = c1 - exp(-sl*beta*(c1-ai))
799 72392 : ctecwf = cfw*p5*(log(hdraft*ocnrufi)/log(zref*ocnrufi))**c2*scw
800 72392 : Cdn_ocn_floe = ctecwf * hdraft / lfloe
801 72392 : Cdn_ocn_floe = max(min(Cdn_ocn_floe,cwmax),c0)
802 : endif
803 :
804 : !------------------------------------------------------------
805 : ! Total drag coefficient (atmo)
806 : !------------------------------------------------------------
807 :
808 72392 : Cdn_atm = Cdn_atm_skin + Cdn_atm_floe + Cdn_atm_pond + Cdn_atm_rdg
809 72392 : Cdn_atm = min(Cdn_atm,camax)
810 :
811 : !------------------------------------------------------------
812 : ! Total drag coefficient (ocean)
813 : !------------------------------------------------------------
814 :
815 72392 : Cdn_ocn = Cdn_ocn_skin + Cdn_ocn_floe + Cdn_ocn_keel
816 72392 : Cdn_ocn = min(Cdn_ocn,cwmax)
817 :
818 : endif
819 :
820 96528 : end subroutine neutral_drag_coeffs
821 :
822 : !=======================================================================
823 : !autodocument_start icepack_atm_boundary
824 : !
825 :
826 10750054 : subroutine icepack_atm_boundary(sfctype, &
827 : Tsf, potT, &
828 : uatm, vatm, &
829 : wind, zlvl, &
830 : Qa, rhoa, &
831 : strx, stry, &
832 : Tref, Qref, &
833 : delt, delq, &
834 : lhcoef, shcoef, &
835 : Cdn_atm, &
836 : Cdn_atm_ratio_n, &
837 10750054 : Qa_iso, Qref_iso, &
838 : uvel, vvel, &
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
846 : potT , & ! air potential temperature (K)
847 : uatm , & ! x-direction wind speed (m/s)
848 : vatm , & ! y-direction wind speed (m/s)
849 : wind , & ! wind speed (m/s)
850 : zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m)
851 : Qa , & ! specific humidity (kg/kg)
852 : rhoa ! air density (kg/m^3)
853 :
854 : real (kind=dbl_kind), intent(inout) :: &
855 : Cdn_atm , & ! neutral drag coefficient
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)
860 : stry ! y surface stress (N)
861 :
862 : real (kind=dbl_kind), intent(inout) :: &
863 : Tref , & ! reference height temperature (K)
864 : Qref , & ! reference height specific humidity (kg/kg)
865 : delt , & ! potential T difference (K)
866 : delq , & ! humidity difference (kg/kg)
867 : shcoef , & ! transfer coefficient for sensible heat
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)
878 : vvel , & ! y-direction ice speed (m/s)
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 10750054 : if (sfctype == 'ice') then
902 8777938 : if (icepack_chkoptargflag(first_call_ice)) then
903 83 : if (tr_iso) then
904 2 : 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 10750054 : l_uvel = c0
914 10750054 : l_vvel = c0
915 10750054 : l_Uref = c0
916 10750054 : if (present(uvel)) then
917 8777938 : l_uvel = uvel
918 : endif
919 10750054 : if (present(vvel)) then
920 8777938 : l_vvel = vvel
921 : endif
922 :
923 10750054 : Cdn_atm_ratio_n = c1
924 :
925 10750054 : if (trim(atmbndy) == 'constant') then
926 : call atmo_boundary_const (sfctype, calc_strair, &
927 : uatm, vatm, &
928 : wind, rhoa, &
929 : strx, stry, &
930 : Tsf, potT, &
931 : Qa, &
932 : delt, delq, &
933 588714 : lhcoef, shcoef )
934 588714 : if (icepack_warnings_aborted(subname)) return
935 : else
936 : call atmo_boundary_layer (sfctype, &
937 : calc_strair, formdrag, &
938 : Tsf, potT, &
939 : uatm, vatm, &
940 : wind, zlvl, &
941 : Qa, rhoa, &
942 : strx, stry, &
943 : Tref, Qref, &
944 : delt, delq, &
945 : lhcoef, shcoef, &
946 : Cdn_atm, &
947 : Cdn_atm_ratio_n, &
948 : Qa_iso=Qa_iso, &
949 : Qref_iso=Qref_iso, &
950 : uvel=l_uvel, vvel=l_vvel,&
951 10161340 : Uref=l_Uref, zlvs=zlvs )
952 10161340 : if (icepack_warnings_aborted(subname)) return
953 : endif ! atmbndy
954 :
955 10750054 : if (present(Uref)) then
956 8777938 : Uref = l_Uref
957 : endif
958 :
959 10750054 : if (sfctype == 'ice') then
960 8777938 : first_call_ice = .false.
961 : endif
962 :
963 : end subroutine icepack_atm_boundary
964 :
965 : !=======================================================================
966 :
967 50625655 : function compute_stability_parameter(zlvl , thva , &
968 : ustar, tstar, &
969 : qstar, Qa) &
970 : result(hol)
971 :
972 : real (kind=dbl_kind), intent(in) :: &
973 : zlvl , & ! atm level height (m)
974 : thva , & ! virtual temperature (K)
975 : ustar , & ! turbulent scale for momentum
976 : tstar , & ! turbulent scale for temperature
977 : qstar , & ! turbulent scale for humidity
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 &
987 : + qstar/(c1/zvir+Qa)) &
988 50625655 : / ustar**2
989 50625655 : hol = sign( min(abs(hol),c10), hol)
990 :
991 50625655 : end function compute_stability_parameter
992 :
993 : !=======================================================================
994 :
995 101251310 : 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
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
1011 : psi_unstable ! unstable stability funcion at hol
1012 :
1013 : character(len=*),parameter :: subname='(compute_stability_function)'
1014 :
1015 101251310 : stable = p5 + sign(p5 , hol)
1016 :
1017 : psi_stable = -(0.7_dbl_kind*hol &
1018 : + 0.75_dbl_kind*(hol-14.3_dbl_kind) &
1019 101251310 : * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
1020 :
1021 101251310 : if(trim(qty) == 'momentum') then
1022 50625655 : psi_unstable = psi_momentum_unstable(hol)
1023 50625655 : elseif(trim(qty) == 'scalar') then
1024 50625655 : 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 101251310 : psi = psi_stable*stable + (c1 - stable)*psi_unstable
1031 :
1032 101251310 : end subroutine compute_stability_function
1033 :
1034 : !------------------------------------------------------------
1035 : ! Define functions
1036 : !------------------------------------------------------------
1037 :
1038 : !=======================================================================
1039 :
1040 50625655 : 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 50625655 : xd = capital_X(hol)
1047 :
1048 : psi_momentum_unstable = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
1049 50625655 : - c2*atan(xd) + pih
1050 :
1051 50625655 : end function psi_momentum_unstable
1052 :
1053 : !=======================================================================
1054 :
1055 60786995 : 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 60786995 : xd = capital_X(hol)
1062 :
1063 60786995 : psi_scalar_unstable = c2 * log((c1 + xd*xd)/c2)
1064 :
1065 60786995 : end function psi_scalar_unstable
1066 :
1067 : !=======================================================================
1068 :
1069 111412650 : real(kind=dbl_kind) function capital_X(hol)
1070 :
1071 : real(kind=dbl_kind), intent(in) :: hol
1072 :
1073 111412650 : capital_X = sqrt(max(sqrt(abs(c1 - c16*hol)) , c1))
1074 :
1075 111412650 : end function capital_X
1076 :
1077 : !=======================================================================
1078 :
1079 : end module icepack_atmo
1080 :
1081 : !=======================================================================
|
