Line data Source code
1 : !=======================================================================
2 : !
3 : ! The albedo and absorbed/transmitted flux parameterizations for
4 : ! snow over ice, bare ice and ponded ice.
5 : !
6 : ! Presently, two methods are included:
7 : ! (1) CCSM3
8 : ! (2) Delta-Eddington
9 : ! as two distinct routines.
10 : ! Either can be called from the ice driver.
11 : !
12 : ! The Delta-Eddington method is described here:
13 : !
14 : ! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple
15 : ! Scattering Parameterization for Solar Radiation in the Sea Ice
16 : ! Component of the Community Climate System Model, NCAR Technical
17 : ! Note NCAR/TN-472+STR February 2007
18 : !
19 : ! name: originally ice_albedo
20 : !
21 : ! authors: Bruce P. Briegleb, NCAR
22 : ! Elizabeth C. Hunke and William H. Lipscomb, LANL
23 : ! 2005, WHL: Moved absorbed_solar from icepack_therm_vertical to this
24 : ! module and changed name from ice_albedo
25 : ! 2006, WHL: Added Delta Eddington routines from Bruce Briegleb
26 : ! 2006, ECH: Changed data statements in Delta Eddington routines (no
27 : ! longer hardwired)
28 : ! Converted to free source form (F90)
29 : ! 2007, BPB: Completely updated Delta-Eddington code, so that:
30 : ! (1) multiple snow layers enabled (i.e. nslyr > 1)
31 : ! (2) included SSL for snow surface absorption
32 : ! (3) added Sswabs for internal snow layer absorption
33 : ! (4) variable sea ice layers allowed (i.e. not hardwired)
34 : ! (5) updated all inherent optical properties
35 : ! (6) included algae absorption for sea ice lowest layer
36 : ! (7) very complete internal documentation included
37 : ! 2007, ECH: Improved efficiency
38 : ! 2008, BPB: Added aerosols to Delta Eddington code
39 : ! 2013, ECH: merged with NCAR version, cleaned up
40 :
41 : module icepack_shortwave
42 :
43 : use icepack_kinds
44 : use icepack_parameters, only: c0, c1, c1p5, c2, c3, c4, c10
45 : use icepack_parameters, only: p01, p1, p15, p25, p5, p75, puny
46 : use icepack_parameters, only: albocn, Timelt, snowpatch, awtvdr, awtidr, awtvdf, awtidf
47 : use icepack_parameters, only: kappav, hs_min, rhofresh, rhos, nspint
48 : use icepack_parameters, only: hi_ssl, hs_ssl, min_bgc, sk_l
49 : use icepack_parameters, only: z_tracers, skl_bgc, calc_tsfc, shortwave, kalg, heat_capacity
50 : use icepack_parameters, only: r_ice, r_pnd, r_snw, dt_mlt, rsnw_mlt, hs0, hs1, hp1
51 : use icepack_parameters, only: pndaspect, albedo_type, albicev, albicei, albsnowv, albsnowi, ahmax
52 : use icepack_tracers, only: ntrcr, nbtrcr_sw
53 : use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
54 : use icepack_tracers, only: tr_bgc_N, tr_aero
55 : use icepack_tracers, only: nt_bgc_N, nt_zaero, tr_bgc_N
56 : use icepack_tracers, only: tr_zaero, nlt_chl_sw, nlt_zaero_sw
57 : use icepack_tracers, only: n_algae, n_aero, n_zaero
58 : use icepack_warnings, only: warnstr, icepack_warnings_add
59 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
60 :
61 : use icepack_zbgc_shared,only: R_chl2N, F_abs_chl
62 : use icepack_zbgc_shared,only: remap_zbgc
63 : use icepack_orbital, only: compute_coszen
64 :
65 :
66 : implicit none
67 :
68 : private
69 : public :: run_dEdd, &
70 : shortwave_ccsm3, &
71 : compute_shortwave_trcr, &
72 : icepack_prep_radiation, &
73 : icepack_step_radiation
74 :
75 : real (kind=dbl_kind), parameter :: &
76 : hpmin = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
77 : hp0 = 0.200_dbl_kind ! pond depth below which transition to bare ice
78 :
79 : real (kind=dbl_kind), parameter :: &
80 : exp_argmax = c10 ! maximum argument of exponential
81 :
82 : !=======================================================================
83 :
84 : contains
85 :
86 : !=======================================================================
87 : !
88 : ! Driver for basic solar radiation from CCSM3. Albedos and absorbed solar.
89 :
90 289620 : subroutine shortwave_ccsm3 (aicen, vicen, &
91 144810 : vsnon, Tsfcn, &
92 : swvdr, swvdf, &
93 : swidr, swidf, &
94 : heat_capacity, &
95 52566 : albedo_type, &
96 : albicev, albicei, &
97 : albsnowv, albsnowi, &
98 : ahmax, &
99 144810 : alvdrn, alidrn, &
100 144810 : alvdfn, alidfn, &
101 144810 : fswsfc, fswint, &
102 144810 : fswthru, fswpenl, &
103 289620 : Iswabs, SSwabs, &
104 144810 : albin, albsn, &
105 : coszen, ncat, &
106 : nilyr)
107 :
108 : integer (kind=int_kind), intent(in) :: &
109 : nilyr , & ! number of ice layers
110 : ncat ! number of ice thickness categories
111 :
112 : real (kind=dbl_kind), dimension (:), intent(in) :: &
113 : aicen , & ! concentration of ice per category
114 : vicen , & ! volume of ice per category
115 : vsnon , & ! volume of ice per category
116 : Tsfcn ! surface temperature
117 :
118 : real (kind=dbl_kind), intent(in) :: &
119 : swvdr , & ! sw down, visible, direct (W/m^2)
120 : swvdf , & ! sw down, visible, diffuse (W/m^2)
121 : swidr , & ! sw down, near IR, direct (W/m^2)
122 : swidf ! sw down, near IR, diffuse (W/m^2)
123 :
124 : ! baseline albedos for ccsm3 shortwave, set in namelist
125 : real (kind=dbl_kind), intent(in) :: &
126 : albicev , & ! visible ice albedo for h > ahmax
127 : albicei , & ! near-ir ice albedo for h > ahmax
128 : albsnowv, & ! cold snow albedo, visible
129 : albsnowi, & ! cold snow albedo, near IR
130 : ahmax ! thickness above which ice albedo is constant (m)
131 :
132 : logical(kind=log_kind), intent(in) :: &
133 : heat_capacity! if true, ice has nonzero heat capacity
134 :
135 : character (len=char_len), intent(in) :: &
136 : albedo_type ! albedo parameterization, 'ccsm3' or 'constant'
137 :
138 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
139 : alvdrn , & ! visible, direct, avg (fraction)
140 : alidrn , & ! near-ir, direct, avg (fraction)
141 : alvdfn , & ! visible, diffuse, avg (fraction)
142 : alidfn , & ! near-ir, diffuse, avg (fraction)
143 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
144 : fswint , & ! SW absorbed in ice interior, below surface (W m-2)
145 : fswthru , & ! SW through ice to ocean (W m-2)
146 : albin , & ! bare ice albedo
147 : albsn ! snow albedo
148 :
149 : real (kind=dbl_kind), intent(inout) :: &
150 : coszen ! cosine(zenith angle)
151 :
152 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
153 : fswpenl , & ! SW entering ice layers (W m-2)
154 : Iswabs , & ! SW absorbed in particular layer (W m-2)
155 : Sswabs ! SW absorbed in particular layer (W m-2)
156 :
157 : ! local variables
158 :
159 : integer (kind=int_kind) :: &
160 : n ! thickness category index
161 :
162 : ! ice and snow albedo for each category
163 :
164 : real (kind=dbl_kind) :: &
165 52566 : alvdrni, & ! visible, direct, ice (fraction)
166 52566 : alidrni, & ! near-ir, direct, ice (fraction)
167 52566 : alvdfni, & ! visible, diffuse, ice (fraction)
168 52566 : alidfni, & ! near-ir, diffuse, ice (fraction)
169 52566 : alvdrns, & ! visible, direct, snow (fraction)
170 52566 : alidrns, & ! near-ir, direct, snow (fraction)
171 52566 : alvdfns, & ! visible, diffuse, snow (fraction)
172 52566 : alidfns ! near-ir, diffuse, snow (fraction)
173 :
174 : character(len=*),parameter :: subname='(shortwave_ccsm3)'
175 :
176 : !-----------------------------------------------------------------
177 : ! Solar radiation: albedo and absorbed shortwave
178 : !-----------------------------------------------------------------
179 :
180 : ! For basic shortwave, set coszen to a constant between 0 and 1.
181 144810 : coszen = p5 ! sun above the horizon
182 :
183 868860 : do n = 1, ncat
184 :
185 1448100 : Sswabs(:,n) = c0
186 :
187 724050 : alvdrni = albocn
188 724050 : alidrni = albocn
189 724050 : alvdfni = albocn
190 724050 : alidfni = albocn
191 :
192 724050 : alvdrns = albocn
193 724050 : alidrns = albocn
194 724050 : alvdfns = albocn
195 724050 : alidfns = albocn
196 :
197 724050 : alvdrn(n) = albocn
198 724050 : alidrn(n) = albocn
199 724050 : alvdfn(n) = albocn
200 724050 : alidfn(n) = albocn
201 :
202 724050 : albin(n) = c0
203 724050 : albsn(n) = c0
204 :
205 724050 : fswsfc(n) = c0
206 724050 : fswint(n) = c0
207 724050 : fswthru(n) = c0
208 4344300 : fswpenl(:,n) = c0
209 3620250 : Iswabs (:,n) = c0
210 :
211 868860 : if (aicen(n) > puny) then
212 :
213 : !-----------------------------------------------------------------
214 : ! Compute albedos for ice and snow.
215 : !-----------------------------------------------------------------
216 :
217 591299 : if (trim(albedo_type) == 'constant') then
218 :
219 81020 : call constant_albedos (aicen(n), &
220 81020 : vsnon(n), &
221 81020 : Tsfcn(n), &
222 : alvdrni, alidrni, &
223 : alvdfni, alidfni, &
224 : alvdrns, alidrns, &
225 : alvdfns, alidfns, &
226 81020 : alvdrn(n), &
227 81020 : alidrn(n), &
228 81020 : alvdfn(n), &
229 81020 : alidfn(n), &
230 81020 : albin(n), &
231 230918 : albsn(n))
232 230918 : if (icepack_warnings_aborted(subname)) return
233 :
234 360381 : elseif (trim(albedo_type) == 'ccsm3') then
235 :
236 130625 : call compute_albedos (aicen(n), &
237 130625 : vicen(n), &
238 130625 : vsnon(n), &
239 130625 : Tsfcn(n), &
240 : albicev, albicei, &
241 : albsnowv, albsnowi, &
242 : ahmax, &
243 : alvdrni, alidrni, &
244 : alvdfni, alidfni, &
245 : alvdrns, alidrns, &
246 : alvdfns, alidfns, &
247 130625 : alvdrn(n), &
248 130625 : alidrn(n), &
249 130625 : alvdfn(n), &
250 130625 : alidfn(n), &
251 130625 : albin(n), &
252 360381 : albsn(n))
253 360381 : if (icepack_warnings_aborted(subname)) return
254 :
255 : else
256 :
257 0 : call icepack_warnings_add(subname//' ERROR: albedo_type '//trim(albedo_type)//' unknown')
258 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
259 0 : return
260 :
261 : endif
262 :
263 : !-----------------------------------------------------------------
264 : ! Compute solar radiation absorbed in ice and penetrating to ocean.
265 : !-----------------------------------------------------------------
266 :
267 : call absorbed_solar (heat_capacity, &
268 : nilyr, &
269 211645 : aicen(n), &
270 211645 : vicen(n), &
271 211645 : vsnon(n), &
272 : swvdr, swvdf, &
273 : swidr, swidf, &
274 : alvdrni, alvdfni, &
275 : alidrni, alidfni, &
276 : alvdrns, alvdfns, &
277 : alidrns, alidfns, &
278 211645 : fswsfc(n), &
279 211645 : fswint(n), &
280 211645 : fswthru(n), &
281 0 : fswpenl(:,n), &
282 591299 : Iswabs(:,n))
283 591299 : if (icepack_warnings_aborted(subname)) return
284 :
285 : endif ! aicen > puny
286 :
287 : enddo ! ncat
288 :
289 : end subroutine shortwave_ccsm3
290 :
291 : !=======================================================================
292 : !
293 : ! Compute albedos for each thickness category
294 :
295 360381 : subroutine compute_albedos (aicen, vicen, &
296 : vsnon, Tsfcn, &
297 : albicev, albicei, &
298 : albsnowv, albsnowi, &
299 : ahmax, &
300 : alvdrni, alidrni, &
301 : alvdfni, alidfni, &
302 : alvdrns, alidrns, &
303 : alvdfns, alidfns, &
304 : alvdrn, alidrn, &
305 : alvdfn, alidfn, &
306 : albin, albsn)
307 :
308 : real (kind=dbl_kind), intent(in) :: &
309 : aicen , & ! concentration of ice per category
310 : vicen , & ! volume of ice per category
311 : vsnon , & ! volume of ice per category
312 : Tsfcn ! surface temperature
313 :
314 : ! baseline albedos for ccsm3 shortwave, set in namelist
315 : real (kind=dbl_kind), intent(in) :: &
316 : albicev , & ! visible ice albedo for h > ahmax
317 : albicei , & ! near-ir ice albedo for h > ahmax
318 : albsnowv, & ! cold snow albedo, visible
319 : albsnowi, & ! cold snow albedo, near IR
320 : ahmax ! thickness above which ice albedo is constant (m)
321 :
322 : real (kind=dbl_kind), intent(out) :: &
323 : alvdrni , & ! visible, direct, ice (fraction)
324 : alidrni , & ! near-ir, direct, ice (fraction)
325 : alvdfni , & ! visible, diffuse, ice (fraction)
326 : alidfni , & ! near-ir, diffuse, ice (fraction)
327 : alvdrns , & ! visible, direct, snow (fraction)
328 : alidrns , & ! near-ir, direct, snow (fraction)
329 : alvdfns , & ! visible, diffuse, snow (fraction)
330 : alidfns , & ! near-ir, diffuse, snow (fraction)
331 : alvdrn , & ! visible, direct, avg (fraction)
332 : alidrn , & ! near-ir, direct, avg (fraction)
333 : alvdfn , & ! visible, diffuse, avg (fraction)
334 : alidfn , & ! near-ir, diffuse, avg (fraction)
335 : albin , & ! bare ice
336 : albsn ! snow
337 :
338 : ! local variables
339 :
340 : real (kind=dbl_kind), parameter :: &
341 : dT_melt = c1 , & ! change in temp to give dalb_mlt
342 : ! albedo change
343 : dalb_mlt = -0.075_dbl_kind, & ! albedo change per dT_melt change
344 : ! in temp for ice
345 : dalb_mltv = -p1 , & ! albedo vis change per dT_melt change
346 : ! in temp for snow
347 : dalb_mlti = -p15 ! albedo nir change per dT_melt change
348 : ! in temp for snow
349 :
350 : real (kind=dbl_kind) :: &
351 130625 : hi , & ! ice thickness (m)
352 130625 : hs , & ! snow thickness (m)
353 130625 : albo, & ! effective ocean albedo, function of ice thickness
354 130625 : fh , & ! piecewise linear function of thickness
355 130625 : fT , & ! piecewise linear function of surface temperature
356 130625 : dTs , & ! difference of Tsfc and Timelt
357 130625 : fhtan,& ! factor used in albedo dependence on ice thickness
358 130625 : asnow ! fractional area of snow cover
359 :
360 : character(len=*),parameter :: subname='(compute_albedos)'
361 :
362 360381 : fhtan = atan(ahmax*c4)
363 :
364 : !-----------------------------------------------------------------
365 : ! Compute albedo for each thickness category.
366 : !-----------------------------------------------------------------
367 :
368 360381 : hi = vicen / aicen
369 360381 : hs = vsnon / aicen
370 :
371 : ! bare ice, thickness dependence
372 360381 : fh = min(atan(hi*c4)/fhtan,c1)
373 360381 : albo = albocn*(c1-fh)
374 360381 : alvdfni = albicev*fh + albo
375 360381 : alidfni = albicei*fh + albo
376 :
377 : ! bare ice, temperature dependence
378 360381 : dTs = Timelt - Tsfcn
379 360381 : fT = min(dTs/dT_melt-c1,c0)
380 360381 : alvdfni = alvdfni - dalb_mlt*fT
381 360381 : alidfni = alidfni - dalb_mlt*fT
382 :
383 : ! avoid negative albedos for thin, bare, melting ice
384 360381 : alvdfni = max (alvdfni, albocn)
385 360381 : alidfni = max (alidfni, albocn)
386 :
387 360381 : if (hs > puny) then
388 :
389 341291 : alvdfns = albsnowv
390 341291 : alidfns = albsnowi
391 :
392 : ! snow on ice, temperature dependence
393 341291 : alvdfns = alvdfns - dalb_mltv*fT
394 341291 : alidfns = alidfns - dalb_mlti*fT
395 :
396 : endif ! hs > puny
397 :
398 : ! direct albedos (same as diffuse for now)
399 360381 : alvdrni = alvdfni
400 360381 : alidrni = alidfni
401 360381 : alvdrns = alvdfns
402 360381 : alidrns = alidfns
403 :
404 : ! fractional area of snow cover
405 360381 : if (hs > puny) then
406 341291 : asnow = hs / (hs + snowpatch)
407 : else
408 19090 : asnow = c0
409 : endif
410 :
411 : ! combine ice and snow albedos (for coupler)
412 : alvdfn = alvdfni*(c1-asnow) + &
413 360381 : alvdfns*asnow
414 : alidfn = alidfni*(c1-asnow) + &
415 360381 : alidfns*asnow
416 : alvdrn = alvdrni*(c1-asnow) + &
417 360381 : alvdrns*asnow
418 : alidrn = alidrni*(c1-asnow) + &
419 360381 : alidrns*asnow
420 :
421 : ! save ice and snow albedos (for history)
422 : albin = awtvdr*alvdrni + awtidr*alidrni &
423 360381 : + awtvdf*alvdfni + awtidf*alidfni
424 : albsn = awtvdr*alvdrns + awtidr*alidrns &
425 360381 : + awtvdf*alvdfns + awtidf*alidfns
426 :
427 360381 : end subroutine compute_albedos
428 :
429 : !=======================================================================
430 : !
431 : ! Compute albedos for each thickness category
432 :
433 230918 : subroutine constant_albedos (aicen, &
434 : vsnon, Tsfcn, &
435 : alvdrni, alidrni, &
436 : alvdfni, alidfni, &
437 : alvdrns, alidrns, &
438 : alvdfns, alidfns, &
439 : alvdrn, alidrn, &
440 : alvdfn, alidfn, &
441 : albin, albsn)
442 :
443 : real (kind=dbl_kind), intent(in) :: &
444 : aicen , & ! concentration of ice per category
445 : vsnon , & ! volume of ice per category
446 : Tsfcn ! surface temperature
447 :
448 : real (kind=dbl_kind), intent(out) :: &
449 : alvdrni , & ! visible, direct, ice (fraction)
450 : alidrni , & ! near-ir, direct, ice (fraction)
451 : alvdfni , & ! visible, diffuse, ice (fraction)
452 : alidfni , & ! near-ir, diffuse, ice (fraction)
453 : alvdrns , & ! visible, direct, snow (fraction)
454 : alidrns , & ! near-ir, direct, snow (fraction)
455 : alvdfns , & ! visible, diffuse, snow (fraction)
456 : alidfns , & ! near-ir, diffuse, snow (fraction)
457 : alvdrn , & ! visible, direct, avg (fraction)
458 : alidrn , & ! near-ir, direct, avg (fraction)
459 : alvdfn , & ! visible, diffuse, avg (fraction)
460 : alidfn , & ! near-ir, diffuse, avg (fraction)
461 : albin , & ! bare ice
462 : albsn ! snow
463 :
464 : ! local variables
465 :
466 : real (kind=dbl_kind), parameter :: &
467 : warmice = 0.68_dbl_kind, &
468 : coldice = 0.70_dbl_kind, &
469 : warmsnow = 0.77_dbl_kind, &
470 : coldsnow = 0.81_dbl_kind
471 :
472 : real (kind=dbl_kind) :: &
473 81020 : hs ! snow thickness (m)
474 :
475 : character(len=*),parameter :: subname='(constant_albedos)'
476 :
477 : !-----------------------------------------------------------------
478 : ! Compute albedo for each thickness category.
479 : !-----------------------------------------------------------------
480 :
481 230918 : hs = vsnon / aicen
482 :
483 230918 : if (hs > puny) then
484 : ! snow, temperature dependence
485 95564 : if (Tsfcn >= -c2*puny) then
486 0 : alvdfn = warmsnow
487 0 : alidfn = warmsnow
488 : else
489 95564 : alvdfn = coldsnow
490 95564 : alidfn = coldsnow
491 : endif
492 : else ! hs < puny
493 : ! bare ice, temperature dependence
494 135354 : if (Tsfcn >= -c2*puny) then
495 0 : alvdfn = warmice
496 0 : alidfn = warmice
497 : else
498 135354 : alvdfn = coldice
499 135354 : alidfn = coldice
500 : endif
501 : endif ! hs > puny
502 :
503 : ! direct albedos (same as diffuse for now)
504 230918 : alvdrn = alvdfn
505 230918 : alidrn = alidfn
506 :
507 230918 : alvdrni = alvdrn
508 230918 : alidrni = alidrn
509 230918 : alvdrns = alvdrn
510 230918 : alidrns = alidrn
511 230918 : alvdfni = alvdfn
512 230918 : alidfni = alidfn
513 230918 : alvdfns = alvdfn
514 230918 : alidfns = alidfn
515 :
516 : ! save ice and snow albedos (for history)
517 : albin = awtvdr*alvdrni + awtidr*alidrni &
518 230918 : + awtvdf*alvdfni + awtidf*alidfni
519 : albsn = awtvdr*alvdrns + awtidr*alidrns &
520 230918 : + awtvdf*alvdfns + awtidf*alidfns
521 :
522 230918 : end subroutine constant_albedos
523 :
524 : !=======================================================================
525 : !
526 : ! Compute solar radiation absorbed in ice and penetrating to ocean
527 : !
528 : ! authors William H. Lipscomb, LANL
529 : ! C. M. Bitz, UW
530 :
531 591299 : subroutine absorbed_solar (heat_capacity, &
532 : nilyr, aicen, &
533 : vicen, vsnon, &
534 : swvdr, swvdf, &
535 : swidr, swidf, &
536 : alvdrni, alvdfni, &
537 : alidrni, alidfni, &
538 : alvdrns, alvdfns, &
539 : alidrns, alidfns, &
540 : fswsfc, fswint, &
541 591299 : fswthru, fswpenl, &
542 591299 : Iswabs)
543 :
544 : logical(kind=log_kind), intent(in) :: &
545 : heat_capacity ! if true, ice has nonzero heat capacity
546 :
547 : integer (kind=int_kind), intent(in) :: &
548 : nilyr ! number of ice layers
549 :
550 : real (kind=dbl_kind), intent(in) :: &
551 : aicen , & ! fractional ice area
552 : vicen , & ! ice volume
553 : vsnon , & ! snow volume
554 : swvdr , & ! sw down, visible, direct (W/m^2)
555 : swvdf , & ! sw down, visible, diffuse (W/m^2)
556 : swidr , & ! sw down, near IR, direct (W/m^2)
557 : swidf , & ! sw down, near IR, diffuse (W/m^2)
558 : alvdrni , & ! visible, direct albedo,ice
559 : alidrni , & ! near-ir, direct albedo,ice
560 : alvdfni , & ! visible, diffuse albedo,ice
561 : alidfni , & ! near-ir, diffuse albedo,ice
562 : alvdrns , & ! visible, direct albedo, snow
563 : alidrns , & ! near-ir, direct albedo, snow
564 : alvdfns , & ! visible, diffuse albedo, snow
565 : alidfns ! near-ir, diffuse albedo, snow
566 :
567 : real (kind=dbl_kind), intent(out):: &
568 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
569 : fswint , & ! SW absorbed in ice interior, below surface (W m-2)
570 : fswthru ! SW through ice to ocean (W m-2)
571 :
572 : real (kind=dbl_kind), dimension (:), intent(out) :: &
573 : Iswabs , & ! SW absorbed in particular layer (W m-2)
574 : fswpenl ! visible SW entering ice layers (W m-2)
575 :
576 : ! local variables
577 :
578 : real (kind=dbl_kind), parameter :: &
579 : i0vis = 0.70_dbl_kind ! fraction of penetrating solar rad (visible)
580 :
581 : integer (kind=int_kind) :: &
582 : k ! ice layer index
583 :
584 : real (kind=dbl_kind) :: &
585 211645 : fswpen , & ! SW penetrating beneath surface (W m-2)
586 211645 : trantop , & ! transmitted frac of penetrating SW at layer top
587 211645 : tranbot ! transmitted frac of penetrating SW at layer bot
588 :
589 : real (kind=dbl_kind) :: &
590 211645 : swabs , & ! net SW down at surface (W m-2)
591 211645 : swabsv , & ! swabs in vis (wvlngth < 700nm) (W/m^2)
592 211645 : swabsi , & ! swabs in nir (wvlngth > 700nm) (W/m^2)
593 211645 : fswpenvdr , & ! penetrating SW, vis direct
594 211645 : fswpenvdf , & ! penetrating SW, vis diffuse
595 211645 : hi , & ! ice thickness (m)
596 211645 : hs , & ! snow thickness (m)
597 211645 : hilyr , & ! ice layer thickness
598 211645 : asnow ! fractional area of snow cover
599 :
600 : character(len=*),parameter :: subname='(absorbed_solar)'
601 :
602 : !-----------------------------------------------------------------
603 : ! Initialize
604 : !-----------------------------------------------------------------
605 :
606 591299 : trantop = c0
607 591299 : tranbot = c0
608 :
609 591299 : hs = vsnon / aicen
610 :
611 : !-----------------------------------------------------------------
612 : ! Fractional snow cover
613 : !-----------------------------------------------------------------
614 591299 : if (hs > puny) then
615 436855 : asnow = hs / (hs + snowpatch)
616 : else
617 154444 : asnow = c0
618 : endif
619 :
620 : !-----------------------------------------------------------------
621 : ! Shortwave flux absorbed at surface, absorbed internally,
622 : ! and penetrating to mixed layer.
623 : ! This parameterization assumes that all IR is absorbed at the
624 : ! surface; only visible is absorbed in the ice interior or
625 : ! transmitted to the ocean.
626 : !-----------------------------------------------------------------
627 :
628 : swabsv = swvdr * ( (c1-alvdrni)*(c1-asnow) &
629 : + (c1-alvdrns)*asnow ) &
630 : + swvdf * ( (c1-alvdfni)*(c1-asnow) &
631 591299 : + (c1-alvdfns)*asnow )
632 :
633 : swabsi = swidr * ( (c1-alidrni)*(c1-asnow) &
634 : + (c1-alidrns)*asnow ) &
635 : + swidf * ( (c1-alidfni)*(c1-asnow) &
636 591299 : + (c1-alidfns)*asnow )
637 :
638 591299 : swabs = swabsv + swabsi
639 :
640 591299 : fswpenvdr = swvdr * (c1-alvdrni) * (c1-asnow) * i0vis
641 591299 : fswpenvdf = swvdf * (c1-alvdfni) * (c1-asnow) * i0vis
642 :
643 : ! no penetrating radiation in near IR
644 : ! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir
645 : ! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir
646 :
647 591299 : fswpen = fswpenvdr + fswpenvdf
648 :
649 591299 : fswsfc = swabs - fswpen
650 :
651 591299 : trantop = c1 ! transmittance at top of ice
652 :
653 : !-----------------------------------------------------------------
654 : ! penetrating SW absorbed in each ice layer
655 : !-----------------------------------------------------------------
656 :
657 3344884 : do k = 1, nilyr
658 :
659 2753585 : hi = vicen / aicen
660 2753585 : hilyr = hi / real(nilyr,kind=dbl_kind)
661 :
662 2753585 : tranbot = exp (-kappav * hilyr * real(k,kind=dbl_kind))
663 2753585 : Iswabs(k) = fswpen * (trantop-tranbot)
664 :
665 : ! bottom of layer k = top of layer k+1
666 2753585 : trantop = tranbot
667 :
668 : ! bgc layer model
669 3344884 : if (k == 1) then ! surface flux
670 591299 : fswpenl(k) = fswpen
671 591299 : fswpenl(k+1) = fswpen * tranbot
672 : else
673 2162286 : fswpenl(k+1) = fswpen * tranbot
674 : endif
675 : enddo ! nilyr
676 :
677 : ! SW penetrating thru ice into ocean
678 591299 : fswthru = fswpen * tranbot
679 :
680 : ! SW absorbed in ice interior
681 591299 : fswint = fswpen - fswthru
682 :
683 : !----------------------------------------------------------------
684 : ! if zero-layer model (no heat capacity), no SW is absorbed in ice
685 : ! interior, so add to surface absorption
686 : !----------------------------------------------------------------
687 :
688 591299 : if (.not. heat_capacity) then
689 :
690 : ! SW absorbed at snow/ice surface
691 230918 : fswsfc = fswsfc + fswint
692 :
693 : ! SW absorbed in ice interior (nilyr = 1)
694 230918 : fswint = c0
695 230918 : Iswabs(1) = c0
696 :
697 : endif ! heat_capacity
698 :
699 591299 : end subroutine absorbed_solar
700 :
701 : ! End ccsm3 shortwave method
702 : !=======================================================================
703 : ! Begin Delta-Eddington shortwave method
704 :
705 : ! Compute initial data for Delta-Eddington method, specifically,
706 : ! the approximate exponential look-up table.
707 : !
708 : ! author: Bruce P. Briegleb, NCAR
709 : ! 2011 ECH modified for melt pond tracers
710 : ! 2013 ECH merged with NCAR version
711 :
712 882945 : subroutine run_dEdd(dt, ncat, &
713 : dEdd_algae, &
714 : nilyr, nslyr, &
715 1765890 : aicen, vicen, &
716 882945 : vsnon, Tsfcn, &
717 1765890 : alvln, apndn, &
718 1765890 : hpndn, ipndn, &
719 882945 : aeron, kalg, &
720 882945 : trcrn_bgcsw, &
721 : heat_capacity, &
722 : tlat, tlon, &
723 0 : calendar_type, &
724 : days_per_year, &
725 : nextsw_cday, yday, &
726 : sec, R_ice, &
727 : R_pnd, R_snw, &
728 : dT_mlt, rsnw_mlt, &
729 : hs0, hs1, hp1, &
730 : pndaspect, &
731 882945 : kaer_tab, waer_tab, &
732 882945 : gaer_tab, &
733 882945 : kaer_bc_tab, &
734 882945 : waer_bc_tab, &
735 882945 : gaer_bc_tab, &
736 882945 : bcenh, &
737 : modal_aero, &
738 : swvdr, swvdf, &
739 : swidr, swidf, &
740 : coszen, fsnow, &
741 882945 : alvdrn, alvdfn, &
742 882945 : alidrn, alidfn, &
743 882945 : fswsfcn, fswintn, &
744 1765890 : fswthrun, fswpenln, &
745 882945 : Sswabsn, Iswabsn, &
746 882945 : albicen, albsnon, &
747 1765890 : albpndn, apeffn, &
748 882945 : snowfracn, &
749 882945 : dhsn, ffracn, &
750 : l_print_point, &
751 : initonly)
752 :
753 : integer (kind=int_kind), intent(in) :: &
754 : ncat , & ! number of ice thickness categories
755 : nilyr , & ! number of ice layers
756 : nslyr ! number of snow layers
757 :
758 : logical(kind=log_kind), intent(in) :: &
759 : heat_capacity,& ! if true, ice has nonzero heat capacity
760 : dEdd_algae, & ! .true. use prognostic chla in dEdd
761 : modal_aero ! .true. use modal aerosol treatment
762 :
763 : ! dEdd tuning parameters, set in namelist
764 : real (kind=dbl_kind), intent(in) :: &
765 : R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
766 : R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
767 : R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
768 : dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
769 : rsnw_mlt, & ! maximum melting snow grain radius (10^-6 m)
770 : hs0 , & ! snow depth for transition to bare sea ice (m)
771 : pndaspect, & ! ratio of pond depth to pond fraction
772 : hs1 , & ! tapering parameter for snow on pond ice
773 : hp1 , & ! critical parameter for pond ice thickness
774 : kalg ! algae absorption coefficient
775 :
776 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
777 : kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
778 : waer_tab, & ! aerosol single scatter albedo (fraction)
779 : gaer_tab ! aerosol asymmetry parameter (cos(theta))
780 :
781 : real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
782 : kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
783 : waer_bc_tab, & ! aerosol single scatter albedo (fraction)
784 : gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
785 :
786 : real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
787 : bcenh ! BC absorption enhancement factor
788 :
789 : character (len=char_len), intent(in) :: &
790 : calendar_type ! differentiates Gregorian from other calendars
791 :
792 : integer (kind=int_kind), intent(in) :: &
793 : days_per_year, & ! number of days in one year
794 : sec ! elapsed seconds into date
795 :
796 : real (kind=dbl_kind), intent(in) :: &
797 : nextsw_cday , & ! julian day of next shortwave calculation
798 : yday ! day of the year
799 :
800 : real(kind=dbl_kind), intent(in) :: &
801 : dt, & ! time step (s)
802 : tlat, & ! latitude of temp pts (radians)
803 : tlon, & ! longitude of temp pts (radians)
804 : swvdr, & ! sw down, visible, direct (W/m^2)
805 : swvdf, & ! sw down, visible, diffuse (W/m^2)
806 : swidr, & ! sw down, near IR, direct (W/m^2)
807 : swidf, & ! sw down, near IR, diffuse (W/m^2)
808 : fsnow ! snowfall rate (kg/m^2 s)
809 :
810 : real(kind=dbl_kind), dimension(:), intent(in) :: &
811 : aicen, & ! concentration of ice
812 : vicen, & ! volume per unit area of ice (m)
813 : vsnon, & ! volume per unit area of snow (m)
814 : Tsfcn, & ! surface temperature (deg C)
815 : alvln, & ! level-ice area fraction
816 : apndn, & ! pond area fraction
817 : hpndn, & ! pond depth (m)
818 : ipndn ! pond refrozen lid thickness (m)
819 :
820 : real(kind=dbl_kind), dimension(:,:), intent(in) :: &
821 : aeron, & ! aerosols (kg/m^3)
822 : trcrn_bgcsw ! zaerosols (kg/m^3) + chlorophyll on shorthwave grid
823 :
824 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
825 : ffracn,& ! fraction of fsurfn used to melt ipond
826 : dhsn ! depth difference for snow on sea ice and pond ice
827 :
828 : real(kind=dbl_kind), intent(inout) :: &
829 : coszen ! cosine solar zenith angle, < 0 for sun below horizon
830 :
831 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
832 : alvdrn, & ! visible direct albedo (fraction)
833 : alvdfn, & ! near-ir direct albedo (fraction)
834 : alidrn, & ! visible diffuse albedo (fraction)
835 : alidfn, & ! near-ir diffuse albedo (fraction)
836 : fswsfcn, & ! SW absorbed at ice/snow surface (W m-2)
837 : fswintn, & ! SW absorbed in ice interior, below surface (W m-2)
838 : fswthrun, & ! SW through ice to ocean (W/m^2)
839 : albicen, & ! albedo bare ice
840 : albsnon, & ! albedo snow
841 : albpndn, & ! albedo pond
842 : apeffn, & ! effective pond area used for radiation calculation
843 : snowfracn ! snow fraction on each category used for radiation
844 :
845 : real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
846 : Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
847 : Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
848 : fswpenln ! visible SW entering ice layers (W m-2)
849 :
850 : logical (kind=log_kind), intent(in) :: &
851 : l_print_point
852 :
853 : logical (kind=log_kind), optional :: &
854 : initonly ! flag to indicate init only, default is false
855 :
856 : ! local temporary variables
857 :
858 : ! other local variables
859 : ! snow variables for Delta-Eddington shortwave
860 : real (kind=dbl_kind) :: &
861 303198 : fsn , & ! snow horizontal fraction
862 303198 : hsn ! snow depth (m)
863 :
864 : real (kind=dbl_kind), dimension (nslyr) :: &
865 2477418 : rhosnwn , & ! snow density (kg/m3)
866 1897671 : rsnwn ! snow grain radius (micrometers)
867 :
868 : ! pond variables for Delta-Eddington shortwave
869 : real (kind=dbl_kind) :: &
870 303198 : fpn , & ! pond fraction of ice cover
871 303198 : hpn ! actual pond depth (m)
872 :
873 : integer (kind=int_kind) :: &
874 : n ! thickness category index
875 :
876 : real (kind=dbl_kind) :: &
877 303198 : ipn , & ! refrozen pond ice thickness (m), mean over ice fraction
878 303198 : hp , & ! pond depth
879 303198 : hs , & ! snow depth
880 303198 : asnow , & ! fractional area of snow cover
881 303198 : rp , & ! volume fraction of retained melt water to total liquid content
882 303198 : hmx , & ! maximum available snow infiltration equivalent depth
883 303198 : dhs , & ! local difference in snow depth on sea ice and pond ice
884 303198 : spn , & ! snow depth on refrozen pond (m)
885 303198 : tmp ! 0 or 1
886 :
887 : logical (kind=log_kind) :: &
888 : linitonly ! local initonly value
889 :
890 : character(len=*),parameter :: subname='(run_dEdd)'
891 :
892 882945 : linitonly = .false.
893 882945 : if (present(initonly)) then
894 882945 : linitonly = initonly
895 : endif
896 :
897 : ! cosine of the zenith angle
898 : #ifdef CESMCOUPLED
899 : call compute_coszen (tlat, tlon, &
900 : yday, sec, coszen, &
901 : days_per_year, nextsw_cday, calendar_type)
902 : #else
903 : call compute_coszen (tlat, tlon, &
904 882945 : yday, sec, coszen)
905 : #endif
906 882945 : if (icepack_warnings_aborted(subname)) return
907 :
908 5008050 : do n = 1, ncat
909 :
910 : ! note that rhoswn, rsnw, fp, hp and Sswabs ARE NOT dimensioned with ncat
911 : ! BPB 19 Dec 2006
912 :
913 : ! set snow properties
914 4125105 : fsn = c0
915 4125105 : hsn = c0
916 9698310 : rhosnwn(:) = c0
917 9698310 : rsnwn(:) = c0
918 4125105 : apeffn(n) = c0 ! for history
919 4125105 : snowfracn(n) = c0 ! for history
920 :
921 5008050 : if (aicen(n) > puny) then
922 :
923 : call shortwave_dEdd_set_snow(nslyr, R_snw, &
924 : dT_mlt, rsnw_mlt, &
925 1275244 : aicen(n), vsnon(n), &
926 1275244 : Tsfcn(n), fsn, &
927 : hs0, hsn, &
928 3770075 : rhosnwn, rsnwn)
929 3770075 : if (icepack_warnings_aborted(subname)) return
930 :
931 : ! set pond properties
932 3770075 : if (tr_pond_cesm) then
933 : ! fraction of ice area
934 229795 : fpn = apndn(n)
935 : ! pond depth over fraction fpn
936 229795 : hpn = hpndn(n)
937 : ! snow infiltration
938 229795 : if (hsn >= hs_min .and. hs0 > puny) then
939 0 : asnow = min(hsn/hs0, c1) ! delta-Eddington formulation
940 0 : fpn = (c1 - asnow) * fpn
941 0 : hpn = pndaspect * fpn
942 : endif
943 : ! Zero out fraction of thin ponds for radiation only
944 229795 : if (hpn < hpmin) fpn = c0
945 229795 : fsn = min(fsn, c1-fpn)
946 229795 : apeffn(n) = fpn ! for history
947 3540280 : elseif (tr_pond_lvl) then
948 3010800 : fpn = c0 ! fraction of ice covered in pond
949 3010800 : hpn = c0 ! pond depth over fpn
950 : ! refrozen pond lid thickness avg over ice
951 : ! allow snow to cover pond ice
952 3010800 : ipn = alvln(n) * apndn(n) * ipndn(n)
953 3010800 : dhs = dhsn(n) ! snow depth difference, sea ice - pond
954 : if (.not. linitonly .and. ipn > puny .and. &
955 3010800 : dhs < puny .and. fsnow*dt > hs_min) &
956 10898 : dhs = hsn - fsnow*dt ! initialize dhs>0
957 3010800 : spn = hsn - dhs ! snow depth on pond ice
958 3010800 : if (.not. linitonly .and. ipn*spn < puny) dhs = c0
959 3010800 : dhsn(n) = dhs ! save: constant until reset to 0
960 :
961 : ! not using ipn assumes that lid ice is perfectly clear
962 : ! if (ipn <= 0.3_dbl_kind) then
963 :
964 : ! fraction of ice area
965 3010800 : fpn = apndn(n) * alvln(n)
966 : ! pond depth over fraction fpn
967 3010800 : hpn = hpndn(n)
968 :
969 : ! reduce effective pond area absorbing surface heat flux
970 : ! due to flux already having been used to melt pond ice
971 3010800 : fpn = (c1 - ffracn(n)) * fpn
972 :
973 : ! taper pond area with snow on pond ice
974 3010800 : if (dhs > puny .and. spn >= puny .and. hs1 > puny) then
975 18591 : asnow = min(spn/hs1, c1)
976 18591 : fpn = (c1 - asnow) * fpn
977 : endif
978 :
979 : ! infiltrate snow
980 3010800 : hp = hpn
981 3010800 : if (hp > puny) then
982 1433285 : hs = hsn
983 1433285 : rp = rhofresh*hp/(rhofresh*hp + rhos*hs)
984 1433285 : if (rp < p15) then
985 1146071 : fpn = c0
986 1146071 : hpn = c0
987 : else
988 287214 : hmx = hs*(rhofresh - rhos)/rhofresh
989 287214 : tmp = max(c0, sign(c1, hp-hmx)) ! 1 if hp>=hmx, else 0
990 : hp = (rhofresh*hp + rhos*hs*tmp) &
991 287214 : / (rhofresh - rhos*(c1-tmp))
992 287214 : hsn = hs - hp*fpn*(c1-tmp)
993 287214 : hpn = hp * tmp
994 287214 : fpn = fpn * tmp
995 : endif
996 : endif ! hp > puny
997 :
998 : ! Zero out fraction of thin ponds for radiation only
999 3010800 : if (hpn < hpmin) fpn = c0
1000 3010800 : fsn = min(fsn, c1-fpn)
1001 :
1002 : ! endif ! masking by lid ice
1003 3010800 : apeffn(n) = fpn ! for history
1004 :
1005 529480 : elseif (tr_pond_topo) then
1006 : ! Lid effective if thicker than hp1
1007 457081 : if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
1008 148413 : fpn = apndn(n)
1009 : else
1010 308668 : fpn = c0
1011 : endif
1012 457081 : if (apndn(n) > puny) then
1013 184923 : hpn = hpndn(n)
1014 : else
1015 272158 : fpn = c0
1016 272158 : hpn = c0
1017 : endif
1018 :
1019 : ! Zero out fraction of thin ponds for radiation only
1020 457081 : if (hpn < hpmin) fpn = c0
1021 :
1022 : ! If ponds are present snow fraction reduced to
1023 : ! non-ponded part dEdd scheme
1024 457081 : fsn = min(fsn, c1-fpn)
1025 :
1026 457081 : apeffn(n) = fpn
1027 : else
1028 72399 : fpn = c0
1029 72399 : hpn = c0
1030 26280 : call shortwave_dEdd_set_pond(Tsfcn(n), &
1031 : fsn, fpn, &
1032 72399 : hpn)
1033 72399 : if (icepack_warnings_aborted(subname)) return
1034 :
1035 72399 : apeffn(n) = fpn ! for history
1036 72399 : fpn = c0
1037 72399 : hpn = c0
1038 : endif ! pond type
1039 :
1040 3770075 : snowfracn(n) = fsn ! for history
1041 :
1042 : call shortwave_dEdd(dEdd_algae, &
1043 : nslyr, nilyr, &
1044 : coszen, heat_capacity, &
1045 1275244 : aicen(n), vicen(n), &
1046 : hsn, fsn, &
1047 0 : rhosnwn, rsnwn, &
1048 : fpn, hpn, &
1049 0 : aeron(:,n), &
1050 : R_ice, R_pnd, &
1051 0 : kaer_tab, waer_tab, &
1052 0 : gaer_tab, &
1053 0 : kaer_bc_tab, &
1054 0 : waer_bc_tab, &
1055 0 : gaer_bc_tab, &
1056 0 : bcenh, modal_aero, &
1057 : kalg, &
1058 : swvdr, swvdf, &
1059 : swidr, swidf, &
1060 1275244 : alvdrn(n), alvdfn(n), &
1061 1275244 : alidrn(n), alidfn(n), &
1062 1275244 : fswsfcn(n), fswintn(n), &
1063 1275244 : fswthrun(n), &
1064 0 : Sswabsn(:,n), &
1065 0 : Iswabsn(:,n), &
1066 1275244 : albicen(n), &
1067 1275244 : albsnon(n), albpndn(n), &
1068 0 : fswpenln(:,n), &
1069 0 : trcrn_bgcsw(:,n), &
1070 6320563 : l_print_point)
1071 3770075 : if (icepack_warnings_aborted(subname)) return
1072 :
1073 : endif ! aicen > puny
1074 :
1075 : enddo ! ncat
1076 :
1077 : end subroutine run_dEdd
1078 :
1079 : !=======================================================================
1080 : !
1081 : ! Compute snow/bare ice/ponded ice shortwave albedos, absorbed and transmitted
1082 : ! flux using the Delta-Eddington solar radiation method as described in:
1083 : !
1084 : ! A Delta-Eddington Multiple Scattering Parameterization for Solar Radiation
1085 : ! in the Sea Ice Component of the Community Climate System Model
1086 : ! B.P.Briegleb and B.Light NCAR/TN-472+STR February 2007
1087 : !
1088 : ! Compute shortwave albedos and fluxes for three surface types:
1089 : ! snow over ice, bare ice and ponded ice.
1090 : !
1091 : ! Albedos and fluxes are output for later use by thermodynamic routines.
1092 : ! Invokes three calls to compute_dEdd, which sets inherent optical properties
1093 : ! appropriate for the surface type. Within compute_dEdd, a call to solution_dEdd
1094 : ! evaluates the Delta-Eddington solution. The final albedos and fluxes are then
1095 : ! evaluated in compute_dEdd. Albedos and fluxes are transferred to output in
1096 : ! this routine.
1097 : !
1098 : ! NOTE regarding albedo diagnostics: This method yields zero albedo values
1099 : ! if there is no incoming solar and thus the albedo diagnostics are masked
1100 : ! out when the sun is below the horizon. To estimate albedo from the history
1101 : ! output (post-processing), compute ice albedo using
1102 : ! (1 - albedo)*swdn = swabs. -ECH
1103 : !
1104 : ! author: Bruce P. Briegleb, NCAR
1105 : ! 2013: E Hunke merged with NCAR version
1106 : !
1107 3770075 : subroutine shortwave_dEdd (dEdd_algae, &
1108 : nslyr, nilyr, &
1109 : coszen, heat_capacity,&
1110 : aice, vice, &
1111 : hs, fs, &
1112 3770075 : rhosnw, rsnw, &
1113 : fp, hp, &
1114 3770075 : aero, &
1115 : R_ice, R_pnd, &
1116 3770075 : kaer_tab, waer_tab, &
1117 3770075 : gaer_tab, &
1118 3770075 : kaer_bc_tab, &
1119 3770075 : waer_bc_tab, &
1120 3770075 : gaer_bc_tab, &
1121 3770075 : bcenh, modal_aero, &
1122 : kalg, &
1123 : swvdr, swvdf, &
1124 : swidr, swidf, &
1125 : alvdr, alvdf, &
1126 : alidr, alidf, &
1127 : fswsfc, fswint, &
1128 3770075 : fswthru, Sswabs, &
1129 3770075 : Iswabs, albice, &
1130 : albsno, albpnd, &
1131 7540150 : fswpenl, zbio, &
1132 : l_print_point)
1133 :
1134 : integer (kind=int_kind), intent(in) :: &
1135 : nilyr , & ! number of ice layers
1136 : nslyr ! number of snow layers
1137 :
1138 : logical (kind=log_kind), intent(in) :: &
1139 : heat_capacity, & ! if true, ice has nonzero heat capacity
1140 : dEdd_algae, & ! .true. use prognostic chla in dEdd
1141 : modal_aero ! .true. use modal aerosol treatment
1142 :
1143 : real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
1144 : kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
1145 : waer_bc_tab, & ! aerosol single scatter albedo (fraction)
1146 : gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
1147 :
1148 : real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
1149 : bcenh ! BC absorption enhancement factor
1150 :
1151 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
1152 : kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
1153 : waer_tab, & ! aerosol single scatter albedo (fraction)
1154 : gaer_tab ! aerosol asymmetry parameter (cos(theta))
1155 :
1156 : real (kind=dbl_kind), intent(in) :: &
1157 : kalg , & ! algae absorption coefficient
1158 : R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
1159 : R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
1160 : aice , & ! concentration of ice
1161 : vice , & ! volume of ice
1162 : hs , & ! snow depth
1163 : fs ! horizontal coverage of snow
1164 :
1165 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1166 : rhosnw , & ! density in snow layer (kg/m3)
1167 : rsnw , & ! grain radius in snow layer (m)
1168 : aero , & ! aerosol tracers
1169 : zbio ! shortwave tracers (zaero+chla)
1170 :
1171 : real (kind=dbl_kind), intent(in) :: &
1172 : fp , & ! pond fractional coverage (0 to 1)
1173 : hp , & ! pond depth (m)
1174 : swvdr , & ! sw down, visible, direct (W/m^2)
1175 : swvdf , & ! sw down, visible, diffuse (W/m^2)
1176 : swidr , & ! sw down, near IR, direct (W/m^2)
1177 : swidf ! sw down, near IR, diffuse (W/m^2)
1178 :
1179 : real (kind=dbl_kind), intent(inout) :: &
1180 : coszen , & ! cosine of solar zenith angle
1181 : alvdr , & ! visible, direct, albedo (fraction)
1182 : alvdf , & ! visible, diffuse, albedo (fraction)
1183 : alidr , & ! near-ir, direct, albedo (fraction)
1184 : alidf , & ! near-ir, diffuse, albedo (fraction)
1185 : fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
1186 : fswint , & ! SW interior absorption (below surface, above ocean,W m-2)
1187 : fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2)
1188 :
1189 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1190 : fswpenl , & ! visible SW entering ice layers (W m-2)
1191 : Sswabs , & ! SW absorbed in snow layer (W m-2)
1192 : Iswabs ! SW absorbed in ice layer (W m-2)
1193 :
1194 : real (kind=dbl_kind), intent(out) :: &
1195 : albice , & ! bare ice albedo, for history
1196 : albsno , & ! snow albedo, for history
1197 : albpnd ! pond albedo, for history
1198 :
1199 : logical (kind=log_kind) , intent(in) :: &
1200 : l_print_point
1201 :
1202 : ! local variables
1203 :
1204 : real (kind=dbl_kind) :: &
1205 1275244 : netsw , & ! net shortwave
1206 1275244 : fnidr , & ! fraction of direct to total down surface flux in nir
1207 1275244 : hstmp , & ! snow thickness (set to 0 for bare ice case)
1208 1275244 : hi , & ! ice thickness (all sea ice layers, m)
1209 1275244 : fi ! snow/bare ice fractional coverage (0 to 1)
1210 :
1211 : real (kind=dbl_kind), dimension (4*n_aero) :: &
1212 9920831 : aero_mp ! aerosol mass path in kg/m2
1213 :
1214 : integer (kind=int_kind) :: &
1215 : srftyp ! surface type over ice: (0=air, 1=snow, 2=pond)
1216 :
1217 : integer (kind=int_kind) :: &
1218 : k , & ! level index
1219 : na , & ! aerosol index
1220 : klev , & ! number of radiation layers - 1
1221 : klevp ! number of radiation interfaces - 1
1222 : ! (0 layer is included also)
1223 :
1224 : real (kind=dbl_kind) :: &
1225 1275244 : vsno ! volume of snow
1226 :
1227 : real (kind=dbl_kind) :: &
1228 1275244 : swdn , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j)
1229 1275244 : swab , & ! fswsfc(i,j)+fswint(i,j)+fswthru(i,j)
1230 1275244 : swalb ! (1.-swab/(swdn+.0001))
1231 :
1232 : ! for history
1233 : real (kind=dbl_kind) :: &
1234 1275244 : avdrl , & ! visible, direct, albedo (fraction)
1235 1275244 : avdfl , & ! visible, diffuse, albedo (fraction)
1236 1275244 : aidrl , & ! near-ir, direct, albedo (fraction)
1237 1275244 : aidfl ! near-ir, diffuse, albedo (fraction)
1238 :
1239 : character(len=*),parameter :: subname='(shortwave_dEdd)'
1240 :
1241 : !-----------------------------------------------------------------------
1242 :
1243 3770075 : klev = nslyr + nilyr + 1 ! number of radiation layers - 1
1244 3770075 : klevp = klev + 1 ! number of radiation interfaces - 1
1245 : ! (0 layer is included also)
1246 :
1247 : ! zero storage albedos and fluxes for accumulation over surface types:
1248 3770075 : hstmp = c0
1249 3770075 : hi = c0
1250 3770075 : fi = c0
1251 3770075 : alvdr = c0
1252 3770075 : alvdf = c0
1253 3770075 : alidr = c0
1254 3770075 : alidf = c0
1255 3770075 : avdrl = c0
1256 3770075 : avdfl = c0
1257 3770075 : aidrl = c0
1258 3770075 : aidfl = c0
1259 3770075 : fswsfc = c0
1260 3770075 : fswint = c0
1261 3770075 : fswthru = c0
1262 : ! compute fraction of nir down direct to total over all points:
1263 3770075 : fnidr = c0
1264 3770075 : if( swidr + swidf > puny ) then
1265 2909240 : fnidr = swidr/(swidr+swidf)
1266 : endif
1267 3770075 : albice = c0
1268 3770075 : albsno = c0
1269 3770075 : albpnd = c0
1270 33496281 : fswpenl(:) = c0
1271 8449182 : Sswabs(:) = c0
1272 29726206 : Iswabs(:) = c0
1273 :
1274 : ! compute aerosol mass path
1275 :
1276 18018659 : aero_mp(:) = c0
1277 3770075 : if( tr_aero ) then
1278 : ! check 4 layers for each aerosol, a snow SSL, snow below SSL,
1279 : ! sea ice SSL, and sea ice below SSL, in that order.
1280 229795 : if (size(aero) < 4*n_aero) then
1281 0 : call icepack_warnings_add(subname//' ERROR: size(aero) too small')
1282 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1283 0 : return
1284 : endif
1285 459590 : do na = 1, 4*n_aero, 4
1286 229795 : vsno = hs * aice
1287 229795 : netsw = swvdr + swidr + swvdf + swidf
1288 459590 : if (netsw > puny) then ! sun above horizon
1289 181316 : aero_mp(na ) = aero(na )*vsno
1290 181316 : aero_mp(na+1) = aero(na+1)*vsno
1291 181316 : aero_mp(na+2) = aero(na+2)*vice
1292 181316 : aero_mp(na+3) = aero(na+3)*vice
1293 : endif ! aice > 0 and netsw > 0
1294 : enddo ! na
1295 : endif ! if aerosols
1296 :
1297 : ! compute shortwave radiation accounting for snow/ice (both snow over
1298 : ! ice and bare ice) and ponded ice (if any):
1299 :
1300 : ! sea ice points with sun above horizon
1301 3770075 : netsw = swvdr + swidr + swvdf + swidf
1302 3770075 : if (netsw > puny) then ! sun above horizon
1303 2909240 : coszen = max(puny,coszen)
1304 : ! evaluate sea ice thickness and fraction
1305 2909240 : hi = vice / aice
1306 2909240 : fi = c1 - fs - fp
1307 : ! bare sea ice points
1308 2909240 : if(fi > c0) then
1309 : ! calculate bare sea ice
1310 :
1311 829017 : srftyp = 0
1312 : call compute_dEdd(nilyr, nslyr, klev, klevp, &
1313 0 : zbio, dEdd_algae, &
1314 : heat_capacity, fnidr, coszen, &
1315 : R_ice, R_pnd, &
1316 0 : kaer_tab, waer_tab, gaer_tab, &
1317 0 : kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
1318 0 : bcenh, modal_aero, kalg, &
1319 : swvdr, swvdf, swidr, swidf, srftyp, &
1320 0 : hstmp, rhosnw, rsnw, hi, hp, &
1321 0 : fi, aero_mp, avdrl, avdfl, &
1322 : aidrl, aidfl, &
1323 : fswsfc, fswint, &
1324 0 : fswthru, Sswabs, &
1325 829017 : Iswabs, fswpenl)
1326 829017 : if (icepack_warnings_aborted(subname)) return
1327 :
1328 829017 : alvdr = alvdr + avdrl *fi
1329 829017 : alvdf = alvdf + avdfl *fi
1330 829017 : alidr = alidr + aidrl *fi
1331 829017 : alidf = alidf + aidfl *fi
1332 : ! for history
1333 : albice = albice &
1334 : + awtvdr*avdrl + awtidr*aidrl &
1335 829017 : + awtvdf*avdfl + awtidf*aidfl
1336 : endif
1337 : endif
1338 :
1339 : ! sea ice points with sun above horizon
1340 3770075 : netsw = swvdr + swidr + swvdf + swidf
1341 3770075 : if (netsw > puny) then ! sun above horizon
1342 2909240 : coszen = max(puny,coszen)
1343 : ! snow-covered sea ice points
1344 2909240 : if(fs > c0) then
1345 : ! calculate snow covered sea ice
1346 :
1347 2092288 : srftyp = 1
1348 : call compute_dEdd(nilyr, nslyr, klev, klevp, &
1349 0 : zbio, dEdd_algae, &
1350 : heat_capacity, fnidr, coszen, &
1351 : R_ice, R_pnd, &
1352 0 : kaer_tab, waer_tab, gaer_tab, &
1353 0 : kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
1354 0 : bcenh, modal_aero, kalg, &
1355 : swvdr, swvdf, swidr, swidf, srftyp, &
1356 0 : hs, rhosnw, rsnw, hi, hp, &
1357 0 : fs, aero_mp, avdrl, avdfl, &
1358 : aidrl, aidfl, &
1359 : fswsfc, fswint, &
1360 0 : fswthru, Sswabs, &
1361 2092288 : Iswabs, fswpenl)
1362 2092288 : if (icepack_warnings_aborted(subname)) return
1363 :
1364 2092288 : alvdr = alvdr + avdrl *fs
1365 2092288 : alvdf = alvdf + avdfl *fs
1366 2092288 : alidr = alidr + aidrl *fs
1367 2092288 : alidf = alidf + aidfl *fs
1368 : ! for history
1369 : albsno = albsno &
1370 : + awtvdr*avdrl + awtidr*aidrl &
1371 2092288 : + awtvdf*avdfl + awtidf*aidfl
1372 : endif
1373 : endif
1374 :
1375 3770075 : hi = c0
1376 :
1377 : ! sea ice points with sun above horizon
1378 3770075 : netsw = swvdr + swidr + swvdf + swidf
1379 3770075 : if (netsw > puny) then ! sun above horizon
1380 2909240 : coszen = max(puny,coszen)
1381 2909240 : hi = vice / aice
1382 : ! if nonzero pond fraction and sufficient pond depth
1383 : ! if( fp > puny .and. hp > hpmin ) then
1384 2909240 : if (fp > puny) then
1385 :
1386 : ! calculate ponded ice
1387 :
1388 297565 : srftyp = 2
1389 : call compute_dEdd(nilyr, nslyr, klev, klevp, &
1390 0 : zbio, dEdd_algae, &
1391 : heat_capacity, fnidr, coszen, &
1392 : R_ice, R_pnd, &
1393 0 : kaer_tab, waer_tab, gaer_tab, &
1394 0 : kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
1395 0 : bcenh, modal_aero, kalg, &
1396 : swvdr, swvdf, swidr, swidf, srftyp, &
1397 0 : hs, rhosnw, rsnw, hi, hp, &
1398 0 : fp, aero_mp, avdrl, avdfl, &
1399 : aidrl, aidfl, &
1400 : fswsfc, fswint, &
1401 0 : fswthru, Sswabs, &
1402 297565 : Iswabs, fswpenl)
1403 297565 : if (icepack_warnings_aborted(subname)) return
1404 :
1405 297565 : alvdr = alvdr + avdrl *fp
1406 297565 : alvdf = alvdf + avdfl *fp
1407 297565 : alidr = alidr + aidrl *fp
1408 297565 : alidf = alidf + aidfl *fp
1409 : ! for history
1410 : albpnd = albpnd &
1411 : + awtvdr*avdrl + awtidr*aidrl &
1412 297565 : + awtvdf*avdfl + awtidf*aidfl
1413 : endif
1414 : endif
1415 :
1416 : ! if no incoming shortwave, set albedos to 1
1417 3770075 : netsw = swvdr + swidr + swvdf + swidf
1418 3770075 : if (netsw <= puny) then ! sun above horizon
1419 860835 : alvdr = c1
1420 860835 : alvdf = c1
1421 860835 : alidr = c1
1422 860835 : alidf = c1
1423 : endif
1424 :
1425 3770075 : if (l_print_point .and. netsw > puny) then
1426 :
1427 0 : write(warnstr,*) subname, ' printing point'
1428 0 : call icepack_warnings_add(warnstr)
1429 0 : write(warnstr,*) subname, ' coszen = ', &
1430 0 : coszen
1431 0 : call icepack_warnings_add(warnstr)
1432 0 : write(warnstr,*) subname, ' swvdr swvdf = ', &
1433 0 : swvdr,swvdf
1434 0 : call icepack_warnings_add(warnstr)
1435 0 : write(warnstr,*) subname, ' swidr swidf = ', &
1436 0 : swidr,swidf
1437 0 : call icepack_warnings_add(warnstr)
1438 0 : write(warnstr,*) subname, ' aice = ', &
1439 0 : aice
1440 0 : call icepack_warnings_add(warnstr)
1441 0 : write(warnstr,*) subname, ' hs = ', &
1442 0 : hs
1443 0 : call icepack_warnings_add(warnstr)
1444 0 : write(warnstr,*) subname, ' hp = ', &
1445 0 : hp
1446 0 : call icepack_warnings_add(warnstr)
1447 0 : write(warnstr,*) subname, ' fs = ', &
1448 0 : fs
1449 0 : call icepack_warnings_add(warnstr)
1450 0 : write(warnstr,*) subname, ' fi = ', &
1451 0 : fi
1452 0 : call icepack_warnings_add(warnstr)
1453 0 : write(warnstr,*) subname, ' fp = ', &
1454 0 : fp
1455 0 : call icepack_warnings_add(warnstr)
1456 0 : write(warnstr,*) subname, ' hi = ', &
1457 0 : hi
1458 0 : call icepack_warnings_add(warnstr)
1459 0 : write(warnstr,*) subname, ' alvdr alvdf = ', &
1460 0 : alvdr,alvdf
1461 0 : call icepack_warnings_add(warnstr)
1462 0 : write(warnstr,*) subname, ' alidr alidf = ', &
1463 0 : alidr,alidf
1464 0 : call icepack_warnings_add(warnstr)
1465 0 : write(warnstr,*) subname, ' fswsfc fswint fswthru = ', &
1466 0 : fswsfc,fswint,fswthru
1467 0 : call icepack_warnings_add(warnstr)
1468 0 : swdn = swvdr+swvdf+swidr+swidf
1469 0 : swab = fswsfc+fswint+fswthru
1470 0 : swalb = (1.-swab/(swdn+.0001))
1471 0 : write(warnstr,*) subname, ' swdn swab swalb = ',swdn,swab,swalb
1472 0 : do k = 1, nslyr
1473 0 : write(warnstr,*) subname, ' snow layer k = ', k, &
1474 0 : ' rhosnw = ', &
1475 0 : rhosnw(k), &
1476 0 : ' rsnw = ', &
1477 0 : rsnw(k)
1478 0 : call icepack_warnings_add(warnstr)
1479 : enddo
1480 0 : do k = 1, nslyr
1481 0 : write(warnstr,*) subname, ' snow layer k = ', k, &
1482 0 : ' Sswabs(k) = ', Sswabs(k)
1483 0 : call icepack_warnings_add(warnstr)
1484 : enddo
1485 0 : do k = 1, nilyr
1486 0 : write(warnstr,*) subname, ' sea ice layer k = ', k, &
1487 0 : ' Iswabs(k) = ', Iswabs(k)
1488 0 : call icepack_warnings_add(warnstr)
1489 : enddo
1490 :
1491 : endif ! l_print_point .and. coszen > .01
1492 :
1493 : end subroutine shortwave_dEdd
1494 :
1495 : !=======================================================================
1496 : !
1497 : ! Evaluate snow/ice/ponded ice inherent optical properties (IOPs), and
1498 : ! then calculate the multiple scattering solution by calling solution_dEdd.
1499 : !
1500 : ! author: Bruce P. Briegleb, NCAR
1501 : ! 2013: E Hunke merged with NCAR version
1502 :
1503 3218870 : subroutine compute_dEdd (nilyr, nslyr, klev, klevp, &
1504 3218870 : zbio, dEdd_algae, &
1505 : heat_capacity, fnidr, coszen, &
1506 : R_ice, R_pnd, &
1507 3218870 : kaer_tab, waer_tab, gaer_tab, &
1508 3218870 : kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
1509 3218870 : bcenh, modal_aero, kalg, &
1510 : swvdr, swvdf, swidr, swidf, srftyp, &
1511 6437740 : hs, rhosnw, rsnw, hi, hp, &
1512 3218870 : fi, aero_mp, alvdr, alvdf, &
1513 : alidr, alidf, &
1514 : fswsfc, fswint, &
1515 3218870 : fswthru, Sswabs, &
1516 3218870 : Iswabs, fswpenl)
1517 :
1518 : integer (kind=int_kind), intent(in) :: &
1519 : nilyr , & ! number of ice layers
1520 : nslyr , & ! number of snow layers
1521 : klev , & ! number of radiation layers - 1
1522 : klevp ! number of radiation interfaces - 1
1523 : ! (0 layer is included also)
1524 :
1525 : logical (kind=log_kind), intent(in) :: &
1526 : heat_capacity,& ! if true, ice has nonzero heat capacity
1527 : dEdd_algae, & ! .true. use prognostic chla in dEdd
1528 : modal_aero ! .true. use modal aerosol treatment
1529 :
1530 : real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
1531 : kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
1532 : waer_bc_tab, & ! aerosol single scatter albedo (fraction)
1533 : gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
1534 :
1535 : real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
1536 : bcenh ! BC absorption enhancement factor
1537 :
1538 : ! dEdd tuning parameters, set in namelist
1539 : real (kind=dbl_kind), intent(in) :: &
1540 : R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
1541 : R_pnd ! ponded ice tuning parameter; +1 > 1sig increase in albedo
1542 :
1543 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
1544 : kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
1545 : waer_tab, & ! aerosol single scatter albedo (fraction)
1546 : gaer_tab ! aerosol asymmetry parameter (cos(theta))
1547 :
1548 : real (kind=dbl_kind), intent(in) :: &
1549 : kalg , & ! algae absorption coefficient
1550 : fnidr , & ! fraction of direct to total down flux in nir
1551 : coszen , & ! cosine solar zenith angle
1552 : swvdr , & ! shortwave down at surface, visible, direct (W/m^2)
1553 : swvdf , & ! shortwave down at surface, visible, diffuse (W/m^2)
1554 : swidr , & ! shortwave down at surface, near IR, direct (W/m^2)
1555 : swidf ! shortwave down at surface, near IR, diffuse (W/m^2)
1556 :
1557 : integer (kind=int_kind), intent(in) :: &
1558 : srftyp ! surface type over ice: (0=air, 1=snow, 2=pond)
1559 :
1560 : real (kind=dbl_kind), intent(in) :: &
1561 : hs ! snow thickness (m)
1562 :
1563 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1564 : rhosnw , & ! snow density in snow layer (kg/m3)
1565 : rsnw , & ! snow grain radius in snow layer (m)
1566 : zbio , & ! zaerosol + chla shortwave tracers kg/m^3
1567 : aero_mp ! aerosol mass path in kg/m2
1568 :
1569 : real (kind=dbl_kind), intent(in) :: &
1570 : hi , & ! ice thickness (m)
1571 : hp , & ! pond depth (m)
1572 : fi ! snow/bare ice fractional coverage (0 to 1)
1573 :
1574 : real (kind=dbl_kind), intent(inout) :: &
1575 : alvdr , & ! visible, direct, albedo (fraction)
1576 : alvdf , & ! visible, diffuse, albedo (fraction)
1577 : alidr , & ! near-ir, direct, albedo (fraction)
1578 : alidf , & ! near-ir, diffuse, albedo (fraction)
1579 : fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
1580 : fswint , & ! SW interior absorption (below surface, above ocean,W m-2)
1581 : fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2)
1582 :
1583 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1584 : fswpenl , & ! visible SW entering ice layers (W m-2)
1585 : Sswabs , & ! SW absorbed in snow layer (W m-2)
1586 : Iswabs ! SW absorbed in ice layer (W m-2)
1587 :
1588 : !-----------------------------------------------------------------------
1589 : !
1590 : ! Set up optical property profiles, based on snow, sea ice and ponded
1591 : ! ice IOPs from:
1592 : !
1593 : ! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple
1594 : ! Scattering Parameterization for Solar Radiation in the Sea Ice
1595 : ! Component of the Community Climate System Model, NCAR Technical
1596 : ! Note NCAR/TN-472+STR February 2007
1597 : !
1598 : ! Computes column Delta-Eddington radiation solution for specific
1599 : ! surface type: either snow over sea ice, bare sea ice, or ponded sea ice.
1600 : !
1601 : ! Divides solar spectrum into 3 intervals: 0.2-0.7, 0.7-1.19, and
1602 : ! 1.19-5.0 micro-meters. The latter two are added (using an assumed
1603 : ! partition of incident shortwave in the 0.7-5.0 micro-meter band between
1604 : ! the 0.7-1.19 and 1.19-5.0 micro-meter band) to give the final output
1605 : ! of 0.2-0.7 visible and 0.7-5.0 near-infrared albedos and fluxes.
1606 : !
1607 : ! Specifies vertical layer optical properties based on input snow depth,
1608 : ! density and grain radius, along with ice and pond depths, then computes
1609 : ! layer by layer Delta-Eddington reflectivity, transmissivity and combines
1610 : ! layers (done by calling routine solution_dEdd). Finally, surface albedos
1611 : ! and internal fluxes/flux divergences are evaluated.
1612 : !
1613 : ! Description of the level and layer index conventions. This is
1614 : ! for the standard case of one snow layer and four sea ice layers.
1615 : !
1616 : ! Please read the following; otherwise, there is 99.9% chance you
1617 : ! will be confused about indices at some point in time........ :)
1618 : !
1619 : ! CICE4.0 snow treatment has one snow layer above the sea ice. This
1620 : ! snow layer has finite heat capacity, so that surface absorption must
1621 : ! be distinguished from internal. The Delta-Eddington solar radiation
1622 : ! thus adds extra surface scattering layers to both snow and sea ice.
1623 : ! Note that in the following, we assume a fixed vertical layer structure
1624 : ! for the radiation calculation. In other words, we always have the
1625 : ! structure shown below for one snow and four sea ice layers, but for
1626 : ! ponded ice the pond fills "snow" layer 1 over the sea ice, and for
1627 : ! bare sea ice the top layers over sea ice are treated as transparent air.
1628 : !
1629 : ! SSL = surface scattering layer for either snow or sea ice
1630 : ! DL = drained layer for sea ice immediately under sea ice SSL
1631 : ! INT = interior layers for sea ice below the drained layer.
1632 : !
1633 : ! Notice that the radiation level starts with 0 at the top. Thus,
1634 : ! the total number radiation layers is klev+1, where klev is the
1635 : ! sum of nslyr, the number of CCSM snow layers, and nilyr, the
1636 : ! number of CCSM sea ice layers, plus the sea ice SSL:
1637 : ! klev = 1 + nslyr + nilyr
1638 : !
1639 : ! For the standard case illustrated below, nslyr=1, nilyr=4,
1640 : ! and klev=6, with the number of layer interfaces klevp=klev+1.
1641 : ! Layer interfaces are the surfaces on which reflectivities,
1642 : ! transmissivities and fluxes are evaluated.
1643 : !
1644 : ! CCSM3 Sea Ice Model Delta-Eddington Solar Radiation
1645 : ! Layers and Interfaces
1646 : ! Layer Index Interface Index
1647 : ! --------------------- --------------------- 0
1648 : ! 0 \\\ snow SSL \\\
1649 : ! snow layer 1 --------------------- 1
1650 : ! 1 rest of snow layer
1651 : ! +++++++++++++++++++++ +++++++++++++++++++++ 2
1652 : ! 2 \\\ sea ice SSL \\\
1653 : ! sea ice layer 1 --------------------- 3
1654 : ! 3 sea ice DL
1655 : ! --------------------- --------------------- 4
1656 : !
1657 : ! sea ice layer 2 4 sea ice INT
1658 : !
1659 : ! --------------------- --------------------- 5
1660 : !
1661 : ! sea ice layer 3 5 sea ice INT
1662 : !
1663 : ! --------------------- --------------------- 6
1664 : !
1665 : ! sea ice layer 4 6 sea ice INT
1666 : !
1667 : ! --------------------- --------------------- 7
1668 : !
1669 : ! When snow lies over sea ice, the radiation absorbed in the
1670 : ! snow SSL is used for surface heating, and that in the rest
1671 : ! of the snow layer for its internal heating. For sea ice in
1672 : ! this case, all of the radiant heat absorbed in both the
1673 : ! sea ice SSL and the DL are used for sea ice layer 1 heating.
1674 : !
1675 : ! When pond lies over sea ice, and for bare sea ice, all of the
1676 : ! radiant heat absorbed within and above the sea ice SSL is used
1677 : ! for surface heating, and that absorbed in the sea ice DL is
1678 : ! used for sea ice layer 1 heating.
1679 : !
1680 : ! Basically, vertical profiles of the layer extinction optical depth (tau),
1681 : ! single scattering albedo (w0) and asymmetry parameter (g) are required over
1682 : ! the klev+1 layers, where klev+1 = 2 + nslyr + nilyr. All of the surface type
1683 : ! information and snow/ice iop properties are evaulated in this routine, so
1684 : ! the tau,w0,g profiles can be passed to solution_dEdd for multiple scattering
1685 : ! evaluation. Snow, bare ice and ponded ice iops are contained in data arrays
1686 : ! in this routine.
1687 : !
1688 : !-----------------------------------------------------------------------
1689 :
1690 : ! local variables
1691 :
1692 : integer (kind=int_kind) :: &
1693 : k , & ! level index
1694 : ns , & ! spectral index
1695 : nr , & ! index for grain radius tables
1696 : ki , & ! index for internal absorption
1697 : km , & ! k starting index for snow, sea ice internal absorption
1698 : kp , & ! k+1 or k+2 index for snow, sea ice internal absorption
1699 : ksrf , & ! level index for surface absorption
1700 : ksnow , & ! level index for snow density and grain size
1701 : kii ! level starting index for sea ice (nslyr+1)
1702 :
1703 : integer (kind=int_kind), parameter :: &
1704 : nmbrad = 32 ! number of snow grain radii in tables
1705 :
1706 : real (kind=dbl_kind) :: &
1707 1142733 : avdr , & ! visible albedo, direct (fraction)
1708 1142733 : avdf , & ! visible albedo, diffuse (fraction)
1709 1142733 : aidr , & ! near-ir albedo, direct (fraction)
1710 1142733 : aidf ! near-ir albedo, diffuse (fraction)
1711 :
1712 : real (kind=dbl_kind) :: &
1713 1142733 : fsfc , & ! shortwave absorbed at snow/bare ice/ponded ice surface (W m-2)
1714 1142733 : fint , & ! shortwave absorbed in interior (W m-2)
1715 1142733 : fthru ! shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
1716 :
1717 : real (kind=dbl_kind), dimension(nslyr) :: &
1718 6870812 : Sabs ! shortwave absorbed in snow layer (W m-2)
1719 :
1720 : real (kind=dbl_kind), dimension(nilyr) :: &
1721 13294138 : Iabs ! shortwave absorbed in ice layer (W m-2)
1722 :
1723 : real (kind=dbl_kind), dimension(nilyr+1) :: &
1724 14436871 : fthrul ! shortwave through to ice layers (W m-2)
1725 :
1726 : real (kind=dbl_kind), dimension (nspint) :: &
1727 4570932 : wghtns ! spectral weights
1728 :
1729 : real (kind=dbl_kind), parameter :: &
1730 : cp67 = 0.67_dbl_kind , & ! nir band weight parameter
1731 : cp78 = 0.78_dbl_kind , & ! nir band weight parameter
1732 : cp01 = 0.01_dbl_kind ! for ocean visible albedo
1733 :
1734 : real (kind=dbl_kind), dimension (0:klev) :: &
1735 17155409 : tau , & ! layer extinction optical depth
1736 17155409 : w0 , & ! layer single scattering albedo
1737 17155409 : g ! layer asymmetry parameter
1738 :
1739 : ! following arrays are defined at model interfaces; 0 is the top of the
1740 : ! layer above the sea ice; klevp is the sea ice/ocean interface.
1741 : real (kind=dbl_kind), dimension (0:klevp) :: &
1742 18298142 : trndir , & ! solar beam down transmission from top
1743 18298142 : trntdr , & ! total transmission to direct beam for layers above
1744 18298142 : trndif , & ! diffuse transmission to diffuse beam for layers above
1745 18298142 : rupdir , & ! reflectivity to direct radiation for layers below
1746 18298142 : rupdif , & ! reflectivity to diffuse radiation for layers below
1747 18298142 : rdndif ! reflectivity to diffuse radiation for layers above
1748 :
1749 : real (kind=dbl_kind), dimension (0:klevp) :: &
1750 18298142 : dfdir , & ! down-up flux at interface due to direct beam at top surface
1751 20583608 : dfdif ! down-up flux at interface due to diffuse beam at top surface
1752 :
1753 : real (kind=dbl_kind) :: &
1754 1142733 : refk , & ! interface k multiple scattering term
1755 1142733 : delr , & ! snow grain radius interpolation parameter
1756 : ! inherent optical properties (iop) for snow
1757 1142733 : Qs , & ! Snow extinction efficiency
1758 1142733 : ks , & ! Snow extinction coefficient (/m)
1759 1142733 : ws , & ! Snow single scattering albedo
1760 1142733 : gs ! Snow asymmetry parameter
1761 :
1762 : real (kind=dbl_kind), dimension(nslyr) :: &
1763 6870812 : frsnw ! snow grain radius in snow layer * adjustment factor (m)
1764 :
1765 : ! actual used ice and ponded ice IOPs, allowing for tuning
1766 : ! modifications of the above "_mn" value
1767 : real (kind=dbl_kind), dimension (nspint) :: &
1768 4570932 : ki_ssl , & ! Surface-scattering-layer ice extinction coefficient (/m)
1769 4570932 : wi_ssl , & ! Surface-scattering-layer ice single scattering albedo
1770 4570932 : gi_ssl , & ! Surface-scattering-layer ice asymmetry parameter
1771 4570932 : ki_dl , & ! Drained-layer ice extinction coefficient (/m)
1772 4570932 : wi_dl , & ! Drained-layer ice single scattering albedo
1773 4570932 : gi_dl , & ! Drained-layer ice asymmetry parameter
1774 4570932 : ki_int , & ! Interior-layer ice extinction coefficient (/m)
1775 4570932 : wi_int , & ! Interior-layer ice single scattering albedo
1776 4570932 : gi_int , & ! Interior-layer ice asymmetry parameter
1777 4570932 : ki_p_ssl , & ! Ice under pond srf scat layer extinction coefficient (/m)
1778 4570932 : wi_p_ssl , & ! Ice under pond srf scat layer single scattering albedo
1779 4570932 : gi_p_ssl , & ! Ice under pond srf scat layer asymmetry parameter
1780 4570932 : ki_p_int , & ! Ice under pond extinction coefficient (/m)
1781 4570932 : wi_p_int , & ! Ice under pond single scattering albedo
1782 4570932 : gi_p_int ! Ice under pond asymmetry parameter
1783 :
1784 : real (kind=dbl_kind), dimension(0:klev) :: &
1785 17155409 : dzk ! layer thickness
1786 :
1787 : real (kind=dbl_kind) :: &
1788 1142733 : dz , & ! snow, sea ice or pond water layer thickness
1789 1142733 : dz_ssl , & ! snow or sea ice surface scattering layer thickness
1790 1142733 : fs ! scaling factor to reduce (nilyr<4) or increase (nilyr>4) DL
1791 : ! extinction coefficient to maintain DL optical depth constant
1792 : ! with changing number of sea ice layers, to approximately
1793 : ! conserve computed albedo for constant physical depth of sea
1794 : ! ice when the number of sea ice layers vary
1795 : real (kind=dbl_kind) :: &
1796 1142733 : sig , & ! scattering coefficient for tuning
1797 1142733 : kabs , & ! absorption coefficient for tuning
1798 1142733 : sigp ! modified scattering coefficient for tuning
1799 :
1800 : real (kind=dbl_kind), dimension(nspint, 0:klev) :: &
1801 52736615 : kabs_chl , & ! absorption coefficient for chlorophyll (/m)
1802 52736615 : tzaer , & ! total aerosol extinction optical depth
1803 52736615 : wzaer , & ! total aerosol single scatter albedo
1804 52736615 : gzaer ! total aerosol asymmetry parameter
1805 :
1806 : real (kind=dbl_kind) :: &
1807 1142733 : albodr , & ! spectral ocean albedo to direct rad
1808 1142733 : albodf ! spectral ocean albedo to diffuse rad
1809 :
1810 : ! for melt pond transition to bare sea ice for small pond depths
1811 : real (kind=dbl_kind) :: &
1812 1142733 : sig_i , & ! ice scattering coefficient (/m)
1813 1142733 : sig_p , & ! pond scattering coefficient (/m)
1814 1142733 : kext ! weighted extinction coefficient (/m)
1815 :
1816 : ! aerosol optical properties from Mark Flanner, 26 June 2008
1817 : ! order assumed: hydrophobic black carbon, hydrophilic black carbon,
1818 : ! four dust aerosols by particle size range:
1819 : ! dust1(.05-0.5 micron), dust2(0.5-1.25 micron),
1820 : ! dust3(1.25-2.5 micron), dust4(2.5-5.0 micron)
1821 : ! spectral bands same as snow/sea ice: (0.3-0.7 micron, 0.7-1.19 micron
1822 : ! and 1.19-5.0 micron in wavelength)
1823 :
1824 : integer (kind=int_kind) :: &
1825 : na , n ! aerosol index
1826 :
1827 : real (kind=dbl_kind) :: &
1828 1142733 : taer , & ! total aerosol extinction optical depth
1829 1142733 : waer , & ! total aerosol single scatter albedo
1830 1142733 : gaer , & ! total aerosol asymmetry parameter
1831 1142733 : swdr , & ! shortwave down at surface, direct (W/m^2)
1832 1142733 : swdf , & ! shortwave down at surface, diffuse (W/m^2)
1833 1142733 : rnilyr , & ! real(nilyr)
1834 1142733 : rnslyr , & ! real(nslyr)
1835 1142733 : rns , & ! real(ns)
1836 2285466 : tmp_0, tmp_ks, tmp_kl ! temp variables
1837 :
1838 : integer(kind=int_kind), dimension(0:klev) :: &
1839 6437740 : k_bcini , & !
1840 6437740 : k_bcins , &
1841 4361603 : k_bcexs
1842 : real(kind=dbl_kind):: &
1843 1142733 : tmp_gs, tmp1 ! temp variables
1844 :
1845 : ! snow grain radii (micro-meters) for table
1846 : real (kind=dbl_kind), dimension(nmbrad), parameter :: &
1847 : rsnw_tab = (/ & ! snow grain radius for each table entry (micro-meters)
1848 : 5._dbl_kind, 7._dbl_kind, 10._dbl_kind, 15._dbl_kind, &
1849 : 20._dbl_kind, 30._dbl_kind, 40._dbl_kind, 50._dbl_kind, &
1850 : 65._dbl_kind, 80._dbl_kind, 100._dbl_kind, 120._dbl_kind, &
1851 : 140._dbl_kind, 170._dbl_kind, 200._dbl_kind, 240._dbl_kind, &
1852 : 290._dbl_kind, 350._dbl_kind, 420._dbl_kind, 500._dbl_kind, &
1853 : 570._dbl_kind, 660._dbl_kind, 760._dbl_kind, 870._dbl_kind, &
1854 : 1000._dbl_kind, 1100._dbl_kind, 1250._dbl_kind, 1400._dbl_kind, &
1855 : 1600._dbl_kind, 1800._dbl_kind, 2000._dbl_kind, 2500._dbl_kind/)
1856 :
1857 : ! snow extinction efficiency (unitless)
1858 : real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
1859 : Qs_tab = reshape((/ &
1860 : 2.131798_dbl_kind, 2.187756_dbl_kind, 2.267358_dbl_kind, &
1861 : 2.104499_dbl_kind, 2.148345_dbl_kind, 2.236078_dbl_kind, &
1862 : 2.081580_dbl_kind, 2.116885_dbl_kind, 2.175067_dbl_kind, &
1863 : 2.062595_dbl_kind, 2.088937_dbl_kind, 2.130242_dbl_kind, &
1864 : 2.051403_dbl_kind, 2.072422_dbl_kind, 2.106610_dbl_kind, &
1865 : 2.039223_dbl_kind, 2.055389_dbl_kind, 2.080586_dbl_kind, &
1866 : 2.032383_dbl_kind, 2.045751_dbl_kind, 2.066394_dbl_kind, &
1867 : 2.027920_dbl_kind, 2.039388_dbl_kind, 2.057224_dbl_kind, &
1868 : 2.023444_dbl_kind, 2.033137_dbl_kind, 2.048055_dbl_kind, &
1869 : 2.020412_dbl_kind, 2.028840_dbl_kind, 2.041874_dbl_kind, &
1870 : 2.017608_dbl_kind, 2.024863_dbl_kind, 2.036046_dbl_kind, &
1871 : 2.015592_dbl_kind, 2.022021_dbl_kind, 2.031954_dbl_kind, &
1872 : 2.014083_dbl_kind, 2.019887_dbl_kind, 2.028853_dbl_kind, &
1873 : 2.012368_dbl_kind, 2.017471_dbl_kind, 2.025353_dbl_kind, &
1874 : 2.011092_dbl_kind, 2.015675_dbl_kind, 2.022759_dbl_kind, &
1875 : 2.009837_dbl_kind, 2.013897_dbl_kind, 2.020168_dbl_kind, &
1876 : 2.008668_dbl_kind, 2.012252_dbl_kind, 2.017781_dbl_kind, &
1877 : 2.007627_dbl_kind, 2.010813_dbl_kind, 2.015678_dbl_kind, &
1878 : 2.006764_dbl_kind, 2.009577_dbl_kind, 2.013880_dbl_kind, &
1879 : 2.006037_dbl_kind, 2.008520_dbl_kind, 2.012382_dbl_kind, &
1880 : 2.005528_dbl_kind, 2.007807_dbl_kind, 2.011307_dbl_kind, &
1881 : 2.005025_dbl_kind, 2.007079_dbl_kind, 2.010280_dbl_kind, &
1882 : 2.004562_dbl_kind, 2.006440_dbl_kind, 2.009333_dbl_kind, &
1883 : 2.004155_dbl_kind, 2.005898_dbl_kind, 2.008523_dbl_kind, &
1884 : 2.003794_dbl_kind, 2.005379_dbl_kind, 2.007795_dbl_kind, &
1885 : 2.003555_dbl_kind, 2.005041_dbl_kind, 2.007329_dbl_kind, &
1886 : 2.003264_dbl_kind, 2.004624_dbl_kind, 2.006729_dbl_kind, &
1887 : 2.003037_dbl_kind, 2.004291_dbl_kind, 2.006230_dbl_kind, &
1888 : 2.002776_dbl_kind, 2.003929_dbl_kind, 2.005700_dbl_kind, &
1889 : 2.002590_dbl_kind, 2.003627_dbl_kind, 2.005276_dbl_kind, &
1890 : 2.002395_dbl_kind, 2.003391_dbl_kind, 2.004904_dbl_kind, &
1891 : 2.002071_dbl_kind, 2.002922_dbl_kind, 2.004241_dbl_kind/), &
1892 : (/nspint,nmbrad/))
1893 :
1894 : ! snow single scattering albedo (unitless)
1895 : real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
1896 : ws_tab = reshape((/ &
1897 : 0.9999994_dbl_kind, 0.9999673_dbl_kind, 0.9954589_dbl_kind, &
1898 : 0.9999992_dbl_kind, 0.9999547_dbl_kind, 0.9938576_dbl_kind, &
1899 : 0.9999990_dbl_kind, 0.9999382_dbl_kind, 0.9917989_dbl_kind, &
1900 : 0.9999985_dbl_kind, 0.9999123_dbl_kind, 0.9889724_dbl_kind, &
1901 : 0.9999979_dbl_kind, 0.9998844_dbl_kind, 0.9866190_dbl_kind, &
1902 : 0.9999970_dbl_kind, 0.9998317_dbl_kind, 0.9823021_dbl_kind, &
1903 : 0.9999960_dbl_kind, 0.9997800_dbl_kind, 0.9785269_dbl_kind, &
1904 : 0.9999951_dbl_kind, 0.9997288_dbl_kind, 0.9751601_dbl_kind, &
1905 : 0.9999936_dbl_kind, 0.9996531_dbl_kind, 0.9706974_dbl_kind, &
1906 : 0.9999922_dbl_kind, 0.9995783_dbl_kind, 0.9667577_dbl_kind, &
1907 : 0.9999903_dbl_kind, 0.9994798_dbl_kind, 0.9621007_dbl_kind, &
1908 : 0.9999885_dbl_kind, 0.9993825_dbl_kind, 0.9579541_dbl_kind, &
1909 : 0.9999866_dbl_kind, 0.9992862_dbl_kind, 0.9541924_dbl_kind, &
1910 : 0.9999838_dbl_kind, 0.9991434_dbl_kind, 0.9490959_dbl_kind, &
1911 : 0.9999810_dbl_kind, 0.9990025_dbl_kind, 0.9444940_dbl_kind, &
1912 : 0.9999772_dbl_kind, 0.9988171_dbl_kind, 0.9389141_dbl_kind, &
1913 : 0.9999726_dbl_kind, 0.9985890_dbl_kind, 0.9325819_dbl_kind, &
1914 : 0.9999670_dbl_kind, 0.9983199_dbl_kind, 0.9256405_dbl_kind, &
1915 : 0.9999605_dbl_kind, 0.9980117_dbl_kind, 0.9181533_dbl_kind, &
1916 : 0.9999530_dbl_kind, 0.9976663_dbl_kind, 0.9101540_dbl_kind, &
1917 : 0.9999465_dbl_kind, 0.9973693_dbl_kind, 0.9035031_dbl_kind, &
1918 : 0.9999382_dbl_kind, 0.9969939_dbl_kind, 0.8953134_dbl_kind, &
1919 : 0.9999289_dbl_kind, 0.9965848_dbl_kind, 0.8865789_dbl_kind, &
1920 : 0.9999188_dbl_kind, 0.9961434_dbl_kind, 0.8773350_dbl_kind, &
1921 : 0.9999068_dbl_kind, 0.9956323_dbl_kind, 0.8668233_dbl_kind, &
1922 : 0.9998975_dbl_kind, 0.9952464_dbl_kind, 0.8589990_dbl_kind, &
1923 : 0.9998837_dbl_kind, 0.9946782_dbl_kind, 0.8476493_dbl_kind, &
1924 : 0.9998699_dbl_kind, 0.9941218_dbl_kind, 0.8367318_dbl_kind, &
1925 : 0.9998515_dbl_kind, 0.9933966_dbl_kind, 0.8227881_dbl_kind, &
1926 : 0.9998332_dbl_kind, 0.9926888_dbl_kind, 0.8095131_dbl_kind, &
1927 : 0.9998148_dbl_kind, 0.9919968_dbl_kind, 0.7968620_dbl_kind, &
1928 : 0.9997691_dbl_kind, 0.9903277_dbl_kind, 0.7677887_dbl_kind/), &
1929 : (/nspint,nmbrad/))
1930 :
1931 : ! snow asymmetry parameter (unitless)
1932 : real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
1933 : gs_tab = reshape((/ &
1934 : 0.859913_dbl_kind, 0.848003_dbl_kind, 0.824415_dbl_kind, &
1935 : 0.867130_dbl_kind, 0.858150_dbl_kind, 0.848445_dbl_kind, &
1936 : 0.873381_dbl_kind, 0.867221_dbl_kind, 0.861714_dbl_kind, &
1937 : 0.878368_dbl_kind, 0.874879_dbl_kind, 0.874036_dbl_kind, &
1938 : 0.881462_dbl_kind, 0.879661_dbl_kind, 0.881299_dbl_kind, &
1939 : 0.884361_dbl_kind, 0.883903_dbl_kind, 0.890184_dbl_kind, &
1940 : 0.885937_dbl_kind, 0.886256_dbl_kind, 0.895393_dbl_kind, &
1941 : 0.886931_dbl_kind, 0.887769_dbl_kind, 0.899072_dbl_kind, &
1942 : 0.887894_dbl_kind, 0.889255_dbl_kind, 0.903285_dbl_kind, &
1943 : 0.888515_dbl_kind, 0.890236_dbl_kind, 0.906588_dbl_kind, &
1944 : 0.889073_dbl_kind, 0.891127_dbl_kind, 0.910152_dbl_kind, &
1945 : 0.889452_dbl_kind, 0.891750_dbl_kind, 0.913100_dbl_kind, &
1946 : 0.889730_dbl_kind, 0.892213_dbl_kind, 0.915621_dbl_kind, &
1947 : 0.890026_dbl_kind, 0.892723_dbl_kind, 0.918831_dbl_kind, &
1948 : 0.890238_dbl_kind, 0.893099_dbl_kind, 0.921540_dbl_kind, &
1949 : 0.890441_dbl_kind, 0.893474_dbl_kind, 0.924581_dbl_kind, &
1950 : 0.890618_dbl_kind, 0.893816_dbl_kind, 0.927701_dbl_kind, &
1951 : 0.890762_dbl_kind, 0.894123_dbl_kind, 0.930737_dbl_kind, &
1952 : 0.890881_dbl_kind, 0.894397_dbl_kind, 0.933568_dbl_kind, &
1953 : 0.890975_dbl_kind, 0.894645_dbl_kind, 0.936148_dbl_kind, &
1954 : 0.891035_dbl_kind, 0.894822_dbl_kind, 0.937989_dbl_kind, &
1955 : 0.891097_dbl_kind, 0.895020_dbl_kind, 0.939949_dbl_kind, &
1956 : 0.891147_dbl_kind, 0.895212_dbl_kind, 0.941727_dbl_kind, &
1957 : 0.891189_dbl_kind, 0.895399_dbl_kind, 0.943339_dbl_kind, &
1958 : 0.891225_dbl_kind, 0.895601_dbl_kind, 0.944915_dbl_kind, &
1959 : 0.891248_dbl_kind, 0.895745_dbl_kind, 0.945950_dbl_kind, &
1960 : 0.891277_dbl_kind, 0.895951_dbl_kind, 0.947288_dbl_kind, &
1961 : 0.891299_dbl_kind, 0.896142_dbl_kind, 0.948438_dbl_kind, &
1962 : 0.891323_dbl_kind, 0.896388_dbl_kind, 0.949762_dbl_kind, &
1963 : 0.891340_dbl_kind, 0.896623_dbl_kind, 0.950916_dbl_kind, &
1964 : 0.891356_dbl_kind, 0.896851_dbl_kind, 0.951945_dbl_kind, &
1965 : 0.891386_dbl_kind, 0.897399_dbl_kind, 0.954156_dbl_kind/), &
1966 : (/nspint,nmbrad/))
1967 :
1968 : ! inherent optical property (iop) arrays for ice and ponded ice
1969 : ! mn = specified mean (or base) value
1970 : ! ki = extinction coefficient (/m)
1971 : ! wi = single scattering albedo
1972 : ! gi = asymmetry parameter
1973 :
1974 : ! ice surface scattering layer (ssl) iops
1975 : real (kind=dbl_kind), dimension (nspint), parameter :: &
1976 : ki_ssl_mn = (/ 1000.1_dbl_kind, 1003.7_dbl_kind, 7042._dbl_kind/), &
1977 : wi_ssl_mn = (/ .9999_dbl_kind, .9963_dbl_kind, .9088_dbl_kind/), &
1978 : gi_ssl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind/)
1979 :
1980 : ! ice drained layer (dl) iops
1981 : real (kind=dbl_kind), dimension (nspint), parameter :: &
1982 : ki_dl_mn = (/ 100.2_dbl_kind, 107.7_dbl_kind, 1309._dbl_kind /), &
1983 : wi_dl_mn = (/ .9980_dbl_kind, .9287_dbl_kind, .0305_dbl_kind /), &
1984 : gi_dl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
1985 :
1986 : ! ice interior layer (int) iops
1987 : real (kind=dbl_kind), dimension (nspint), parameter :: &
1988 : ki_int_mn = (/ 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind /), &
1989 : wi_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind /), &
1990 : gi_int_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
1991 :
1992 : ! ponded ice surface scattering layer (ssl) iops
1993 : real (kind=dbl_kind), dimension (nspint), parameter :: &
1994 : ki_p_ssl_mn = (/ 70.2_dbl_kind, 77.7_dbl_kind, 1309._dbl_kind/), &
1995 : wi_p_ssl_mn = (/ .9972_dbl_kind, .9009_dbl_kind, .0305_dbl_kind/), &
1996 : gi_p_ssl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
1997 :
1998 : ! ponded ice interior layer (int) iops
1999 : real (kind=dbl_kind), dimension (nspint), parameter :: &
2000 : ki_p_int_mn = (/ 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind/), &
2001 : wi_p_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind/), &
2002 : gi_p_int_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
2003 :
2004 : ! inherent optical property (iop) arrays for pond water and underlying ocean
2005 : ! kw = Pond water extinction coefficient (/m)
2006 : ! ww = Pond water single scattering albedo
2007 : ! gw = Pond water asymmetry parameter
2008 : real (kind=dbl_kind), dimension (nspint), parameter :: &
2009 : kw = (/ 0.20_dbl_kind, 12.0_dbl_kind, 729._dbl_kind /), &
2010 : ww = (/ 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /), &
2011 : gw = (/ 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /)
2012 :
2013 : real (kind=dbl_kind), parameter :: &
2014 : rhoi = 917.0_dbl_kind,& ! pure ice mass density (kg/m3)
2015 : fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max
2016 : fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min
2017 : ! tuning parameters
2018 : ! ice and pond scat coeff fractional change for +- one-sigma in albedo
2019 : fp_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for + stn dev in alb
2020 : fm_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for - stn dev in alb
2021 : fp_pnd = 2.00_dbl_kind, & ! ponded ice fraction of scat coeff for + stn dev in alb
2022 : fm_pnd = 0.50_dbl_kind ! ponded ice fraction of scat coeff for - stn dev in alb
2023 :
2024 : real (kind=dbl_kind), parameter :: & !chla-specific absorption coefficient
2025 : kchl_tab = 0.01 !0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1
2026 : ! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014
2027 : ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011)
2028 : ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014)
2029 : ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting)
2030 : !chlorophyll mass extinction cross section (m^2/mg chla)
2031 :
2032 : character(len=*),parameter :: subname='(compute_dEdd)'
2033 :
2034 : !-----------------------------------------------------------------------
2035 : ! Initialize and tune bare ice/ponded ice iops
2036 :
2037 36638086 : k_bcini(:) = c0
2038 36638086 : k_bcins(:) = c0
2039 36638086 : k_bcexs(:) = c0
2040 :
2041 3218870 : rnilyr = c1/real(nilyr,kind=dbl_kind)
2042 3218870 : rnslyr = c1/real(nslyr,kind=dbl_kind)
2043 3218870 : kii = nslyr + 1
2044 :
2045 : ! initialize albedos and fluxes to 0
2046 28969830 : fthrul = c0
2047 25750960 : Iabs = c0
2048 136895734 : kabs_chl(:,:) = c0
2049 136895734 : tzaer(:,:) = c0
2050 136895734 : wzaer(:,:) = c0
2051 136895734 : gzaer(:,:) = c0
2052 :
2053 3218870 : avdr = c0
2054 3218870 : avdf = c0
2055 3218870 : aidr = c0
2056 3218870 : aidf = c0
2057 3218870 : fsfc = c0
2058 3218870 : fint = c0
2059 3218870 : fthru = c0
2060 :
2061 : ! spectral weights
2062 : ! weights 2 (0.7-1.19 micro-meters) and 3 (1.19-5.0 micro-meters)
2063 : ! are chosen based on 1D calculations using ratio of direct to total
2064 : ! near-infrared solar (0.7-5.0 micro-meter) which indicates clear/cloudy
2065 : ! conditions: more cloud, the less 1.19-5.0 relative to the
2066 : ! 0.7-1.19 micro-meter due to cloud absorption.
2067 3218870 : wghtns(1) = c1
2068 3218870 : wghtns(2) = cp67 + (cp78-cp67)*(c1-fnidr)
2069 : ! wghtns(3) = cp33 + (cp22-cp33)*(c1-fnidr)
2070 3218870 : wghtns(3) = c1 - wghtns(2)
2071 :
2072 : ! find snow grain adjustment factor, dependent upon clear/overcast sky
2073 : ! estimate. comparisons with SNICAR show better agreement with DE when
2074 : ! this factor is included (clear sky near 1 and overcast near 0.8 give
2075 : ! best agreement). Multiply by rnsw here for efficiency.
2076 7668256 : do k = 1, nslyr
2077 4449386 : frsnw(k) = (fr_max*fnidr + fr_min*(c1-fnidr))*rsnw(k)
2078 7668256 : Sabs(k) = c0
2079 : enddo
2080 :
2081 : ! layer thicknesses
2082 : ! snow
2083 3218870 : dz = hs*rnslyr
2084 : ! for small enough snow thickness, ssl thickness half of top snow layer
2085 : !ech: note this is highly resolution dependent!
2086 3218870 : dzk(0) = min(hs_ssl, dz/c2)
2087 3218870 : dzk(1) = dz - dzk(0)
2088 3218870 : if (nslyr > 1) then
2089 1538145 : do k = 2, nslyr
2090 1538145 : dzk(k) = dz
2091 : enddo
2092 : endif
2093 :
2094 : ! ice
2095 3218870 : dz = hi*rnilyr
2096 : ! empirical reduction in sea ice ssl thickness for ice thinner than 1.5m;
2097 : ! factor of 30 gives best albedo comparison with limited observations
2098 3218870 : dz_ssl = hi_ssl
2099 : !ech: note hardwired parameters
2100 : ! if( hi < 1.5_dbl_kind ) dz_ssl = hi/30._dbl_kind
2101 3218870 : dz_ssl = min(hi_ssl, hi/30._dbl_kind)
2102 : ! set sea ice ssl thickness to half top layer if sea ice thin enough
2103 : !ech: note this is highly resolution dependent!
2104 3218870 : dz_ssl = min(dz_ssl, dz/c2)
2105 :
2106 3218870 : dzk(kii) = dz_ssl
2107 3218870 : dzk(kii+1) = dz - dz_ssl
2108 3218870 : if (kii+2 <= klev) then
2109 22532090 : do k = kii+2, klev
2110 22532090 : dzk(k) = dz
2111 : enddo
2112 : endif
2113 :
2114 : ! adjust sea ice iops with tuning parameters; tune only the
2115 : ! scattering coefficient by factors of R_ice, R_pnd, where
2116 : ! R values of +1 correspond approximately to +1 sigma changes in albedo, and
2117 : ! R values of -1 correspond approximately to -1 sigma changes in albedo
2118 : ! Note: the albedo change becomes non-linear for R values > +1 or < -1
2119 3218870 : if( R_ice >= c0 ) then
2120 12875480 : do ns = 1, nspint
2121 9656610 : sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fp_ice*R_ice)
2122 9656610 : ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
2123 9656610 : wi_ssl(ns) = sigp/ki_ssl(ns)
2124 9656610 : gi_ssl(ns) = gi_ssl_mn(ns)
2125 :
2126 9656610 : sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fp_ice*R_ice)
2127 9656610 : ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
2128 9656610 : wi_dl(ns) = sigp/ki_dl(ns)
2129 9656610 : gi_dl(ns) = gi_dl_mn(ns)
2130 :
2131 9656610 : sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fp_ice*R_ice)
2132 9656610 : ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
2133 9656610 : wi_int(ns) = sigp/ki_int(ns)
2134 12875480 : gi_int(ns) = gi_int_mn(ns)
2135 : enddo
2136 : else !if( R_ice < c0 ) then
2137 0 : do ns = 1, nspint
2138 0 : sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fm_ice*R_ice)
2139 0 : sigp = max(sigp, c0)
2140 0 : ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
2141 0 : wi_ssl(ns) = sigp/ki_ssl(ns)
2142 0 : gi_ssl(ns) = gi_ssl_mn(ns)
2143 :
2144 0 : sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fm_ice*R_ice)
2145 0 : sigp = max(sigp, c0)
2146 0 : ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
2147 0 : wi_dl(ns) = sigp/ki_dl(ns)
2148 0 : gi_dl(ns) = gi_dl_mn(ns)
2149 :
2150 0 : sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fm_ice*R_ice)
2151 0 : sigp = max(sigp, c0)
2152 0 : ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
2153 0 : wi_int(ns) = sigp/ki_int(ns)
2154 0 : gi_int(ns) = gi_int_mn(ns)
2155 : enddo
2156 : endif ! adjust ice iops
2157 :
2158 : ! adjust ponded ice iops with tuning parameters
2159 3218870 : if( R_pnd >= c0 ) then
2160 12875480 : do ns = 1, nspint
2161 9656610 : sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fp_pnd*R_pnd)
2162 9656610 : ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
2163 9656610 : wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
2164 9656610 : gi_p_ssl(ns) = gi_p_ssl_mn(ns)
2165 :
2166 9656610 : sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fp_pnd*R_pnd)
2167 9656610 : ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
2168 9656610 : wi_p_int(ns) = sigp/ki_p_int(ns)
2169 12875480 : gi_p_int(ns) = gi_p_int_mn(ns)
2170 : enddo
2171 : else !if( R_pnd < c0 ) then
2172 0 : do ns = 1, nspint
2173 0 : sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fm_pnd*R_pnd)
2174 0 : sigp = max(sigp, c0)
2175 0 : ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
2176 0 : wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
2177 0 : gi_p_ssl(ns) = gi_p_ssl_mn(ns)
2178 :
2179 0 : sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fm_pnd*R_pnd)
2180 0 : sigp = max(sigp, c0)
2181 0 : ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
2182 0 : wi_p_int(ns) = sigp/ki_p_int(ns)
2183 0 : gi_p_int(ns) = gi_p_int_mn(ns)
2184 : enddo
2185 : endif ! adjust ponded ice iops
2186 :
2187 : ! use srftyp to determine interface index of surface absorption
2188 3218870 : if (srftyp == 1) then
2189 : ! snow covered sea ice
2190 2092288 : ksrf = 1
2191 : else
2192 : ! bare sea ice or ponded ice
2193 1126582 : ksrf = nslyr + 2
2194 : endif
2195 :
2196 3218870 : if (tr_bgc_N .and. dEdd_algae) then ! compute kabs_chl for chlorophyll
2197 0 : do k = 0, klev
2198 0 : kabs_chl(1,k) = kchl_tab*zbio(nlt_chl_sw+k)
2199 : enddo
2200 : else
2201 3218870 : k = klev
2202 3218870 : kabs_chl(1,k) = kalg*(0.50_dbl_kind/dzk(k))
2203 : endif ! kabs_chl
2204 :
2205 : !mgf++
2206 3218870 : if (modal_aero) then
2207 0 : do k=0,klev
2208 0 : if (k < nslyr+1) then ! define indices for snow layer
2209 : ! use top rsnw, rhosnw for snow ssl and rest of top layer
2210 0 : ksnow = k - min(k-1,0)
2211 0 : tmp_gs = frsnw(ksnow)
2212 :
2213 : ! get grain size index:
2214 : ! works for 25 < snw_rds < 1625 um:
2215 0 : if (tmp_gs < 125) then
2216 0 : tmp1 = tmp_gs/50
2217 0 : k_bcini(k) = nint(tmp1)
2218 0 : elseif (tmp_gs < 175) then
2219 0 : k_bcini(k) = 2
2220 : else
2221 0 : tmp1 = (tmp_gs/250)+2
2222 0 : k_bcini(k) = nint(tmp1)
2223 : endif
2224 : else ! use the largest snow grain size for ice
2225 0 : k_bcini(k) = 8
2226 : endif
2227 : ! Set index corresponding to BC effective radius. Here,
2228 : ! asssume constant BC effective radius of 100nm
2229 : ! (corresponding to index 2)
2230 0 : k_bcins(k) = 2
2231 0 : k_bcexs(k) = 2
2232 :
2233 : ! check bounds:
2234 0 : if (k_bcini(k) < 1) k_bcini(k) = 1
2235 0 : if (k_bcini(k) > 8) k_bcini(k) = 8
2236 0 : if (k_bcins(k) < 1) k_bcins(k) = 1
2237 0 : if (k_bcins(k) > 10) k_bcins(k) = 10
2238 0 : if (k_bcexs(k) < 1) k_bcexs(k) = 1
2239 0 : if (k_bcexs(k) > 10) k_bcexs(k) = 10
2240 :
2241 : ! print ice radius index:
2242 : ! write(warnstr,*) subname, "MGFICE2:k, ice index= ",k, k_bcini(k)
2243 : ! call icepack_warnings_add(warnstr)
2244 : enddo ! k
2245 :
2246 0 : if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll
2247 0 : do n = 1,n_zaero
2248 0 : if (n == 1) then ! interstitial BC
2249 0 : do k = 0, klev
2250 0 : do ns = 1,nspint ! not weighted by aice
2251 0 : tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
2252 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2253 0 : wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
2254 0 : waer_bc_tab(ns,k_bcexs(k))* &
2255 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2256 0 : gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
2257 0 : waer_bc_tab(ns,k_bcexs(k))* &
2258 0 : gaer_bc_tab(ns,k_bcexs(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
2259 : enddo ! nspint
2260 : enddo
2261 0 : elseif (n==2) then ! within-ice BC
2262 0 : do k = 0, klev
2263 0 : do ns = 1,nspint
2264 0 : tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k)) * &
2265 0 : bcenh(ns,k_bcins(k),k_bcini(k))* &
2266 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2267 0 : wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
2268 0 : waer_bc_tab(ns,k_bcins(k))* &
2269 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2270 0 : gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
2271 0 : waer_bc_tab(ns,k_bcins(k))* &
2272 0 : gaer_bc_tab(ns,k_bcins(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
2273 : enddo ! nspint
2274 : enddo
2275 : else ! dust
2276 0 : do k = 0, klev
2277 0 : do ns = 1,nspint ! not weighted by aice
2278 0 : tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
2279 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2280 0 : wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
2281 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2282 0 : gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
2283 0 : gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
2284 : enddo ! nspint
2285 : enddo
2286 : endif !(n=1)
2287 : enddo ! n_zaero
2288 : endif ! tr_zaero and dEdd_algae
2289 :
2290 : else ! Bulk aerosol treatment
2291 3218870 : if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll
2292 0 : do n = 1,n_zaero ! multiply by aice?
2293 0 : do k = 0, klev
2294 0 : do ns = 1,nspint ! not weighted by aice
2295 0 : tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
2296 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2297 0 : wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
2298 0 : zbio(nlt_zaero_sw(n)+k)*dzk(k)
2299 0 : gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
2300 0 : gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
2301 : enddo ! nspint
2302 : enddo
2303 : enddo
2304 : endif !tr_zaero
2305 :
2306 : endif ! modal_aero
2307 :
2308 : !-----------------------------------------------------------------------
2309 :
2310 : ! begin spectral loop
2311 12875480 : do ns = 1, nspint
2312 :
2313 : ! set optical properties of air/snow/pond overlying sea ice
2314 : ! air
2315 9656610 : if( srftyp == 0 ) then
2316 9804729 : do k=0,nslyr
2317 7317678 : tau(k) = c0
2318 7317678 : w0(k) = c0
2319 9804729 : g(k) = c0
2320 : enddo
2321 : ! snow
2322 7169559 : else if( srftyp == 1 ) then
2323 : ! interpolate snow iops using input snow grain radius,
2324 : ! snow density and tabular data
2325 19214112 : do k=0,nslyr
2326 : ! use top rsnw, rhosnw for snow ssl and rest of top layer
2327 12937248 : ksnow = k - min(k-1,0)
2328 : ! find snow iops using input snow density and snow grain radius:
2329 12937248 : if( frsnw(ksnow) < rsnw_tab(1) ) then
2330 0 : Qs = Qs_tab(ns,1)
2331 0 : ws = ws_tab(ns,1)
2332 0 : gs = gs_tab(ns,1)
2333 12937248 : else if( frsnw(ksnow) >= rsnw_tab(nmbrad) ) then
2334 0 : Qs = Qs_tab(ns,nmbrad)
2335 0 : ws = ws_tab(ns,nmbrad)
2336 0 : gs = gs_tab(ns,nmbrad)
2337 : else
2338 : ! linear interpolation in rsnw
2339 413991936 : do nr=2,nmbrad
2340 539686812 : if( rsnw_tab(nr-1) <= frsnw(ksnow) .and. &
2341 151569372 : frsnw(ksnow) < rsnw_tab(nr)) then
2342 8944008 : delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / &
2343 17409252 : (rsnw_tab(nr) - rsnw_tab(nr-1))
2344 4472004 : Qs = Qs_tab(ns,nr-1)*(c1-delr) + &
2345 17409252 : Qs_tab(ns,nr)*delr
2346 4472004 : ws = ws_tab(ns,nr-1)*(c1-delr) + &
2347 17409252 : ws_tab(ns,nr)*delr
2348 4472004 : gs = gs_tab(ns,nr-1)*(c1-delr) + &
2349 17409252 : gs_tab(ns,nr)*delr
2350 : endif
2351 : enddo ! nr
2352 : endif
2353 4472004 : ks = Qs*((rhosnw(ksnow)/rhoi)*3._dbl_kind / &
2354 12937248 : (4._dbl_kind*frsnw(ksnow)*1.0e-6_dbl_kind))
2355 :
2356 12937248 : tau(k) = (ks + kabs_chl(ns,k))*dzk(k)
2357 12937248 : w0(k) = ks/(ks + kabs_chl(ns,k)) *ws
2358 19214112 : g(k) = gs
2359 : enddo ! k
2360 :
2361 :
2362 : ! aerosol in snow
2363 6276864 : if (tr_zaero .and. dEdd_algae) then
2364 0 : do k = 0,nslyr
2365 0 : gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
2366 0 : wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
2367 0 : g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
2368 0 : (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
2369 0 : w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
2370 0 : (tau(k) + tzaer(ns,k))
2371 0 : tau(k) = tau(k) + tzaer(ns,k)
2372 : enddo
2373 6276864 : elseif (tr_aero) then
2374 440610 : k = 0 ! snow SSL
2375 440610 : taer = c0
2376 440610 : waer = c0
2377 440610 : gaer = c0
2378 :
2379 881220 : do na=1,4*n_aero,4
2380 : ! mgf++
2381 881220 : if (modal_aero) then
2382 0 : if (na == 1) then
2383 : !interstitial BC
2384 : taer = taer + &
2385 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))
2386 : waer = waer + &
2387 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
2388 0 : waer_bc_tab(ns,k_bcexs(k))
2389 : gaer = gaer + &
2390 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
2391 0 : waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
2392 0 : elseif (na == 5)then
2393 : !within-ice BC
2394 : taer = taer + &
2395 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
2396 0 : bcenh(ns,k_bcins(k),k_bcini(k))
2397 : waer = waer + &
2398 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
2399 0 : waer_bc_tab(ns,k_bcins(k))
2400 : gaer = gaer + &
2401 0 : aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
2402 0 : waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
2403 : else
2404 : ! other species (dust)
2405 : taer = taer + &
2406 0 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
2407 : waer = waer + &
2408 0 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
2409 0 : waer_tab(ns,(1+(na-1)/4))
2410 : gaer = gaer + &
2411 0 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
2412 0 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2413 : endif
2414 : else
2415 : taer = taer + &
2416 440610 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
2417 : waer = waer + &
2418 0 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
2419 440610 : waer_tab(ns,(1+(na-1)/4))
2420 : gaer = gaer + &
2421 0 : aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
2422 440610 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2423 : endif !modal_aero
2424 : !mgf--
2425 : enddo ! na
2426 440610 : gaer = gaer/(waer+puny)
2427 440610 : waer = waer/(taer+puny)
2428 :
2429 881220 : do k=1,nslyr
2430 440610 : taer = c0
2431 440610 : waer = c0
2432 440610 : gaer = c0
2433 881220 : do na=1,4*n_aero,4
2434 881220 : if (modal_aero) then
2435 : !mgf++
2436 0 : if (na==1) then
2437 : ! interstitial BC
2438 : taer = taer + &
2439 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))
2440 : waer = waer + &
2441 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
2442 0 : waer_bc_tab(ns,k_bcexs(k))
2443 : gaer = gaer + &
2444 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
2445 0 : waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
2446 0 : elseif (na==5) then
2447 : ! within-ice BC
2448 : taer = taer + &
2449 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))*&
2450 0 : bcenh(ns,k_bcins(k),k_bcini(k))
2451 : waer = waer + &
2452 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
2453 0 : waer_bc_tab(ns,k_bcins(k))
2454 : gaer = gaer + &
2455 0 : (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
2456 0 : waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
2457 :
2458 : else
2459 : ! other species (dust)
2460 : taer = taer + &
2461 0 : (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))
2462 : waer = waer + &
2463 0 : (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
2464 0 : waer_tab(ns,(1+(na-1)/4))
2465 : gaer = gaer + &
2466 0 : (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
2467 0 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2468 : endif !(na==1)
2469 :
2470 : else
2471 : taer = taer + &
2472 440610 : (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))
2473 : waer = waer + &
2474 0 : (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
2475 440610 : waer_tab(ns,(1+(na-1)/4))
2476 : gaer = gaer + &
2477 0 : (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
2478 440610 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2479 : endif ! modal_aero
2480 : !mgf--
2481 : enddo ! na
2482 440610 : gaer = gaer/(waer+puny)
2483 440610 : waer = waer/(taer+puny)
2484 0 : g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
2485 440610 : (w0(k)*tau(k) + waer*taer)
2486 0 : w0(k) = (w0(k)*tau(k) + waer*taer) / &
2487 440610 : (tau(k) + taer)
2488 881220 : tau(k) = tau(k) + taer
2489 : enddo ! k
2490 : endif ! tr_aero
2491 :
2492 : ! pond
2493 : else !if( srftyp == 2 ) then
2494 : ! pond water layers evenly spaced
2495 892695 : dz = hp/(c1/rnslyr+c1)
2496 3642537 : do k=0,nslyr
2497 2749842 : tau(k) = kw(ns)*dz
2498 2749842 : w0(k) = ww(ns)
2499 3642537 : g(k) = gw(ns)
2500 : ! no aerosol in pond
2501 : enddo ! k
2502 : endif ! srftyp
2503 :
2504 : ! set optical properties of sea ice
2505 :
2506 : ! bare or snow-covered sea ice layers
2507 9656610 : if( srftyp <= 1 ) then
2508 : ! ssl
2509 8763915 : k = kii
2510 8763915 : tau(k) = (ki_ssl(ns)+kabs_chl(ns,k))*dzk(k)
2511 8763915 : w0(k) = ki_ssl(ns)/(ki_ssl(ns) + kabs_chl(ns,k))*wi_ssl(ns)
2512 8763915 : g(k) = gi_ssl(ns)
2513 : ! dl
2514 8763915 : k = kii + 1
2515 : ! scale dz for dl relative to 4 even-layer-thickness 1.5m case
2516 8763915 : fs = p25/rnilyr
2517 8763915 : tau(k) = (ki_dl(ns) + kabs_chl(ns,k)) *dzk(k)*fs
2518 8763915 : w0(k) = ki_dl(ns)/(ki_dl(ns) + kabs_chl(ns,k)) *wi_dl(ns)
2519 8763915 : g(k) = gi_dl(ns)
2520 : ! int above lowest layer
2521 8763915 : if (kii+2 <= klev-1) then
2522 52583490 : do k = kii+2, klev-1
2523 43819575 : tau(k) = (ki_int(ns) + kabs_chl(ns,k))*dzk(k)
2524 43819575 : w0(k) = ki_int(ns)/(ki_int(ns) + kabs_chl(ns,k)) *wi_int(ns)
2525 52583490 : g(k) = gi_int(ns)
2526 : enddo
2527 : endif
2528 : ! lowest layer
2529 8763915 : k = klev
2530 : ! add algae to lowest sea ice layer, visible only:
2531 8763915 : kabs = ki_int(ns)*(c1-wi_int(ns))
2532 8763915 : if( ns == 1 ) then
2533 : ! total layer absorption optical depth fixed at value
2534 : ! of kalg*0.50m, independent of actual layer thickness
2535 2921305 : kabs = kabs + kabs_chl(ns,k)
2536 : endif
2537 8763915 : sig = ki_int(ns)*wi_int(ns)
2538 8763915 : tau(k) = (kabs+sig)*dzk(k)
2539 8763915 : w0(k) = (sig/(sig+kabs))
2540 8763915 : g(k) = gi_int(ns)
2541 : ! aerosol in sea ice
2542 8763915 : if (tr_zaero .and. dEdd_algae) then
2543 0 : do k = kii, klev
2544 0 : gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
2545 0 : wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
2546 0 : g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
2547 0 : (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
2548 0 : w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
2549 0 : (tau(k) + tzaer(ns,k))
2550 0 : tau(k) = tau(k) + tzaer(ns,k)
2551 : enddo
2552 8763915 : elseif (tr_aero) then
2553 551421 : k = kii ! sea ice SSL
2554 551421 : taer = c0
2555 551421 : waer = c0
2556 551421 : gaer = c0
2557 1102842 : do na=1,4*n_aero,4
2558 : !mgf++
2559 1102842 : if (modal_aero) then
2560 0 : if (na==1) then
2561 : ! interstitial BC
2562 : taer = taer + &
2563 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))
2564 : waer = waer + &
2565 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
2566 0 : waer_bc_tab(ns,k_bcexs(k))
2567 : gaer = gaer + &
2568 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
2569 0 : waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
2570 0 : elseif (na==5) then
2571 : ! within-ice BC
2572 : taer = taer + &
2573 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
2574 0 : bcenh(ns,k_bcins(k),k_bcini(k))
2575 : waer = waer + &
2576 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
2577 0 : waer_bc_tab(ns,k_bcins(k))
2578 : gaer = gaer + &
2579 0 : aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
2580 0 : waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
2581 : else
2582 : ! other species (dust)
2583 : taer = taer + &
2584 0 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
2585 : waer = waer + &
2586 0 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
2587 0 : waer_tab(ns,(1+(na-1)/4))
2588 : gaer = gaer + &
2589 0 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
2590 0 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2591 : endif
2592 : else !bulk
2593 : taer = taer + &
2594 551421 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
2595 : waer = waer + &
2596 0 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
2597 551421 : waer_tab(ns,(1+(na-1)/4))
2598 : gaer = gaer + &
2599 0 : aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
2600 551421 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2601 : endif ! modal_aero
2602 : !mgf--
2603 : enddo ! na
2604 :
2605 551421 : gaer = gaer/(waer+puny)
2606 551421 : waer = waer/(taer+puny)
2607 0 : g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
2608 551421 : (w0(k)*tau(k) + waer*taer)
2609 0 : w0(k) = (w0(k)*tau(k) + waer*taer) / &
2610 551421 : (tau(k) + taer)
2611 551421 : tau(k) = tau(k) + taer
2612 4411368 : do k = kii+1, klev
2613 3859947 : taer = c0
2614 3859947 : waer = c0
2615 3859947 : gaer = c0
2616 7719894 : do na=1,4*n_aero,4
2617 : !mgf++
2618 7719894 : if (modal_aero) then
2619 0 : if (na==1) then
2620 : ! interstitial BC
2621 : taer = taer + &
2622 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))
2623 : waer = waer + &
2624 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
2625 0 : waer_bc_tab(ns,k_bcexs(k))
2626 : gaer = gaer + &
2627 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
2628 0 : waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
2629 0 : elseif (na==5) then
2630 : ! within-ice BC
2631 : taer = taer + &
2632 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
2633 0 : bcenh(ns,k_bcins(k),k_bcini(k))
2634 : waer = waer + &
2635 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
2636 0 : waer_bc_tab(ns,k_bcins(k))
2637 : gaer = gaer + &
2638 0 : (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
2639 0 : waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
2640 :
2641 : else
2642 : ! other species (dust)
2643 : taer = taer + &
2644 0 : (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))
2645 : waer = waer + &
2646 0 : (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
2647 0 : waer_tab(ns,(1+(na-1)/4))
2648 : gaer = gaer + &
2649 0 : (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
2650 0 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2651 : endif
2652 : else !bulk
2653 :
2654 : taer = taer + &
2655 3859947 : (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))
2656 : waer = waer + &
2657 0 : (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
2658 3859947 : waer_tab(ns,(1+(na-1)/4))
2659 : gaer = gaer + &
2660 0 : (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
2661 3859947 : waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
2662 : endif ! modal_aero
2663 : !mgf--
2664 : enddo ! na
2665 3859947 : gaer = gaer/(waer+puny)
2666 3859947 : waer = waer/(taer+puny)
2667 0 : g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
2668 3859947 : (w0(k)*tau(k) + waer*taer)
2669 0 : w0(k) = (w0(k)*tau(k) + waer*taer) / &
2670 3859947 : (tau(k) + taer)
2671 4411368 : tau(k) = tau(k) + taer
2672 : enddo ! k
2673 : endif ! tr_aero
2674 :
2675 : ! sea ice layers under ponds
2676 : else !if( srftyp == 2 ) then
2677 892695 : k = kii
2678 892695 : tau(k) = ki_p_ssl(ns)*dzk(k)
2679 892695 : w0(k) = wi_p_ssl(ns)
2680 892695 : g(k) = gi_p_ssl(ns)
2681 892695 : k = kii + 1
2682 892695 : tau(k) = ki_p_int(ns)*dzk(k)
2683 892695 : w0(k) = wi_p_int(ns)
2684 892695 : g(k) = gi_p_int(ns)
2685 892695 : if (kii+2 <= klev) then
2686 6248865 : do k = kii+2, klev
2687 5356170 : tau(k) = ki_p_int(ns)*dzk(k)
2688 5356170 : w0(k) = wi_p_int(ns)
2689 6248865 : g(k) = gi_p_int(ns)
2690 : enddo ! k
2691 : endif
2692 : ! adjust pond iops if pond depth within specified range
2693 892695 : if( hpmin <= hp .and. hp <= hp0 ) then
2694 651600 : k = kii
2695 651600 : sig_i = ki_ssl(ns)*wi_ssl(ns)
2696 651600 : sig_p = ki_p_ssl(ns)*wi_p_ssl(ns)
2697 651600 : sig = sig_i + (sig_p-sig_i)*(hp/hp0)
2698 651600 : kext = sig + ki_p_ssl(ns)*(c1-wi_p_ssl(ns))
2699 651600 : tau(k) = kext*dzk(k)
2700 651600 : w0(k) = sig/kext
2701 651600 : g(k) = gi_p_int(ns)
2702 651600 : k = kii + 1
2703 : ! scale dz for dl relative to 4 even-layer-thickness 1.5m case
2704 651600 : fs = p25/rnilyr
2705 651600 : sig_i = ki_dl(ns)*wi_dl(ns)*fs
2706 651600 : sig_p = ki_p_int(ns)*wi_p_int(ns)
2707 651600 : sig = sig_i + (sig_p-sig_i)*(hp/hp0)
2708 651600 : kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
2709 651600 : tau(k) = kext*dzk(k)
2710 651600 : w0(k) = sig/kext
2711 651600 : g(k) = gi_p_int(ns)
2712 651600 : if (kii+2 <= klev) then
2713 4561200 : do k = kii+2, klev
2714 3909600 : sig_i = ki_int(ns)*wi_int(ns)
2715 3909600 : sig_p = ki_p_int(ns)*wi_p_int(ns)
2716 3909600 : sig = sig_i + (sig_p-sig_i)*(hp/hp0)
2717 3909600 : kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
2718 3909600 : tau(k) = kext*dzk(k)
2719 3909600 : w0(k) = sig/kext
2720 4561200 : g(k) = gi_p_int(ns)
2721 : enddo ! k
2722 : endif
2723 : endif ! small pond depth transition to bare sea ice
2724 : endif ! srftyp
2725 :
2726 : ! set reflectivities for ocean underlying sea ice
2727 9656610 : rns = real(ns-1, kind=dbl_kind)
2728 9656610 : albodr = cp01 * (c1 - min(rns, c1))
2729 9656610 : albodf = cp01 * (c1 - min(rns, c1))
2730 :
2731 : ! layer input properties now completely specified: tau, w0, g,
2732 : ! albodr, albodf; now compute the Delta-Eddington solution
2733 : ! reflectivities and transmissivities for each layer; then,
2734 : ! combine the layers going downwards accounting for multiple
2735 : ! scattering between layers, and finally start from the
2736 : ! underlying ocean and combine successive layers upwards to
2737 : ! the surface; see comments in solution_dEdd for more details.
2738 :
2739 : call solution_dEdd &
2740 : (coszen, srftyp, klev, klevp, nslyr, &
2741 : tau, w0, g, albodr, albodf, &
2742 : trndir, trntdr, trndif, rupdir, rupdif, &
2743 9656610 : rdndif)
2744 9656610 : if (icepack_warnings_aborted(subname)) return
2745 :
2746 : ! the interface reflectivities and transmissivities required
2747 : ! to evaluate interface fluxes are returned from solution_dEdd;
2748 : ! now compute up and down fluxes for each interface, using the
2749 : ! combined layer properties at each interface:
2750 : !
2751 : ! layers interface
2752 : !
2753 : ! --------------------- k
2754 : ! k
2755 : ! ---------------------
2756 :
2757 119570868 : do k = 0, klevp
2758 : ! interface scattering
2759 109914258 : refk = c1/(c1 - rdndif(k)*rupdif(k))
2760 : ! dir tran ref from below times interface scattering, plus diff
2761 : ! tran and ref from below times interface scattering
2762 : ! fdirup(k) = (trndir(k)*rupdir(k) + &
2763 : ! (trntdr(k)-trndir(k)) &
2764 : ! *rupdif(k))*refk
2765 : ! dir tran plus total diff trans times interface scattering plus
2766 : ! dir tran with up dir ref and down dif ref times interface scattering
2767 : ! fdirdn(k) = trndir(k) + (trntdr(k) &
2768 : ! - trndir(k) + trndir(k) &
2769 : ! *rupdir(k)*rdndif(k))*refk
2770 : ! diffuse tran ref from below times interface scattering
2771 : ! fdifup(k) = trndif(k)*rupdif(k)*refk
2772 : ! diffuse tran times interface scattering
2773 : ! fdifdn(k) = trndif(k)*refk
2774 :
2775 : ! dfdir = fdirdn - fdirup
2776 39009405 : dfdir(k) = trndir(k) &
2777 39009405 : + (trntdr(k)-trndir(k)) * (c1 - rupdif(k)) * refk &
2778 109914258 : - trndir(k)*rupdir(k) * (c1 - rdndif(k)) * refk
2779 109914258 : if (dfdir(k) < puny) dfdir(k) = c0 !echmod necessary?
2780 : ! dfdif = fdifdn - fdifup
2781 109914258 : dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk
2782 119570868 : if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary?
2783 : enddo ! k
2784 :
2785 : ! calculate final surface albedos and fluxes-
2786 : ! all absorbed flux above ksrf is included in surface absorption
2787 :
2788 12875480 : if( ns == 1) then ! visible
2789 :
2790 3218870 : swdr = swvdr
2791 3218870 : swdf = swvdf
2792 3218870 : avdr = rupdir(0)
2793 3218870 : avdf = rupdif(0)
2794 :
2795 3218870 : tmp_0 = dfdir(0 )*swdr + dfdif(0 )*swdf
2796 3218870 : tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
2797 3218870 : tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
2798 :
2799 : ! for layer biology: save visible only
2800 28969830 : do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering
2801 28969830 : fthrul(k-nslyr-1) = dfdir(k)*swdr + dfdif(k)*swdf
2802 : enddo
2803 :
2804 3218870 : fsfc = fsfc + tmp_0 - tmp_ks
2805 3218870 : fint = fint + tmp_ks - tmp_kl
2806 3218870 : fthru = fthru + tmp_kl
2807 :
2808 : ! if snow covered ice, set snow internal absorption; else, Sabs=0
2809 3218870 : if( srftyp == 1 ) then
2810 2092288 : ki = 0
2811 4312416 : do k=1,nslyr
2812 : ! skip snow SSL, since SSL absorption included in the surface
2813 : ! absorption fsfc above
2814 2220128 : km = k
2815 2220128 : kp = km + 1
2816 2220128 : ki = ki + 1
2817 777294 : Sabs(ki) = Sabs(ki) &
2818 777294 : + dfdir(km)*swdr + dfdif(km)*swdf &
2819 4312416 : - (dfdir(kp)*swdr + dfdif(kp)*swdf)
2820 : enddo ! k
2821 : endif
2822 :
2823 : ! complex indexing to insure proper absorptions for sea ice
2824 3218870 : ki = 0
2825 25750960 : do k=nslyr+2,nslyr+1+nilyr
2826 : ! for bare ice, DL absorption for sea ice layer 1
2827 22532090 : km = k
2828 22532090 : kp = km + 1
2829 : ! modify for top sea ice layer for snow over sea ice
2830 22532090 : if( srftyp == 1 ) then
2831 : ! must add SSL and DL absorption for sea ice layer 1
2832 14646016 : if( k == nslyr+2 ) then
2833 2092288 : km = k - 1
2834 2092288 : kp = km + 2
2835 : endif
2836 : endif
2837 22532090 : ki = ki + 1
2838 7999131 : Iabs(ki) = Iabs(ki) &
2839 7999131 : + dfdir(km)*swdr + dfdif(km)*swdf &
2840 25750960 : - (dfdir(kp)*swdr + dfdif(kp)*swdf)
2841 : enddo ! k
2842 :
2843 : else !if(ns > 1) then ! near IR
2844 :
2845 6437740 : swdr = swidr
2846 6437740 : swdf = swidf
2847 :
2848 : ! let fr1 = alb_1*swd*wght1 and fr2 = alb_2*swd*wght2 be the ns=2,3
2849 : ! reflected fluxes respectively, where alb_1, alb_2 are the band
2850 : ! albedos, swd = nir incident shortwave flux, and wght1, wght2 are
2851 : ! the 2,3 band weights. thus, the total reflected flux is:
2852 : ! fr = fr1 + fr2 = alb_1*swd*wght1 + alb_2*swd*wght2 hence, the
2853 : ! 2,3 nir band albedo is alb = fr/swd = alb_1*wght1 + alb_2*wght2
2854 :
2855 6437740 : aidr = aidr + rupdir(0)*wghtns(ns)
2856 6437740 : aidf = aidf + rupdif(0)*wghtns(ns)
2857 :
2858 6437740 : tmp_0 = dfdir(0 )*swdr + dfdif(0 )*swdf
2859 6437740 : tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
2860 6437740 : tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
2861 :
2862 6437740 : tmp_0 = tmp_0 * wghtns(ns)
2863 6437740 : tmp_ks = tmp_ks * wghtns(ns)
2864 6437740 : tmp_kl = tmp_kl * wghtns(ns)
2865 :
2866 6437740 : fsfc = fsfc + tmp_0 - tmp_ks
2867 6437740 : fint = fint + tmp_ks - tmp_kl
2868 6437740 : fthru = fthru + tmp_kl
2869 :
2870 : ! if snow covered ice, set snow internal absorption; else, Sabs=0
2871 6437740 : if( srftyp == 1 ) then
2872 4184576 : ki = 0
2873 8624832 : do k=1,nslyr
2874 : ! skip snow SSL, since SSL absorption included in the surface
2875 : ! absorption fsfc above
2876 4440256 : km = k
2877 4440256 : kp = km + 1
2878 4440256 : ki = ki + 1
2879 1554588 : Sabs(ki) = Sabs(ki) &
2880 1554588 : + (dfdir(km)*swdr + dfdif(km)*swdf &
2881 1554588 : - (dfdir(kp)*swdr + dfdif(kp)*swdf)) &
2882 8624832 : * wghtns(ns)
2883 : enddo ! k
2884 : endif
2885 :
2886 : ! complex indexing to insure proper absorptions for sea ice
2887 6437740 : ki = 0
2888 51501920 : do k=nslyr+2,nslyr+1+nilyr
2889 : ! for bare ice, DL absorption for sea ice layer 1
2890 45064180 : km = k
2891 45064180 : kp = km + 1
2892 : ! modify for top sea ice layer for snow over sea ice
2893 45064180 : if( srftyp == 1 ) then
2894 : ! must add SSL and DL absorption for sea ice layer 1
2895 29292032 : if( k == nslyr+2 ) then
2896 4184576 : km = k - 1
2897 4184576 : kp = km + 2
2898 : endif
2899 : endif
2900 45064180 : ki = ki + 1
2901 15998262 : Iabs(ki) = Iabs(ki) &
2902 15998262 : + (dfdir(km)*swdr + dfdif(km)*swdf &
2903 15998262 : - (dfdir(kp)*swdr + dfdif(kp)*swdf)) &
2904 51501920 : * wghtns(ns)
2905 : enddo ! k
2906 :
2907 : endif ! ns = 1, ns > 1
2908 :
2909 : enddo ! end spectral loop ns
2910 :
2911 : ! accumulate fluxes over bare sea ice
2912 3218870 : alvdr = avdr
2913 3218870 : alvdf = avdf
2914 3218870 : alidr = aidr
2915 3218870 : alidf = aidf
2916 3218870 : fswsfc = fswsfc + fsfc *fi
2917 3218870 : fswint = fswint + fint *fi
2918 3218870 : fswthru = fswthru + fthru*fi
2919 :
2920 7668256 : do k = 1, nslyr
2921 7668256 : Sswabs(k) = Sswabs(k) + Sabs(k)*fi
2922 : enddo ! k
2923 :
2924 25750960 : do k = 1, nilyr
2925 22532090 : Iswabs(k) = Iswabs(k) + Iabs(k)*fi
2926 :
2927 : ! bgc layer
2928 22532090 : fswpenl(k) = fswpenl(k) + fthrul(k)* fi
2929 25750960 : if (k == nilyr) then
2930 3218870 : fswpenl(k+1) = fswpenl(k+1) + fthrul(k+1)*fi
2931 : endif
2932 : enddo ! k
2933 :
2934 : !----------------------------------------------------------------
2935 : ! if ice has zero heat capacity, no SW can be absorbed
2936 : ! in the ice/snow interior, so add to surface absorption.
2937 : ! Note: nilyr = nslyr = 1 for this case
2938 : !----------------------------------------------------------------
2939 :
2940 3218870 : if (.not. heat_capacity) then
2941 :
2942 : ! SW absorbed at snow/ice surface
2943 0 : fswsfc = fswsfc + Iswabs(1) + Sswabs(1)
2944 :
2945 : ! SW absorbed in ice interior
2946 0 : fswint = c0
2947 0 : Iswabs(1) = c0
2948 0 : Sswabs(1) = c0
2949 :
2950 : endif ! heat_capacity
2951 :
2952 : end subroutine compute_dEdd
2953 :
2954 : !=======================================================================
2955 : !
2956 : ! Given input vertical profiles of optical properties, evaluate the
2957 : ! monochromatic Delta-Eddington solution.
2958 : !
2959 : ! author: Bruce P. Briegleb, NCAR
2960 : ! 2013: E Hunke merged with NCAR version
2961 9656610 : subroutine solution_dEdd &
2962 : (coszen, srftyp, klev, klevp, nslyr, &
2963 9656610 : tau, w0, g, albodr, albodf, &
2964 9656610 : trndir, trntdr, trndif, rupdir, rupdif, &
2965 9656610 : rdndif)
2966 :
2967 : real (kind=dbl_kind), intent(in) :: &
2968 : coszen ! cosine solar zenith angle
2969 :
2970 : integer (kind=int_kind), intent(in) :: &
2971 : srftyp , & ! surface type over ice: (0=air, 1=snow, 2=pond)
2972 : klev , & ! number of radiation layers - 1
2973 : klevp , & ! number of radiation interfaces - 1
2974 : ! (0 layer is included also)
2975 : nslyr ! number of snow layers
2976 :
2977 : real (kind=dbl_kind), dimension(0:klev), intent(in) :: &
2978 : tau , & ! layer extinction optical depth
2979 : w0 , & ! layer single scattering albedo
2980 : g ! layer asymmetry parameter
2981 :
2982 : real (kind=dbl_kind), intent(in) :: &
2983 : albodr , & ! ocean albedo to direct rad
2984 : albodf ! ocean albedo to diffuse rad
2985 :
2986 : ! following arrays are defined at model interfaces; 0 is the top of the
2987 : ! layer above the sea ice; klevp is the sea ice/ocean interface.
2988 : real (kind=dbl_kind), dimension (0:klevp), intent(out) :: &
2989 : trndir , & ! solar beam down transmission from top
2990 : trntdr , & ! total transmission to direct beam for layers above
2991 : trndif , & ! diffuse transmission to diffuse beam for layers above
2992 : rupdir , & ! reflectivity to direct radiation for layers below
2993 : rupdif , & ! reflectivity to diffuse radiation for layers below
2994 : rdndif ! reflectivity to diffuse radiation for layers above
2995 :
2996 : !-----------------------------------------------------------------------
2997 : !
2998 : ! Delta-Eddington solution for snow/air/pond over sea ice
2999 : !
3000 : ! Generic solution for a snow/air/pond input column of klev+1 layers,
3001 : ! with srftyp determining at what interface fresnel refraction occurs.
3002 : !
3003 : ! Computes layer reflectivities and transmissivities, from the top down
3004 : ! to the lowest interface using the Delta-Eddington solutions for each
3005 : ! layer; combines layers from top down to lowest interface, and from the
3006 : ! lowest interface (underlying ocean) up to the top of the column.
3007 : !
3008 : ! Note that layer diffuse reflectivity and transmissivity are computed
3009 : ! by integrating the direct over several gaussian angles. This is
3010 : ! because the diffuse reflectivity expression sometimes is negative,
3011 : ! but the direct reflectivity is always well-behaved. We assume isotropic
3012 : ! radiation in the upward and downward hemispheres for this integration.
3013 : !
3014 : ! Assumes monochromatic (spectrally uniform) properties across a band
3015 : ! for the input optical parameters.
3016 : !
3017 : ! If total transmission of the direct beam to the interface above a particular
3018 : ! layer is less than trmin, then no further Delta-Eddington solutions are
3019 : ! evaluated for layers below.
3020 : !
3021 : ! The following describes how refraction is handled in the calculation.
3022 : !
3023 : ! First, we assume that radiation is refracted when entering either
3024 : ! sea ice at the base of the surface scattering layer, or water (i.e. melt
3025 : ! pond); we assume that radiation does not refract when entering snow, nor
3026 : ! upon entering sea ice from a melt pond, nor upon entering the underlying
3027 : ! ocean from sea ice.
3028 : !
3029 : ! To handle refraction, we define a "fresnel" layer, which physically
3030 : ! is of neglible thickness and is non-absorbing, which can be combined to
3031 : ! any sea ice layer or top of melt pond. The fresnel layer accounts for
3032 : ! refraction of direct beam and associated reflection and transmission for
3033 : ! solar radiation. A fresnel layer is combined with the top of a melt pond
3034 : ! or to the surface scattering layer of sea ice if no melt pond lies over it.
3035 : !
3036 : ! Some caution must be exercised for the fresnel layer, because any layer
3037 : ! to which it is combined is no longer a homogeneous layer, as are all other
3038 : ! individual layers. For all other layers for example, the direct and diffuse
3039 : ! reflectivities/transmissivities (R/T) are the same for radiation above or
3040 : ! below the layer. This is the meaning of homogeneous! But for the fresnel
3041 : ! layer this is not so. Thus, the R/T for this layer must be distinguished
3042 : ! for radiation above from that from radiation below. For generality, we
3043 : ! treat all layers to be combined as inhomogeneous.
3044 : !
3045 : !-----------------------------------------------------------------------
3046 :
3047 : ! local variables
3048 :
3049 : integer (kind=int_kind) :: &
3050 : kfrsnl ! radiation interface index for fresnel layer
3051 :
3052 : ! following variables are defined for each layer; 0 refers to the top
3053 : ! layer. In general we must distinguish directions above and below in
3054 : ! the diffuse reflectivity and transmissivity, as layers are not assumed
3055 : ! to be homogeneous (apart from the single layer Delta-Edd solutions);
3056 : ! the direct is always from above.
3057 : real (kind=dbl_kind), dimension (0:klev) :: &
3058 51466227 : rdir , & ! layer reflectivity to direct radiation
3059 51466227 : rdif_a , & ! layer reflectivity to diffuse radiation from above
3060 51466227 : rdif_b , & ! layer reflectivity to diffuse radiation from below
3061 51466227 : tdir , & ! layer transmission to direct radiation (solar beam + diffuse)
3062 51466227 : tdif_a , & ! layer transmission to diffuse radiation from above
3063 51466227 : tdif_b , & ! layer transmission to diffuse radiation from below
3064 51466227 : trnlay ! solar beam transm for layer (direct beam only)
3065 :
3066 : integer (kind=int_kind) :: &
3067 : k ! level index
3068 :
3069 : real (kind=dbl_kind), parameter :: &
3070 : trmin = 0.001_dbl_kind ! minimum total transmission allowed
3071 : ! total transmission is that due to the direct beam; i.e. it includes
3072 : ! both the directly transmitted solar beam and the diffuse downwards
3073 : ! transmitted radiation resulting from scattering out of the direct beam
3074 : real (kind=dbl_kind) :: &
3075 3428199 : tautot , & ! layer optical depth
3076 3428199 : wtot , & ! layer single scattering albedo
3077 3428199 : gtot , & ! layer asymmetry parameter
3078 3428199 : ftot , & ! layer forward scattering fraction
3079 3428199 : ts , & ! layer scaled extinction optical depth
3080 3428199 : ws , & ! layer scaled single scattering albedo
3081 3428199 : gs , & ! layer scaled asymmetry parameter
3082 3428199 : rintfc , & ! reflection (multiple) at an interface
3083 3428199 : refkp1 , & ! interface multiple scattering for k+1
3084 3428199 : refkm1 , & ! interface multiple scattering for k-1
3085 3428199 : tdrrdir , & ! direct tran times layer direct ref
3086 3428199 : tdndif ! total down diffuse = tot tran - direct tran
3087 :
3088 : ! perpendicular and parallel relative to plane of incidence and scattering
3089 : real (kind=dbl_kind) :: &
3090 3428199 : R1 , & ! perpendicular polarization reflection amplitude
3091 3428199 : R2 , & ! parallel polarization reflection amplitude
3092 3428199 : T1 , & ! perpendicular polarization transmission amplitude
3093 3428199 : T2 , & ! parallel polarization transmission amplitude
3094 3428199 : Rf_dir_a , & ! fresnel reflection to direct radiation
3095 3428199 : Tf_dir_a , & ! fresnel transmission to direct radiation
3096 3428199 : Rf_dif_a , & ! fresnel reflection to diff radiation from above
3097 3428199 : Rf_dif_b , & ! fresnel reflection to diff radiation from below
3098 3428199 : Tf_dif_a , & ! fresnel transmission to diff radiation from above
3099 3428199 : Tf_dif_b ! fresnel transmission to diff radiation from below
3100 :
3101 : ! refractive index for sea ice, water; pre-computed, band-independent,
3102 : ! diffuse fresnel reflectivities
3103 : real (kind=dbl_kind), parameter :: &
3104 : refindx = 1.310_dbl_kind , & ! refractive index of sea ice (water also)
3105 : cp063 = 0.063_dbl_kind , & ! diffuse fresnel reflectivity from above
3106 : cp455 = 0.455_dbl_kind ! diffuse fresnel reflectivity from below
3107 :
3108 : real (kind=dbl_kind) :: &
3109 3428199 : mu0 , & ! cosine solar zenith angle incident
3110 3428199 : mu0nij ! cosine solar zenith angle in medium below fresnel level
3111 :
3112 : real (kind=dbl_kind) :: &
3113 3428199 : mu0n ! cosine solar zenith angle in medium
3114 :
3115 : real (kind=dbl_kind) :: &
3116 3428199 : alp , & ! temporary for alpha
3117 3428199 : gam , & ! temporary for agamm
3118 3428199 : lm , & ! temporary for el
3119 3428199 : mu , & ! temporary for gauspt
3120 3428199 : ne , & ! temporary for n
3121 3428199 : ue , & ! temporary for u
3122 3428199 : extins , & ! extinction
3123 3428199 : amg , & ! alp - gam
3124 3428199 : apg ! alp + gam
3125 :
3126 : integer (kind=int_kind), parameter :: &
3127 : ngmax = 8 ! number of gaussian angles in hemisphere
3128 :
3129 : real (kind=dbl_kind), dimension (ngmax), parameter :: &
3130 : gauspt & ! gaussian angles (radians)
3131 : = (/ .9894009_dbl_kind, .9445750_dbl_kind, &
3132 : .8656312_dbl_kind, .7554044_dbl_kind, &
3133 : .6178762_dbl_kind, .4580168_dbl_kind, &
3134 : .2816036_dbl_kind, .0950125_dbl_kind/), &
3135 : gauswt & ! gaussian weights
3136 : = (/ .0271525_dbl_kind, .0622535_dbl_kind, &
3137 : .0951585_dbl_kind, .1246290_dbl_kind, &
3138 : .1495960_dbl_kind, .1691565_dbl_kind, &
3139 : .1826034_dbl_kind, .1894506_dbl_kind/)
3140 :
3141 : integer (kind=int_kind) :: &
3142 : ng ! gaussian integration index
3143 :
3144 : real (kind=dbl_kind) :: &
3145 3428199 : gwt , & ! gaussian weight
3146 3428199 : swt , & ! sum of weights
3147 3428199 : trn , & ! layer transmission
3148 3428199 : rdr , & ! rdir for gaussian integration
3149 3428199 : tdr , & ! tdir for gaussian integration
3150 3428199 : smr , & ! accumulator for rdif gaussian integration
3151 3428199 : smt ! accumulator for tdif gaussian integration
3152 :
3153 : real (kind=dbl_kind) :: &
3154 3428199 : exp_min ! minimum exponential value
3155 :
3156 : character(len=*),parameter :: subname='(solution_dEdd)'
3157 :
3158 : !-----------------------------------------------------------------------
3159 :
3160 119570868 : do k = 0, klevp
3161 109914258 : trndir(k) = c0
3162 109914258 : trntdr(k) = c0
3163 109914258 : trndif(k) = c0
3164 109914258 : rupdir(k) = c0
3165 109914258 : rupdif(k) = c0
3166 119570868 : rdndif(k) = c0
3167 : enddo
3168 :
3169 : ! initialize top interface of top layer
3170 9656610 : trndir(0) = c1
3171 9656610 : trntdr(0) = c1
3172 9656610 : trndif(0) = c1
3173 9656610 : rdndif(0) = c0
3174 :
3175 : ! mu0 is cosine solar zenith angle above the fresnel level; make
3176 : ! sure mu0 is large enough for stable and meaningful radiation
3177 : ! solution: .01 is like sun just touching horizon with its lower edge
3178 9656610 : mu0 = max(coszen,p01)
3179 :
3180 : ! mu0n is cosine solar zenith angle used to compute the layer
3181 : ! Delta-Eddington solution; it is initially computed to be the
3182 : ! value below the fresnel level, i.e. the cosine solar zenith
3183 : ! angle below the fresnel level for the refracted solar beam:
3184 9656610 : mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
3185 :
3186 : ! compute level of fresnel refraction
3187 : ! if ponded sea ice, fresnel level is the top of the pond.
3188 9656610 : kfrsnl = 0
3189 : ! if snow over sea ice or bare sea ice, fresnel level is
3190 : ! at base of sea ice SSL (and top of the sea ice DL); the
3191 : ! snow SSL counts for one, then the number of snow layers,
3192 : ! then the sea ice SSL which also counts for one:
3193 9656610 : if( srftyp < 2 ) kfrsnl = nslyr + 2
3194 :
3195 : ! proceed down one layer at a time; if the total transmission to
3196 : ! the interface just above a given layer is less than trmin, then no
3197 : ! Delta-Eddington computation for that layer is done.
3198 :
3199 : ! begin main level loop
3200 109914258 : do k = 0, klev
3201 :
3202 : ! initialize all layer apparent optical properties to 0
3203 100257648 : rdir (k) = c0
3204 100257648 : rdif_a(k) = c0
3205 100257648 : rdif_b(k) = c0
3206 100257648 : tdir (k) = c0
3207 100257648 : tdif_a(k) = c0
3208 100257648 : tdif_b(k) = c0
3209 100257648 : trnlay(k) = c0
3210 :
3211 : ! compute next layer Delta-eddington solution only if total transmission
3212 : ! of radiation to the interface just above the layer exceeds trmin.
3213 :
3214 100257648 : if (trntdr(k) > trmin ) then
3215 :
3216 : ! calculation over layers with penetrating radiation
3217 :
3218 52146940 : tautot = tau(k)
3219 52146940 : wtot = w0(k)
3220 52146940 : gtot = g(k)
3221 52146940 : ftot = gtot*gtot
3222 :
3223 52146940 : ts = taus(wtot,ftot,tautot)
3224 52146940 : ws = omgs(wtot,ftot)
3225 52146940 : gs = asys(gtot,ftot)
3226 52146940 : lm = el(ws,gs)
3227 52146940 : ue = u(ws,gs,lm)
3228 :
3229 52146940 : mu0n = mu0nij
3230 : ! if level k is above fresnel level and the cell is non-pond, use the
3231 : ! non-refracted beam instead
3232 52146940 : if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
3233 :
3234 52146940 : exp_min = min(exp_argmax,lm*ts)
3235 52146940 : extins = exp(-exp_min)
3236 52146940 : ne = n(ue,extins)
3237 :
3238 : ! first calculation of rdif, tdif using Delta-Eddington formulas
3239 : ! rdif_a(k) = (ue+c1)*(ue-c1)*(c1/extins - extins)/ne
3240 52146940 : rdif_a(k) = (ue**2-c1)*(c1/extins - extins)/ne
3241 52146940 : tdif_a(k) = c4*ue/ne
3242 :
3243 : ! evaluate rdir,tdir for direct beam
3244 52146940 : exp_min = min(exp_argmax,ts/mu0n)
3245 52146940 : trnlay(k) = exp(-exp_min)
3246 52146940 : alp = alpha(ws,mu0n,gs,lm)
3247 52146940 : gam = agamm(ws,mu0n,gs,lm)
3248 52146940 : apg = alp + gam
3249 52146940 : amg = alp - gam
3250 52146940 : rdir(k) = apg*rdif_a(k) + amg*(tdif_a(k)*trnlay(k) - c1)
3251 52146940 : tdir(k) = apg*tdif_a(k) + (amg* rdif_a(k)-apg+c1)*trnlay(k)
3252 :
3253 : ! recalculate rdif,tdif using direct angular integration over rdir,tdir,
3254 : ! since Delta-Eddington rdif formula is not well-behaved (it is usually
3255 : ! biased low and can even be negative); use ngmax angles and gaussian
3256 : ! integration for most accuracy:
3257 52146940 : R1 = rdif_a(k) ! use R1 as temporary
3258 52146940 : T1 = tdif_a(k) ! use T1 as temporary
3259 52146940 : swt = c0
3260 52146940 : smr = c0
3261 52146940 : smt = c0
3262 469322460 : do ng=1,ngmax
3263 417175520 : mu = gauspt(ng)
3264 417175520 : gwt = gauswt(ng)
3265 417175520 : swt = swt + mu*gwt
3266 417175520 : exp_min = min(exp_argmax,ts/mu)
3267 417175520 : trn = exp(-exp_min)
3268 417175520 : alp = alpha(ws,mu,gs,lm)
3269 417175520 : gam = agamm(ws,mu,gs,lm)
3270 417175520 : apg = alp + gam
3271 417175520 : amg = alp - gam
3272 417175520 : rdr = apg*R1 + amg*T1*trn - amg
3273 417175520 : tdr = apg*T1 + amg*R1*trn - apg*trn + trn
3274 417175520 : smr = smr + mu*rdr*gwt
3275 469322460 : smt = smt + mu*tdr*gwt
3276 : enddo ! ng
3277 52146940 : rdif_a(k) = smr/swt
3278 52146940 : tdif_a(k) = smt/swt
3279 :
3280 : ! homogeneous layer
3281 52146940 : rdif_b(k) = rdif_a(k)
3282 52146940 : tdif_b(k) = tdif_a(k)
3283 :
3284 : ! add fresnel layer to top of desired layer if either
3285 : ! air or snow overlies ice; we ignore refraction in ice
3286 : ! if a melt pond overlies it:
3287 :
3288 52146940 : if( k == kfrsnl ) then
3289 : ! compute fresnel reflection and transmission amplitudes
3290 : ! for two polarizations: 1=perpendicular and 2=parallel to
3291 : ! the plane containing incident, reflected and refracted rays.
3292 : R1 = (mu0 - refindx*mu0n) / &
3293 5612695 : (mu0 + refindx*mu0n)
3294 : R2 = (refindx*mu0 - mu0n) / &
3295 5612695 : (refindx*mu0 + mu0n)
3296 : T1 = c2*mu0 / &
3297 5612695 : (mu0 + refindx*mu0n)
3298 : T2 = c2*mu0 / &
3299 5612695 : (refindx*mu0 + mu0n)
3300 :
3301 : ! unpolarized light for direct beam
3302 5612695 : Rf_dir_a = p5 * (R1*R1 + R2*R2)
3303 5612695 : Tf_dir_a = p5 * (T1*T1 + T2*T2)*refindx*mu0n/mu0
3304 :
3305 : ! precalculated diffuse reflectivities and transmissivities
3306 : ! for incident radiation above and below fresnel layer, using
3307 : ! the direct albedos and accounting for complete internal
3308 : ! reflection from below; precalculated because high order
3309 : ! number of gaussian points (~256) is required for convergence:
3310 :
3311 : ! above
3312 5612695 : Rf_dif_a = cp063
3313 5612695 : Tf_dif_a = c1 - Rf_dif_a
3314 : ! below
3315 5612695 : Rf_dif_b = cp455
3316 5612695 : Tf_dif_b = c1 - Rf_dif_b
3317 :
3318 : ! the k = kfrsnl layer properties are updated to combined
3319 : ! the fresnel (refractive) layer, always taken to be above
3320 : ! the present layer k (i.e. be the top interface):
3321 :
3322 5612695 : rintfc = c1 / (c1-Rf_dif_b*rdif_a(k))
3323 2073335 : tdir(k) = Tf_dir_a*tdir(k) + &
3324 2073335 : Tf_dir_a*rdir(k) * &
3325 5612695 : Rf_dif_b*rintfc*tdif_a(k)
3326 2073335 : rdir(k) = Rf_dir_a + &
3327 2073335 : Tf_dir_a*rdir(k) * &
3328 5612695 : rintfc*Tf_dif_b
3329 2073335 : rdif_a(k) = Rf_dif_a + &
3330 2073335 : Tf_dif_a*rdif_a(k) * &
3331 5612695 : rintfc*Tf_dif_b
3332 2073335 : rdif_b(k) = rdif_b(k) + &
3333 2073335 : tdif_b(k)*Rf_dif_b * &
3334 5612695 : rintfc*tdif_a(k)
3335 5612695 : tdif_a(k) = tdif_a(k)*rintfc*Tf_dif_a
3336 5612695 : tdif_b(k) = tdif_b(k)*rintfc*Tf_dif_b
3337 :
3338 : ! update trnlay to include fresnel transmission
3339 5612695 : trnlay(k) = Tf_dir_a*trnlay(k)
3340 :
3341 : endif ! k = kfrsnl
3342 :
3343 : endif ! trntdr(k) > trmin
3344 :
3345 : ! initialize current layer properties to zero; only if total
3346 : ! transmission to the top interface of the current layer exceeds the
3347 : ! minimum, will these values be computed below:
3348 : ! Calculate the solar beam transmission, total transmission, and
3349 : ! reflectivity for diffuse radiation from below at interface k,
3350 : ! the top of the current layer k:
3351 : !
3352 : ! layers interface
3353 : !
3354 : ! --------------------- k-1
3355 : ! k-1
3356 : ! --------------------- k
3357 : ! k
3358 : ! ---------------------
3359 : ! For k = klevp
3360 : ! note that we ignore refraction between sea ice and underlying ocean:
3361 : !
3362 : ! layers interface
3363 : !
3364 : ! --------------------- k-1
3365 : ! k-1
3366 : ! --------------------- k
3367 : ! \\\\\\\ ocean \\\\\\\
3368 :
3369 100257648 : trndir(k+1) = trndir(k)*trnlay(k)
3370 100257648 : refkm1 = c1/(c1 - rdndif(k)*rdif_a(k))
3371 100257648 : tdrrdir = trndir(k)*rdir(k)
3372 100257648 : tdndif = trntdr(k) - trndir(k)
3373 35581206 : trntdr(k+1) = trndir(k)*tdir(k) + &
3374 135838854 : (tdndif + tdrrdir*rdndif(k))*refkm1*tdif_a(k)
3375 35581206 : rdndif(k+1) = rdif_b(k) + &
3376 135838854 : (tdif_b(k)*rdndif(k)*refkm1*tdif_a(k))
3377 109914258 : trndif(k+1) = trndif(k)*refkm1*tdif_a(k)
3378 :
3379 : enddo ! k end main level loop
3380 :
3381 : ! compute reflectivity to direct and diffuse radiation for layers
3382 : ! below by adding succesive layers starting from the underlying
3383 : ! ocean and working upwards:
3384 : !
3385 : ! layers interface
3386 : !
3387 : ! --------------------- k
3388 : ! k
3389 : ! --------------------- k+1
3390 : ! k+1
3391 : ! ---------------------
3392 :
3393 9656610 : rupdir(klevp) = albodr
3394 9656610 : rupdif(klevp) = albodf
3395 :
3396 109914258 : do k=klev,0,-1
3397 : ! interface scattering
3398 100257648 : refkp1 = c1/( c1 - rdif_b(k)*rupdif(k+1))
3399 : ! dir from top layer plus exp tran ref from lower layer, interface
3400 : ! scattered and tran thru top layer from below, plus diff tran ref
3401 : ! from lower layer with interface scattering tran thru top from below
3402 35581206 : rupdir(k) = rdir(k) &
3403 71162412 : + ( trnlay(k) *rupdir(k+1) &
3404 135838854 : + (tdir(k)-trnlay(k))*rupdif(k+1))*refkp1*tdif_b(k)
3405 : ! dif from top layer from above, plus dif tran upwards reflected and
3406 : ! interface scattered which tran top from below
3407 109914258 : rupdif(k) = rdif_a(k) + tdif_a(k)*rupdif(k+1)*refkp1*tdif_b(k)
3408 : enddo ! k
3409 :
3410 9656610 : end subroutine solution_dEdd
3411 :
3412 : !=======================================================================
3413 : !
3414 : ! Set snow horizontal coverage, density and grain radius diagnostically
3415 : ! for the Delta-Eddington solar radiation method.
3416 : !
3417 : ! author: Bruce P. Briegleb, NCAR
3418 : ! 2013: E Hunke merged with NCAR version
3419 :
3420 3770075 : subroutine shortwave_dEdd_set_snow(nslyr, R_snw, &
3421 : dT_mlt, rsnw_mlt, &
3422 : aice, vsno, &
3423 : Tsfc, fs, &
3424 : hs0, hs, &
3425 3770075 : rhosnw, rsnw)
3426 :
3427 : integer (kind=int_kind), intent(in) :: &
3428 : nslyr ! number of snow layers
3429 :
3430 : real (kind=dbl_kind), intent(in) :: &
3431 : R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
3432 : dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
3433 : rsnw_mlt ! maximum melting snow grain radius (10^-6 m)
3434 :
3435 : real (kind=dbl_kind), intent(in) :: &
3436 : aice , & ! concentration of ice
3437 : vsno , & ! volume of snow
3438 : Tsfc , & ! surface temperature
3439 : hs0 ! snow depth for transition to bare sea ice (m)
3440 :
3441 : real (kind=dbl_kind), intent(inout) :: &
3442 : fs , & ! horizontal coverage of snow
3443 : hs ! snow depth
3444 :
3445 : real (kind=dbl_kind), dimension (:), intent(out) :: &
3446 : rhosnw , & ! density in snow layer (kg/m3)
3447 : rsnw ! grain radius in snow layer (micro-meters)
3448 :
3449 : ! local variables
3450 :
3451 : integer (kind=int_kind) :: &
3452 : ks ! snow vertical index
3453 :
3454 : real (kind=dbl_kind) :: &
3455 1275244 : fT , & ! piecewise linear function of surface temperature
3456 1275244 : dTs , & ! difference of Tsfc and Timelt
3457 1275244 : rsnw_nm ! actual used nonmelt snow grain radius (micro-meters)
3458 :
3459 : real (kind=dbl_kind), parameter :: &
3460 : ! units for the following are 1.e-6 m (micro-meters)
3461 : rsnw_fresh = 100._dbl_kind, & ! freshly-fallen snow grain radius
3462 : rsnw_nonmelt = 500._dbl_kind, & ! nonmelt snow grain radius
3463 : rsnw_sig = 250._dbl_kind ! assumed sigma for snow grain radius
3464 :
3465 : character(len=*),parameter :: subname='(shortwave_dEdd_set_snow)'
3466 :
3467 : !-----------------------------------------------------------------------
3468 :
3469 : ! set snow horizontal fraction
3470 3770075 : hs = vsno / aice
3471 :
3472 3770075 : if (hs >= hs_min) then
3473 2808328 : fs = c1
3474 2808328 : if (hs0 > puny) fs = min(hs/hs0, c1)
3475 : endif
3476 :
3477 : ! bare ice, temperature dependence
3478 3770075 : dTs = Timelt - Tsfc
3479 3770075 : fT = -min(dTs/dT_mlt-c1,c0)
3480 : ! tune nonmelt snow grain radius if desired: note that
3481 : ! the sign is negative so that if R_snw is 1, then the
3482 : ! snow grain radius is reduced and thus albedo increased.
3483 3770075 : rsnw_nm = rsnw_nonmelt - R_snw*rsnw_sig
3484 3770075 : rsnw_nm = max(rsnw_nm, rsnw_fresh)
3485 3770075 : rsnw_nm = min(rsnw_nm, rsnw_mlt)
3486 8449182 : do ks = 1, nslyr
3487 : ! snow density ccsm3 constant value
3488 4679107 : rhosnw(ks) = rhos
3489 : ! snow grain radius between rsnw_nonmelt and rsnw_mlt
3490 4679107 : rsnw(ks) = rsnw_nm + (rsnw_mlt-rsnw_nm)*fT
3491 4679107 : rsnw(ks) = max(rsnw(ks), rsnw_fresh)
3492 8449182 : rsnw(ks) = min(rsnw(ks), rsnw_mlt)
3493 : enddo ! ks
3494 :
3495 3770075 : end subroutine shortwave_dEdd_set_snow
3496 :
3497 : !=======================================================================
3498 : !
3499 : ! Set pond fraction and depth diagnostically for
3500 : ! the Delta-Eddington solar radiation method.
3501 : !
3502 : ! author: Bruce P. Briegleb, NCAR
3503 : ! 2013: E Hunke merged with NCAR version
3504 :
3505 72399 : subroutine shortwave_dEdd_set_pond(Tsfc, &
3506 : fs, fp, &
3507 : hp)
3508 :
3509 : real (kind=dbl_kind), intent(in) :: &
3510 : Tsfc , & ! surface temperature
3511 : fs ! horizontal coverage of snow
3512 :
3513 : real (kind=dbl_kind), intent(out) :: &
3514 : fp , & ! pond fractional coverage (0 to 1)
3515 : hp ! pond depth (m)
3516 :
3517 : ! local variables
3518 :
3519 : real (kind=dbl_kind) :: &
3520 26280 : fT , & ! piecewise linear function of surface temperature
3521 26280 : dTs ! difference of Tsfc and Timelt
3522 :
3523 : real (kind=dbl_kind), parameter :: &
3524 : dT_pnd = c1 ! change in temp for pond fraction and depth
3525 :
3526 : character(len=*),parameter :: subname='(shortwave_dEdd_set_pond)'
3527 :
3528 : !-----------------------------------------------------------------------
3529 :
3530 : ! bare ice, temperature dependence
3531 72399 : dTs = Timelt - Tsfc
3532 72399 : fT = -min(dTs/dT_pnd-c1,c0)
3533 : ! pond
3534 72399 : fp = 0.3_dbl_kind*fT*(c1-fs)
3535 72399 : hp = 0.3_dbl_kind*fT*(c1-fs)
3536 :
3537 72399 : end subroutine shortwave_dEdd_set_pond
3538 :
3539 : ! End Delta-Eddington shortwave method
3540 :
3541 : !=======================================================================
3542 : !
3543 : ! authors Nicole Jeffery, LANL
3544 :
3545 0 : subroutine compute_shortwave_trcr(nslyr, &
3546 0 : bgcN, zaero, &
3547 0 : trcrn_bgcsw, &
3548 0 : sw_grid, hin, &
3549 : hbri, &
3550 : nilyr, nblyr, &
3551 0 : i_grid, &
3552 : skl_bgc, z_tracers )
3553 :
3554 : integer (kind=int_kind), intent(in) :: &
3555 : nslyr ! number of snow layers
3556 :
3557 : integer (kind=int_kind), intent(in) :: &
3558 : nblyr , & ! number of bio layers
3559 : nilyr ! number of ice layers
3560 :
3561 : real (kind=dbl_kind), dimension (:), intent(in) :: &
3562 : bgcN , & ! Nit tracer
3563 : zaero ! zaero tracer
3564 :
3565 : real (kind=dbl_kind), dimension (:), intent(out):: &
3566 : trcrn_bgcsw ! ice on shortwave grid tracers
3567 :
3568 : real (kind=dbl_kind), dimension (:), intent(in) :: &
3569 : sw_grid , & !
3570 : i_grid ! CICE bio grid
3571 :
3572 : real(kind=dbl_kind), intent(in) :: &
3573 : hin , & ! CICE ice thickness
3574 : hbri ! brine height
3575 :
3576 : logical (kind=log_kind), intent(in) :: &
3577 : skl_bgc , & ! skeletal layer bgc
3578 : z_tracers ! zbgc
3579 :
3580 : ! local variables
3581 :
3582 : integer (kind=int_kind) :: k, n, nn
3583 :
3584 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
3585 0 : trtmp0, & ! temporary, remapped tracers
3586 0 : trtmp
3587 :
3588 : real (kind=dbl_kind), dimension (nilyr+1):: &
3589 0 : icegrid ! correct for large ice surface layers
3590 :
3591 : real (kind=dbl_kind):: &
3592 0 : top_conc ! 1% (min_bgc) of surface concentration
3593 : ! when hin > hbri: just used in sw calculation
3594 :
3595 : character(len=*),parameter :: subname='(compute_shortwave_trcr)'
3596 :
3597 : !-----------------------------------------------------------------
3598 : ! Compute aerosols and algal chlorophyll on shortwave grid
3599 : !-----------------------------------------------------------------
3600 :
3601 0 : trtmp0(:) = c0
3602 0 : trtmp(:) = c0
3603 0 : trcrn_bgcsw(:) = c0
3604 :
3605 0 : do k = 1,nilyr+1
3606 0 : icegrid(k) = sw_grid(k)
3607 : enddo
3608 0 : if (sw_grid(1)*hin*c2 > hi_ssl) then
3609 0 : icegrid(1) = hi_ssl/c2/hin
3610 : endif
3611 :
3612 0 : if (z_tracers) then
3613 0 : if (tr_bgc_N) then
3614 0 : if (size(bgcN) < n_algae*(nblyr+3)) then
3615 0 : call icepack_warnings_add(subname//' ERROR: size(bgcN) too small')
3616 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
3617 0 : return
3618 : endif
3619 :
3620 0 : do k = 1, nblyr+1
3621 0 : do n = 1, n_algae
3622 0 : trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) + &
3623 0 : R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1 + k-1)
3624 : enddo ! n
3625 : enddo ! k
3626 :
3627 0 : top_conc = trtmp0(nt_bgc_N(1))*min_bgc
3628 : call remap_zbgc (nilyr+1, &
3629 : nt_bgc_N(1), &
3630 0 : trtmp0(1:ntrcr ), &
3631 0 : trtmp (1:ntrcr+2), &
3632 : 1, nblyr+1, &
3633 : hin, hbri, &
3634 0 : icegrid(1:nilyr+1), &
3635 0 : i_grid(1:nblyr+1), top_conc )
3636 0 : if (icepack_warnings_aborted(subname)) return
3637 :
3638 0 : do k = 1, nilyr+1
3639 0 : trcrn_bgcsw(nlt_chl_sw+nslyr+k) = trtmp(nt_bgc_N(1) + k-1)
3640 : enddo ! k
3641 :
3642 0 : do n = 1, n_algae ! snow contribution
3643 0 : trcrn_bgcsw(nlt_chl_sw)= trcrn_bgcsw(nlt_chl_sw) &
3644 0 : + R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+1)
3645 : ! snow surface layer
3646 0 : trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) = &
3647 0 : trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) &
3648 0 : + R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+2)
3649 : ! only 1 snow layer in zaero
3650 : enddo ! n
3651 : endif ! tr_bgc_N
3652 :
3653 0 : if (tr_zaero) then
3654 0 : if (size(zaero) < n_zaero*(nblyr+3)) then
3655 0 : call icepack_warnings_add(subname//' ERROR: size(zaero) too small')
3656 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
3657 0 : return
3658 : endif
3659 :
3660 0 : do n = 1, n_zaero
3661 :
3662 0 : trtmp0(:) = c0
3663 0 : trtmp(:) = c0
3664 :
3665 0 : do k = 1, nblyr+1
3666 0 : trtmp0(nt_zaero(n) + k-1) = zaero(nt_zaero(n)-nt_zaero(1)+1+k-1)
3667 : enddo
3668 :
3669 0 : top_conc = trtmp0(nt_zaero(n))*min_bgc
3670 : call remap_zbgc (nilyr+1, &
3671 0 : nt_zaero(n), &
3672 0 : trtmp0(1:ntrcr ), &
3673 0 : trtmp (1:ntrcr+2), &
3674 : 1, nblyr+1, &
3675 : hin, hbri, &
3676 0 : icegrid(1:nilyr+1), &
3677 0 : i_grid(1:nblyr+1), top_conc )
3678 0 : if (icepack_warnings_aborted(subname)) return
3679 :
3680 0 : do k = 1,nilyr+1
3681 0 : trcrn_bgcsw(nlt_zaero_sw(n)+nslyr+k) = trtmp(nt_zaero(n) + k-1)
3682 : enddo
3683 0 : trcrn_bgcsw(nlt_zaero_sw(n))= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+1) !snow ssl
3684 0 : trcrn_bgcsw(nlt_zaero_sw(n)+1:nlt_zaero_sw(n)+nslyr)= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+2)
3685 : enddo ! n
3686 : endif ! tr_zaero
3687 0 : elseif (skl_bgc) then
3688 :
3689 0 : do nn = 1,n_algae
3690 0 : trcrn_bgcsw(nbtrcr_sw) = trcrn_bgcsw(nbtrcr_sw) &
3691 0 : + F_abs_chl(nn)*R_chl2N(nn) &
3692 0 : * bgcN(nt_bgc_N(nn)-nt_bgc_N(1)+1)*sk_l/hin &
3693 0 : * real(nilyr,kind=dbl_kind)
3694 : enddo
3695 :
3696 : endif
3697 : end subroutine compute_shortwave_trcr
3698 :
3699 : !=======================================================================
3700 : !autodocument_start icepack_prep_radiation
3701 : ! Scales radiation fields computed on the previous time step.
3702 : !
3703 : ! authors: Elizabeth Hunke, LANL
3704 :
3705 1370160 : subroutine icepack_prep_radiation (ncat, nilyr, nslyr, &
3706 1370160 : aice, aicen, &
3707 : swvdr, swvdf, &
3708 : swidr, swidf, &
3709 : alvdr_ai, alvdf_ai, &
3710 : alidr_ai, alidf_ai, &
3711 : scale_factor, &
3712 1370160 : fswsfcn, fswintn, &
3713 1370160 : fswthrun, fswpenln, &
3714 1370160 : Sswabsn, Iswabsn)
3715 :
3716 : integer (kind=int_kind), intent(in) :: &
3717 : ncat , & ! number of ice thickness categories
3718 : nilyr , & ! number of ice layers
3719 : nslyr ! number of snow layers
3720 :
3721 : real (kind=dbl_kind), intent(in) :: &
3722 : aice , & ! ice area fraction
3723 : swvdr , & ! sw down, visible, direct (W/m^2)
3724 : swvdf , & ! sw down, visible, diffuse (W/m^2)
3725 : swidr , & ! sw down, near IR, direct (W/m^2)
3726 : swidf , & ! sw down, near IR, diffuse (W/m^2)
3727 : ! grid-box-mean albedos aggregated over categories (if calc_Tsfc)
3728 : alvdr_ai , & ! visible, direct (fraction)
3729 : alidr_ai , & ! near-ir, direct (fraction)
3730 : alvdf_ai , & ! visible, diffuse (fraction)
3731 : alidf_ai ! near-ir, diffuse (fraction)
3732 :
3733 : real (kind=dbl_kind), dimension(:), intent(in) :: &
3734 : aicen ! ice area fraction in each category
3735 :
3736 : real (kind=dbl_kind), intent(inout) :: &
3737 : scale_factor ! shortwave scaling factor, ratio new:old
3738 :
3739 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
3740 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
3741 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
3742 : fswthrun ! SW through ice to ocean (W/m^2)
3743 :
3744 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
3745 : fswpenln , & ! visible SW entering ice layers (W m-2)
3746 : Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
3747 : Sswabsn ! SW radiation absorbed in snow layers (W m-2)
3748 :
3749 : !autodocument_end
3750 :
3751 : ! local variables
3752 :
3753 : integer (kind=int_kind) :: &
3754 : k , & ! vertical index
3755 : n ! thickness category index
3756 :
3757 474288 : real (kind=dbl_kind) :: netsw
3758 :
3759 : character(len=*),parameter :: subname='(icepack_prep_radiation)'
3760 :
3761 : !-----------------------------------------------------------------
3762 : ! Compute netsw scaling factor (new netsw / old netsw)
3763 : !-----------------------------------------------------------------
3764 :
3765 1370160 : if (aice > c0 .and. scale_factor > puny) then
3766 : netsw = swvdr*(c1 - alvdr_ai) &
3767 : + swvdf*(c1 - alvdf_ai) &
3768 : + swidr*(c1 - alidr_ai) &
3769 741605 : + swidf*(c1 - alidf_ai)
3770 741605 : scale_factor = netsw / scale_factor
3771 : else
3772 628555 : scale_factor = c1
3773 : endif
3774 :
3775 7834848 : do n = 1, ncat
3776 :
3777 7834848 : if (aicen(n) > puny) then
3778 :
3779 : !-----------------------------------------------------------------
3780 : ! Scale absorbed solar radiation for change in net shortwave
3781 : !-----------------------------------------------------------------
3782 :
3783 4360794 : fswsfcn(n) = scale_factor*fswsfcn (n)
3784 4360794 : fswintn(n) = scale_factor*fswintn (n)
3785 4360794 : fswthrun(n) = scale_factor*fswthrun(n)
3786 37427508 : do k = 1,nilyr+1
3787 37427508 : fswpenln(k,n) = scale_factor*fswpenln(k,n)
3788 : enddo !k
3789 9630500 : do k=1,nslyr
3790 9630500 : Sswabsn(k,n) = scale_factor*Sswabsn(k,n)
3791 : enddo
3792 33066714 : do k=1,nilyr
3793 33066714 : Iswabsn(k,n) = scale_factor*Iswabsn(k,n)
3794 : enddo
3795 :
3796 : endif
3797 : enddo ! ncat
3798 :
3799 1370160 : end subroutine icepack_prep_radiation
3800 :
3801 : !=======================================================================
3802 : !autodocument_start icepack_step_radiation
3803 : ! Computes radiation fields
3804 : !
3805 : ! authors: William H. Lipscomb, LANL
3806 : ! David Bailey, NCAR
3807 : ! Elizabeth C. Hunke, LANL
3808 :
3809 1027755 : subroutine icepack_step_radiation (dt, ncat, &
3810 : nblyr, &
3811 : nilyr, nslyr, &
3812 : dEdd_algae, &
3813 1027755 : swgrid, igrid, &
3814 1027755 : fbri, &
3815 1027755 : aicen, vicen, &
3816 1027755 : vsnon, Tsfcn, &
3817 1027755 : alvln, apndn, &
3818 1027755 : hpndn, ipndn, &
3819 1027755 : aeron, &
3820 1027755 : bgcNn, zaeron, &
3821 1027755 : trcrn_bgcsw, &
3822 : TLAT, TLON, &
3823 0 : calendar_type, &
3824 : days_per_year, &
3825 : nextsw_cday, &
3826 : yday, sec, &
3827 1027755 : kaer_tab, waer_tab, &
3828 1027755 : gaer_tab, &
3829 1027755 : kaer_bc_tab, &
3830 1027755 : waer_bc_tab, &
3831 1027755 : gaer_bc_tab, &
3832 1027755 : bcenh, &
3833 : modal_aero, &
3834 : swvdr, swvdf, &
3835 : swidr, swidf, &
3836 : coszen, fsnow, &
3837 2055510 : alvdrn, alvdfn, &
3838 2055510 : alidrn, alidfn, &
3839 1027755 : fswsfcn, fswintn, &
3840 1027755 : fswthrun, fswpenln, &
3841 1027755 : Sswabsn, Iswabsn, &
3842 1027755 : albicen, albsnon, &
3843 1027755 : albpndn, apeffn, &
3844 1027755 : snowfracn, &
3845 1027755 : dhsn, ffracn, &
3846 : l_print_point, &
3847 : initonly)
3848 :
3849 : integer (kind=int_kind), intent(in) :: &
3850 : ncat , & ! number of ice thickness categories
3851 : nilyr , & ! number of ice layers
3852 : nslyr , & ! number of snow layers
3853 : nblyr ! number of bgc layers
3854 :
3855 : real (kind=dbl_kind), intent(in) :: &
3856 : dt , & ! time step (s)
3857 : swvdr , & ! sw down, visible, direct (W/m^2)
3858 : swvdf , & ! sw down, visible, diffuse (W/m^2)
3859 : swidr , & ! sw down, near IR, direct (W/m^2)
3860 : swidf , & ! sw down, near IR, diffuse (W/m^2)
3861 : fsnow , & ! snowfall rate (kg/m^2 s)
3862 : TLAT, TLON ! latitude and longitude (radian)
3863 :
3864 : character (len=char_len), intent(in) :: &
3865 : calendar_type ! differentiates Gregorian from other calendars
3866 :
3867 : integer (kind=int_kind), intent(in) :: &
3868 : days_per_year, & ! number of days in one year
3869 : sec ! elapsed seconds into date
3870 :
3871 : real (kind=dbl_kind), intent(in) :: &
3872 : nextsw_cday , & ! julian day of next shortwave calculation
3873 : yday ! day of the year
3874 :
3875 : real (kind=dbl_kind), intent(inout) :: &
3876 : coszen ! cosine solar zenith angle, < 0 for sun below horizon
3877 :
3878 : real (kind=dbl_kind), dimension (:), intent(in) :: &
3879 : igrid ! biology vertical interface points
3880 :
3881 : real (kind=dbl_kind), dimension (:), intent(in) :: &
3882 : swgrid ! grid for ice tracers used in dEdd scheme
3883 :
3884 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
3885 : kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
3886 : waer_tab, & ! aerosol single scatter albedo (fraction)
3887 : gaer_tab ! aerosol asymmetry parameter (cos(theta))
3888 :
3889 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
3890 : kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
3891 : waer_bc_tab, & ! aerosol single scatter albedo (fraction)
3892 : gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
3893 :
3894 : real (kind=dbl_kind), dimension(:,:,:), intent(in) :: &
3895 : bcenh
3896 :
3897 : real (kind=dbl_kind), dimension(:), intent(in) :: &
3898 : aicen , & ! ice area fraction in each category
3899 : vicen , & ! ice volume in each category (m)
3900 : vsnon , & ! snow volume in each category (m)
3901 : Tsfcn , & ! surface temperature (deg C)
3902 : alvln , & ! level-ice area fraction
3903 : apndn , & ! pond area fraction
3904 : hpndn , & ! pond depth (m)
3905 : ipndn , & ! pond refrozen lid thickness (m)
3906 : fbri ! brine fraction
3907 :
3908 : real(kind=dbl_kind), dimension(:,:), intent(in) :: &
3909 : aeron , & ! aerosols (kg/m^3)
3910 : bgcNn , & ! bgc Nit tracers
3911 : zaeron ! bgcz aero tracers
3912 :
3913 : real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
3914 : trcrn_bgcsw ! zaerosols (kg/m^3) and chla (mg/m^3)
3915 :
3916 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
3917 : alvdrn , & ! visible, direct albedo (fraction)
3918 : alidrn , & ! near-ir, direct (fraction)
3919 : alvdfn , & ! visible, diffuse (fraction)
3920 : alidfn , & ! near-ir, diffuse (fraction)
3921 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
3922 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
3923 : fswthrun , & ! SW through ice to ocean (W/m^2)
3924 : snowfracn , & ! snow fraction on each category
3925 : dhsn , & ! depth difference for snow on sea ice and pond ice
3926 : ffracn , & ! fraction of fsurfn used to melt ipond
3927 : ! albedo components for history
3928 : albicen , & ! bare ice
3929 : albsnon , & ! snow
3930 : albpndn , & ! pond
3931 : apeffn ! effective pond area used for radiation calculation
3932 :
3933 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
3934 : fswpenln , & ! visible SW entering ice layers (W m-2)
3935 : Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
3936 : Sswabsn ! SW radiation absorbed in snow layers (W m-2)
3937 :
3938 : logical (kind=log_kind), intent(in) :: &
3939 : l_print_point, & ! flag for printing diagnostics
3940 : dEdd_algae , & ! .true. use prognostic chla in dEdd
3941 : modal_aero ! .true. use modal aerosol optical treatment
3942 :
3943 : logical (kind=log_kind), optional :: &
3944 : initonly ! flag to indicate init only, default is false
3945 :
3946 : !autodocument_end
3947 :
3948 : ! local variables
3949 :
3950 : integer (kind=int_kind) :: &
3951 : n ! thickness category index
3952 :
3953 : logical (kind=log_kind) :: &
3954 : linitonly ! local flag for initonly
3955 :
3956 : real(kind=dbl_kind) :: &
3957 355764 : hin, & ! Ice thickness (m)
3958 355764 : hbri ! brine thickness (m)
3959 :
3960 : character(len=*),parameter :: subname='(icepack_step_radiation)'
3961 :
3962 1027755 : hin = c0
3963 1027755 : hbri = c0
3964 1027755 : linitonly = .false.
3965 1027755 : if (present(initonly)) then
3966 135 : linitonly = initonly
3967 : endif
3968 :
3969 : ! Initialize
3970 5876910 : do n = 1, ncat
3971 4849155 : alvdrn (n) = c0
3972 4849155 : alidrn (n) = c0
3973 4849155 : alvdfn (n) = c0
3974 4849155 : alidfn (n) = c0
3975 4849155 : fswsfcn (n) = c0
3976 4849155 : fswintn (n) = c0
3977 5876910 : fswthrun(n) = c0
3978 : enddo ! ncat
3979 42063570 : fswpenln (:,:) = c0
3980 37214415 : Iswabsn (:,:) = c0
3981 12174165 : Sswabsn (:,:) = c0
3982 5876910 : trcrn_bgcsw(:,:) = c0
3983 :
3984 : ! Interpolate z-shortwave tracers to shortwave grid
3985 1027755 : if (dEdd_algae) then
3986 0 : do n = 1, ncat
3987 0 : if (aicen(n) .gt. puny) then
3988 0 : hin = vicen(n)/aicen(n)
3989 0 : hbri= fbri(n)*hin
3990 : call compute_shortwave_trcr(nslyr, &
3991 0 : bgcNn(:,n), &
3992 0 : zaeron(:,n), &
3993 0 : trcrn_bgcsw(:,n), &
3994 0 : swgrid, hin, &
3995 : hbri, &
3996 : nilyr, nblyr, &
3997 0 : igrid, &
3998 0 : skl_bgc, z_tracers )
3999 0 : if (icepack_warnings_aborted(subname)) return
4000 : endif
4001 : enddo
4002 : endif
4003 :
4004 1027755 : if (calc_Tsfc) then
4005 1027755 : if (trim(shortwave) == 'dEdd') then ! delta Eddington
4006 :
4007 : call run_dEdd(dt, ncat, &
4008 : dEdd_algae, &
4009 : nilyr, nslyr, &
4010 0 : aicen, vicen, &
4011 0 : vsnon, Tsfcn, &
4012 0 : alvln, apndn, &
4013 0 : hpndn, ipndn, &
4014 0 : aeron, kalg, &
4015 0 : trcrn_bgcsw, &
4016 : heat_capacity, &
4017 : TLAT, TLON, &
4018 : calendar_type,days_per_year, &
4019 : nextsw_cday, yday, &
4020 : sec, R_ice, &
4021 : R_pnd, R_snw, &
4022 : dT_mlt, rsnw_mlt, &
4023 : hs0, hs1, &
4024 : hp1, pndaspect, &
4025 0 : kaer_tab, waer_tab, &
4026 0 : gaer_tab, &
4027 0 : kaer_bc_tab, &
4028 0 : waer_bc_tab, &
4029 0 : gaer_bc_tab, &
4030 0 : bcenh, &
4031 : modal_aero, &
4032 : swvdr, swvdf, &
4033 : swidr, swidf, &
4034 : coszen, fsnow, &
4035 0 : alvdrn, alvdfn, &
4036 0 : alidrn, alidfn, &
4037 0 : fswsfcn, fswintn, &
4038 0 : fswthrun, fswpenln, &
4039 0 : Sswabsn, Iswabsn, &
4040 0 : albicen, albsnon, &
4041 0 : albpndn, apeffn, &
4042 0 : snowfracn, &
4043 0 : dhsn, ffracn, &
4044 : l_print_point, &
4045 882945 : linitonly)
4046 882945 : if (icepack_warnings_aborted(subname)) return
4047 :
4048 144810 : elseif (trim(shortwave) == 'ccsm3') then
4049 :
4050 0 : call shortwave_ccsm3(aicen, vicen, &
4051 0 : vsnon, &
4052 0 : Tsfcn, &
4053 : swvdr, swvdf, &
4054 : swidr, swidf, &
4055 : heat_capacity, &
4056 : albedo_type, &
4057 : albicev, albicei, &
4058 : albsnowv, albsnowi, &
4059 : ahmax, &
4060 0 : alvdrn, alidrn, &
4061 0 : alvdfn, alidfn, &
4062 0 : fswsfcn, fswintn, &
4063 0 : fswthrun, &
4064 0 : fswpenln, &
4065 0 : Iswabsn, &
4066 0 : Sswabsn, &
4067 0 : albicen, albsnon, &
4068 : coszen, ncat, &
4069 144810 : nilyr)
4070 144810 : if (icepack_warnings_aborted(subname)) return
4071 :
4072 : else
4073 :
4074 0 : call icepack_warnings_add(subname//' ERROR: shortwave '//trim(shortwave)//' unknown')
4075 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
4076 0 : return
4077 :
4078 : endif ! shortwave
4079 :
4080 : else ! .not. calc_Tsfc
4081 :
4082 : ! Calculate effective pond area for HadGEM
4083 :
4084 0 : if (tr_pond_topo) then
4085 0 : do n = 1, ncat
4086 0 : apeffn(n) = c0
4087 0 : if (aicen(n) > puny) then
4088 : ! Lid effective if thicker than hp1
4089 0 : if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
4090 0 : apeffn(n) = apndn(n)
4091 : else
4092 0 : apeffn(n) = c0
4093 : endif
4094 0 : if (apndn(n) < puny) apeffn(n) = c0
4095 : endif
4096 : enddo ! ncat
4097 :
4098 : endif ! tr_pond_topo
4099 :
4100 : ! Initialize for safety
4101 0 : do n = 1, ncat
4102 0 : alvdrn(n) = c0
4103 0 : alidrn(n) = c0
4104 0 : alvdfn(n) = c0
4105 0 : alidfn(n) = c0
4106 0 : fswsfcn(n) = c0
4107 0 : fswintn(n) = c0
4108 0 : fswthrun(n) = c0
4109 : enddo ! ncat
4110 0 : Iswabsn(:,:) = c0
4111 0 : Sswabsn(:,:) = c0
4112 :
4113 : endif ! calc_Tsfc
4114 :
4115 : end subroutine icepack_step_radiation
4116 :
4117 : ! Delta-Eddington solution expressions
4118 :
4119 : !=======================================================================
4120 :
4121 469322460 : real(kind=dbl_kind) function alpha(w,uu,gg,e)
4122 :
4123 : real(kind=dbl_kind), intent(in) :: w, uu, gg, e
4124 :
4125 469322460 : alpha = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
4126 :
4127 469322460 : end function alpha
4128 :
4129 : !=======================================================================
4130 :
4131 469322460 : real(kind=dbl_kind) function agamm(w,uu,gg,e)
4132 :
4133 : real(kind=dbl_kind), intent(in) :: w, uu, gg, e
4134 :
4135 469322460 : agamm = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
4136 :
4137 469322460 : end function agamm
4138 :
4139 : !=======================================================================
4140 :
4141 52146940 : real(kind=dbl_kind) function n(uu,et)
4142 :
4143 : real(kind=dbl_kind), intent(in) :: uu, et
4144 :
4145 52146940 : n = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
4146 :
4147 52146940 : end function n
4148 :
4149 : !=======================================================================
4150 :
4151 52146940 : real(kind=dbl_kind) function u(w,gg,e)
4152 :
4153 : real(kind=dbl_kind), intent(in) :: w, gg, e
4154 :
4155 52146940 : u = c1p5*(c1 - w*gg)/e
4156 :
4157 52146940 : end function u
4158 :
4159 : !=======================================================================
4160 :
4161 52146940 : real(kind=dbl_kind) function el(w,gg)
4162 :
4163 : real(kind=dbl_kind), intent(in) :: w, gg
4164 :
4165 52146940 : el = sqrt(c3*(c1-w)*(c1 - w*gg))
4166 :
4167 52146940 : end function el
4168 :
4169 : !=======================================================================
4170 :
4171 52146940 : real(kind=dbl_kind) function taus(w,f,t)
4172 :
4173 : real(kind=dbl_kind), intent(in) :: w, f, t
4174 :
4175 52146940 : taus = (c1 - w*f)*t
4176 :
4177 52146940 : end function taus
4178 :
4179 : !=======================================================================
4180 :
4181 52146940 : real(kind=dbl_kind) function omgs(w,f)
4182 :
4183 : real(kind=dbl_kind), intent(in) :: w, f
4184 :
4185 52146940 : omgs = (c1 - f)*w/(c1 - w*f)
4186 :
4187 52146940 : end function omgs
4188 :
4189 : !=======================================================================
4190 :
4191 52146940 : real(kind=dbl_kind) function asys(gg,f)
4192 :
4193 : real(kind=dbl_kind), intent(in) :: gg, f
4194 :
4195 52146940 : asys = (gg - f)/(c1 - f)
4196 :
4197 52146940 : end function asys
4198 :
4199 : !=======================================================================
4200 :
4201 : end module icepack_shortwave
4202 :
4203 : !=======================================================================
|
