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