Line data Source code
1 : !=======================================================================
2 : !
3 : ! snow redistribution and metamorphism
4 : !
5 : ! authors Elizabeth Hunke, LANL
6 : ! Nicole Jeffery, LANL
7 : !
8 : module icepack_snow
9 :
10 : use icepack_kinds
11 : use icepack_parameters, only: puny, p1, p5, c0, c1, c4, c10, c100, pi
12 : use icepack_parameters, only: rhos, rhow, rhoi, rhofresh, snwgrain
13 : use icepack_parameters, only: snwlvlfac, Tffresh, cp_ice, Lfresh
14 : use icepack_parameters, only: snwredist, rsnw_fall, rsnw_tmax, rhosnew
15 : use icepack_parameters, only: rhosmin, rhosmax, windmin, drhosdwind
16 : use icepack_parameters, only: isnw_T, isnw_Tgrd, isnw_rhos
17 : use icepack_parameters, only: snowage_rhos, snowage_Tgrd, snowage_T
18 : use icepack_parameters, only: snowage_tau, snowage_kappa, snowage_drdt0
19 : use icepack_parameters, only: snw_aging_table, use_smliq_pnd
20 : use icepack_tracers, only: ncat, nilyr, nslyr
21 :
22 : use icepack_therm_shared, only: icepack_ice_temperature
23 : use icepack_therm_shared, only: adjust_enthalpy
24 :
25 : use icepack_warnings, only: icepack_warnings_add, icepack_warnings_setabort
26 : use icepack_warnings, only: icepack_warnings_aborted
27 :
28 : implicit none
29 : private
30 :
31 : public :: icepack_step_snow, drain_snow, icepack_init_snow
32 :
33 : real (kind=dbl_kind), parameter, public :: &
34 : S_r = 0.033_dbl_kind, & ! irreducible saturation (Anderson 1976)
35 : S_wet= 4.22e5_dbl_kind ! wet metamorphism parameter (um^3/s)
36 : ! = 1.e18 * 4.22e-13 (Oleson 2010)
37 :
38 : real (kind=dbl_kind) :: &
39 : min_rhos, & ! snowtable axis data, assumes linear data
40 : del_rhos, &
41 : min_Tgrd, &
42 : del_Tgrd, &
43 : min_T , &
44 : del_T
45 :
46 : logical (kind=log_kind) :: &
47 : lin_rhos = .false., & ! flag to specify whether the snowtable dimensions are linear
48 : lin_Tgrd = .false., & ! this will allow quick lookup
49 : lin_T = .false.
50 :
51 : !=======================================================================
52 :
53 : contains
54 :
55 : !=======================================================================
56 : !autodocument_start icepack_init_snow
57 : ! Updates snow tracers
58 : !
59 : ! authors: Elizabeth C. Hunke, LANL
60 : ! Nicole Jeffery, LANL
61 :
62 12 : subroutine icepack_init_snow
63 :
64 : !autodocument_end
65 :
66 : ! local variables
67 :
68 : integer (kind=int_kind) :: n
69 :
70 : character (len=*),parameter :: subname='(icepack_init_snow)'
71 :
72 : !-----------------------------------------------------------------
73 : ! Snow metamorphism lookup table
74 : !-----------------------------------------------------------------
75 :
76 : ! if snw_aging_table = 'snicar'
77 : ! best-fit parameters are read from a table (netcdf)
78 : ! snowage_tau, snowage_kappa, snowage_drdt0
79 : ! 11 temperatures from 223.15 to 273.15 K by 5.0
80 : ! 31 temperature gradients from 0 to 300 K/m by 10
81 : ! 8 snow densities from 50 to 400 kg/m3 by 50
82 :
83 : ! if snw_aging_table = 'test'
84 : ! for testing Icepack without netcdf,
85 : ! use a subsampled, hard-coded table
86 : ! 5 temperatures from 243.15 by 5.0 K
87 : ! 5 temperature gradients from 0 to 40 K/m by 10
88 : ! 1 snow density at 300 kg/m3
89 :
90 : ! if snw_aging_table = 'file'
91 : ! all data has to be passed into icepack_parameters
92 :
93 12 : if (snwgrain) then
94 11 : if (trim(snw_aging_table) == 'snicar') then ! read netcdf file
95 0 : isnw_rhos = 8 ! maxiumum snow density index
96 0 : isnw_Tgrd = 31 ! maxiumum temperature gradient index
97 0 : isnw_T = 11 ! maxiumum temperature index
98 0 : min_rhos = 50.0_dbl_kind ! snowtable dimension data
99 0 : del_rhos = 50.0_dbl_kind
100 0 : lin_rhos = .true.
101 0 : min_Tgrd = 0.0_dbl_kind
102 0 : del_Tgrd = 10.0_dbl_kind
103 0 : lin_Tgrd = .true.
104 0 : min_T = 223.15_dbl_kind
105 0 : del_T = 5.0_dbl_kind
106 0 : lin_T = .true.
107 : ! check if these are already allocated/set, if so, make sure size is OK
108 0 : if (allocated(snowage_tau)) then
109 : if (size(snowage_tau,dim=1) /= isnw_rhos .or. &
110 0 : size(snowage_tau,dim=2) /= isnw_Tgrd .or. &
111 : size(snowage_tau,dim=3) /= isnw_T ) then
112 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
113 0 : call icepack_warnings_add(subname//'ERROR: snowage_tau size snw_aging_table=snicar')
114 0 : return
115 : endif
116 : else
117 0 : allocate (snowage_tau (isnw_rhos,isnw_Tgrd,isnw_T))
118 : endif
119 :
120 0 : if (allocated(snowage_kappa)) then
121 : if (size(snowage_kappa,dim=1) /= isnw_rhos .or. &
122 0 : size(snowage_kappa,dim=2) /= isnw_Tgrd .or. &
123 : size(snowage_kappa,dim=3) /= isnw_T ) then
124 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
125 0 : call icepack_warnings_add(subname//'ERROR: snowage_kappa size snw_aging_table=snicar')
126 0 : return
127 : endif
128 : else
129 0 : allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
130 : endif
131 :
132 0 : if (allocated(snowage_drdt0)) then
133 : if (size(snowage_drdt0,dim=1) /= isnw_rhos .or. &
134 0 : size(snowage_drdt0,dim=2) /= isnw_Tgrd .or. &
135 : size(snowage_drdt0,dim=3) /= isnw_T ) then
136 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
137 0 : call icepack_warnings_add(subname//'ERROR: snowage_drdt0 size snw_aging_table=snicar')
138 0 : return
139 : endif
140 : else
141 0 : allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
142 : endif
143 :
144 0 : if (allocated(snowage_rhos)) deallocate(snowage_rhos)
145 0 : if (allocated(snowage_Tgrd)) deallocate(snowage_Tgrd)
146 0 : if (allocated(snowage_T)) deallocate(snowage_T)
147 0 : allocate (snowage_rhos (isnw_rhos))
148 0 : allocate (snowage_Tgrd (isnw_Tgrd))
149 0 : allocate (snowage_T (isnw_T))
150 0 : do n = 1, isnw_rhos
151 0 : snowage_rhos(n) = min_rhos + real((n-1),dbl_kind)*del_rhos
152 : enddo
153 0 : do n = 1, isnw_Tgrd
154 0 : snowage_Tgrd(n) = min_Tgrd + real((n-1),dbl_kind)*del_Tgrd
155 : enddo
156 0 : do n = 1, isnw_T
157 0 : snowage_T(n) = min_T + real((n-1),dbl_kind)*del_T
158 : enddo
159 :
160 11 : elseif (trim(snw_aging_table) == 'file') then
161 0 : isnw_rhos = -1 ! maxiumum snow density index
162 0 : isnw_Tgrd = -1 ! maxiumum temperature gradient index
163 0 : isnw_T = -1 ! maxiumum temperature index
164 :
165 11 : elseif (trim(snw_aging_table) == 'test') then
166 11 : isnw_rhos = 1 ! maxiumum snow density index
167 11 : isnw_Tgrd = 5 ! maxiumum temperature gradient index
168 11 : isnw_T = 5 ! maxiumum temperature index
169 11 : min_rhos = 300.0_dbl_kind ! snowtable dimension data
170 11 : del_rhos = 50.0_dbl_kind
171 11 : lin_rhos = .true.
172 11 : min_Tgrd = 0.0_dbl_kind
173 11 : del_Tgrd = 10.0_dbl_kind
174 11 : lin_Tgrd = .true.
175 11 : min_T = 243.15_dbl_kind
176 11 : del_T = 5.0_dbl_kind
177 11 : lin_T = .true.
178 11 : if (allocated(snowage_tau)) deallocate(snowage_tau)
179 11 : if (allocated(snowage_kappa)) deallocate(snowage_kappa)
180 11 : if (allocated(snowage_drdt0)) deallocate(snowage_drdt0)
181 11 : if (allocated(snowage_rhos)) deallocate(snowage_rhos)
182 11 : if (allocated(snowage_Tgrd)) deallocate(snowage_Tgrd)
183 11 : if (allocated(snowage_T)) deallocate(snowage_T)
184 11 : allocate (snowage_tau (isnw_rhos,isnw_Tgrd,isnw_T))
185 11 : allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
186 11 : allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
187 11 : allocate (snowage_rhos (isnw_rhos))
188 11 : allocate (snowage_Tgrd (isnw_Tgrd))
189 11 : allocate (snowage_T (isnw_T))
190 22 : do n = 1, isnw_rhos
191 22 : snowage_rhos(n) = min_rhos + real((n-1),dbl_kind)*del_rhos
192 : enddo
193 66 : do n = 1, isnw_Tgrd
194 66 : snowage_Tgrd(n) = min_Tgrd + real((n-1),dbl_kind)*del_Tgrd
195 : enddo
196 66 : do n = 1, isnw_T
197 66 : snowage_T(n) = min_T + real((n-1),dbl_kind)*del_T
198 : enddo
199 :
200 : ! Subset of dry snow aging parameters
201 : snowage_tau = reshape((/ &
202 : 3.34947394_dbl_kind, 4.02124159_dbl_kind, 4.03328223_dbl_kind, &
203 : 3.02686921_dbl_kind, 2.14125851_dbl_kind, 3.97008737_dbl_kind, &
204 : 4.72725821_dbl_kind, 3.65313459_dbl_kind, 2.41198936_dbl_kind, &
205 : 2.53065623e-1_dbl_kind, 4.60286630_dbl_kind, 4.99721440_dbl_kind, &
206 : 3.29168191_dbl_kind, 2.66426779e-1_dbl_kind, 9.15830596e-5_dbl_kind, &
207 : 5.33186128_dbl_kind, 4.90833452_dbl_kind, 1.55269141_dbl_kind, &
208 : 1.31225526e-3_dbl_kind, 9.36078196e-4_dbl_kind, 6.25428631_dbl_kind, &
209 : 5.04394794_dbl_kind, 2.92857366e-3_dbl_kind, 9.01488751e-3_dbl_kind, &
210 : 1.19037046e-2_dbl_kind/), &
211 44 : (/isnw_rhos,isnw_Tgrd,isnw_T/))
212 :
213 : snowage_kappa = reshape((/ &
214 : 0.60824438, 0.56442972, 0.5527807, 0.64299537, 0.77672359, &
215 : 0.57105932, 0.52791041, 0.59868076, 0.7487191, 1.57946877, &
216 : 0.54236508, 0.52458285, 0.65520877, 1.52356017, 4.37789838, &
217 : 0.51449138, 0.54494334, 0.91628508, 3.28847035, 3.64418487, &
218 : 0.48538708, 0.55386601, 2.81247103, 2.72445522, 2.8230216/), &
219 660 : (/isnw_rhos,isnw_Tgrd,isnw_T/))
220 :
221 : snowage_drdt0 = reshape((/ &
222 : 1.26597871, 1.26602416, 1.26613263, 1.26620414, 1.26629424, &
223 : 1.92418877, 1.92430063, 1.92445964, 1.92451557, 1.92469806, &
224 : 2.79086547, 2.79147315, 2.79137562, 2.79150846, 2.79216439, &
225 : 3.85605846, 3.85668001, 3.85844559, 3.86073682, 3.863199, &
226 : 5.0861907, 5.08765668, 5.09200195, 5.09665276, 5.10079895/), &
227 660 : (/isnw_rhos,isnw_Tgrd,isnw_T/))
228 : else
229 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
230 0 : call icepack_warnings_add(subname//'ERROR: snw_aging_table value')
231 0 : return
232 : endif
233 : endif
234 :
235 : end subroutine icepack_init_snow
236 :
237 : !=======================================================================
238 : !autodocument_start icepack_step_snow
239 : ! Updates snow tracers
240 : !
241 : ! authors: Elizabeth C. Hunke, LANL
242 : ! Nicole Jeffery, LANL
243 :
244 339264 : subroutine icepack_step_snow(dt, &
245 : wind, aice, &
246 339264 : aicen, vicen, &
247 678528 : vsnon, Tsfc, &
248 339264 : zqin1, zSin1, &
249 339264 : zqsn, &
250 339264 : alvl, vlvl, &
251 339264 : smice, smliq, &
252 678528 : rsnw, rhos_cmpn, &
253 : fresh, fhocn, &
254 : fsloss, fsnow)
255 :
256 : real (kind=dbl_kind), intent(in) :: &
257 : dt , & ! time step
258 : wind , & ! wind speed (m/s)
259 : fsnow , & ! snowfall rate (kg m-2 s-1)
260 : aice ! ice area fraction
261 :
262 : real (kind=dbl_kind), dimension(:), intent(in) :: &
263 : aicen, & ! ice area fraction
264 : vicen, & ! ice volume (m)
265 : Tsfc , & ! surface temperature (C)
266 : zqin1, & ! ice upper layer enthalpy
267 : zSin1, & ! ice upper layer salinity
268 : alvl, & ! level ice area tracer
269 : vlvl ! level ice volume tracer
270 :
271 : real (kind=dbl_kind), intent(inout) :: &
272 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
273 : fhocn , & ! net heat flux to ocean (W/m^2)
274 : fsloss ! rate of snow loss to leads (kg/m^2/s)
275 :
276 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
277 : vsnon ! snow volume (m)
278 :
279 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
280 : zqsn , & ! snow enthalpy (J/m^3)
281 : smice , & ! tracer for mass of ice in snow (kg/m^3)
282 : smliq , & ! tracer for mass of liquid in snow (kg/m^3)
283 : rsnw , & ! snow grain radius (10^-6 m)
284 : rhos_cmpn ! effective snow density: compaction (kg/m^3)
285 :
286 : !autodocument_end
287 :
288 : ! local variables
289 :
290 : integer (kind=int_kind) :: k, n
291 :
292 : real (kind=dbl_kind), dimension(ncat) :: &
293 678528 : zTin1, & ! ice upper layer temperature (C)
294 678528 : hsn , & ! snow thickness (m)
295 339264 : hin ! ice thickness
296 :
297 : real (kind=dbl_kind) :: &
298 : vsno, & ! snow volume (m)
299 : tmp1, tmp2
300 :
301 : character (len=*),parameter :: subname='(icepack_step_snow)'
302 :
303 : !-----------------------------------------------------------------
304 : ! Initialize effective snow density (compaction) for new snow
305 : !-----------------------------------------------------------------
306 :
307 339264 : if (snwredist(1:3) == 'ITD') then
308 1825344 : do n = 1, ncat
309 9430944 : do k = 1, nslyr
310 9126720 : if (rhos_cmpn(k,n) < rhosmin) rhos_cmpn(k,n) = rhosnew
311 : enddo
312 : enddo
313 : endif
314 :
315 : !-----------------------------------------------------------------
316 : ! Redistribute snow based on wind
317 : !-----------------------------------------------------------------
318 :
319 339264 : vsno = c0
320 2035584 : do n = 1, ncat
321 2035584 : vsno = vsno + vsnon(n)
322 : enddo
323 339264 : tmp1 = rhos*vsno + fresh*dt
324 :
325 339264 : if (snwredist(1:3) == 'ITD' .and. aice > puny) then
326 : call snow_redist(dt, &
327 : wind, aicen(:), &
328 : vicen(:), vsnon(:), &
329 : zqsn(:,:), &
330 : alvl(:), vlvl(:), &
331 : fresh, fhocn, &
332 : fsloss, rhos_cmpn, &
333 217083 : fsnow)
334 217083 : if (icepack_warnings_aborted(subname)) return
335 : endif
336 :
337 339264 : vsno = c0
338 2035584 : do n = 1, ncat
339 2035584 : vsno = vsno + vsnon(n)
340 : enddo
341 339264 : tmp2 = rhos*vsno + fresh*dt
342 :
343 : ! check conservation
344 339264 : if (abs(tmp1-tmp2)>puny) then
345 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
346 0 : call icepack_warnings_add(subname//'ERROR: snow redistribution')
347 : endif
348 :
349 : !-----------------------------------------------------------------
350 : ! Adjust snow grain radius
351 : !-----------------------------------------------------------------
352 :
353 339264 : if (snwgrain) then
354 1825344 : do n = 1, ncat
355 1521120 : zTin1(n) = c0
356 1521120 : hsn (n) = c0
357 1521120 : hin (n) = c0
358 1825344 : if (aicen(n) > puny) then
359 885331 : zTin1(n) = icepack_ice_temperature(zqin1(n), zSin1(n))
360 885331 : hsn(n) = vsnon(n)/aicen(n)
361 885331 : hin(n) = vicen(n)/aicen(n)
362 : endif
363 : enddo
364 :
365 : call update_snow_radius (dt, &
366 : rsnw, hin, &
367 : Tsfc, zTin1, &
368 : hsn, zqsn, &
369 304224 : smice, smliq)
370 304224 : if (icepack_warnings_aborted(subname)) return
371 : endif
372 :
373 : end subroutine icepack_step_snow
374 :
375 : !=======================================================================
376 :
377 : ! Snow redistribution by wind, based on O. Lecomte Ph.D. (2014).
378 : ! The original formulation:
379 : ! Snow in suspension depends on wind speed, density and the standard
380 : ! deviation of the ice thickness distribution. Snow is redistributed
381 : ! among ice categories proportionally to the category areas.
382 : !
383 : ! Namelist option snwredist = 'ITDrdg' modifies the approach to use
384 : ! the level and ridged ice tracers:
385 : ! As above, but use the standard deviation of the level and ridged
386 : ! ice thickness distribution for snow in suspension, and redistribute
387 : ! based on ridged ice area.
388 : !
389 : ! convention:
390 : ! volume, mass and energy include factor of ain
391 : ! thickness does not
392 :
393 217083 : subroutine snow_redist(dt, wind, ain, vin, vsn, zqsn, &
394 217083 : alvl, vlvl, fresh, fhocn, fsloss, rhos_cmpn, fsnow)
395 :
396 : real (kind=dbl_kind), intent(in) :: &
397 : dt , & ! time step (s)
398 : wind , & ! wind speed (m/s)
399 : fsnow ! snowfall rate (kg m-2 s-1)
400 :
401 : real (kind=dbl_kind), dimension(:), intent(in) :: &
402 : ain , & ! ice area fraction
403 : vin , & ! ice volume (m)
404 : alvl , & ! level ice area tracer
405 : vlvl ! level ice volume tracer
406 :
407 : real (kind=dbl_kind), intent(inout) :: &
408 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
409 : fhocn , & ! net heat flux to ocean (W/m^2)
410 : fsloss ! rate of snow loss to leads (kg/m^2/s)
411 :
412 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
413 : vsn ! snow volume (m)
414 :
415 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
416 : zqsn , & ! snow enthalpy (J/m^3)
417 : rhos_cmpn ! effective snow density: compaction (kg/m^3)
418 :
419 : ! local variables
420 :
421 : integer (kind=int_kind) :: &
422 : n , & ! category index
423 : k ! layer index
424 :
425 : integer (kind=int_kind), dimension(ncat) :: &
426 434166 : klyr ! layer index
427 :
428 : real (kind=dbl_kind), parameter :: &
429 : refsd = c1 , & ! standard deviation reference
430 : gamma = 1.e-5_dbl_kind ! tuning coefficient
431 :
432 : real (kind=dbl_kind) :: &
433 : Vseas , & ! critical seasonal wind speed (m/s)
434 : ITDsd , & ! standard deviation of ITD
435 : flost , & ! fraction of snow lost in leads
436 : alost , & ! effective lead area for snow lost in leads
437 : suma , & ! sum of ice area over categories
438 : sumv , & ! sum of ice volume over categories (m)
439 : summ , & ! sum of snow mass over categories (kg/m^2)
440 : sumq , & ! sum of snow enthalpy over categories (kg/m^2)
441 : msusp , & ! potential mass of snow in suspension (kg/m^2)
442 : msnw_susp , & ! mass of snow in suspension (kg/m^2)
443 : esnw_susp , & ! energy of snow in suspension (J/m^2)
444 : asnw_lvl , & ! mass of snow redeposited on level ice (kg/m^2)
445 : e_redeptmp, & ! redeposited energy (J/m^2)
446 : dhsn , & ! change in snow depth (m)
447 : dmp , & ! mass difference in previous layer (kg/m^2)
448 : hslyr , & ! snow layer thickness (m)
449 : hslab , & ! new snow thickness (m)
450 : drhos , & ! change in snow density due to compaction (kg/m^3)
451 : mlost , & ! mass of suspended snow lost in leads (kg/m^2)
452 : elost , & ! energy of suspended snow lost in leads (J/m^2)
453 : de , & ! change in energy (J/m^2)
454 : al, ar , & ! areas of level and ridged ice
455 : hlvl, hrdg, & ! thicknesses of level and ridged ice
456 : tmp1, tmp2, & ! temporary values
457 : tmp3, tmp4, & ! temporary values
458 : tmp5 , & ! temporary values
459 : work ! temporary value
460 :
461 : real (kind=dbl_kind), dimension(ncat) :: &
462 434166 : sfac , & ! temporary for snwlvlfac
463 217083 : ardg , & ! ridged ice area tracer
464 434166 : m_erosion , & ! eroded mass (kg/m^2)
465 434166 : e_erosion , & ! eroded energy (J/m^2)
466 434166 : m_redep , & ! redeposited mass (kg/m^2)
467 434166 : e_redep , & ! redeposited energy (J/m^2)
468 434166 : vsn_init , & ! initial volume (m)
469 434166 : esn_init , & ! initial energy (J/m^2)
470 434166 : esn_final , & ! final energy (J/m^2)
471 434166 : atmp , & ! temporary variable for ain, for debugging convenience
472 434166 : hin , & ! ice thickness (m)
473 434166 : hsn , & ! snow depth (m)
474 434166 : hsn_new ! new snow depth (m)
475 :
476 : real (kind=dbl_kind), dimension (nslyr) :: &
477 434166 : dzs ! snow layer thickness after redistribution (m)
478 :
479 : real (kind=dbl_kind), dimension (nslyr+1) :: &
480 434166 : zs1 , & ! depth of snow layer boundaries (m)
481 217083 : zs2 ! adjusted depths, with equal hslyr (m)
482 :
483 : character (len=*),parameter :: subname='(snow_redist)'
484 :
485 : !-----------------------------------------------------------------
486 : ! Conservation checks
487 : !-----------------------------------------------------------------
488 :
489 217083 : tmp1 = c0
490 217083 : tmp3 = c0
491 1302498 : do n = 1, ncat
492 : ! mass conservation check
493 1085415 : tmp1 = tmp1 + vsn(n)
494 1085415 : vsn_init(n) = vsn(n)
495 1085415 : esn_init(n) = c0
496 : ! energy conservation check
497 6729573 : do k = 1, nslyr
498 5427075 : tmp3 = tmp3 + vsn(n)*zqsn(k,n)/nslyr
499 6512490 : esn_init(n) = esn_init(n) + vsn(n)*zqsn(k,n)/nslyr
500 : enddo
501 : enddo
502 :
503 : !-----------------------------------------------------------------
504 : ! category thickness and sums
505 : !-----------------------------------------------------------------
506 :
507 1302498 : hin(:) = c0
508 1302498 : hsn(:) = c0
509 217083 : suma = c0
510 217083 : sumv = c0
511 1302498 : do n = 1, ncat
512 1085415 : atmp(n) = ain(n)
513 1085415 : if (atmp(n) > puny) then
514 887689 : hin(n) = vin(n)/atmp(n)
515 887689 : hsn(n) = vsn(n)/atmp(n)
516 : endif
517 1085415 : hsn_new(n) = hsn(n)
518 1085415 : suma = suma + atmp(n)
519 1085415 : sumv = sumv + vin(n)
520 : ! maintain positive definite enthalpy
521 6729573 : do k = 1, nslyr
522 6512490 : zqsn(k,n) = min(zqsn(k,n) + Lfresh*rhos, c0)
523 : enddo
524 : enddo ! ncat
525 :
526 : !-----------------------------------------------------------------
527 : ! standard deviation of ice thickness distribution
528 : !-----------------------------------------------------------------
529 :
530 217083 : work = c0
531 217083 : asnw_lvl = c0
532 217083 : if (trim(snwredist) == 'ITDrdg') then ! use level and ridged ice
533 1302498 : do n = 1, ncat
534 1085415 : ardg(n) = c1 - alvl(n) ! ridged ice tracer
535 1085415 : al = alvl(n) * atmp(n) ! level
536 1085415 : ar = ardg(n) * atmp(n) ! ridged
537 1085415 : hlvl = c0
538 1085415 : hrdg = c0
539 1085415 : if (al > puny) hlvl = vin(n)*vlvl(n)/al
540 1085415 : if (ar > puny) hrdg = vin(n)*(c1-vlvl(n))/ar
541 1085415 : work = work + al*(hlvl - sumv)**2 + ar*(hrdg - sumv)**2
542 :
543 : ! for redeposition of snow on level ice
544 1085415 : sfac(n) = snwlvlfac
545 1085415 : if (ardg(n) > c0) sfac(n) = min(snwlvlfac, alvl(n)/ardg(n))
546 1302498 : asnw_lvl = asnw_lvl + al - sfac(n)*ar
547 : enddo
548 217083 : asnw_lvl = asnw_lvl/suma
549 : ! else ! snwredist = 'ITDsd' ! use standard ITD
550 : ! do n = 1, ncat
551 : ! work = work + atmp(n)*(hin(n) - sumv)**2
552 : ! enddo
553 : endif
554 217083 : ITDsd = sqrt(work)
555 :
556 : !-----------------------------------------------------------------
557 : ! fraction of suspended snow lost in leads
558 : !-----------------------------------------------------------------
559 :
560 217083 : flost = (c1 - suma) * exp(-ITDsd/refsd)
561 : ! flost = c0 ! echmod for testing
562 217083 : alost = c1 - suma * (c1-flost)
563 :
564 : !-----------------------------------------------------------------
565 : ! suspended snow
566 : !-----------------------------------------------------------------
567 :
568 217083 : msusp = c0
569 1302498 : do n = 1, ncat
570 : ! critical seasonal wind speed needed to compact snow to density rhos
571 1085415 : Vseas = (rhos_cmpn(1,n) - 44.6_dbl_kind)/174.0_dbl_kind ! use top layer
572 1085415 : Vseas = max(Vseas, c0)
573 : ! maximum mass per unit area of snow in suspension (kg/m^2)
574 1085415 : if (ITDsd > puny) &
575 : msusp = msusp + atmp(n)*gamma*dt*max(wind-Vseas,c0) &
576 1302378 : * (rhosmax-rhos_cmpn(1,n))/(rhosmax*ITDsd)
577 : enddo
578 :
579 : !-----------------------------------------------------------------
580 : ! erosion
581 : !-----------------------------------------------------------------
582 :
583 217083 : msnw_susp = c0
584 217083 : esnw_susp = c0
585 1302498 : klyr(:) = 1
586 1302498 : do n = 1, ncat
587 1085415 : m_erosion(n) = c0 ! mass
588 1085415 : e_erosion(n) = c0 ! energy
589 1302498 : if (atmp(n) > puny) then
590 887689 : m_erosion(n) = min(msusp, rhos*vsn(n))
591 887689 : if (m_erosion(n) > puny) then
592 828533 : summ = c0
593 828533 : dmp = m_erosion(n)
594 4971198 : do k = 1, nslyr
595 4971198 : if (dmp > c0) then
596 3059450 : dhsn = min(hsn(n)/nslyr, dmp/(rhos*atmp(n)))
597 3059450 : msnw_susp = msnw_susp + dhsn*rhos*atmp(n) ! total mass in suspension
598 3059450 : hsn_new(n) = hsn_new(n) - dhsn
599 3059450 : e_erosion(n) = e_erosion(n) + dhsn*zqsn(k,n)*atmp(n)
600 3059450 : klyr(n) = k ! number of affected layers
601 3059450 : summ = summ + rhos*vsn(n)/nslyr ! mass, partial sum
602 3059450 : dmp = max(m_erosion(n) - summ, c0)
603 : endif ! dmp
604 : enddo
605 828533 : esnw_susp = esnw_susp + e_erosion(n) ! total energy in suspension
606 : endif
607 : endif
608 : enddo
609 :
610 : !-----------------------------------------------------------------
611 : ! redeposition
612 : !-----------------------------------------------------------------
613 :
614 1302498 : do n = 1, ncat
615 1085415 : if (trim(snwredist) == 'ITDrdg') then ! use level and ridged ice
616 1085415 : work = atmp(n)*(c1-flost)*(ardg(n)*(c1+sfac(n)) + asnw_lvl)
617 : else ! use standard ITD
618 0 : work = atmp(n)*(c1-flost)
619 : endif
620 1085415 : m_redep(n) = msnw_susp*work ! mass
621 1085415 : e_redep(n) = c0
622 1085415 : e_redeptmp = esnw_susp*work ! energy
623 :
624 : ! change in snow depth
625 1085415 : dhsn = c0
626 1085415 : if (atmp(n) > puny) then
627 887689 : dhsn = m_redep(n) / (rhos*atmp(n))
628 :
629 887689 : if (abs(dhsn) > c0) then
630 :
631 832612 : e_redep(n) = e_redeptmp
632 832612 : vsn(n) = (hsn_new(n)+dhsn)*atmp(n)
633 :
634 : ! change in snow energy
635 832612 : de = e_redeptmp / klyr(n)
636 : ! spread among affected layers
637 832612 : sumq = c0
638 3896141 : do k = 1, klyr(n)
639 : zqsn(k,n) = (atmp(n)*hsn_new(n)*zqsn(k,n) + de) &
640 3063529 : / (vsn(n)) ! factor of nslyr cancels out
641 :
642 3896141 : if (zqsn(k,n) > c0) then
643 7 : sumq = sumq + zqsn(k,n)
644 7 : zqsn(k,n) = c0
645 : endif
646 :
647 : enddo ! klyr
648 832612 : zqsn(klyr(n),n) = min(zqsn(klyr(n),n) + sumq, c0) ! may lose energy here
649 :
650 : !-----------------------------------------------------------------
651 : ! Conserving energy, compute the enthalpy of the new equal layers
652 : !-----------------------------------------------------------------
653 :
654 832612 : if (nslyr > 1) then
655 :
656 4995672 : dzs(:) = hsn(n) / real(nslyr,kind=dbl_kind) ! old layer thickness
657 3896141 : do k = 1, klyr(n)
658 3896141 : dzs(k) = dzs(k) + dhsn / klyr(n) ! old layer thickness (updated)
659 : enddo
660 832612 : hsn_new(n) = hsn_new(n) + dhsn
661 832612 : hslyr = hsn_new(n) / real(nslyr,kind=dbl_kind) ! new layer thickness
662 :
663 832612 : zs1(1) = c0
664 832612 : zs1(1+nslyr) = hsn_new(n)
665 :
666 832612 : zs2(1) = c0
667 832612 : zs2(1+nslyr) = hsn_new(n)
668 :
669 4163060 : do k = 1, nslyr-1
670 3330448 : zs1(k+1) = zs1(k) + dzs(k) ! old layer depths (unequal thickness)
671 4163060 : zs2(k+1) = zs2(k) + hslyr ! new layer depths (equal thickness)
672 : enddo
673 :
674 : call adjust_enthalpy (nslyr, &
675 : zs1(:), zs2(:), &
676 : hslyr , hsn_new(n), &
677 832612 : zqsn(:,n))
678 832612 : if (icepack_warnings_aborted(subname)) return
679 : else
680 0 : hsn_new(n) = hsn_new(n) + dhsn
681 : endif ! nslyr > 1
682 : endif ! |dhsn| > puny
683 : endif ! ain > puny
684 :
685 : ! maintain positive definite enthalpy
686 6729573 : do k = 1, nslyr
687 6512490 : zqsn(k,n) = zqsn(k,n) - Lfresh*rhos
688 : enddo
689 : enddo ! ncat
690 :
691 : !-----------------------------------------------------------------
692 : ! mass of suspended snow lost in leads
693 : !-----------------------------------------------------------------
694 217083 : mlost = msnw_susp*alost
695 217083 : fsloss = fsloss + mlost / dt
696 :
697 : !-----------------------------------------------------------------
698 : ! mass conservation check
699 : !-----------------------------------------------------------------
700 :
701 217083 : tmp2 = c0
702 1302498 : do n = 1, ncat
703 1302498 : tmp2 = tmp2 + vsn(n)
704 : enddo
705 :
706 217083 : if (tmp2 > tmp1) then ! correct roundoff error
707 24 : vsn(:) = vsn(:) * tmp1/tmp2
708 4 : tmp2 = c0
709 24 : do n = 1, ncat
710 24 : tmp2 = tmp2 + vsn(n)
711 : enddo
712 : endif
713 :
714 217083 : if (tmp2 < tmp1) fresh = fresh + rhos*(tmp1-tmp2)/dt
715 :
716 217083 : tmp2 = tmp2 + (mlost/rhos)
717 :
718 217083 : if (abs(tmp1-tmp2) > puny) then
719 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
720 0 : call icepack_warnings_add(subname//'ERROR: snow redistribution mass conservation error')
721 : ! write(warning,*)'mass conservation error in snow_redist', tmp1, tmp2
722 : ! write(warning,*)'klyr',klyr
723 : ! write(warning,*)'ain',atmp(:)
724 : ! write(warning,*)'vsn final',vsn(:)
725 : ! write(warning,*)'vsn init',vsn_init(:)
726 : ! write(warning,*)'rhos*vsn init',rhos*vsn_init(:)
727 : ! write(warning,*)'m_erosion',m_erosion(:)
728 : ! write(warning,*)'m_redep',m_redep(:)
729 : ! write(warning,*)'mlost',mlost
730 : ! write(warning,*)'v_erosion',m_erosion(:)/rhos
731 : ! write(warning,*)'v_redep',m_redep(:)/rhos
732 : ! write(warning,*)'v lost',mlost/rhos
733 : ! write(warning,*)'hsn',hsn(:)
734 : ! write(warning,*)'hsn_new',hsn_new(:)
735 : ! write(warning,*)'vsn_new',hsn_new(:)*atmp(:)
736 : ! write(warning,*)'lost',suma,flost,alost,msnw_susp
737 : endif
738 :
739 : !-----------------------------------------------------------------
740 : ! energy conservation check
741 : !-----------------------------------------------------------------
742 :
743 217083 : tmp4 = c0
744 217083 : tmp5 = c0
745 1302498 : esn_final(:) = c0
746 1302498 : do n = 1, ncat
747 6512490 : do k = 1, nslyr
748 5427075 : tmp4 = tmp4 + vsn(n)*zqsn(k,n)/nslyr
749 6512490 : esn_final(n) = esn_final(n) + vsn(n)*zqsn(k,n)/nslyr
750 : enddo
751 1302498 : tmp5 = tmp5 - e_erosion(n) + e_redep(n)
752 : enddo
753 217083 : tmp5 = tmp5 + esnw_susp*alost
754 :
755 : !-----------------------------------------------------------------
756 : ! energy of suspended snow lost in leads
757 : !-----------------------------------------------------------------
758 217083 : elost = tmp3 - tmp4
759 217083 : fhocn = fhocn + elost / dt
760 :
761 217083 : if (abs(tmp5) > nslyr*Lfresh*puny) then
762 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
763 0 : call icepack_warnings_add(subname//'ERROR: snow redistribution energy conservation error')
764 : ! write(warning,*)'energy conservation error in snow_redist', tmp3, tmp4, tmp5
765 : ! write(warning,*)'klyr',klyr
766 : ! write(warning,*)'ain',atmp(:)
767 : ! write(warning,*)'vsn final',vsn(:)
768 : ! write(warning,*)'vsn init',vsn_init(:)
769 : ! write(warning,*)'rhos*vsn init',rhos*vsn_init(:)
770 : ! write(warning,*)'m_erosion',m_erosion(:)
771 : ! write(warning,*)'m_redep',m_redep(:)
772 : ! write(warning,*)'mlost',mlost
773 : ! write(warning,*)'v_erosion',m_erosion(:)/rhos
774 : ! write(warning,*)'v_redep',m_redep(:)/rhos
775 : ! write(warning,*)'v lost',mlost/rhos
776 : ! write(warning,*)'hsn',hsn(:)
777 : ! write(warning,*)'hsn_new',hsn_new(:)
778 : ! write(warning,*)'vsn_new',hsn_new(:)*atmp(:)
779 : ! write(warning,*)'lost',suma,flost,alost,msnw_susp
780 : ! write(warning,*)'tmp3(1)', (vsn(1)*zqsn(k,1)/nslyr,k=1,nslyr)
781 : ! write(warning,*)'esn init',esn_init(:)
782 : ! write(warning,*)'esn final',esn_final(:)
783 : ! write(warning,*)'e_erosion',e_erosion(:)
784 : ! write(warning,*)'e_redep',e_redep(:)
785 : ! write(warning,*)'elost',elost,esnw_susp*alost,Lfresh*mlost
786 : ! write(warning,*)'esnw_susp',esnw_susp
787 : endif
788 :
789 : !-----------------------------------------------------------------
790 : ! wind compaction
791 : !-----------------------------------------------------------------
792 :
793 1302498 : do n = 1, ncat
794 1302498 : if (vsn(n) > puny) then
795 : ! compact freshly fallen or redistributed snow
796 863880 : drhos = drhosdwind * max(wind - windmin, c0)
797 863880 : hslab = c0
798 863880 : if (fsnow > c0) &
799 833455 : hslab = max(min(fsnow*dt/(rhos+drhos), hsn_new(n)-hsn(n)), c0)
800 863880 : hslyr = hsn_new(n) / real(nslyr,kind=dbl_kind)
801 5183280 : do k = 1, nslyr
802 4319400 : work = hslab - hslyr * real(k-1,kind=dbl_kind)
803 4319400 : work = max(c0, min(hslyr, work))
804 4319400 : rhos_cmpn(k,n) = rhos_cmpn(k,n) + drhos*work/hslyr
805 5183280 : rhos_cmpn(k,n) = min(rhos_cmpn(k,n), rhosmax)
806 : enddo
807 : endif
808 : enddo
809 :
810 : end subroutine snow_redist
811 :
812 : !=======================================================================
813 :
814 : ! Snow grain metamorphism
815 :
816 304224 : subroutine update_snow_radius (dt, rsnw, hin, &
817 304224 : Tsfc, zTin, hsn, zqsn, smice, smliq)
818 :
819 : real (kind=dbl_kind), intent(in) :: &
820 : dt ! time step
821 :
822 : real (kind=dbl_kind), dimension(ncat), intent(in) :: &
823 : zTin , & ! surface ice temperature (oC)
824 : Tsfc , & ! surface temperature (oC)
825 : hin , & ! ice thickness (m)
826 : hsn ! snow thickness (m)
827 :
828 : real (kind=dbl_kind), dimension(nslyr,ncat), intent(in) :: &
829 : zqsn ! enthalpy of snow (J m-3)
830 :
831 : real (kind=dbl_kind), dimension(nslyr,ncat), intent(inout) :: &
832 : rsnw ! snow grain radius
833 :
834 : real (kind=dbl_kind), dimension(nslyr,ncat), intent(inout) :: &
835 : smice , & ! tracer for mass of ice in snow (kg/m^3)
836 : smliq ! tracer for mass of liquid in snow (kg/m^3)
837 :
838 : ! local temporary variables
839 :
840 : integer (kind=int_kind) :: k, n
841 :
842 : real (kind=dbl_kind), dimension(nslyr) :: &
843 608448 : drsnw_wet, & ! wet metamorphism (10^-6 m)
844 304224 : drsnw_dry ! dry (temperature gradient) metamorphism (10^-6 m)
845 :
846 : character (len=*),parameter :: subname='(update_snow_radius)'
847 :
848 1825344 : do n = 1, ncat
849 :
850 1825344 : if (hsn(n) > puny .and. hin(n) > puny) then
851 :
852 5183388 : drsnw_dry(:) = c0
853 5183388 : drsnw_wet(:) = c0
854 :
855 : !-----------------------------------------------------------------
856 : ! dry metamorphism
857 : !-----------------------------------------------------------------
858 : call snow_dry_metamorph (dt, rsnw(:,n), &
859 : drsnw_dry, zqsn(:,n), Tsfc(n), &
860 863898 : zTin(n), hsn(n), hin(n))
861 863898 : if (icepack_warnings_aborted(subname)) return
862 :
863 : !-----------------------------------------------------------------
864 : ! wet metamorphism
865 : !-----------------------------------------------------------------
866 5183388 : do k = 1,nslyr
867 : call snow_wet_metamorph (dt, drsnw_wet(k), rsnw(k,n), &
868 4319490 : smice(k,n), smliq(k,n))
869 4319490 : if (icepack_warnings_aborted(subname)) return
870 5183388 : rsnw(k,n) = min(rsnw_tmax, rsnw(k,n) + drsnw_dry(k) + drsnw_wet(k))
871 : enddo
872 :
873 : else ! hsn or hin < puny
874 3943332 : do k = 1,nslyr
875 : ! rsnw_fall < rsnw < rsnw_tmax
876 3286110 : rsnw (k,n) = max(rsnw_fall, min(rsnw_tmax, rsnw(k,n)))
877 3286110 : smice(k,n) = rhos
878 3943332 : smliq(k,n) = c0
879 : enddo
880 : endif ! hsn, hin
881 : enddo
882 :
883 : end subroutine update_snow_radius
884 :
885 : !=======================================================================
886 :
887 : ! Snow grain metamorphism
888 :
889 863898 : subroutine snow_dry_metamorph (dt, rsnw, drsnw_dry, zqsn, &
890 : Tsfc, zTin1, hsn, hin)
891 :
892 : ! Vapor redistribution: Method is to retrieve 3 best-fit parameters that
893 : ! depend on snow temperature, temperature gradient, and density,
894 : ! that are derived from the microphysical model described in:
895 : ! Flanner and Zender (2006), Linking snowpack microphysics and albedo
896 : ! evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834.
897 : ! The parametric equation has the form:
898 : ! dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), where:
899 : ! r is the effective radius,
900 : ! tau and kappa are best-fit parameters,
901 : ! drdt_0 is the initial rate of change of effective radius, and
902 : ! dr_fresh is the difference between the current and fresh snow states
903 : ! (r_current - r_fresh).
904 :
905 : real (kind=dbl_kind), intent(in) :: &
906 : dt ! time step (s)
907 :
908 : real (kind=dbl_kind), dimension(nslyr), intent(in) :: &
909 : rsnw, & ! snow grain radius (10^-6 m)
910 : zqsn ! snow enthalpy (J m-3)
911 :
912 : real (kind=dbl_kind), dimension(nslyr), intent(inout) :: &
913 : drsnw_dry ! change due to snow aging (10^-6 m)
914 :
915 : real (kind=dbl_kind), intent(in) :: &
916 : Tsfc, & ! surface temperature (C)
917 : zTin1, & ! top ice layer temperature (C)
918 : hsn, & ! snow thickness (m)
919 : hin ! ice thickness (m)
920 :
921 : ! local temporary variables
922 :
923 : integer (kind=int_kind) :: k
924 :
925 : integer (kind=int_kind) :: &
926 : T_idx, & ! temperature index
927 : Tgrd_idx, & ! temperature gradient index
928 : rhos_idx ! density index
929 :
930 : real (kind=dbl_kind), dimension(nslyr):: &
931 2591694 : zdTdz, & ! temperature gradient (K/s)
932 863898 : zTsn ! snow temperature (C)
933 :
934 : real (kind=dbl_kind) :: &
935 : bst_tau, & ! snow aging parameter retrieved from lookup table [hour]
936 : bst_kappa, & ! snow aging parameter retrieved from lookup table [unitless]
937 : bst_drdt0, & ! snow aging parameter retrieved from lookup table [um hr-1]
938 : dr_fresh, & ! change in snow radius from fresh (10^-6 m)
939 : dzs, & ! snow layer thickness (m)
940 : dzi, & ! ice layer thickness (m)
941 : dz ! dzs + dzi (m)
942 :
943 : logical (kind=log_kind), save :: &
944 : first_call = .true. ! first call flag
945 :
946 : character (char_len) :: &
947 : string ! generic string for writing messages
948 :
949 : character (len=*),parameter :: subname='(snow_dry_metamorph)'
950 :
951 : !-----------------------------------------------------------------
952 : ! On the first call, check that the table is setup properly
953 : ! Check sizes of 1D and 3D data
954 : ! Check that the 1D data is regularly spaced and set min, del, and lin values
955 : ! for each 1D variable. This info will be used later for the table lookup.
956 : !-----------------------------------------------------------------
957 :
958 863898 : if (first_call) then
959 : if (isnw_rhos < 1 .or. isnw_Tgrd < 1 .or. isnw_T < 1 .or. &
960 : size(snowage_rhos) /= isnw_rhos .or. &
961 : size(snowage_Tgrd) /= isnw_Tgrd .or. &
962 : size(snowage_T) /= isnw_T .or. &
963 : size(snowage_tau) /= isnw_rhos*isnw_Tgrd*isnw_T .or. &
964 110 : size(snowage_kappa) /= isnw_rhos*isnw_Tgrd*isnw_T .or. &
965 : size(snowage_drdt0) /= isnw_rhos*isnw_Tgrd*isnw_T) then
966 0 : write(string,'(a,3i4)') subname//' snowtable size1 = ',isnw_rhos, isnw_Tgrd, isnw_T
967 0 : call icepack_warnings_add(string)
968 0 : write(string,'(a,3i4)') subname//' snowtable size2 = ',size(snowage_rhos),size(snowage_Tgrd),size(snowage_T)
969 0 : call icepack_warnings_add(string)
970 0 : write(string,'(a,3i9)') subname//' snowtable size3 = ',size(snowage_tau),size(snowage_kappa),size(snowage_drdt0)
971 0 : call icepack_warnings_add(string)
972 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
973 0 : call icepack_warnings_add(subname//'ERROR: arrays sizes error')
974 0 : return
975 : endif
976 11 : call snowtable_check_dimension(snowage_rhos, min_rhos, del_rhos, lin_rhos)
977 11 : if (icepack_warnings_aborted(subname)) return
978 11 : call snowtable_check_dimension(snowage_Tgrd, min_Tgrd, del_Tgrd, lin_Tgrd)
979 11 : if (icepack_warnings_aborted(subname)) return
980 11 : call snowtable_check_dimension(snowage_T , min_T , del_T , lin_T )
981 11 : if (icepack_warnings_aborted(subname)) return
982 : endif
983 :
984 : !-----------------------------------------------------------------
985 : ! Needed for variable snow density not currently modeled
986 : ! calculate density based on liquid and ice content of snow
987 : !-----------------------------------------------------------------
988 :
989 5183388 : drsnw_dry(:) = c0
990 5183388 : zTsn(:) = c0
991 5183388 : zdTdz(:) = c0
992 :
993 863898 : dzs = hsn/real(nslyr,kind=dbl_kind)
994 863898 : dzi = hin/real(nilyr,kind=dbl_kind)
995 863898 : dz = dzs + dzi
996 :
997 863898 : zTsn(1) = (Lfresh + zqsn(1)/rhos)/cp_ice
998 863898 : if (nslyr == 1) then
999 : zdTdz(1) = min(c10*isnw_Tgrd, &
1000 : !ech refactored abs((zTsn(1)*dzi+zTin1*dzs)/(dzs+dzi+puny) - Tsfc)/(hsn+puny))
1001 0 : abs(zTsn(1)*dzi + zTin1*dzs - Tsfc*dz)/(dz*hsn+puny))
1002 : else
1003 4319490 : do k = 2, nslyr
1004 3455592 : zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
1005 3455592 : if (k == 2) then
1006 863898 : zdTdz(k-1) = abs((zTsn(k-1)+zTsn(k))*p5 - Tsfc)/(dzs+puny)
1007 : else
1008 2591694 : zdTdz(k-1) = abs (zTsn(k-2)-zTsn(k))*p5 /(dzs+puny)
1009 : endif
1010 4319490 : zdTdz(k-1) = min(c10*isnw_Tgrd,zdTdz(k-1))
1011 : enddo
1012 : !ech refactored zdTdz(nslyr) = abs((zTsn(nslyr)*dzi + zTin1*dzs)/(dzs + dzi+puny) &
1013 : !ech refactored - (zTsn(nslyr) + zTsn(nslyr-1))*p5) / (dzs+puny)
1014 : zdTdz(nslyr) = abs((zTsn(nslyr)*dzi + zTin1*dzs) &
1015 863898 : - (zTsn(nslyr) + zTsn(nslyr-1))*p5*dz) / (dz*dzs+puny)
1016 863898 : zdTdz(nslyr) = min(c10*isnw_Tgrd, zdTdz(nslyr))
1017 : endif
1018 :
1019 5183388 : do k = 1, nslyr
1020 :
1021 : !-----------------------------------------------------------------
1022 : ! best-fit table indices:
1023 : ! Will abort if 1D data is not regularly spaced (lin_* flag must be true)
1024 : ! Leave option for doing a search into non regularly spaced 1D data in the future
1025 : ! via a binary search or similar.
1026 : !-----------------------------------------------------------------
1027 :
1028 4319490 : if (lin_rhos) then
1029 4319490 : rhos_idx = nint( (rhos - min_rhos) / del_rhos, kind=int_kind) + 1
1030 : else
1031 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1032 0 : call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for rhos not supported yet')
1033 0 : return
1034 : endif
1035 :
1036 4319490 : if (lin_Tgrd) then
1037 4319490 : Tgrd_idx = nint( (zdTdz(k) - min_Tgrd) / del_Tgrd, kind=int_kind) + 1
1038 : else
1039 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1040 0 : call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for Tgrd not supported yet')
1041 0 : return
1042 : endif
1043 :
1044 4319490 : if (lin_T) then
1045 4319490 : T_idx = nint(abs(zTsn(k)+ Tffresh - min_T ) / del_T , kind=int_kind) + 1
1046 : else
1047 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1048 0 : call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for T not supported yet')
1049 0 : return
1050 : endif
1051 :
1052 : ! boundary check:
1053 4319490 : rhos_idx = min(isnw_rhos, max(1,rhos_idx))
1054 4319490 : Tgrd_idx = min(isnw_Tgrd, max(1,Tgrd_idx))
1055 4319490 : T_idx = min(isnw_T , max(1,T_idx ))
1056 :
1057 4319490 : bst_tau = snowage_tau (rhos_idx,Tgrd_idx,T_idx)
1058 4319490 : bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx)
1059 4319490 : bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx)
1060 :
1061 : ! change in snow effective radius, using best-fit parameters
1062 4319490 : dr_fresh = max(c0,rsnw(k)-rsnw_fall)
1063 : drsnw_dry(k) = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1/bst_kappa))&
1064 5183388 : * (dt/3600.0_dbl_kind)
1065 : enddo
1066 :
1067 863898 : first_call = .false.
1068 :
1069 : end subroutine snow_dry_metamorph
1070 :
1071 : !=======================================================================
1072 :
1073 : ! Snow grain metamorphism
1074 :
1075 4319490 : subroutine snow_wet_metamorph (dt, dr_wet, rsnw, smice, smliq)
1076 : !
1077 : ! Liquid water redistribution: Apply the grain growth function from:
1078 : ! Brun, E. (1989), Investigation of wet-snow metamorphism in respect of
1079 : ! liquid-water content, Annals of Glaciology, 13, 22-26.
1080 : ! There are two parameters that describe the grain growth rate as
1081 : ! a function of snow liquid water content (LWC). The "LWC=0" parameter
1082 : ! is zeroed here because we are accounting for dry snowing with a
1083 : ! different representation
1084 : !
1085 : real (kind=dbl_kind), intent(in) :: &
1086 : dt ! time step
1087 :
1088 : real (kind=dbl_kind), intent(in) :: &
1089 : rsnw , & ! snow grain radius (10^-6 m)
1090 : smice, & ! snow ice density (kg/m^3)
1091 : smliq ! snow liquid density (kg/m^3)
1092 :
1093 : real (kind=dbl_kind), intent(inout) :: &
1094 : dr_wet
1095 :
1096 : real (kind=dbl_kind) :: &
1097 : fliq ! liquid mass fraction
1098 :
1099 : character (len=*),parameter :: subname='(snow_wet_metamorph)'
1100 :
1101 4319490 : dr_wet = c0
1102 4319490 : fliq = c1
1103 4319490 : if (smice + smliq > c0 .and. rsnw > c0) then
1104 4319490 : fliq = min(smliq/(smice + smliq),p1)
1105 4319490 : dr_wet = S_wet * fliq**3*dt/(c4*pi*rsnw**2)
1106 : endif
1107 :
1108 4319490 : end subroutine snow_wet_metamorph
1109 :
1110 : !=======================================================================
1111 :
1112 : ! Analyze 1D array for regular spacing for snow table lookup
1113 : ! and set the min, del, and lin values.
1114 : ! Tolerance for regular spacing set at 1.0e-8 * typical array value
1115 :
1116 33 : subroutine snowtable_check_dimension(array, min_fld, del_fld, lin_fld)
1117 :
1118 : real (kind=dbl_kind), intent(in), dimension(:) :: &
1119 : array ! array to check
1120 :
1121 : real (kind=dbl_kind), intent(inout) :: &
1122 : min_fld, & ! min value if linear
1123 : del_fld ! delta value if linear
1124 :
1125 : logical (kind=log_kind), intent(inout) :: &
1126 : lin_fld ! linear flag
1127 :
1128 : ! local temporary variables
1129 :
1130 : integer (kind=int_kind) :: n, asize
1131 :
1132 : real (kind=dbl_kind) :: tolerance ! tolerance for linear checking
1133 :
1134 : character (len=*),parameter :: subname='(snowtable_check_dimension)'
1135 :
1136 33 : asize = size(array)
1137 :
1138 33 : if (asize == 1) then
1139 11 : min_fld = array(1)
1140 11 : del_fld = array(1) ! arbitrary
1141 11 : lin_fld = .true.
1142 : else
1143 22 : lin_fld = .true.
1144 22 : min_fld = array(1)
1145 22 : del_fld = array(2) - array(1)
1146 22 : tolerance = 1.0e-08_dbl_kind * max(abs(array(1)),abs(array(2))) ! relative to typical value
1147 88 : do n = 3,asize
1148 88 : if (abs(array(n) - array(n-1) - del_fld) > tolerance) lin_fld = .false.
1149 : enddo
1150 : endif
1151 :
1152 33 : end subroutine snowtable_check_dimension
1153 :
1154 : !=======================================================================
1155 :
1156 : ! Conversions between ice mass, liquid water mass in snow
1157 :
1158 1521120 : subroutine drain_snow (vsnon, aicen, &
1159 1521120 : massice, massliq, meltsliq)
1160 :
1161 : real (kind=dbl_kind), intent(in) :: &
1162 : vsnon, & ! snow volume (m)
1163 : aicen ! aice area fraction
1164 :
1165 : real (kind=dbl_kind), intent(inout) :: &
1166 : meltsliq ! total liquid content (kg/m^2)
1167 :
1168 : real (kind=dbl_kind), dimension(nslyr), intent(in) :: &
1169 : massice ! mass of ice in snow (kg/m^2)
1170 :
1171 : real (kind=dbl_kind), dimension(nslyr), intent(inout) :: &
1172 : massliq ! mass of liquid in snow (kg/m^2)
1173 :
1174 : ! local temporary variables
1175 :
1176 : integer (kind=int_kind) :: k
1177 :
1178 : real (kind=dbl_kind) :: &
1179 : hslyr, & ! snow layer thickness (m)
1180 : hsn, & ! snow thickness (m)
1181 : sliq ! snow liquid content (kg/m^2)
1182 :
1183 : real (kind=dbl_kind), dimension(nslyr) :: &
1184 3042240 : dlin , & ! liquid mass into the layer from above (kg/m^2)
1185 3042240 : dlout , & ! liquid mass out of the layer (kg/m^2)
1186 3042240 : phi_liq , & ! volumetric liquid fraction
1187 3042240 : phi_ice ! volumetric ice fraction
1188 :
1189 : character (len=*), parameter :: subname='(drain_snow)'
1190 :
1191 1521120 : hsn = c0
1192 1521120 : sliq = c0
1193 1521120 : if (aicen > c0) hsn = vsnon/aicen
1194 1521120 : if (hsn > puny) then
1195 5153766 : dlin (:) = c0
1196 5153766 : dlout(:) = c0
1197 858961 : hslyr = hsn / real(nslyr,kind=dbl_kind)
1198 5153766 : do k = 1, nslyr
1199 4294805 : massliq(k) = massliq(k) + dlin(k) ! add liquid in from layer above
1200 4294805 : phi_ice(k) = min(c1, massice(k) / (rhoi *hslyr))
1201 4294805 : phi_liq(k) = massliq(k) / (rhofresh*hslyr)
1202 4294805 : dlout(k) = max(c0, (phi_liq(k) - S_r*(c1-phi_ice(k))) * rhofresh * hslyr)
1203 4294805 : massliq(k) = massliq(k) - dlout(k)
1204 5153766 : if (k < nslyr) then
1205 3435844 : dlin(k+1) = dlout(k)
1206 : else
1207 858961 : sliq = dlout(nslyr) ! this (re)initializes meltsliq
1208 : endif
1209 : enddo
1210 : else
1211 662159 : sliq = meltsliq ! computed in thickness_changes
1212 : endif
1213 1521120 : if (use_smliq_pnd) meltsliq = sliq
1214 :
1215 1521120 : end subroutine drain_snow
1216 :
1217 : !=======================================================================
1218 :
1219 : end module icepack_snow
1220 :
1221 : !=======================================================================
|
