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