SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
snow_module Module Reference

Functions/Subroutines

subroutine snow_cal_meltheat (snowUse, tstep, tau_r, SnowDensMax, lvS_J_kg, lv_J_kg, tstep_real, RadMeltFact, TempMeltFact, SnowAlbMax, SnowDensMin, Temp_C, Precip, PrecipLimit, PrecipLimitAlb, nsh_real, sfr, Tsurf_ind, Tsurf_ind_snow, state_id, qn1_ind_snow, kup_ind_snow, SnowWater, deltaQi, alb1, SnowPack_in, SnowFrac_in, SnowAlb_in, SnowDens_in, SnowfallCum_in, SnowPack_out, SnowFrac_out, SnowAlb_out, SnowDens_out, SnowfallCum_out, mwh, Qm, QmFreez, QmRain, veg_fr, snowCalcSwitch, Qm_melt, Qm_freezState, Qm_rain, FreezMelt, FreezState, FreezStateVol, rainOnSnow, SnowDepth, mw_ind, dataOutLineSnow)
 
subroutine meltheat (lvS_J_kg, lv_J_kg, tstep_real, RadMeltFact, TempMeltFact, SnowAlbMax, SnowDensMin, Temp_C, Precip, PrecipLimit, PrecipLimitAlb, nsh_real, waterdens, sfr, Tsurf_ind, state_id, qn1_ind_snow, SnowWater, deltaQi, SnowPack, SnowFrac, SnowAlb, SnowDens, SnowfallCum, mwh, fwh, Qm, QmFreez, QmRain, snowCalcSwitch, Qm_melt, Qm_freezState, Qm_rain, FreezMelt, FreezState, FreezStateVol, rainOnSnow, SnowDepth, mw_ind)
 
subroutine snowcalc (tstep, imin, it, dectime, is, EvapMethod, CRWmin, CRWmax, nsh_real, lvS_J_kg, avdens, avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, avcp, sIce_hPa, PervFraction, vegfraction, addimpervious, vpd_hPa, qn_e, s_hPa, ResistSurf, RA, rb, tlv, snowdensmin, SnowProf_24hr, precip, PipeCapacity, RunoffToWater, addVeg, SnowLimPaved, SnowLimBldg, FlowChange, drain, WetThresh, stateOld, mw_ind, SoilStoreCap, rainonsnow, freezmelt, freezstate, freezstatevol, Qm_Melt, Qm_rain, Tsurf_ind, sfr, dayofWeek_id, StoreDrainPrm, SnowPackLimit, AddWater, addwaterrunoff, soilstore_id, SnowPack, SurplusEvap, SnowFrac, SnowWater, iceFrac, SnowDens, runoffAGimpervious, runoffAGveg, surplusWaterBody, rss_nsurf, runoffSnow, runoff, runoffSoil, chang, changSnow, SnowToSurf, state_id, ev_snow, SnowDepth, SnowRemoval, swe, ev, chSnow_tot, ev_tot, qe_tot, runoff_tot, surf_chang_tot, runoffPipes, mwstore, runoffwaterbody)
 
subroutine evap_suews_snow (Qm, QP, lvS_J_kg, avdens, avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, tstep, avcp, sIce_hPa, dectime, ev_snow, tlv_sub)
 
subroutine snowrem (is, PavSurf, BldgSurf, nsurf, SnowFrac, sfr, SnowPack, SnowRemoval, SnowLimPaved, SnowLimBldg)
 
real(kind(1d0)) function snowdepletioncurve (is, swe, sweD)
 
subroutine veg_fr_snow (sfr, SnowFrac, veg_fr)
 
subroutine snowupdate (tstep, Temp_C, tau_a, tau_f, tau_r, SnowDensMax, SnowDensMin, SnowAlbMax, SnowAlbMin, SnowPack_prev, SnowAlb_prev, SnowDens_prev, SnowAlb_next, SnowDens_next)
 
real(kind(1d0)) function update_snow_albedo (tstep, SnowPack_prev, SnowAlb_prev, Temp_C, tau_a, tau_f, SnowAlbMax, SnowAlbMin)
 
real(kind(1d0)) function, dimension(nsurf) update_snow_dens (tstep, SnowPack_prev, SnowDens_prev, tau_r, SnowDensMax, SnowDensMin)
 

Function/Subroutine Documentation

◆ evap_suews_snow()

subroutine snow_module::evap_suews_snow ( real(kind(1d0)), intent(in)  Qm,
real(kind(1d0)), intent(in)  QP,
real(kind(1d0)), intent(in)  lvS_J_kg,
real(kind(1d0)), intent(in)  avdens,
real(kind(1d0)), intent(in)  avRh,
real(kind(1d0)), intent(in)  Press_hPa,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  RAsnow,
real(kind(1d0)), intent(in)  psyc_hPa,
integer, intent(in)  tstep,
real(kind(1d0)), intent(in)  avcp,
real(kind(1d0)), intent(in)  sIce_hPa,
real(kind(1d0)), intent(in)  dectime,
real(kind(1d0)), intent(out)  ev_snow,
real(kind(1d0)), intent(out)  tlv_sub 
)

Definition at line 1156 of file suews_phys_snow.f95.

References meteo::sat_vap_pressice().

Referenced by snowcalc().

1156 
1157  USE meteo, ONLY: sat_vap_pressice
1158 
1159  IMPLICIT NONE
1160 
1161  !INPUT
1162  REAL(KIND(1d0)), INTENT(in)::qm !melt heat,
1163  REAL(KIND(1d0)), INTENT(in)::qp !advect. heat
1164  REAL(KIND(1d0)), INTENT(in)::lvs_j_kg !latent heat of sublimation
1165  REAL(KIND(1d0)), INTENT(in)::avdens !air density
1166  REAL(KIND(1d0)), INTENT(in)::avrh !relative humidity
1167  REAL(KIND(1d0)), INTENT(in)::press_hpa !air pressure
1168  REAL(KIND(1d0)), INTENT(in)::temp_c !air temperature
1169  REAL(KIND(1d0)), INTENT(in)::rasnow !aerodyn res snow
1170  REAL(KIND(1d0)), INTENT(in)::psyc_hpa !psychometric constant
1171  REAL(KIND(1d0)), INTENT(in)::avcp !spec. heat,
1172  REAL(KIND(1d0)), INTENT(in)::sice_hpa !satured curve on snow
1173  REAL(KIND(1d0)), INTENT(in)::dectime
1174  INTEGER, INTENT(in):: tstep
1175 
1176  REAL(KIND(1d0)), INTENT(out)::ev_snow !Evaporation
1177  REAL(KIND(1d0)), INTENT(out)::tlv_sub !Latent heat for sublimation
1178 
1179  !OTHER VARIABLES
1180  REAL(KIND(1d0))::e_snow, & !PM equation obe line
1181  sae_snow, & !s * (Available energy)
1182  qe_snow, & !Latent heat flux
1183  vdrcice, & !Vapour pressure deficit
1184  esice_hpa, & !Saturation vapor pressure over ice
1185  eaice_hpa, & !Vapour pressure
1186  tstep_real !timestep as real
1187 
1188  INTEGER::from = 1
1189  !-----------------------------------------------------
1190 
1191  tstep_real = REAL(tstep, kind(1d0))
1192 
1193  sae_snow = sice_hpa*(qp - qm) !Calculate the driving parameter in calculation of evaporation. Järvi et al. (2015)
1194 
1195  esice_hpa = sat_vap_pressice(temp_c, press_hpa, from, dectime) !Saturation vapor pressure over ice
1196  eaice_hpa = avrh/100*esice_hpa !Vapour pressure of water
1197  vdrcice = (esice_hpa - eaice_hpa)*avdens*avcp !Vapour pressure deficit
1198  tlv_sub = lvs_j_kg/tstep_real !Latent heat for sublimation
1199  e_snow = sae_snow + vdrcice/rasnow !PM equation
1200  qe_snow = e_snow/(sice_hpa + psyc_hpa) !Latent heat (W/m^2)
1201  ev_snow = qe_snow/tlv_sub !Evaporation (in mm)
1202 
1203  RETURN
1204 
real(kind(1d0)) lvs_j_kg
real(kind(1d0)) qm
real(kind(1d0)) avcp
real(kind(1d0)) psyc_hpa
real(kind(1d0)) function sat_vap_pressice(Temp_c, PRESS_hPa, from, dectime)
real(kind(1d0)) sice_hpa
real(kind(1d0)) dectime
real(kind(1d0)) avdens
Here is the call graph for this function:
Here is the caller graph for this function:

◆ meltheat()

subroutine snow_module::meltheat ( real(kind(1d0)), intent(in)  lvS_J_kg,
real(kind(1d0)), intent(in)  lv_J_kg,
real(kind(1d0)), intent(in)  tstep_real,
real(kind(1d0)), intent(in)  RadMeltFact,
real(kind(1d0)), intent(in)  TempMeltFact,
real(kind(1d0)), intent(in)  SnowAlbMax,
real(kind(1d0)), intent(in)  SnowDensMin,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  Precip,
real(kind(1d0)), intent(in)  PrecipLimit,
real(kind(1d0)), intent(in)  PrecipLimitAlb,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), intent(in)  waterdens,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), dimension(nsurf), intent(in)  Tsurf_ind,
real(kind(1d0)), dimension(nsurf), intent(in)  state_id,
real(kind(1d0)), dimension(nsurf), intent(in)  qn1_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowWater,
real(kind(1d0)), dimension(nsurf), intent(in)  deltaQi,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowPack,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowFrac,
real(kind(1d0)), intent(inout)  SnowAlb,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowDens,
real(kind(1d0)), intent(inout)  SnowfallCum,
real(kind(1d0)), intent(out)  mwh,
real(kind(1d0)), intent(out)  fwh,
real(kind(1d0)), intent(out)  Qm,
real(kind(1d0)), intent(out)  QmFreez,
real(kind(1d0)), intent(out)  QmRain,
integer, dimension(nsurf), intent(out)  snowCalcSwitch,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_melt,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_freezState,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_rain,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezMelt,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezState,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezStateVol,
real(kind(1d0)), dimension(nsurf), intent(out)  rainOnSnow,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowDepth,
real(kind(1d0)), dimension(nsurf), intent(out)  mw_ind 
)

Definition at line 235 of file suews_phys_snow.f95.

References allocatearray::bldgsurf, allocatearray::nsurf, allocatearray::pavsurf, and allocatearray::watersurf.

Referenced by snow_cal_meltheat().

235 
236  IMPLICIT NONE
237 
238  !These are input to the module
239  ! INTEGER, INTENT(in)::bldgsurf
240  ! INTEGER, INTENT(in)::nsurf
241  ! INTEGER, INTENT(in)::PavSurf
242  ! INTEGER, INTENT(in)::WaterSurf
243 
244  REAL(KIND(1d0)), INTENT(in)::lvs_j_kg
245  REAL(KIND(1d0)), INTENT(in)::lv_j_kg
246  REAL(KIND(1d0)), INTENT(in)::tstep_real
247  REAL(KIND(1d0)), INTENT(in)::radmeltfact
248  REAL(KIND(1d0)), INTENT(in)::tempmeltfact
249  REAL(KIND(1d0)), INTENT(in)::snowalbmax
250  REAL(KIND(1d0)), INTENT(in)::snowdensmin
251  REAL(KIND(1d0)), INTENT(in)::temp_c
252  REAL(KIND(1d0)), INTENT(in)::precip
253  REAL(KIND(1d0)), INTENT(in)::preciplimit
254  REAL(KIND(1d0)), INTENT(in)::preciplimitalb
255  REAL(KIND(1d0)), INTENT(in)::nsh_real
256  REAL(KIND(1d0)), INTENT(in)::waterdens
257 
258  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
259  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::tsurf_ind
260  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::state_id
261  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::qn1_ind_snow
262  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowwater
263  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::deltaqi
264 
265  !Input and output as this is updated in this subroutine
266  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowpack
267  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowfrac
268  REAL(KIND(1d0)), INTENT(inout)::snowalb
269  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowdens
270  REAL(KIND(1d0)), INTENT(inout)::snowfallcum
271 
272  !Output:
273  REAL(KIND(1d0)), INTENT(out)::mwh
274  REAL(KIND(1d0)), INTENT(out)::fwh
275  REAL(KIND(1d0)), INTENT(out)::qm
276  REAL(KIND(1d0)), INTENT(out)::qmfreez
277  REAL(KIND(1d0)), INTENT(out)::qmrain
278 
279  !Output, dimension nsurf
280  INTEGER, DIMENSION(nsurf), INTENT(out)::snowcalcswitch
281 
282  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_melt
283  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_freezstate
284  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_rain
285  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezmelt
286  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezstate
287  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezstatevol
288  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::rainonsnow
289  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowdepth
290  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::mw_ind
291 
292  ! local variables:
293  REAL(KIND(1d0))::adjmeltfact
294  REAL(KIND(1d0))::watfreeze
295 
296  REAL(KIND(1d0)), PARAMETER::cw = 4190 !,ci=2090 !Specific heat capacity of water
297 
298  INTEGER::is, xx
299 
300  !Initialize snow variables
301  mwh = 0 !Initialize snow melt and heat related to snowmelt
302  fwh = 0
303  qm = 0
304  qmfreez = 0
305  qmrain = 0
306  snowcalcswitch = 0
307  qm_melt = 0
308  qm_freezstate = 0
309  qm_rain = 0
310  freezmelt = 0
311  freezstate = 0
312  freezstatevol = 0
313  rainonsnow = 0
314  snowdepth = 0
315  mw_ind = 0
316 
317  !===dummy calculations===
318  xx = bldgsurf
319  xx = pavsurf
320  !===dummy calculations end===
321 
322  !=========================================================================================
323  DO is = 1, nsurf !Go each surface type through
324  IF (sfr(is) /= 0) THEN !If surface type existing,
325 
326  IF (snowpack(is) > 0) THEN !If SnowPack existing, calculate meltwater related water flows
327 
328  snowdepth(is) = (snowpack(is)/1000)*waterdens/snowdens(is) !Snow depth in m
329 
330  !Calculate meltwater related water flows with hourly degree-day method.
331 
332  !These are for snow melting
333  IF (temp_c >= 0) THEN
334  IF (qn1_ind_snow(is) < 0) THEN
335  mw_ind(is) = tempmeltfact*temp_c !(mm C−1 h−1)*(C) = in mm h-1
336  ELSE
337  mw_ind(is) = radmeltfact*(qn1_ind_snow(is)) !(mm m2 W−1 h−1)*(W m-2)= mm h-1 ??
338  ENDIF
339 
340  ELSE !Freezing equations
341  adjmeltfact = 1 !Relationship between the temperature melt and freezing factors
342  mw_ind(is) = tempmeltfact*temp_c*adjmeltfact ! in mm h-1
343  ENDIF
344  !Previous equation give the hourly values, divide these with the timestep number
345  mw_ind(is) = mw_ind(is)/nsh_real
346 
347  IF (mw_ind(is) > snowpack(is)) mw_ind(is) = snowpack(is)!Limited by the previous timestep SnowPack
348 
349  !-----------------------------------------------------
350  ! Heat consumed to snowmelt/refreezing within Tstep.
351  ! Converted from mm nsh-1 to mm nsh-1 and to m s-1
352  qm_melt(is) = waterdens*((mw_ind(is)/tstep_real)/1000)*(lvs_j_kg - lv_j_kg)
353 
354  !If melt is negative this means freezing water in the SnowPack
355  IF (mw_ind(is) < 0) THEN
356 
357  freezmelt(is) = -mw_ind(is) !Save this to variable FreezMelt
358  mw_ind(is) = 0
359 
360  !Freezing water cannot exceed meltwater store
361  IF (freezmelt(is) > snowwater(is)) freezmelt(is) = snowwater(is)
362 
363  !Recalculate melt related energy
364  qm_melt(is) = waterdens*((-freezmelt(is)/tstep_real)/1000)*(lvs_j_kg - lv_j_kg)
365  ENDIF
366 
367  !-----------------------------------------------------
368  ! If air temperature is above zero, precipitation causes advective heat to the
369  ! SnowPack. Eq (23) in Sun et al., 1999
370  ! Calculation done at resolution of the model timestep
371  IF (temp_c >= preciplimit .AND. precip > 0) THEN
372  qm_rain(is) = waterdens*cw*(temp_c - preciplimit)*(precip*0.001/tstep_real) !in W m-2
373  IF (qm_rain(is) < 0) THEN !Can only be positive
374  qm_rain(is) = 0
375  ELSE
376  rainonsnow(is) = precip !Save information on the rain on snow event
377  ENDIF
378  ENDIF
379 
380  ENDIF !End if SnowPack
381 
382  !=================================================================
383 
384  !Freeze surface water state_id if cold enough.
385  IF (tsurf_ind(is) < 0 .AND. state_id(is) > 0) THEN
386 
387  snowcalcswitch(is) = 1 !If water on ground this forms ice and snow calculations are made
388 
389  !Other surfaces than water treated first
390  IF (is /= watersurf) THEN
391 
392  !FreezState(is) = state_id(is)
393  !Previously all state_id could freeze in 5-min timestep. Now we calculate how much water
394  !can freeze in a timestep based on the same temperature freezing fraction.
395  freezstate(is) = -tempmeltfact*tsurf_ind(is)/nsh_real
396 
397  !The amount of freezing water cannot be greater than the surface state_id
398  IF (freezstate(is) > state_id(is)) freezstate(is) = state_id(is)
399 
400  IF (snowpack(is) == 0 .OR. snowfrac(is) == 0) THEN !SnowPack forms
401  freezstatevol(is) = freezstate(is)
402  ELSE !There is snow already on ground
403  freezstatevol(is) = freezstate(is)*(1 - snowfrac(is))/snowfrac(is)
404  ENDIF
405 
406  ! If the amount of freezing water is very small and there is state_id left to the ground
407  ! no freezing of water will take place
408  IF (freezstatevol(is) < 0.00000000001 .AND. freezstate(is) < state_id(is)) THEN
409  freezstate(is) = 0
410  freezstatevol(is) = 0
411  ENDIF
412 
413  !Calculate the heat exchange in W m-2
414  qm_freezstate(is) = -waterdens*(freezstate(is)/tstep_real/1000)*(lvs_j_kg - lv_j_kg)
415 
416  !Water surface separately
417  ELSE
418  !Calculate average value how much water can freeze above the water areas
419  !Equation is -hA(T-T0) = rhoV(Cp+dT +Lf) in 5-min timestep
420  !h=convective heat trasnfer,A, area of water,rwo water density,V volume, dT temperature difference
421  !before and end of the 5-min period. dT equals zero, h=100 and when multiplied with Area, the equation
422  !simplyfies to the this. LJ 14 July 2015
423  watfreeze = 100*(0 - temp_c)/(waterdens*(lvs_j_kg - lv_j_kg))
424  freezstate(is) = watfreeze
425  qm_freezstate(is) = -waterdens*(watfreeze/tstep_real/1000)*(lvs_j_kg - lv_j_kg)
426  ENDIF
427 
428  ENDIF
429 
430  !======================================================================
431  ! Define if any snowmelt calculations are made: SnowPack existing,
432  ! freezing occuring on ground or from precip
433  IF (is /= watersurf) THEN
434  IF (snowpack(is) > 0 .OR. (precip > 0 .AND. tsurf_ind(is) < 0)) THEN
435  snowcalcswitch(is) = 1
436  ENDIF
437  ELSE !Water surface separately
438  IF (snowpack(watersurf) > 0 .OR. freezstate(watersurf) > 0) THEN
439  snowcalcswitch(watersurf) = 1
440  ENDIF
441  ENDIF
442 
443  !Update snow density of each surface
444  IF (precip > 0 .AND. tsurf_ind(is) < 0 .AND. snowpack(is) > 0) THEN
445  snowdens(is) = snowdens(is)*snowpack(is)/(snowpack(is) + precip) + snowdensmin*precip/(snowpack(is) + precip)
446  ENDIF
447 
448  !Weighted variables for the whole area
449  mwh = mwh + mw_ind(is)*sfr(is)*snowfrac(is) !Snowmelt
450  fwh = fwh + freezmelt(is)*sfr(is)*snowfrac(is) !Freezing water
451  qm = qm + qm_melt(is)*sfr(is)*snowfrac(is) !Energy consumed to the melt/freezing.
452  qmrain = qmrain + qm_rain(is)*sfr(is)*snowfrac(is) !Rain on snow
453  qmfreez = qmfreez + deltaqi(is)*sfr(is)*snowfrac(is) + qm_freezstate(is)*sfr(is)*(1 - snowfrac(is)) !Freezing water
454  ENDIF
455 
456  ENDDO !End surface type
457 
458  !Update snow albedo to its maximum value if precipitation exists
459  IF (precip > 0 .AND. sum(snowpack) > 0 .AND. temp_c < 0) THEN
460 
461  snowfallcum = snowfallcum + precip
462 
463  IF (snowfallcum > preciplimitalb) THEN
464 
466  snowfallcum = 0
467  ENDIF
468  ELSE
469 
470  snowfallcum = 0
471  ENDIF
472 
real(kind(1d0)) waterdens
real(kind(1d0)) snowfallcum
real(kind(1d0)) snowalb
real(kind(1d0)) tempmeltfact
real(kind(1d0)) preciplimit
real(kind(1d0)) lvs_j_kg
real(kind(1d0)) lv_j_kg
real(kind(1d0)) fwh
real(kind(1d0)) qm
real(kind(1d0)) snowdensmin
real(kind(1d0)) adjmeltfact
real(kind(1d0)) snowalbmax
real(kind(1d0)) radmeltfact
real(kind(1d0)) mwh
real(kind(1d0)) qmrain
real(kind(1d0)) qmfreez
real(kind(1d0)) preciplimitalb
Here is the caller graph for this function:

◆ snow_cal_meltheat()

subroutine snow_module::snow_cal_meltheat ( integer, intent(in)  snowUse,
integer, intent(in)  tstep,
real(kind(1d0)), intent(in)  tau_r,
real(kind(1d0)), intent(in)  SnowDensMax,
real(kind(1d0)), intent(in)  lvS_J_kg,
real(kind(1d0)), intent(in)  lv_J_kg,
real(kind(1d0)), intent(in)  tstep_real,
real(kind(1d0)), intent(in)  RadMeltFact,
real(kind(1d0)), intent(in)  TempMeltFact,
real(kind(1d0)), intent(in)  SnowAlbMax,
real(kind(1d0)), intent(in)  SnowDensMin,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  Precip,
real(kind(1d0)), intent(in)  PrecipLimit,
real(kind(1d0)), intent(in)  PrecipLimitAlb,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), dimension(nsurf), intent(in)  Tsurf_ind,
real(kind(1d0)), dimension(nsurf), intent(in)  Tsurf_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(in)  state_id,
real(kind(1d0)), dimension(nsurf), intent(in)  qn1_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(in)  kup_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowWater,
real(kind(1d0)), dimension(nsurf), intent(in)  deltaQi,
real(kind(1d0)), intent(in)  alb1,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowPack_in,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowFrac_in,
real(kind(1d0)), intent(in)  SnowAlb_in,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowDens_in,
real(kind(1d0)), intent(in)  SnowfallCum_in,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowPack_out,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowFrac_out,
real(kind(1d0)), intent(out)  SnowAlb_out,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowDens_out,
real(kind(1d0)), intent(out)  SnowfallCum_out,
real(kind(1d0)), intent(out)  mwh,
real(kind(1d0)), intent(out)  Qm,
real(kind(1d0)), intent(out)  QmFreez,
real(kind(1d0)), intent(out)  QmRain,
real(kind(1d0)), intent(out)  veg_fr,
integer, dimension(nsurf), intent(out)  snowCalcSwitch,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_melt,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_freezState,
real(kind(1d0)), dimension(nsurf), intent(out)  Qm_rain,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezMelt,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezState,
real(kind(1d0)), dimension(nsurf), intent(out)  FreezStateVol,
real(kind(1d0)), dimension(nsurf), intent(out)  rainOnSnow,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowDepth,
real(kind(1d0)), dimension(nsurf), intent(out)  mw_ind,
real(kind(1d0)), dimension(ncolumnsdataoutsnow_notime), intent(out)  dataOutLineSnow 
)

Definition at line 46 of file suews_phys_snow.f95.

References meltheat(), allocatearray::nsurf, update_snow_dens(), and veg_fr_snow().

Referenced by suews_driver::suews_cal_main().

46 
47  IMPLICIT NONE
48  ! INTEGER, PARAMETER::nsurf = 7
49  ! INTEGER, PARAMETER::PavSurf = 1
50  ! INTEGER, PARAMETER::BldgSurf = 2
51  ! INTEGER, PARAMETER::WaterSurf = 7
52  INTEGER, PARAMETER::ncolumnsdataoutsnow_notime = ncolumnsdataoutsnow - 5
53  REAL(KIND(1d0)), PARAMETER::waterdens = 999.8395 !Density of water in 0 cel deg
54 
55  !These are input to the module
56  INTEGER, INTENT(in)::snowuse
57  INTEGER, INTENT(in)::tstep
58  ! INTEGER,INTENT(in)::bldgsurf
59  ! INTEGER,INTENT(in)::nsurf
60  ! INTEGER,INTENT(in)::PavSurf
61  ! INTEGER,INTENT(in)::WaterSurf
62 
63  REAL(KIND(1d0)), INTENT(in)::lvs_j_kg
64  REAL(KIND(1d0)), INTENT(in)::lv_j_kg
65  REAL(KIND(1d0)), INTENT(in)::tstep_real
66  REAL(KIND(1d0)), INTENT(in)::radmeltfact
67  REAL(KIND(1d0)), INTENT(in)::tempmeltfact
68  REAL(KIND(1d0)), INTENT(in)::snowalbmax
69  REAL(KIND(1d0)), INTENT(in)::snowdensmax
70  REAL(KIND(1d0)), INTENT(in)::snowdensmin
71  REAL(KIND(1d0)), INTENT(in)::temp_c
72  REAL(KIND(1d0)), INTENT(in)::precip
73  REAL(KIND(1d0)), INTENT(in)::preciplimit
74  REAL(KIND(1d0)), INTENT(in)::preciplimitalb
75  REAL(KIND(1d0)), INTENT(in)::nsh_real
76  REAL(KIND(1d0)), INTENT(in)::alb1
77  REAL(KIND(1d0)), INTENT(in)::tau_r
78  ! REAL(KIND(1d0)),INTENT(in)::waterdens
79 
80  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
81  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::tsurf_ind
82  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::tsurf_ind_snow
83  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::state_id
84  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::qn1_ind_snow
85  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::kup_ind_snow
86  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowwater
87  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::deltaqi
88 
89  !Input and output as this is updated in this subroutine
90  REAL(KIND(1d0)), INTENT(in)::snowalb_in
91  REAL(KIND(1d0)), INTENT(in)::snowfallcum_in
92  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowpack_in
93  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowfrac_in
94  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowdens_in
95  REAL(KIND(1d0)), INTENT(out)::snowalb_out
96  REAL(KIND(1d0)), INTENT(out)::snowfallcum_out
97  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowpack_out
98  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowfrac_out
99  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowdens_out
100 
101  REAL(KIND(1d0))::snowalb
102  REAL(KIND(1d0))::snowfallcum
103  REAL(KIND(1d0)), DIMENSION(nsurf)::snowpack
104  REAL(KIND(1d0)), DIMENSION(nsurf)::snowfrac
105  REAL(KIND(1d0)), DIMENSION(nsurf)::snowdens
106 
107  !Output:
108  REAL(KIND(1d0)), INTENT(out)::mwh
109  REAL(KIND(1d0))::fwh
110  REAL(KIND(1d0)), INTENT(out)::qm
111  REAL(KIND(1d0)), INTENT(out)::qmfreez
112  REAL(KIND(1d0)), INTENT(out)::qmrain
113 
114  REAL(KIND(1d0)), INTENT(out)::veg_fr
115 
116  INTEGER, DIMENSION(nsurf), INTENT(out)::snowcalcswitch
117 
118  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_melt
119  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_freezstate
120  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::qm_rain
121  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezmelt
122  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezstate
123  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::freezstatevol
124  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::rainonsnow
125  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowdepth
126  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::mw_ind
127 
128  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutSnow_notime), INTENT(out) :: dataoutlinesnow
129 
130  snowalb = snowalb_in
131  snowfallcum = snowfallcum_in
132  snowpack = snowpack_in
133  snowfrac = snowfrac_in
134  snowdens = snowdens_in
135 
136  IF (snowuse == 1) THEN
137  snowdens = update_snow_dens( &
138  tstep, snowfrac_in, snowdens_in, &
140 
141  CALL meltheat( &
142  lvs_j_kg, lv_j_kg, tstep_real, radmeltfact, tempmeltfact, &
143  snowalbmax, snowdensmin, temp_c, precip, preciplimit, preciplimitalb, &
144  nsh_real, waterdens, sfr, tsurf_ind, state_id, qn1_ind_snow, &
145  snowwater, deltaqi, snowpack, snowfrac, snowalb, snowdens, snowfallcum, &
146  mwh, fwh, qm, qmfreez, qmrain, snowcalcswitch, &
147  qm_melt, qm_freezstate, qm_rain, freezmelt, freezstate, freezstatevol, &
148  rainonsnow, snowdepth, mw_ind)
149 
150  ELSE ! no snow calculation
151  mwh = 0
152  fwh = 0
153  qm = 0
154  qmfreez = 0
155  qmrain = 0
156  snowfallcum = 0
157  snowcalcswitch = 0
158  qm_melt = 0
159  qm_freezstate = 0
160  qm_rain = 0
161  freezmelt = 0
162  freezstate = 0
163  freezstatevol = 0
164  rainonsnow = 0
165  snowdepth = 0
166  mw_ind = 0
167  snowfrac = 0
168 
169  END IF
170 
171  ! update veg_fr
172  CALL veg_fr_snow( &
173  sfr, snowfrac, &!input
174  veg_fr)!output
175 
176  snowalb_out = snowalb
177  snowfallcum_out = snowfallcum
178  snowpack_out = snowpack
179  snowfrac_out = snowfrac
180  snowdens_out = snowdens
181 
182  ! pack output into one line
183  dataoutlinesnow = [ &
184  snowpack_out(1:nsurf), mw_ind(1:nsurf), qm_melt(1:nsurf), & !26
185  qm_rain(1:nsurf), qm_freezstate(1:nsurf), snowfrac_out(1:(nsurf - 1)), & !46
186  rainonsnow(1:nsurf), & !53
187  qn1_ind_snow(1:nsurf), kup_ind_snow(1:nsurf), freezmelt(1:nsurf), & !74
188  snowwater(1:nsurf), snowdens_out(1:nsurf), & !88
189  snowdepth(1:nsurf), tsurf_ind_snow(1:nsurf), &
190  alb1]
191  ! dataOutLineSnow=set_nan(dataOutLineSnow)
192 
real(kind(1d0)) tau_r
real(kind(1d0)) waterdens
real(kind(1d0)) snowfallcum
real(kind(1d0)) snowalb
real(kind(1d0)) tempmeltfact
real(kind(1d0)) preciplimit
real(kind(1d0)) lvs_j_kg
real(kind(1d0)) lv_j_kg
real(kind(1d0)) fwh
real(kind(1d0)) qm
real(kind(1d0)) snowdensmin
real(kind(1d0)) snowalbmax
real(kind(1d0)) radmeltfact
real(kind(1d0)) mwh
real(kind(1d0)) qmrain
real(kind(1d0)) qmfreez
real(kind(1d0)) snowdensmax
real(kind(1d0)) preciplimitalb
Here is the call graph for this function:
Here is the caller graph for this function:

◆ snowcalc()

subroutine snow_module::snowcalc ( integer, intent(in)  tstep,
integer, intent(in)  imin,
integer, intent(in)  it,
real(kind(1d0)), intent(in)  dectime,
integer, intent(in)  is,
integer, intent(in)  EvapMethod,
real(kind(1d0)), intent(in)  CRWmin,
real(kind(1d0)), intent(in)  CRWmax,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), intent(in)  lvS_J_kg,
real(kind(1d0)), intent(in)  avdens,
real(kind(1d0)), intent(in)  avRh,
real(kind(1d0)), intent(in)  Press_hPa,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  RAsnow,
real(kind(1d0)), intent(in)  psyc_hPa,
real(kind(1d0)), intent(in)  avcp,
real(kind(1d0)), intent(in)  sIce_hPa,
real(kind(1d0)), intent(in)  PervFraction,
real(kind(1d0)), intent(in)  vegfraction,
real(kind(1d0)), intent(in)  addimpervious,
real(kind(1d0)), intent(in)  vpd_hPa,
real(kind(1d0)), intent(in)  qn_e,
real(kind(1d0)), intent(in)  s_hPa,
real(kind(1d0)), intent(in)  ResistSurf,
real(kind(1d0)), intent(in)  RA,
real(kind(1d0)), intent(in)  rb,
real(kind(1d0)), intent(in)  tlv,
real(kind(1d0)), intent(in)  snowdensmin,
real(kind(1d0)), dimension(0:23, 2), intent(in)  SnowProf_24hr,
real(kind(1d0)), intent(in)  precip,
real(kind(1d0)), intent(in)  PipeCapacity,
real(kind(1d0)), intent(in)  RunoffToWater,
real(kind(1d0)), intent(in)  addVeg,
real(kind(1d0)), intent(in)  SnowLimPaved,
real(kind(1d0)), intent(in)  SnowLimBldg,
real(kind(1d0)), intent(in)  FlowChange,
real(kind(1d0)), dimension(nsurf), intent(in)  drain,
real(kind(1d0)), dimension(nsurf), intent(in)  WetThresh,
real(kind(1d0)), dimension(nsurf), intent(in)  stateOld,
real(kind(1d0)), dimension(nsurf), intent(in)  mw_ind,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap,
real(kind(1d0)), dimension(nsurf), intent(in)  rainonsnow,
real(kind(1d0)), dimension(nsurf), intent(in)  freezmelt,
real(kind(1d0)), dimension(nsurf), intent(in)  freezstate,
real(kind(1d0)), dimension(nsurf), intent(in)  freezstatevol,
real(kind(1d0)), dimension(nsurf), intent(in)  Qm_Melt,
real(kind(1d0)), dimension(nsurf), intent(in)  Qm_rain,
real(kind(1d0)), dimension(nsurf), intent(in)  Tsurf_ind,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
integer, dimension(3), intent(in)  dayofWeek_id,
real(kind(1d0)), dimension(6, nsurf), intent(in)  StoreDrainPrm,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowPackLimit,
real(kind(1d0)), dimension(nsurf), intent(in)  AddWater,
real(kind(1d0)), dimension(nsurf), intent(in)  addwaterrunoff,
real(kind(1d0)), dimension(nsurf), intent(inout)  soilstore_id,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowPack,
real(kind(1d0)), dimension(2), intent(inout)  SurplusEvap,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowFrac,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowWater,
real(kind(1d0)), dimension(nsurf), intent(inout)  iceFrac,
real(kind(1d0)), dimension(nsurf), intent(inout)  SnowDens,
real(kind(1d0)), intent(inout)  runoffAGimpervious,
real(kind(1d0)), intent(inout)  runoffAGveg,
real(kind(1d0)), intent(inout)  surplusWaterBody,
real(kind(1d0)), dimension(nsurf), intent(out)  rss_nsurf,
real(kind(1d0)), dimension(nsurf), intent(out)  runoffSnow,
real(kind(1d0)), dimension(nsurf), intent(out)  runoff,
real(kind(1d0)), dimension(nsurf), intent(out)  runoffSoil,
real(kind(1d0)), dimension(nsurf), intent(out)  chang,
real(kind(1d0)), dimension(nsurf), intent(out)  changSnow,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowToSurf,
real(kind(1d0)), dimension(nsurf), intent(out)  state_id,
real(kind(1d0)), dimension(nsurf), intent(out)  ev_snow,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowDepth,
real(kind(1d0)), dimension(2), intent(out)  SnowRemoval,
real(kind(1d0)), intent(out)  swe,
real(kind(1d0)), intent(out)  ev,
real(kind(1d0)), intent(out)  chSnow_tot,
real(kind(1d0)), intent(out)  ev_tot,
real(kind(1d0)), intent(out)  qe_tot,
real(kind(1d0)), intent(out)  runoff_tot,
real(kind(1d0)), intent(out)  surf_chang_tot,
real(kind(1d0)), intent(out)  runoffPipes,
real(kind(1d0)), intent(out)  mwstore,
real(kind(1d0)), intent(out)  runoffwaterbody 
)

Definition at line 497 of file suews_phys_snow.f95.

References allocatearray::bldgsurf, evap_module::cal_evap(), evap_suews_snow(), allocatearray::nsurf, allocatearray::pavsurf, snowdepletioncurve(), snowrem(), waterdist_module::updateflood(), and allocatearray::watersurf.

Referenced by suews_driver::suews_cal_qe().

497 
498  !Calculation of snow and water balance on 5 min timestep. Treats snowfree and snow covered
499  !areas separately. Weighting is taken into account in the overall values.
500  !Last modified:
501  ! LJ in 6 May 2015 - Modified to run with timestep
502  ! HCW 06 Mar 2015 - Unused variable 'i' removed.
503  ! HCW 26 Jan 2015 - Added weekday/weekend option for snow clearing profiles
504  ! LJ in 24 May 2013
505  !========================================================================
506  USE waterdist_module, ONLY: updateflood
507 
508  IMPLICIT NONE
509  ! INTEGER, PARAMETER::nsurf = 7! number of surface types
510  ! INTEGER, PARAMETER::PavSurf = 1 !New surface classes: Grass = 5th/7 surfaces
511  ! INTEGER, PARAMETER::BldgSurf = 2 !New surface classes: Grass = 5th/7 surfaces
512  ! INTEGER, PARAMETER::ConifSurf = 3 !New surface classes: Grass = 5th/7 surfaces
513  ! INTEGER,PARAMETER::DecidSurf = 4 !New surface classes: Grass = 5th/7 surfaces
514  ! INTEGER,PARAMETER::GrassSurf = 5
515  ! INTEGER, PARAMETER::BSoilSurf = 6!New surface classes: Grass = 5th/7 surfaces
516  ! INTEGER, PARAMETER::WaterSurf = 7
517 
518  INTEGER, PARAMETER::snowfractionchoice = 2 ! this PARAMETER is used all through the model
519  REAL(KIND(1d0)), PARAMETER::waterdens = 999.8395 !Density of water in 0 cel deg
520 
521  ! INTEGER,INTENT(in)::id
522  ! INTEGER,INTENT(in)::nsurf
523  INTEGER, INTENT(in)::tstep
524  INTEGER, INTENT(in)::imin
525  INTEGER, INTENT(in)::it
526  INTEGER, INTENT(in)::is
527 
528  ! INTEGER,INTENT(in)::ConifSurf
529  ! INTEGER,INTENT(in)::BSoilSurf
530  ! INTEGER,INTENT(in)::BldgSurf
531  ! INTEGER,INTENT(in)::PavSurf
532  ! INTEGER,INTENT(in)::WaterSurf
533  INTEGER, INTENT(in)::evapmethod!Evaporation calculated according to Rutter (1) or Shuttleworth (2)
534  INTEGER, DIMENSION(3), INTENT(in) ::dayofweek_id
535 
536  REAL(KIND(1d0)), INTENT(in)::dectime
537  REAL(KIND(1d0)), INTENT(in)::crwmin
538  REAL(KIND(1d0)), INTENT(in)::crwmax
539  REAL(KIND(1d0)), INTENT(in)::nsh_real
540  REAL(KIND(1d0)), INTENT(in)::lvs_j_kg
541  ! REAL(KIND(1d0)), INTENT(in)::lv_j_kg
542  REAL(KIND(1d0)), INTENT(in)::avdens
543  REAL(KIND(1d0)), INTENT(in)::vpd_hpa ! vapour pressure deficit [hPa]
544  REAL(KIND(1d0)), INTENT(in)::qn_e !net available energy for evaporation
545  REAL(KIND(1d0)), INTENT(in)::avrh
546  REAL(KIND(1d0)), INTENT(in)::press_hpa
547  REAL(KIND(1d0)), INTENT(in)::temp_c
548  REAL(KIND(1d0)), INTENT(in)::rasnow
549  REAL(KIND(1d0)), INTENT(in)::psyc_hpa
550  REAL(KIND(1d0)), INTENT(in)::avcp
551  REAL(KIND(1d0)), INTENT(in)::sice_hpa
552  REAL(KIND(1d0)), INTENT(in)::pervfraction
553  REAL(KIND(1d0)), INTENT(in)::vegfraction
554  REAL(KIND(1d0)), INTENT(in)::addimpervious
555  ! REAL(KIND(1d0)),INTENT(in)::numPM
556  REAL(KIND(1d0)), INTENT(in)::s_hpa
557  REAL(KIND(1d0)), INTENT(in)::resistsurf
558  ! REAL(KIND(1d0)),INTENT(in)::sp
559  REAL(KIND(1d0)), INTENT(in)::ra
560  REAL(KIND(1d0)), INTENT(in)::rb
561  REAL(KIND(1d0)), INTENT(in)::tlv
562  REAL(KIND(1d0)), INTENT(in)::snowdensmin
563  REAL(KIND(1d0)), INTENT(in)::precip
564  REAL(KIND(1d0)), INTENT(in)::pipecapacity
565  REAL(KIND(1d0)), INTENT(in)::runofftowater
566  REAL(KIND(1d0)), INTENT(in)::addveg
567  REAL(KIND(1d0)), INTENT(in)::snowlimpaved
568  REAL(KIND(1d0)), INTENT(in)::snowlimbldg
569  REAL(KIND(1d0)), INTENT(in)::flowchange
570 
571  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::drain
572  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::wetthresh
573  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::stateold
574  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::mw_ind
575  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::soilstorecap
576  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::rainonsnow
577  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezmelt
578  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezstate
579  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezstatevol
580  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::qm_melt
581  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::qm_rain
582  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::tsurf_ind
583  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
584  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowpacklimit
585  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::addwater
586  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::addwaterrunoff
587  REAL(KIND(1d0)), DIMENSION(6, nsurf), INTENT(in)::storedrainprm
588  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::snowprof_24hr
589 
590  !Updated status: input and output
591  REAL(KIND(1d0)), INTENT(inout)::runoffagveg
592  REAL(KIND(1d0)), INTENT(inout)::runoffagimpervious
593  REAL(KIND(1d0)), INTENT(inout)::surpluswaterbody
594 
595  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::soilstore_id
596  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowpack
597  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowfrac
598  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowwater
599  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::icefrac
600  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::snowdens
601  REAL(KIND(1d0)), DIMENSION(2), INTENT(inout):: surplusevap
602 
603  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::rss_nsurf
604  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::runoffsnow !Initialize for runoff caused by snowmelting
605  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::runoff
606  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::runoffsoil
607  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::chang
608  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::changsnow
609  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowtosurf
610  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::state_id
611  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowdepth
612  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::ev_snow
613  REAL(KIND(1d0)), DIMENSION(2), INTENT(out)::snowremoval
614 
615  REAL(KIND(1d0)), INTENT(out)::swe
616  REAL(KIND(1d0)), INTENT(out)::ev
617  REAL(KIND(1d0)), INTENT(out)::ev_tot
618  REAL(KIND(1d0)), INTENT(out)::chsnow_tot
619  REAL(KIND(1d0)), INTENT(out)::qe_tot
620  REAL(KIND(1d0)), INTENT(out)::runoff_tot
621  REAL(KIND(1d0)), INTENT(out)::surf_chang_tot
622  REAL(KIND(1d0)), INTENT(out)::runoffpipes
623  REAL(KIND(1d0)), INTENT(out)::mwstore
624  REAL(KIND(1d0)), INTENT(out)::runoffwaterbody
625 
626  REAL(KIND(1d0))::qe
627 
628  ! REAL(KIND(1d0))::Evap_SUEWS_Snow
629  REAL(KIND(1d0))::meltexcess !Excess melt water that needs to leave SnowPack
630  REAL(KIND(1d0))::snowtotinit
631  REAL(KIND(1d0))::evpart
632  REAL(KIND(1d0))::runofftest
633  REAL(KIND(1d0))::snowfracfresh1 !Snow fraction for newly formed SnowPack
634  REAL(KIND(1d0))::snowfracfresh2 !Snow fraction for newly formed SnowPack from state_id only
635  REAL(KIND(1d0))::snowfracold
636  REAL(KIND(1d0))::waterholdcapfrac
637  REAL(KIND(1d0))::fwc !Water holding capacity of snow in mm
638  REAL(KIND(1d0))::tlv_sub
639  ! REAL(KIND(1d0)):: SnowDepletionCurve
640 
641  INTEGER:: iu !1=weekday OR 2=weekend
642  REAL(KIND(1d0)), PARAMETER :: ipthreshold_mmhr = 10 !Threshold for intense precipitation [mm hr-1]
643 
644  REAL(KIND(1d0)), DIMENSION(7)::capstore ! current storage capacity [mm]
645 
646  !========================================================================
647  !Initialize variables for the calculation of water storages and evaporation
648  ev_tot = 0
649  qe_tot = 0
650  runoff_tot = 0
651  surf_chang_tot = 0
652 
653  !swe = 0
654  chsnow_tot = 0
655  !runoffPipes = 0
656  !mwstore = 0
657  !runoffwaterbody = 0
658  rss_nsurf = 0
659  !state_id = 0
660  snowdepth = 0
661  ev_snow = 0
662  !SnowRemoval = 0
663  tlv_sub = 0
664 
665  ! Use weekday or weekend snow clearing profile
666  iu = 1 !Set to 1=weekday
667  IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) iu = 2 !Set to 2=weekend
668 
669  !write(*,*) is
670  runoffsnow(is) = 0 !Initialize for runoff caused by snowmelting
671  runoff(is) = 0
672  runoffsoil(is) = 0
673  chang(is) = 0
674  changsnow(is) = 0
675  runofftest = 0
676  snowtosurf(is) = 0
677  evpart = 0
678  ev = 0
679  snowfracfresh1 = 0
680  snowfracfresh2 = 0
681  snowfracold = 0
682 
683  !Initial SnowPack + meltwater in it
684  snowtotinit = snowpack(is) + snowwater(is)
685 
686  !Calculate water holding capacity (Jin et al. 1999)
687  IF (snowdens(is) >= 200) THEN
689  ELSE
690  waterholdcapfrac = crwmin + (crwmax - crwmin)*(200 - snowdens(is))/200
691  ENDIF
692 
693  !======================================================================
694  ! Calculate evaporation from SnowPack and snow free surfaces (in mm)
695  ! IF (SnowFrac(is)<1) CALL Evap_SUEWS !ev and qe for snow free surface out
696  capstore(is) = storedrainprm(6, is)
697  IF (snowfrac(is) < 1) CALL cal_evap( &
698  evapmethod, state_id(is), wetthresh(is), capstore(is), &!input
699  vpd_hpa, avdens, avcp, qn_e, s_hpa, psyc_hpa, resistsurf, ra, rb, tlv, &
700  rss_nsurf(is), ev, qe) !output
701 
702  IF (snowfrac(is) > 0) THEN
703  call evap_suews_snow(qm_melt(is), qm_rain(is), lvs_j_kg, avdens, avrh, press_hpa, temp_c, rasnow, &
704  psyc_hpa, tstep, avcp, sice_hpa, dectime, ev_snow(is), tlv_sub)
705  ENDIF
706 
707  !If not enough water for evaporation in impervious surfaces,
708  !evaporation is taken from pervious surfaces
709  IF (is > 2) THEN
710  IF (pervfraction /= 0) THEN
711  evpart = (surplusevap(pavsurf)*sfr(pavsurf) + surplusevap(bldgsurf)*sfr(bldgsurf))/pervfraction
712  ENDIF
713  ENDIF
714 
715  !============================================================================
716  !Water surface is treated separately
717  IF (is == watersurf .AND. sfr(watersurf) > 0) GO TO 606
718 
719  !The calculations are divided into 2 main parts
720  ! 1) Surface is fully covered with snow at the beginning of the time step
721  ! 2) Surface is not fully covered with snow but rather part is snow free OR
722  ! surface not orginally covered with snow, but the snow forms at the current timestep
723 
724  !1)------------------------------------------------------------------
725  ! ------------------------------------------------------------------
726  IF (snowpack(is) > 0 .AND. snowfrac(is) == 1) THEN
727 
728  ev_snow(is) = ev_snow(is) + evpart !Evaporation surplus
729 
730  !(Snowfall per interval+freezing of melt water and surface state_id) - (meltwater+evaporation from SnowPack)
731  changsnow(is) = (precip + freezmelt(is)) - (mw_ind(is) + ev_snow(is)) !Calculate change in SnowPack (in mm)
732 
733  !If rain on snow event, add this water to SnowWater
734  IF (rainonsnow(is) > 0) THEN
735  changsnow(is) = changsnow(is) - precip
736  snowwater(is) = snowwater(is) + rainonsnow(is)
737  ENDIF
738 
739  snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
740 
741  !---------If SnowPack exists after the state_id calculations
742  IF (snowpack(is) > 0) THEN
743 
744  !Add melted water to meltstore and freeze water according to freezMelt(is)
745  snowwater(is) = snowwater(is) + mw_ind(is) - freezmelt(is)
746 
747  !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
748  fwc = waterholdcapfrac*snowpack(is)
749 
750  !If FWC is exceeded, excess meltwater (MeltExcess) will leave from the SnowPack
751  IF (snowwater(is) >= fwc) THEN
752  meltexcess = 0 !Initialize the excess meltwater
753  meltexcess = snowwater(is) - fwc !Calculate the exceess water
754  snowwater(is) = fwc !Update the SnowWater to the maximum it can hold
755  runoffsnow(is) = runoffsnow(is) + meltexcess
756  ENDIF
757 
758  !At the end of the hour calculate possible snow removal
759  IF (snowprof_24hr(it, iu) == 1 .AND. is < 3 .AND. (imin == (nsh_real - 1)/nsh_real*60)) &
760  CALL snowrem( &
761  is, pavsurf, bldgsurf, nsurf, &
762  snowfrac, sfr, &
763  snowpack, snowremoval, &
765  !----------If SnowPack is negative, it melts at this timestep
766  ELSEIF (snowpack(is) < 0) THEN
767 
768  !If freezing meltwater inside this timestep, remove it from the SnowWater
769  snowwater(is) = snowwater(is) - freezmelt(is) + mw_ind(is) + snowpack(is)
770  snowpack(is) = 0.0 !Set the snow pack and snow
771  snowfracold = 1
772  snowfrac(is) = 0
773  snowdens(is) = 0
774 
775  IF (snowwater(is) < 0) THEN !Not enough water in the meltwater store,
776  ev_snow(is) = ev_snow(is) + snowwater(is) !QUESTION: evaporation from snow is decreased?
777  IF (ev_snow(is) < 0) ev_snow(is) = 0
778  changsnow(is) = changsnow(is) + snowwater(is)
779  snowwater(is) = 0
780  ELSE
781  chang(is) = snowwater(is) !Meltwater goes to surface state_id as no snow exists anymore
782  state_id(is) = state_id(is) + chang(is)
783  snowwater(is) = 0
784  ENDIF
785  ENDIF !SnowPack negative or positive
786 
787  !2)------Surface not fully covered with snow-------------------------------------------
788  ! ------------------------------------------------------------------------------------
789  ELSEIF (snowfrac(is) < 1) THEN
790 
791  !Snow calculations: SnowPack can either exist or form at the current timestep
792  IF (snowpack(is) > 0) THEN
793  ev_snow(is) = ev_snow(is) + evpart !Evaporation surplus
794 
795  !----SnowPack water balance for the whole surface area. In reality snow depth = SnowPack/SnowFrac(is)
796  !(Snowfall per interval+freezing of melt water and surface state_id) - (meltwater+evaporation from SnowPack)
797  changsnow(is) = (precip + freezmelt(is) + freezstatevol(is)) - (mw_ind(is) + ev_snow(is)) !Calculate change in SnowPack (in mm)
798 
799  !If rain on snow event, add this water to SnowWater
800  IF (rainonsnow(is) > 0) THEN
801  changsnow(is) = changsnow(is) - precip
802  snowwater(is) = snowwater(is) + rainonsnow(is)
803  ENDIF
804  snowpack(is) = snowpack(is) + changsnow(is)
805 
806  !The fraction of snow will update when:
807  !a) Surface state_id is dry but precipitation occurs =1
808  !b) There is both precipitation and all surface state_id freezes =1
809  !c) No precipitation but all state_id freezes at a single timestep =2
810  !d) Part of the surface freezes
811  IF (precip > 0 .AND. freezstate(is) == state_id(is)) THEN !both a) and b)
812  snowfracfresh1 = 1
813  ELSEIF (precip == 0 .AND. freezstate(is) > 0 .AND. freezstate(is) == state_id(is)) THEN
814  snowfracfresh1 = 1
815 
816  !snowFracFresh1=SnowDepletionCurve(is,SnowPack(is),SnowPackLimit(is))
817  !if (snowFracFresh1<0.001) snowFracFresh1=0.001
818  ELSEIF (freezstate(is) > 0 .AND. freezstate(is) < state_id(is)) THEN !This if not all water freezes
819  snowfracfresh1 = 0.95 !Now this fraction set to something close to one. Should be improved in the future at some point
820  !if (is==1)then
821  ! write(*,*) id,it,imin,SnowFrac(is),FreezState(is),state_id(is)
822  ! pause
823  !endif
824  ENDIF
825 
826  !SnowPack can also form at the current timestep (2). If this forms purely from snowfall or/and all water at surface freezes,
827  !the whole surface will be covered with snow. If there is water on ground this snowfall can immediately melt
828  !and in this case the snow fraction is not necessarily 1 but its information is saved to snowFracFresh that
829  !is taken into account in snow fraction after calculation of state_id.
830  ELSEIF (snowpack(is) == 0 .AND. tsurf_ind(is) < 0) THEN
831 
832  !The fraction of snow will get a value of 1 (ie full snow cover):
833  !Surface state_id is dry but precipitation occurs, no precipitation but all state_id freezes at a single timestep,
834  !There is both precipitation and all surface state_id freezes
835  IF ((precip > 0 .AND. state_id(is) == 0) .OR. (precip == 0 .AND. freezstate(is) == state_id(is)) .OR. &
836  (precip > 0 .AND. freezstate(is) == state_id(is))) THEN
837 
838  !ev=ev+EvPart
839  changsnow(is) = precip + freezstatevol(is)
840  snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
841 
842  snowfracfresh1 = 1
843  icefrac(is) = freezstate(is)/(freezstate(is) + precip)
844  snowdens(is) = snowdensmin
845  ENDIF
846 
847  IF (freezstate(is) > 0 .AND. freezstate(is) < state_id(is)) THEN
848 
849  changsnow(is) = precip + freezstatevol(is)
850  snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
851  snowfracfresh2 = 0.95 !Now this fraction set to something close to one. Should be improved in the future at some point
852 
853  !snowFracFresh2=SnowDepletionCurve(is,SnowPack(is),SnowPackLimit(is))
854  !if (snowFracFresh2<0.001) snowFracFresh2=0.001
855  icefrac(is) = 1
856  snowdens(is) = snowdensmin
857  !write(*,*) 2,is,id,it,imin,SnowFrac(is),FreezState(is),state_id(is),state_id(is)+Precip
858  !pause
859 
860  ENDIF
861  ENDIF
862 
863  !---------If SnowPack exists after the state_id calculations
864  IF (snowpack(is) > 0) THEN
865 
866  !Add melted water to meltstore and freeze water according to freezMelt(is)
867  snowwater(is) = snowwater(is) + mw_ind(is) - freezmelt(is)
868 
869  !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
870  fwc = waterholdcapfrac*snowpack(is)
871 
872  !If FWC is exceeded, excess meltwater (MeltExcess) will leave from the SnowPack
873  IF (snowwater(is) >= fwc) THEN
874  meltexcess = 0 !Initialize the excess meltwater
875  meltexcess = snowwater(is) - fwc !Calculate the exceess water
876  snowwater(is) = fwc !Update the SnowWater to the maximum it can hold
877 
878  !If the fraction of snow is greater than 0.8 or if the surface is is buildings,
879  !the excess water will directly go to runoff. Otherwise it will flow to the
880  !snow free area via SnowToSurf(is)
881  IF ((snowfrac(is) > 0.9 .AND. is /= bldgsurf) .OR. (is == bldgsurf)) THEN
882  runoffsnow(is) = runoffsnow(is) + meltexcess
883  ELSE
884  snowtosurf(is) = snowtosurf(is) + meltexcess*snowfrac(is)/(1 - snowfrac(is))
885  ENDIF
886  ENDIF
887 
888  !At the end of the hour calculate possible snow removal
889  IF (snowprof_24hr(it, iu) == 1 .AND. is < 3 .AND. (imin == (nsh_real - 1)/nsh_real*60)) &
890  CALL snowrem( &
891  is, pavsurf, bldgsurf, nsurf, &
892  snowfrac, sfr, &
893  snowpack, snowremoval, &
895 
896  !----------If SnowPack is negative, it melts at this timestep
897  ELSEIF (snowpack(is) < 0) THEN
898 
899  !If freezing meltwater inside this timestep, remove it from the SnowWater
900  snowwater(is) = snowwater(is) - freezmelt(is) + mw_ind(is) + snowpack(is)
901 
902  snowpack(is) = 0.0 !Set the snow pack and snow
903  snowfracfresh1 = 0
904  snowfracfresh2 = 0
905  snowdens(is) = 0
906 
907  IF (snowwater(is) < 0) THEN !Not enough water in the meltwater store,
908  ev_snow(is) = ev_snow(is) + snowwater(is) !QUESTION: evaporation from snow is decreased.?
909  IF (ev_snow(is) < 0) ev_snow(is) = 0
910  changsnow(is) = changsnow(is) + snowwater(is)
911  snowwater(is) = 0
912  ELSE
913  snowtosurf(is) = snowtosurf(is) + snowwater(is)*snowfrac(is)/(1 - snowfrac(is))
914  snowwater(is) = 0
915  ENDIF
916  ENDIF !SnowPack negative or positive
917 
918  !--------
919  !Next the snow free surface (3). Calculations only done if snowfraction is smaller than 1
920  IF ((is == pavsurf .OR. is == bldgsurf) .AND. snowfrac(is) < 1) THEN !Impervious surfaces (paved, buildings)
921 
922  !Surface store update. If precipitation is greater than the threshold, the exceeding water
923  !goes directly to runoff
924  IF (precip > ipthreshold_mmhr/nsh_real) THEN
925  !runoff = runoff + (precipitation+water from the snow surface+water from other surfaces-the thereshold limit)
926  runoff(is) = runoff(is) + (precip + snowtosurf(is) + addwater(is) - ipthreshold_mmhr/nsh_real)
927  chang(is) = ipthreshold_mmhr/nsh_real - (drain(is) + ev + freezstate(is))
928  ELSE
929  !Add precip and water from other surfaces and remove drainage, evap and freezing of state_id
930  chang(is) = precip + snowtosurf(is) + addwater(is) - (drain(is) + ev + freezstate(is))
931  ENDIF
932 
933  state_id(is) = state_id(is) + chang(is) !Change in state_id (for whole surface area areasfr(is))
934 
935  !Add water from impervious grids
936  ! Check sfr/=0 added HCW 08 Dec 2015
937  IF (is == pavsurf .AND. sfr(pavsurf) > 0) state_id(is) = state_id(is) + (addimpervious)/sfr(pavsurf)
938 
939  runoff(is) = runoff(is) + drain(is)*addwaterrunoff(is) !Drainage (not flowing to other surfaces) goes to runoff
940 
941  IF (state_id(is) < 0.0) THEN !Surface state_id cannot be negative
942  surplusevap(is) = abs(state_id(is)) !take evaporation from other surfaces in mm
943  ev = ev - surplusevap(is)
944  state_id(is) = 0.0
945  ENDIF
946 
947  ELSEIF (is >= 3 .AND. snowfrac(is) < 1) THEN ! Pervious surfaces (conif, decid, grass unirr, grass irr)
948 
949  ev = ev + evpart
950 
951  !Change in water stores
952  IF (vegfraction > 0) THEN
953  IF (precip + addveg*(sfr(is)/vegfraction) > (ipthreshold_mmhr/nsh_real)) THEN !if 5min precipitation is larger than 10 mm
954  runoff(is) = runoff(is) + (precip + addveg*(sfr(is)/vegfraction) + &
955  snowtosurf(is) + addwater(is) - (ipthreshold_mmhr/nsh_real))
956  chang(is) = (ipthreshold_mmhr/nsh_real) - (drain(is) + ev + freezstate(is))
957  ELSE
958  chang(is) = precip + addveg*(sfr(is)/vegfraction) + snowtosurf(is) + &
959  addwater(is) - (drain(is) + ev + freezstate(is))
960  ENDIF
961  ELSE
962  chang(is) = precip + snowtosurf(is) + addwater(is) - (drain(is) + ev + freezstate(is))
963  END IF
964 
965  state_id(is) = state_id(is) + chang(is)
966 
967  !Add water in soil store only if ground is not frozen
968  IF (temp_c > 0) THEN
969  soilstore_id(is) = soilstore_id(is) + drain(is)*addwaterrunoff(is)*(1 - snowfrac(is))
970  ELSE
971  runoff(is) = runoff(is) + drain(is)*addwaterrunoff(is)
972  ENDIF
973 
974  !If state_id of the surface is negative, remove water from soilstore
975  IF (state_id(is) < 0.0) THEN
976 
977  IF ((soilstore_id(is) + state_id(is)) >= 0 .AND. temp_c > 0) THEN !If water in soilstore, water is removed
978 
979  soilstore_id(is) = soilstore_id(is) + state_id(is)*(1 - snowfrac(is))
980  state_id(is) = 0.0
981 
982  ELSE !If not water in the soilstore evaporation does not occur
983  chang(is) = chang(is) + state_id(is)
984  ev = ev + state_id(is)
985  state_id(is) = 0.0
986  ENDIF
987  ENDIF !state_id is negative
988 
989  !If soilstorage is full at this point, excess will go to surface runoff
990  IF (soilstore_id(is) > soilstorecap(is)) THEN
991  runofftest = runofftest + (soilstore_id(is) - soilstorecap(is))
992  soilstore_id(is) = soilstorecap(is)
993  ELSEIF (soilstore_id(is) < 0) THEN
994  soilstore_id(is) = 0
995  ENDIF
996 
997  ENDIF !Surface type
998 
999  ENDIF !Surface fraction
1000 
1001  !-------------------------------------------------------------------------------------------------------------------
1002 
1003  !Calculate change in SnowPack and state_id for the respective surface areas
1004  !Here the case where not all surface state_id freezes is handled
1005  IF (snowfracfresh2 > 0) THEN
1006  surf_chang_tot = (state_id(is) - stateold(is))*sfr(is)*(1 - snowfrac(is)) - precip*sfr(is)*(1 - snowfracfresh2)
1007  chsnow_tot = ((snowpack(is) + snowwater(is)) - snowtotinit)*sfr(is)*(1 - snowfrac(is)) - precip*sfr(is)*snowfracfresh2
1008  ELSE
1009  surf_chang_tot = (state_id(is) - stateold(is))*sfr(is)*(1 - snowfrac(is))
1010  chsnow_tot = ((snowpack(is) + snowwater(is)) - snowtotinit)*sfr(is)*max(snowfrac(is), snowfracold)
1011  ENDIF
1012 
1013  !Add evaporation to total
1014  IF (is == bldgsurf .OR. is == pavsurf) THEN
1015  ev_tot = ev*sfr(is)*(1 - snowfrac(is)) + ev_snow(is)*sfr(is)*max(snowfrac(is), snowfracold)
1016  qe_tot = ev_snow(is)*tlv_sub*sfr(is)*snowfrac(is) + ev*tlv*sfr(is)*(1 - snowfrac(is))
1017  ELSE
1018  ev_tot = ev*sfr(is)*(1 - snowfrac(is)) + ev_snow(is)*sfr(is)*max(snowfrac(is), snowfracold)
1019  qe_tot = ev_snow(is)*tlv_sub*sfr(is)*max(snowfrac(is), snowfracold) + ev*tlv*sfr(is)*(1 - snowfrac(is))
1020  ENDIF
1021 
1022  !========RUNOFF=======================
1023 
1024  !Add runoff to pipes
1025  runoffpipes = runoffpipes + runoffsnow(is)*sfr(is)*max(snowfrac(is), snowfracold) + runoff(is)*sfr(is)*(1 - snowfrac(is)) &
1026  + runofftest*sfr(is)
1027  CALL updateflood( &
1028  is, runoff, &! input:
1029  sfr, pipecapacity, runofftowater, &
1030  runoffagimpervious, surpluswaterbody, runoffagveg, runoffpipes)! inout:
1031 
1032  runoff_tot = runoffsnow(is)*sfr(is)*max(snowfrac(is), snowfracold) + runoff(is)*sfr(is)*(1 - snowfrac(is)) &
1033  + runofftest*sfr(is)
1034 
1035  !===Update snow depth, weighted SWE, and Mwstore
1036  IF (snowdens(is) /= 0) THEN
1037  snowdepth(is) = snowpack(is)*waterdens/snowdens(is)
1038  ENDIF
1039 
1040  ! Calculate overall snow water equivalent
1041  swe = swe + snowpack(is)*sfr(is)*max(snowfrac(is), snowfracold)
1042  mwstore = mwstore + snowwater(is)*sfr(is)*max(snowfrac(is), snowfracold)
1043 
1044  !if (id==6.and.it==13.and.imin==20) then!
1045  !if (id==85.and.it==3.and.imin==10) then!
1046  ! if (id==92.and.it==21.and.imin==35) then!
1047  ! write(*,*) ((SnowPack(is)+SnowWater(is))-snowTotInit)*sfr(is)*(1-SnowFrac(is)),&
1048  ! runoff(is)*sfr(is)*(1-SnowFrac(is)),&
1049  ! ev*sfr(is)*(1-SnowFrac(is)),&
1050  ! (state_id(is)-stateOld(is))*sfr(is)*(1-SnowFrac(is)),Precip*sfr(is)
1051  ! write(*,*) changSnow(is),runoff(is),ev,chang(is),runoffTest,FreezState(is) !changSnow(is)-freezMelt(is)
1052  ! write(*,*) is,Precip,runoff_per_tstep,ev_per_tstep,surf_chang_per_tstep,chSnow_per_interval
1053  ! write(*,*) is,Precip-runoff_per_tstep-ev_per_tstep,surf_chang_per_tstep+chSnow_per_interval
1054  ! write(*,*) is,SnowFrac(is),sfr(is),sfr(is)*ev_snow(is)
1055  ! pause
1056  ! endif
1057 
1058  !Only now update the new snow fractions both in the case that snow existing already on ground
1059  !and snow forms at the current timestep
1060  IF (snowfracfresh1 > 0) snowfrac(is) = snowfracfresh1
1061  IF (snowfracfresh2 > 0) snowfrac(is) = snowfracfresh2
1062 
1063  !Calculate new snow fraction here.
1064  !Tässä ongelmana että snow fraction muuttuu vain kun on sulamisvettä ja on vika tunti.
1065  !Tämä ei juuri koskaan toteudu johtuen lämpötilan vuorokausisyklistä
1066  !Kokeile tässä ajaa kahdella tavalla 1) ei tarvita Mw:tä
1067  ! 2) päivitys voi tapahtua millon vain
1068  !if (SnowFractionChoice==2.and.imin==(nsh_real-1)/nsh_real*60) then
1069  IF (snowfractionchoice == 2) THEN
1070  IF (snowpack(is) > 0 .AND. mw_ind(is) > 0) THEN
1071  snowfrac(is) = snowdepletioncurve(is, snowpack(is), snowpacklimit(is))
1072  IF (snowfrac(is) < 0.001) snowfrac(is) = 0.001 !The snow fraction minimum is 1% of the surface
1073  ELSEIF (snowpack(is) == 0) THEN
1074  snowfrac(is) = 0
1075  ENDIF
1076  ENDIF
1077 
1078  RETURN
1079 
1080  !==========================================================================
1081  !WATERBODY is treated separately as state_id always below ice if ice existing
1082  !Calculate change in SnowPack
1083 606 changsnow(watersurf) = (precip + freezmelt(watersurf) + freezstate(watersurf)) - &
1084  (mw_ind(watersurf) + ev_snow(watersurf))
1085 
1086  snowpack(watersurf) = snowpack(watersurf) + changsnow(watersurf) !Update SnowPack
1087  state_id(watersurf) = state_id(watersurf) + flowchange - freezstate(watersurf) !Update state_id below ice
1088 
1089  !If SnowPack exists
1090  IF (snowpack(watersurf) > 0) THEN
1091 
1092  !Add melted water to meltstore and freeze water according to freezMelt(is)
1093  snowwater(watersurf) = snowwater(watersurf) + mw_ind(watersurf) - freezmelt(watersurf)
1094 
1095  !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
1096  fwc = waterholdcapfrac*snowpack(watersurf)
1097 
1098  !If FWC is exceeded, add meltwater to state_id
1099  IF (snowwater(watersurf) >= fwc .AND. temp_c >= 0) THEN
1100  state_id(watersurf) = state_id(watersurf) + (snowwater(watersurf) - fwc)
1101  snowwater(watersurf) = fwc
1102  ENDIF
1103 
1104  !If SnowPack is negative, it melts at this timestep
1105  ELSEIF (snowpack(is) < 0) THEN
1106 
1107  !Add water to the meltwater store
1108  !If freezing meltwater inside this hour, remove it from the SnowWater
1109  snowwater(watersurf) = snowwater(watersurf) - freezmelt(watersurf) &
1110  + mw_ind(watersurf)
1111 
1112  state_id(watersurf) = state_id(watersurf) + snowwater(watersurf) + snowpack(watersurf) !Add meltwater to state_id
1113  snowpack(watersurf) = 0
1114  IF (state_id(watersurf) < 0) ev_snow(watersurf) = ev_snow(watersurf) + state_id(watersurf)
1115 
1116  ENDIF !SnowPack negative or positive
1117 
1118  !Check water state_id separately
1119  IF (state_id(watersurf) > storedrainprm(5, watersurf)) THEN
1120  runoff(watersurf) = runoff(watersurf) + (state_id(watersurf) - storedrainprm(5, watersurf))
1121  state_id(watersurf) = storedrainprm(5, watersurf)
1122  runoffwaterbody = runoffwaterbody + runoff(watersurf)*sfr(watersurf)
1123  ELSE
1124  state_id(watersurf) = state_id(watersurf) + surpluswaterbody
1125 
1126  IF (state_id(watersurf) > storedrainprm(5, watersurf)) THEN
1127  runoffwaterbody = runoffwaterbody + (state_id(watersurf) - storedrainprm(5, watersurf))*sfr(watersurf)
1128  state_id(watersurf) = storedrainprm(5, watersurf)
1129  ENDIF
1130  ENDIF
1131 
1132  !Change state_id of snow and surface
1133  chsnow_tot = ((snowpack(watersurf) + snowwater(watersurf)) - snowtotinit)*sfr(watersurf)
1134  !ch_per_interval=ch_per_interval+(state_id(WaterSurf)-stateOld(WaterSurf))*sfr(WaterSurf)
1135  surf_chang_tot = (state_id(watersurf) - stateold(watersurf))*sfr(watersurf)
1136 
1137  !Evaporation
1138  ev_tot = ev*sfr(watersurf) + ev_snow(watersurf)*sfr(watersurf)
1139  qe_tot = ev_snow(watersurf)*tlv_sub*sfr(watersurf) + ev*tlv*sfr(watersurf)
1140  runoff_tot = runoff(is) !The total runoff from the area
1141 
1142  IF (snowpack(watersurf) > 0) THEN !Fraction only 1 or 0
1143  snowfrac(watersurf) = 1
1144  ELSE
1145  snowfrac(watersurf) = 0
1146  ENDIF
1147 
real(kind(1d0)) tlv
real(kind(1d0)) waterdens
integer snowfractionchoice
real(kind(1d0)) mwstore
real(kind(1d0)) snowlimpaved
real(kind(1d0)) lvs_j_kg
real(kind(1d0)) snowdensmin
real(kind(1d0)) avcp
subroutine updateflood(is, runoff, sfr, PipeCapacity, RunoffToWater, runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes)
integer imin
real(kind(1d0)), dimension(2) snowremoval
real(kind(1d0)), dimension(0:23, 2) snowprof_24hr
real(kind(1d0)) psyc_hpa
real(kind(1d0)) vpd_hpa
real(kind(1d0)) snowlimbldg
real(kind(1d0)) sice_hpa
real(kind(1d0)) s_hpa
real(kind(1d0)) crwmin
integer it
real(kind(1d0)) dectime
real(kind(1d0)) avdens
real(kind(1d0)) waterholdcapfrac
real(kind(1d0)) crwmax
real(kind(1d0)) swe
Here is the call graph for this function:
Here is the caller graph for this function:

◆ snowdepletioncurve()

real(kind(1d0)) function snow_module::snowdepletioncurve ( integer  is,
real(kind(1d0))  swe,
real(kind(1d0))  sweD 
)

Definition at line 1244 of file suews_phys_snow.f95.

References allocatearray::bldgsurf, allocatearray::pavsurf, and allocatearray::watersurf.

Referenced by snowcalc().

1244  !This function calculates surface coverage of snow according to the
1245  !depletion curves in Valeo and Ho (2004).
1246  !INPUT: is Surface type number
1247  ! swe Snow water content
1248  ! sweD Limit for
1249 
1250  USE allocatearray
1251 
1252  IMPLICIT NONE
1253 
1254  INTEGER::is
1255  REAL(KIND(1d0))::asc, swed, swe
1256 
1257  ! initialisation
1258  asc = 1
1259  !Impervious surface
1260  IF (is == pavsurf) THEN
1261 
1262  IF (swe <= swed) THEN !Snow water equivalent below threshold
1263  asc = ((swe/swed))**2
1264  ELSE
1265  asc = 1
1266  ENDIF
1267 
1268  !Bldgs surface
1269  ELSEIF (is == bldgsurf) THEN
1270 
1271  IF (swe <= swed) THEN
1272  IF ((swe/swed) < 0.9) THEN
1273  asc = (swe/swed)*0.5
1274  ELSE
1275  asc = (swe/swed)**8
1276  ENDIF
1277  ELSE
1278  asc = 1
1279  ENDIF
1280  ELSEIF (is == watersurf) THEN
1281  IF (swe > 0) asc = 1
1282 
1283  !Vegetion surfaces
1284  ELSE
1285  IF (swe <= swed) THEN
1286 
1287  asc = 1 - ((1/3.1416)*acos(2*(swe/swed) - 1))**1.7
1288  ELSE
1289  asc = 1
1290  ENDIF
1291 
1292  ENDIF
1293 
1294  !asc=real(int(10000.*asc))/10000 !4 decimal precision
1295 
1296  RETURN
integer, parameter pavsurf
integer, parameter bldgsurf
integer, parameter watersurf
real(kind(1d0)) swe
Here is the caller graph for this function:

◆ snowrem()

subroutine snow_module::snowrem ( integer, intent(in)  is,
integer, intent(in)  PavSurf,
integer, intent(in)  BldgSurf,
integer, intent(in)  nsurf,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowFrac,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowPack,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowRemoval,
real(kind(1d0)), intent(in)  SnowLimPaved,
real(kind(1d0)), intent(in)  SnowLimBldg 
)

Definition at line 1215 of file suews_phys_snow.f95.

Referenced by snowcalc().

1215 
1216  IMPLICIT NONE
1217  INTEGER, INTENT(in) :: is, pavsurf, bldgsurf, nsurf
1218  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) :: snowfrac, sfr
1219  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: snowpack, snowremoval
1220  REAL(KIND(1d0)), INTENT(in) :: snowlimpaved, snowlimbldg
1221  !write(*,*) is, SnowPack(is),SnowLimPaved,SnowLimBldg
1222 
1223  IF (is == pavsurf) THEN
1224  IF (snowpack(pavsurf) > snowlimpaved) THEN
1227  !SnowPack(PavSurf)=SnowPack(PavSurf)/SnowFrac(PavSurf)
1228  ENDIF
1229  ENDIF
1230  IF (is == bldgsurf) THEN
1231  IF (snowpack(bldgsurf) > snowlimbldg) THEN
1234  !SnowPack(BldgSurf)=SnowPack(BldgSurf)/SnowFrac(BldgSurf)
1235  ENDIF
1236  ENDIF
1237  !write(*,*) is, SnowPack(is),SnowLimPaved,SnowLimBldg
1238  !pause
real(kind(1d0)) snowlimpaved
integer, parameter nsurf
real(kind(1d0)), dimension(2) snowremoval
real(kind(1d0)), dimension(nsurf) snowfrac
real(kind(1d0)), dimension(nsurf) snowpack
real(kind(1d0)), dimension(nsurf) sfr
real(kind(1d0)) snowlimbldg
integer, parameter pavsurf
integer, parameter bldgsurf
Here is the caller graph for this function:

◆ snowupdate()

subroutine snow_module::snowupdate ( integer, intent(in)  tstep,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  tau_a,
real(kind(1d0)), intent(in)  tau_f,
real(kind(1d0)), intent(in)  tau_r,
real(kind(1d0)), intent(in)  SnowDensMax,
real(kind(1d0)), intent(in)  SnowDensMin,
real(kind(1d0)), intent(in)  SnowAlbMax,
real(kind(1d0)), intent(in)  SnowAlbMin,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowPack_prev,
real(kind(1d0)), intent(in)  SnowAlb_prev,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowDens_prev,
real(kind(1d0)), intent(out)  SnowAlb_next,
real(kind(1d0)), dimension(nsurf), intent(out)  SnowDens_next 
)

Definition at line 1337 of file suews_phys_snow.f95.

References update_snow_albedo(), and update_snow_dens().

1337 
1338  IMPLICIT NONE
1339 
1340  ! INTEGER, INTENT(in)::nsurf
1341  INTEGER, INTENT(in)::tstep
1342 
1343  REAL(KIND(1D0)), INTENT(in)::temp_c !Air temperature
1344  REAL(KIND(1D0)), INTENT(in)::tau_a
1345  REAL(KIND(1D0)), INTENT(in)::tau_f
1346  REAL(KIND(1D0)), INTENT(in)::tau_r
1347  REAL(KIND(1D0)), INTENT(in)::snowdensmax
1348  REAL(KIND(1D0)), INTENT(in)::snowdensmin
1349  REAL(KIND(1D0)), INTENT(in)::snowalbmax
1350  REAL(KIND(1D0)), INTENT(in)::snowalbmin
1351 
1352  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowpack_prev
1353 
1354  REAL(KIND(1d0)), INTENT(in)::snowalb_prev
1355  REAL(KIND(1d0)), INTENT(out)::snowalb_next
1356 
1357  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowdens_prev
1358  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowdens_next
1359 
1360  ! INTEGER::is
1361  REAL(KIND(1D0))::alb_change !Change in snow albedo
1362  REAL(KIND(1D0))::dens_change !Change in snow density
1363  REAL(KIND(1D0)), parameter:: tau_1 = 24*60*60 !Number of seconds in a day
1364 
1365  !Initialize
1366  alb_change = 0
1367  dens_change = 0
1368 
1369  !==========================================================
1370  !Calculation of snow albedo by Lemonsu et al. 2010
1371  !(org: Verseghy (1991)&Baker et al.(1990))
1372  ! IF (SUM(SnowPack_prev) > 0) THEN !Check if snow on any of the surfaces
1373  ! IF (Temp_C < 0) THEN
1374  ! !alb_change = tau_a*(60*60)/tau_1
1375  ! alb_change = tau_a*(tstep)/tau_1
1376  ! SnowAlb_next = SnowAlb_prev - alb_change
1377  ! ELSE
1378  ! !alb_change = exp(-tau_f*(60*60)/tau_1)
1379  ! alb_change = EXP(-tau_f*(tstep)/tau_1)
1380  ! SnowAlb_next = (SnowAlb_prev - SnowAlbMin)*alb_change + SnowAlbMin
1381  ! ENDIF
1382  ! IF (SnowAlb_next < SnowAlbMin) SnowAlb_next = SnowAlbMin !Albedo cannot be smaller than the min albedo
1383  ! IF (SnowAlb_next > SnowAlbMax) SnowAlb_next = SnowAlbMax !Albedo cannot be larger than the max albedo
1384  ! if (SnowAlb_next < 0) print *, 'SnowAlbMin/max in SnowUpdate', SnowAlbMin, SnowAlbMax, SnowAlb_next
1385  ! ELSE
1386  ! SnowAlb_next = 0
1387  ! ENDIF
1388  ! if (SnowAlb_next < 0) print *, 'SnowAlb in SnowUpdate', SnowAlb_next
1389  snowalb_next = update_snow_albedo( &
1390  tstep, snowpack_prev, snowalb_prev, temp_c, &
1392 
1393  !Update snow density: There is a mistake in Järvi et al. (2014): tau_h should be tau_1
1394  ! DO is = 1, nsurf
1395  ! !If SnowPack existing
1396  ! IF (SnowPack_prev(is) > 0) THEN
1397  ! dens_change = EXP(-tau_r*(tstep)/tau_1)
1398  ! IF (SnowPack_prev(is) > 0) SnowDens_next(is) = (SnowDens_prev(is) - SnowDensMax)*dens_change + SnowDensMax
1399  ! IF (SnowDens_next(is) > SnowDensMax) SnowDens_next(is) = SnowDensMax
1400  ! ELSE
1401  ! SnowDens_next(is) = SnowDensMin
1402  ! ENDIF
1403  ! ENDDO
1404 
1405  snowdens_next = update_snow_dens( &
1406  tstep, snowpack_prev, snowdens_prev, &
1408 
real(kind(1d0)) tau_r
real(kind(1d0)) snowdensmin
real(kind(1d0)) tau_f
real(kind(1d0)) snowalbmax
real(kind(1d0)) tau_a
real(kind(1d0)) snowdensmax
real(kind(1d0)) snowalbmin
Here is the call graph for this function:

◆ update_snow_albedo()

real(kind(1d0)) function snow_module::update_snow_albedo ( integer, intent(in)  tstep,
real(kind(1d0)), dimension(7), intent(in)  SnowPack_prev,
real(kind(1d0)), intent(in)  SnowAlb_prev,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  tau_a,
real(kind(1d0)), intent(in)  tau_f,
real(kind(1d0)), intent(in)  SnowAlbMax,
real(kind(1d0)), intent(in)  SnowAlbMin 
)

Definition at line 1415 of file suews_phys_snow.f95.

Referenced by snowupdate(), and suews_driver::suews_cal_qn().

1415  implicit none
1416  INTEGER, INTENT(in)::tstep
1417 
1418  REAL(KIND(1d0)), DIMENSION(7), INTENT(in)::snowpack_prev
1419  REAL(KIND(1D0)), INTENT(in) ::temp_c
1420  REAL(KIND(1d0)), INTENT(in)::snowalb_prev
1421  REAL(KIND(1D0)), INTENT(in)::tau_a
1422  REAL(KIND(1D0)), INTENT(in)::tau_f
1423  REAL(KIND(1D0)), INTENT(in)::snowalbmax
1424  REAL(KIND(1D0)), INTENT(in)::snowalbmin
1425 
1426  REAL(KIND(1D0)) :: snowalb_next
1427 
1428  REAL(KIND(1D0))::alb_change !Change in snow albedo
1429  REAL(KIND(1D0)), parameter:: tau_1 = 24*60*60 !Number of seconds in a day
1430 
1431  !==========================================================
1432  !Calculation of snow albedo by Lemonsu et al. 2010
1433  !(org: Verseghy (1991)&Baker et al.(1990))
1434  IF (sum(snowpack_prev) > 0) THEN !Check if snow on any of the surfaces
1435  IF (temp_c < 0) THEN
1436  !alb_change = tau_a*(60*60)/tau_1
1437  alb_change = tau_a*(tstep)/tau_1
1438  snowalb_next = snowalb_prev - alb_change
1439  ELSE
1440  !alb_change = exp(-tau_f*(60*60)/tau_1)
1441  alb_change = exp(-tau_f*(tstep)/tau_1)
1442  snowalb_next = (snowalb_prev - snowalbmin)*alb_change + snowalbmin
1443  ENDIF
1444  IF (snowalb_next < snowalbmin) snowalb_next = snowalbmin !Albedo cannot be smaller than the min albedo
1445  IF (snowalb_next > snowalbmax) snowalb_next = snowalbmax !Albedo cannot be larger than the max albedo
1446  if (snowalb_next < 0) print *, 'SnowAlbMin/max in SnowUpdate', snowalbmin, snowalbmax, snowalb_next
1447  ELSE
1448  snowalb_next = 0
1449  ENDIF
1450  if (snowalb_next < 0) print *, 'SnowAlb in SnowUpdate', snowalb_next
1451 
real(kind(1d0)) tau_f
real(kind(1d0)) snowalbmax
real(kind(1d0)) tau_a
real(kind(1d0)) snowalbmin
Here is the caller graph for this function:

◆ update_snow_dens()

real(kind(1d0)) function, dimension(nsurf) snow_module::update_snow_dens ( integer, intent(in)  tstep,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowPack_prev,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowDens_prev,
real(kind(1d0)), intent(in)  tau_r,
real(kind(1d0)), intent(in)  SnowDensMax,
real(kind(1d0)), intent(in)  SnowDensMin 
)

Definition at line 1458 of file suews_phys_snow.f95.

References allocatearray::nsurf.

Referenced by snow_cal_meltheat(), and snowupdate().

1458  implicit none
1459  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowpack_prev
1460  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowdens_prev
1461  REAL(KIND(1D0)), INTENT(in)::snowdensmax
1462  REAL(KIND(1D0)), INTENT(in)::snowdensmin
1463  REAL(KIND(1D0)), INTENT(in)::tau_r
1464  INTEGER, INTENT(in)::tstep
1465 
1466  ! REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowDens_next
1467 
1468  REAL(KIND(1D0)), DIMENSION(nsurf) :: snowdens_next
1469  REAL(KIND(1D0)) :: dens_change
1470 
1471  INTEGER::is
1472 
1473  REAL(KIND(1D0)), parameter:: tau_1 = 24*60*60
1474  ! INTEGER,parameter::nsurf=7
1475 
1476  !Update snow density: There is a mistake in Järvi et al. (2014): tau_h should be tau_1
1477  DO is = 1, nsurf
1478 
1479  !If SnowPack existing
1480  IF (snowpack_prev(is) > 0) THEN
1481  dens_change = exp(-tau_r*(tstep)/tau_1)
1482  IF (snowpack_prev(is) > 0) snowdens_next(is) = (snowdens_prev(is) - snowdensmax)*dens_change + snowdensmax
1483  IF (snowdens_next(is) > snowdensmax) snowdens_next(is) = snowdensmax
1484  ELSE
1485  snowdens_next(is) = snowdensmin
1486  ENDIF
1487  ENDDO
1488 
real(kind(1d0)) tau_r
integer, parameter nsurf
real(kind(1d0)) snowdensmin
real(kind(1d0)) snowdensmax
Here is the caller graph for this function:

◆ veg_fr_snow()

subroutine snow_module::veg_fr_snow ( real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), dimension(nsurf), intent(in)  SnowFrac,
real(kind(1d0)), intent(out)  veg_fr 
)
Parameters
[in]sfrsurface fractions
[in]snowfracsnowy surface fractions [-]
[out]veg_frvegetated surface fractions [-]

Definition at line 1302 of file suews_phys_snow.f95.

Referenced by snow_cal_meltheat().

1302 
1303  IMPLICIT NONE
1304 
1305  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) :: sfr
1306  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) :: snowfrac
1307 
1308  REAL(KIND(1d0)), INTENT(out) :: veg_fr
1309 
1310  veg_fr = dot_product(sfr(3:7), 1 - snowfrac(3:7))
1311 
real(kind(1d0)), dimension(nsurf) snowfrac
real(kind(1d0)), dimension(nsurf) sfr
Here is the caller graph for this function: