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