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