CanopyHydrology_SnowWater_impl.hh 7.57 KB
Newer Older
Ethan Coon's avatar
Ethan Coon committed
1
2
3
#ifndef CANOPY_HYDROLOGY_SNOWWATER_IMPL_HH_
#define CANOPY_HYDROLOGY_SNOWWATER_IMPL_HH_

4
#include <stdio.h>  
5
#include <cmath>     
6
7
8
9
10
11
#include "landunit_varcon.h"    
#include "column_varcon.h" 
#include "clm_varpar.h"         
#include "clm_varctl.h"     

using namespace std;
Pillai, Himanshu's avatar
Pillai, Himanshu committed
12
13
using std::min ;
using std::max ;
14
15
16
17

namespace ELM {

template<typename Array_d>
Pillai, Himanshu's avatar
Pillai, Himanshu committed
18
   NATURE void CanopyHydrology_SnowWater(const double& dtime,
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
    const double& qflx_floodg,
    const int& ltype,
    const int& ctype,
    const bool& urbpoi,
    const bool& do_capsnow,                            
    const int& oldfflag,
    const double& forc_air_temp,
    const double& t_grnd,
    const double& qflx_snow_grnd_col,
    const double& qflx_snow_melt,
    const double& n_melt,
    const double& frac_h2osfc,
    double& snow_depth,
    double& h2osno,
    double& integrated_snow,
    Array_d swe_old,
    Array_d h2osoi_liq,
    Array_d h2osoi_ice,
    Array_d t_soisno,
    Array_d frac_iceold,
    int& snow_level,
    Array_d dz,
    Array_d z,
    Array_d zi,
    int& newnode,
    double& qflx_floodc,
    double& qflx_snow_h2osfc,
    double& frac_sno_eff,
    double& frac_sno)
  {       


  //parameters
    double rpi=4.0e0*atan(1.0e0)  ;
    double tfrz=273.15;
    double zlnd = 0.010;


57
58
59
60
61
  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: swe_old 
  // real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: t_soisno, frac_iceold
  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: dz, z, zi

62
63
64
65
  //local variables 
    double  temp_intsnow, temp_snow_depth, z_avg, fmelt, dz_snowf, snowmelt ;
    double  newsnow, bifall, accum_factor, fsno_new, smr ;
    int j ;
66

67
68
69
70
71
72
  //apply gridcell flood water flux to non-lake columns
    if (ctype != icol_sunwall && ctype != icol_shadewall) {      
     qflx_floodc = qflx_floodg; }
     else{
       qflx_floodc = 0.0;
     }
73

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
  //Determine snow height and snow water

  //Use Alta relationship, Anderson(1976); LaChapelle(1961),
  //U.S.Department of Agriculture Forest Service, Project F,
  //Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.

     qflx_snow_h2osfc = 0.0;
  //set temporary variables prior to updating
     temp_snow_depth=snow_depth;
  //save initial snow content
  // for(j = -nlevsno+1; j < snow_level; j++) {
  //    swe_old[j] = 0.00;
  // }
  // for(j = snow_level+1; j < 0; j++) {
  //    swe_old[j]=h2osoi_liq[j]+h2osoi_ice[j];
  // }

     for(j = 0; j < nlevsno; j++) {
       if(j < nlevsno+snow_level )  swe_old[j]=0.00;
       else  swe_old[j] = h2osoi_liq[j]+h2osoi_ice[j];
     }
95

96
97
98
99
100
101
     if (do_capsnow) {
       dz_snowf = 0.;
       newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
       frac_sno=1.;
       integrated_snow = 5.e2; 
     } else {
102

103
       if (forc_air_temp > tfrz + 2.) {
Pillai, Himanshu's avatar
Pillai, Himanshu committed
104
        bifall=50. + 1.7*pow((17.0),1.5);
105
      } else if (forc_air_temp > tfrz - 15.) {
Pillai, Himanshu's avatar
Pillai, Himanshu committed
106
        bifall=50. + 1.7*pow((forc_air_temp - tfrz + 15.),1.5);
107
      } else {
108
        bifall=50.;
109
      }
110

111
112
   // newsnow is all snow that doesn't fall on h2osfc
      newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
113

114
   // update integrated_snow
Pillai, Himanshu's avatar
Pillai, Himanshu committed
115
   integrated_snow = max(integrated_snow,h2osno) ; //h2osno could be larger due to frost
116

117
118
   // snowmelt from previous time step * dtime
   snowmelt = qflx_snow_melt * dtime;
119

120
121
   // set shape factor for accumulation of snow
   accum_factor=0.1;
122

123
   if (h2osno > 0.0) {
124

125
126
127
128
      //======================  FSCA PARAMETERIZATIONS  ======================
      // fsca parameterization based on *changes* in swe
      // first compute change from melt during previous time step
    if(snowmelt > 0.) {
129

Pillai, Himanshu's avatar
Pillai, Himanshu committed
130
     smr=min(1.,(h2osno)/(integrated_snow));
131

Pillai, Himanshu's avatar
Pillai, Himanshu committed
132
     frac_sno = 1. - pow((acos(min(1.,(2.*smr - 1.)))/rpi),(n_melt)) ;
133

134
   }
135

136
137
138
139
      // update fsca by new snow event, add to previous fsca
   if (newsnow > 0.) {
     fsno_new = 1. - (1. - tanh(accum_factor*newsnow))*(1. - frac_sno);
     frac_sno = fsno_new;
140

141
         // reset integrated_snow after accumulation events
Pillai, Himanshu's avatar
Pillai, Himanshu committed
142
143
     temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*pow((1.0-max(frac_sno,1e-6)),(1.0/n_melt))+1.0))) ;
     integrated_snow = min(1.e8,temp_intsnow) ;
144
   }
145

146
      //====================================================================
147

148
149
150
151
152
153
154
155
156
157
158
      // for subgrid fluxes
   if (subgridflag ==1 && ! urbpoi) {
     if (frac_sno > 0.){
      snow_depth=snow_depth + newsnow/(bifall * frac_sno);
    } else {
      snow_depth=0.;
    }
  } else {
         // for uniform snow cover
   snow_depth=snow_depth+newsnow/bifall;
 }
159

160
161
162
163
      // use original fsca formulation (n&y 07)
 if (oldfflag == 1) { 
         // snow cover fraction in Niu et al. 2007
   if(snow_depth > 0.0)  {
Pillai, Himanshu's avatar
Pillai, Himanshu committed
164
    frac_sno = tanh(snow_depth/(2.5*zlnd*pow((min(800.0,(h2osno+ newsnow)/snow_depth)/100.0),1.0)) ) ;
165
166
  }
  if(h2osno < 1.0)  {
Pillai, Himanshu's avatar
Pillai, Himanshu committed
167
    frac_sno=min(frac_sno,h2osno);
168
169
170
171
172
173
174
175
176
177
178
179
  }
}

   } else { //h2osno == 0
      // initialize frac_sno and snow_depth when no snow present initially
    if (newsnow > 0.) { 
     z_avg = newsnow/bifall;
     fmelt=newsnow;
     frac_sno = tanh(accum_factor*newsnow);

         // make integrated_snow consistent w/ new fsno, h2osno
         integrated_snow = 0. ;//reset prior to adding newsnow below
Pillai, Himanshu's avatar
Pillai, Himanshu committed
180
181
         temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*pow((1.0-max(frac_sno,1e-6)),(1.0/n_melt)))+1.0));
         integrated_snow = min(1.e8,temp_intsnow);
182
183
184
185

         // update snow_depth and h2osno to be consistent with frac_sno, z_avg
         if (subgridflag ==1 && !urbpoi) {
          snow_depth=z_avg/frac_sno;
186
        } else {
187
          snow_depth=newsnow/bifall;
188
        }
189
190
191
192
         // use n&y07 formulation
        if (oldfflag == 1) { 
            // snow cover fraction in Niu et al. 2007
          if(snow_depth > 0.0)  {
Pillai, Himanshu's avatar
Pillai, Himanshu committed
193
           frac_sno = tanh(snow_depth/(2.5*zlnd*pow((min(800.0,newsnow/snow_depth)/100.0),1.0)) );
194
195
196
197
198
199
200
201
         }
       }
     } else {
       z_avg = 0.;
       snow_depth = 0.;
       frac_sno = 0.;
     }
   } // end of h2osno > 0
202

203
204
   // snow directly falling on surface water melts, increases h2osfc
   qflx_snow_h2osfc = frac_h2osfc*qflx_snow_grnd_col;
205

206
207
208
   // update h2osno for new snow
   h2osno = h2osno + newsnow ;
   integrated_snow = integrated_snow + newsnow;
209

210
211
   // update change in snow depth
   dz_snowf = (snow_depth - temp_snow_depth) / dtime;
212
213
214
215
216

  } //end of do_capsnow construct

  // set frac_sno_eff variable
  if (ltype == istsoil || ltype == istcrop) {
217
218
   if (subgridflag ==1) { 
    frac_sno_eff = frac_sno;
219
  } else {
220
    frac_sno_eff = 1.;
221
  }
222
223
224
225
226
227
228
229
230
231
232
233
} else {
 frac_sno_eff = 1.;
}

if (ltype==istwet && t_grnd>tfrz) {
 h2osno=0.;
 snow_depth=0.;
}

  //When the snow accumulation exceeds 10 mm, initialize snow layer
  //Currently, the water temperature for the precipitation is simply set
  //as the surface air temperature
234
  newnode = 0 ; //flag for when snow node will be initialized
235
236
237
238
239
240
241
242
243
244
  if (snow_level == 0 && qflx_snow_grnd_col > 0.00 && frac_sno*snow_depth >= 0.010) {
   newnode = 1;
   snow_level = -1;
   dz[nlevsno-1] = snow_depth ;                    //meter
   z[nlevsno-1] = -0.50*dz[nlevsno-1];
   zi[nlevsno-2] = -dz[nlevsno-1];
   t_soisno[nlevsno-1] = min(tfrz, forc_air_temp) ;   //K
   h2osoi_ice[nlevsno-1] = h2osno ;            //kg/m2
   h2osoi_liq[nlevsno-1] = 0.0  ;               //kg/m2
   frac_iceold[nlevsno-1] = 1.0;
245
246
 }

247
248
249
250
251
252
253
254
255
256
  //The change of ice partial density of surface node due to precipitation.
  //Only ice part of snowfall is added here, the liquid part will be added
  //later.
 if (snow_level < 0 && newnode == 0) {
  h2osoi_ice[nlevsno-1+snow_level+1] = h2osoi_ice[nlevsno-1+snow_level+1]+newsnow;
  dz[nlevsno-1+snow_level+1] = dz[nlevsno-1+snow_level+1]+dz_snowf*dtime;
}
}
}

Ethan Coon's avatar
Ethan Coon committed
257
258

#endif