1 | import numpy as np |
---|
2 | from matplotlib.backends.backend_pdf import PdfPages |
---|
3 | import matplotlib.pyplot as plt |
---|
4 | import matplotlib as mpl |
---|
5 | import matplotlib.cm as cm |
---|
6 | from mpl_toolkits.basemap import Basemap |
---|
7 | from matplotlib.patches import Polygon |
---|
8 | # |
---|
9 | # Compute layer thickness |
---|
10 | # |
---|
11 | def CompThickness(x): |
---|
12 | nbx=len(x) |
---|
13 | dx=0.0 |
---|
14 | for i in range(1,nbx): |
---|
15 | dx = np.append(dx,(x[i-1]+x[i])/2.0) |
---|
16 | # Last layer |
---|
17 | dx=np.append(dx,(x[nbx-1]-dx[nbx-1])*2+dx[nbx-1]) |
---|
18 | # |
---|
19 | tx=[] |
---|
20 | for i in range(len(dx)-1): |
---|
21 | tx = np.append(tx,dx[i+1]-dx[i]) |
---|
22 | # |
---|
23 | return dx,tx |
---|
24 | # |
---|
25 | # Reading booleans |
---|
26 | # |
---|
27 | def str2bool(v): |
---|
28 | return v.lower() in ("yes", "true", "t", "1") |
---|
29 | # |
---|
30 | # Ploting function |
---|
31 | # |
---|
32 | def discretizationplot(wlev, tlev, wdepth_max, tsol_max, overalltitle): |
---|
33 | # |
---|
34 | # Compute thicknesses |
---|
35 | # |
---|
36 | wy,wz=CompThickness(wlev) |
---|
37 | ty,tz=CompThickness(tlev) |
---|
38 | # |
---|
39 | x=np.array([0,1]) |
---|
40 | |
---|
41 | f, (ax1, ax2) = plt.subplots(1, 2, sharey=False) |
---|
42 | |
---|
43 | zmax=max([np.max(wz),np.max(tz)]) |
---|
44 | cNorm = mpl.colors.Normalize(vmin=0, vmax=zmax) |
---|
45 | cmap = plt.get_cmap('Dark2') |
---|
46 | |
---|
47 | ztmp=np.reshape(np.repeat(wz,2), (len(wy)-1, len(x)) ) |
---|
48 | ax1.set_yscale('symlog') |
---|
49 | ax1.yaxis.set_major_formatter(mpl.ticker.ScalarFormatter()) |
---|
50 | ax1.yaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%d')) |
---|
51 | ax1.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: str(int(round(x))))) |
---|
52 | # ax1.set_yticks([-0.02, -0.2, -2]) |
---|
53 | datamap=ax1.pcolor(x,-wy, ztmp, zorder=0, norm=cNorm, cmap=cmap) |
---|
54 | # Add tine lines to show the middle of the layer where SM is positioned |
---|
55 | for i in range(len(wlev)): |
---|
56 | ax1.plot(x,np.array([-wlev[i],-wlev[i]]), color='k', linewidth=0.4) |
---|
57 | # |
---|
58 | ax1.text(0.5, -wdepth_max, "depth_Wmax", fontweight='bold',ha='center', |
---|
59 | va='center') |
---|
60 | ax1.set_xlabel("Min thickness="+str(round(np.min(wz),5))+"[m]") |
---|
61 | ax1.set_title("nslm="+str(len(wlev))) |
---|
62 | # |
---|
63 | ztmp=np.reshape(np.repeat(tz, 2), (len(ty)-1,len(x))) |
---|
64 | datamap = ax2.pcolor(x,-ty, ztmp, zorder=0, norm=cNorm, cmap=cmap) |
---|
65 | # Add tine lines to show the middle of the layer where T is positioned |
---|
66 | for i in range(len(tlev)): |
---|
67 | ax2.plot(x,np.array([-tlev[i],-tlev[i]]), color='k', linewidth=0.4) |
---|
68 | |
---|
69 | ax2.text(0.5, -tsol_max, "depth_max", fontweight='bold',ha='center', |
---|
70 | va='center') |
---|
71 | ax2.set_title("ngrnd="+str(len(tlev))) |
---|
72 | ax2.set_xlabel("Max thickness="+str(round(np.max(tz),5))+"[m]") |
---|
73 | # |
---|
74 | cb=plt.colorbar(datamap, cmap=cmap, shrink=0.5, norm=cNorm) |
---|
75 | cb.set_label("Layer Thickness [m]") |
---|
76 | # |
---|
77 | plt.suptitle(overalltitle) |
---|
78 | # |
---|
79 | # plt.show() |
---|
80 | pdf.savefig() |
---|
81 | # |
---|
82 | # |
---|
83 | # As a reminder, the levels of the hydrology as in the code today |
---|
84 | # |
---|
85 | depth_Wmax=2.0 |
---|
86 | nslm=11 |
---|
87 | zz=np.zeros(nslm, float) |
---|
88 | dz=np.zeros(nslm+1, float) |
---|
89 | thickness=np.zeros(nslm, float) |
---|
90 | depth=np.zeros(nslm+1, float) |
---|
91 | |
---|
92 | for i in range(1,nslm): |
---|
93 | zz[i] = depth_Wmax*((2**i)-1)/((2**(nslm-1))-1) |
---|
94 | dz[i] = zz[i]-zz[i-1] |
---|
95 | |
---|
96 | print "== Old depth of the hydrology for reference ==========================" |
---|
97 | print "zz=", zz |
---|
98 | print "dz=",dz |
---|
99 | print "======================================================================" |
---|
100 | print " " |
---|
101 | print " " |
---|
102 | # |
---|
103 | # |
---|
104 | # New Methode for computing levels with the new parameters to control. |
---|
105 | # |
---|
106 | # All units in meters |
---|
107 | # |
---|
108 | # Maximum depth of the soil thermodynamics |
---|
109 | depth_max = 8.0 |
---|
110 | try : |
---|
111 | depth_max=float(raw_input("Maximum depth of soil : depth_max [8.0m] >")) |
---|
112 | except ValueError: |
---|
113 | depth_max = 8.0 |
---|
114 | # |
---|
115 | # Maximum depth of soil moisture |
---|
116 | depth_Wmax=2.0 |
---|
117 | try : |
---|
118 | depth_Wmax=float(raw_input("Depth of hydrology (smaller than depth_max) : depth_Wmax [2.0m] >")) |
---|
119 | except ValueError: |
---|
120 | depth_Wmax=2.0 |
---|
121 | if ( depth_Wmax > depth_max) : |
---|
122 | print "ERROR : depth_Wmax needs to be smaller than depth_max." |
---|
123 | print "Correction : depth_Wmax set to depth_max" |
---|
124 | depth_Wmax = depth_max |
---|
125 | # |
---|
126 | # Thickness of first Layer |
---|
127 | depth_toptickness = 9.77517107e-04 |
---|
128 | try : |
---|
129 | depth_toptickness=float(raw_input("Depth of top hydrology layer : depth_toptickness [9.77517107e-04 m] >")) |
---|
130 | except ValueError: |
---|
131 | depth_toptickness = 9.77517107e-04 |
---|
132 | # The constant layer thickness starts at : |
---|
133 | depth_cstthickness=0.75 |
---|
134 | try : |
---|
135 | depth_cstthickness=float(raw_input("Depth at which constant layer thickness start (smaller than depth_Wmax/2) : depth_cstthickness [0.75 m] >")) |
---|
136 | except ValueError: |
---|
137 | depth_cstthickness = 0.75 |
---|
138 | if ( depth_cstthickness > depth_Wmax/2.0 ) : |
---|
139 | print "ERROR : depth_cstthickness needs to be smaller than depth_Wmax/2.0" |
---|
140 | print "Correction : Constant thickness disabled" |
---|
141 | depth_cstthickness = depth_Wmax |
---|
142 | # |
---|
143 | # Depth at which we refine the layers as we reach the depth_Wmax |
---|
144 | depth_refinebottom = False |
---|
145 | try : |
---|
146 | depth_refinebottom=str2bool(raw_input("Depth at which depth refinement starts (between depth_cstthickness and depth_Wmax) : depth_refinebottom [False] >")) |
---|
147 | except ValueError: |
---|
148 | depth_refinebottom=False |
---|
149 | # |
---|
150 | # Depth at which we resume geometrical increases |
---|
151 | depth_geom=2.5 |
---|
152 | try : |
---|
153 | depth_geom=float(raw_input("Depth at which geometrical increase starts again (larger than depth_Wmax) : depth_geom [2.5 m] >")) |
---|
154 | except ValueError: |
---|
155 | depth_geom=2.5 |
---|
156 | if ( depth_geom < depth_Wmax ) : |
---|
157 | print "ERROR : depth_geom needs to be larger than depth_Wmax" |
---|
158 | print "Correction : setting depth_geom to depth_Wmax" |
---|
159 | depth_geom=depth_Wmax |
---|
160 | # |
---|
161 | # Reason Below depth_geom |
---|
162 | reason_geom_below=2 |
---|
163 | try : |
---|
164 | reason_geom_below=float(raw_input("Reason of geometrical serie below depth_geom : reason_geom_below [2] >")) |
---|
165 | except ValueError: |
---|
166 | reason_geom_below=2 |
---|
167 | # |
---|
168 | # It is coded FORTRAN style so that translations is easier |
---|
169 | # |
---|
170 | nltmp=500 |
---|
171 | ztmp=np.zeros(nltmp, float) |
---|
172 | dtmp=np.zeros(nltmp+1, float) |
---|
173 | for i in range(1,nltmp): |
---|
174 | if ( ztmp[i-1] < depth_cstthickness ) : |
---|
175 | ztmp[i] = depth_toptickness*2.0*((2**i)-1) |
---|
176 | dtmp[i] = ztmp[i]-ztmp[i-1] |
---|
177 | igeo=i |
---|
178 | zgeoend=ztmp[i] |
---|
179 | dgeoend=dtmp[i] |
---|
180 | else : |
---|
181 | ztmp[i] = ztmp[i-1]+dtmp[i-1] |
---|
182 | dtmp[i] = dtmp[i-1] |
---|
183 | # |
---|
184 | nbrefine=0 |
---|
185 | drefine=np.zeros(nltmp, float) |
---|
186 | if (depth_refinebottom) : |
---|
187 | # |
---|
188 | # Compute parameters for the constant increment interval before refining, |
---|
189 | # If needed ! |
---|
190 | cstint=depth_Wmax-(2.0*zgeoend) |
---|
191 | nbcst=max(int(np.modf(cstint/dgeoend)[1]),0) |
---|
192 | if ( nbcst > 0 ) : |
---|
193 | dcst=cstint/nbcst |
---|
194 | else : |
---|
195 | dcst=dgeoend |
---|
196 | # |
---|
197 | # If we have to add constant increments |
---|
198 | # |
---|
199 | if ( nbcst > 0 ) : |
---|
200 | # Add linear levels |
---|
201 | for i in range(igeo+1,igeo+nbcst+1) : |
---|
202 | ztmp[i] = ztmp[i-1]+dcst |
---|
203 | dtmp[i] = dcst |
---|
204 | # Refine the levels toward the bottom |
---|
205 | for i in range(igeo+nbcst+1,2*igeo+nbcst+1) : |
---|
206 | irev=(2*igeo+nbcst+1)-i |
---|
207 | ztmp[i] = ztmp[i-1]+dtmp[irev] |
---|
208 | dtmp[i] = ztmp[i]-ztmp[i-1] |
---|
209 | drefine[nbrefine]=dtmp[i] |
---|
210 | nbrefine=nbrefine+1 |
---|
211 | # Without constant increments |
---|
212 | else : |
---|
213 | igeo=np.argmin(np.abs(ztmp-(depth_Wmax/2.0))) |
---|
214 | ztmp[igeo]=depth_Wmax/2.0 |
---|
215 | dtmp[igeo] = ztmp[igeo]-ztmp[igeo-1] |
---|
216 | for i in range(igeo+1,2*igeo+1) : |
---|
217 | irev=(2*igeo+1)-i |
---|
218 | ztmp[i] = ztmp[i-1]+dtmp[irev] |
---|
219 | dtmp[i] = ztmp[i]-ztmp[i-1] |
---|
220 | drefine[nbrefine]=dtmp[i] |
---|
221 | nbrefine=nbrefine+1 |
---|
222 | # |
---|
223 | nslm=np.argmin(np.abs(ztmp-depth_Wmax))+1 |
---|
224 | # |
---|
225 | # Extract now the values we need to reach depth_Wmax |
---|
226 | # |
---|
227 | zz=ztmp[0:nslm] |
---|
228 | dz=dtmp[0:nslm+1] |
---|
229 | dz[nslm]=0.0 |
---|
230 | print "== New values ========================================================" |
---|
231 | print "Number of layers in hydrol = ", nslm |
---|
232 | print "zz=", zz |
---|
233 | print "dz=",dz |
---|
234 | # |
---|
235 | # Extension for the thermodynamics below depth_Wmax |
---|
236 | # |
---|
237 | ztmp=np.zeros(nltmp, float) |
---|
238 | ztmp[0]=depth_toptickness/2 |
---|
239 | for i in range(1,nslm): |
---|
240 | ztmp[i]=zz[i] |
---|
241 | # |
---|
242 | # If we have created a refined region at the bottom of the soil moisture. We need to |
---|
243 | # unwinde it for temperature below depth_Wmax |
---|
244 | # |
---|
245 | if ( nbrefine > 1 ) : |
---|
246 | for i in range(nslm,nslm+nbrefine,1): |
---|
247 | iref=(nbrefine-1)-(i-nslm) |
---|
248 | ztmp[i]=ztmp[i-1]+drefine[iref] |
---|
249 | # |
---|
250 | # New nslm |
---|
251 | # |
---|
252 | nslm=nslm+nbrefine |
---|
253 | # |
---|
254 | for i in range(nslm,nltmp): |
---|
255 | if ( ztmp[i-1] < depth_geom ) : |
---|
256 | reason = 1.0 |
---|
257 | else : |
---|
258 | reason = reason_geom_below |
---|
259 | ztmp[i]=ztmp[i-1]+reason*(ztmp[i-1]-ztmp[i-2]) |
---|
260 | # |
---|
261 | ngrnd=np.argmin(np.abs(ztmp-depth_max))+1 |
---|
262 | # |
---|
263 | tdepth=np.zeros(ngrnd, float) |
---|
264 | tdepth=ztmp[0:ngrnd] |
---|
265 | print "Number of layers in thermosoil = ", ngrnd |
---|
266 | print "tdepth=",tdepth |
---|
267 | print "======================================================================" |
---|
268 | # |
---|
269 | # Do the plot |
---|
270 | # |
---|
271 | pdf = PdfPages('VerticalLayers.pdf') |
---|
272 | wdepth=np.copy(zz) |
---|
273 | wdepth[0]=(zz[1]-zz[0])/2 |
---|
274 | discretizationplot(wdepth, tdepth, depth_Wmax, depth_max, "Mixing water and temperature") |
---|
275 | pdf.close() |
---|