CanopyHydrology_SnowWater_impl.hh 7.54 KB
Newer Older
1
#include <stdio.h>  
2
#include <cmath>     
3
4
5
6
7
8
9
10
11
12
13
#include "landunit_varcon.h"    
#include "column_varcon.h" 
#include "clm_varpar.h"         
#include "clm_varctl.h"     
#include "CanopyHydrology_cpp.hh"    

using namespace std;

namespace ELM {

template<typename Array_d>
14
15
16
17
18
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
  void CanopyHydrology_SnowWater(const double& dtime,
    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;


53
54
55
56
57
  // 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

58
59
60
61
  //local variables 
    double  temp_intsnow, temp_snow_depth, z_avg, fmelt, dz_snowf, snowmelt ;
    double  newsnow, bifall, accum_factor, fsno_new, smr ;
    int j ;
62

63
64
65
66
67
68
  //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;
     }
69

70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
  //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];
     }
91

92
93
94
95
96
97
     if (do_capsnow) {
       dz_snowf = 0.;
       newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
       frac_sno=1.;
       integrated_snow = 5.e2; 
     } else {
98

99
       if (forc_air_temp > tfrz + 2.) {
100
        bifall=50. + 1.7*std::pow((17.0),1.5);
101
      } else if (forc_air_temp > tfrz - 15.) {
102
        bifall=50. + 1.7*std::pow((forc_air_temp - tfrz + 15.),1.5);
103
      } else {
104
        bifall=50.;
105
      }
106

107
108
   // newsnow is all snow that doesn't fall on h2osfc
      newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
109

110
111
   // update integrated_snow
   integrated_snow = std::max(integrated_snow,h2osno) ; //h2osno could be larger due to frost
112

113
114
   // snowmelt from previous time step * dtime
   snowmelt = qflx_snow_melt * dtime;
115

116
117
   // set shape factor for accumulation of snow
   accum_factor=0.1;
118

119
   if (h2osno > 0.0) {
120

121
122
123
124
      //======================  FSCA PARAMETERIZATIONS  ======================
      // fsca parameterization based on *changes* in swe
      // first compute change from melt during previous time step
    if(snowmelt > 0.) {
125

126
     smr=std::min(1.,(h2osno)/(integrated_snow));
127

128
     frac_sno = 1. - std::pow((acos(min(1.,(2.*smr - 1.)))/rpi),(n_melt)) ;
129

130
   }
131

132
133
134
135
      // 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;
136

137
138
139
140
         // reset integrated_snow after accumulation events
     temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*std::pow((1.0-std::max(frac_sno,1e-6)),(1.0/n_melt))+1.0))) ;
     integrated_snow = std::min(1.e8,temp_intsnow) ;
   }
141

142
      //====================================================================
143

144
145
146
147
148
149
150
151
152
153
154
      // 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;
 }
155

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

   } 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
         temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*std::pow((1.0-std::max(frac_sno,1e-6)),(1.0/n_melt)))+1.0));
         integrated_snow = std::min(1.e8,temp_intsnow);

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

199
200
   // snow directly falling on surface water melts, increases h2osfc
   qflx_snow_h2osfc = frac_h2osfc*qflx_snow_grnd_col;
201

202
203
204
   // update h2osno for new snow
   h2osno = h2osno + newsnow ;
   integrated_snow = integrated_snow + newsnow;
205

206
207
   // update change in snow depth
   dz_snowf = (snow_depth - temp_snow_depth) / dtime;
208
209
210
211
212

  } //end of do_capsnow construct

  // set frac_sno_eff variable
  if (ltype == istsoil || ltype == istcrop) {
213
214
   if (subgridflag ==1) { 
    frac_sno_eff = frac_sno;
215
  } else {
216
    frac_sno_eff = 1.;
217
  }
218
219
220
221
222
223
224
225
226
227
228
229
} 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
230
  newnode = 0 ; //flag for when snow node will be initialized
231
232
233
234
235
236
237
238
239
240
  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;
241
242
 }

243
244
245
246
247
248
249
250
251
252
  //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;
}
}
}