Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_parameters/vertical_soil.f90 @ 5925

Last change on this file since 5925 was 3224, checked in by josefine.ghattas, 8 years ago

Missing error handeling.

A Jornet-Puig

Property svn:keywords set to Date Revision HeadURL

File size: 19.2 KB

Line
1	MODULE vertical_soil
2	!! $HeadURL$
3	!! $Date$
4	!! $Revision$
5
6	USE defprec
7	USE ioipsl_para
8	USE vertical_soil_var
9
10	IMPLICIT NONE
11
12
13	PUBLIC :: vertical_soil_init
14
15	CONTAINS
16	!! ================================================================================================================================
17	!! SUBROUTINE : vertical_soil_init
18	!!
19	!>\BRIEF
20	!!
21	!! DESCRIPTION : To define the total number of layers
22	!! for hydrology model (nslm),
23	!! and for thermodynamics model ngrnd.
24	!!
25	!! We have 2 type of levels in the hydrology and the thermodynamics :
26	!! (1) the nodes where the scalar variables are computed.
27	!! (2) the intermediate levels: these are the bounds of the layers and thus the
28	!! contact point between layer i and i+1.
29	!!
30	!! Situation in hydrol.f90
31	!! The CWRR model only uses two variables to describe the vertical discretization :
32	!! (1) znh : the depth of the nodes where soil moisture is computed. The nodes are not
33	!! supposed to be in the middle of the layer. For the first and last layers the
34	!! nodes are put onto the intermediate level. At 0 for the first layer and
35	!! zmaxh for the last layer.
36	!! (2) dnh : the distance between two nodes
37	!!
38	!! Situation in thermosoil.f90 :
39	!! For the thermodynamics other variables are needed to describe the vertical structure :
40	!! (1) znt(znh) : the depth of the nodes are the same as in hydrol.f90 except for
41	!! the first and last levels. In older versions than revision XXXX of thermosoil it is
42	!! computed with fz(i+1/2)
43	!! (2) zlt (Zint) (zz_coef in thermosoil) : is the depth of the intermediate levels.
44	!! In revision before XXXX of thermosoil.f90 it's computed by fz(i).
45	!! Since revision XXXX, in the new formulation (vertical.f90) it is computed as
46	!! zint(i) = (znh(i)+znh(i+1))/2
47	!! (3) dlt(dz_tmp) (dz2 in thermosoil): This is the thickness of the layer, i.e the
48	!! distance between the top and bottom of the layer. Thus it is given by
49	!! zint(i+1)-zint(i). It should not be confused with dnh in hydrology which is the
50	!! distance between nodes.
51	!!
52	!! in standard hydrology model (de Rosnay Patricia, 2000; depth_topthickness = 9.77517107e-04, zmaxh=2.0, refinebottom = .FALSE.):
53	!!
54	!! Example for standard depth of the hydrology
55	!! Number of layers in hydrol = 11
56	!! znh = depth of nodes
57	!! [ 0.00000000e+00 1.95503421e-03 5.86510264e-03 1.36852395e-02
58	!! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01
59	!! 4.98533725e-01 9.99022483e-01 2.00000000e+00 ]
60	!! dnh = internode distance
61	!! [ 0. 0.00195503 0.00391007 0.00782014 0.01564027 0.03128055
62	!! 0.06256109 0.12512219 0.25024438 0.50048876 1.00097752 ]
63	!! dlh = soil layer thickness
64	!! [ 0.00097752 0.00293255 0.0058651 0.01173021 0.02346041 0.04692082
65	!! 0.09384164 0.18768328 0.37536657 0.75073314 0.50048876 ]
66	!! hcum = depth of soil layer bottom
67	!! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02
68	!! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01
69	!! 7.48778104e-01 1.49951124e+00 2.00000000e+00 ]
70	!!
71	!!
72	!! In thermal model (an example of: depth_topthickness = 9.77517107e-04, zmaxt=10, depth_geom=10, refinebottom = .FALSE.):
73	!! Number of layers in thermosoil = 19
74	!! The approximate maximal depth for thermosoil = 10.0078201332
75	!! The actual maximal depth for thermosoil = 10.0
76	!! znt=
77	!! [ 4.88758554e-04 1.95503421e-03 5.86510264e-03 1.36852395e-02
78	!! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01
79	!! 4.98533725e-01 9.99022483e-01 1.74975562e+00 2.50048876e+00
80	!! 3.50146627e+00 4.50244379e+00 5.50342131e+00 6.50439882e+00
81	!! 7.50537634e+00 8.50635386e+00 9.50733137e+00 ]
82	!! dlt=
83	!! [ 9.77517107e-04 2.93255132e-03 5.86510264e-03 1.17302053e-02
84	!! 2.34604106e-02 4.69208211e-02 9.38416423e-02 1.87683285e-01
85	!! 3.75366569e-01 7.50733138e-01 5.00488758e-01 1.00097752e+00
86	!! 1.00097752e+00 1.00097752e+00 1.00097752e+00 1.00097752e+00
87	!! 1.00097752e+00 1.00097752e+00 9.93157383e-01 ]
88	!! zlt=
89	!! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02
90	!! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01
91	!! 7.48778104e-01 1.49951124e+00 2.00000000e+00 3.00097752e+00
92	!! 4.00195503e+00 5.00293255e+00 6.00391007e+00 7.00488758e+00
93	!! 8.00586510e+00 9.00684262e+00 1.00000000e+01 ]
94	!!
95	!!
96	!! A simple figure below shows the discretization. ^
97	!! '------' means interface; 'X' means node; '...' means many layers; \| means distance.
98	!! v
99	!! Keep in mind that the nodes are not necessarly in the middle of the layers.
100	!!
101	!! Hydrology Thermodynamics
102	!!
103	!! --------X------- znh(1) ^ ------------------- 1st zlt(1) ^
104	!! \| \|
105	!! \| X znt(1) \|depth_topthickness, dlt(1)
106	!! \|dnh(2) \|
107	!! ----------------- \| ------------------- 2nd zlt(2) v ^
108	!! \| \|
109	!! X znh(2) v ^ X znt(2) \|dlt(2)
110	!! \| \|
111	!! ----------------- \|dnh(3) ------------------- 3rd zlt(3) v
112	!! \|
113	!! X znh(3) v X znt(3)
114	!!
115	!! ----------------- ... ------------------- ...
116	!!
117	!! ... X ...
118	!!
119	!! --------X------- znh(nslm) ------------------- nslm zlt(nslm)
120	!!
121	!! X znt(nslm)
122	!!
123	!! ------------------- nslm+1 zlt(nslm+1)
124	!!
125	!! X znt(nslm+1)
126	!!
127	!! ------------------- ...
128	!!
129	!! X znt(ngrnd)
130	!!
131	!! ------------------- ngrnd zlt(ngrnd)
132	!!
133	!!
134	!! RECENT CHANGE(S): None
135	!!
136	!! MAIN OUTPUT VARIABLE(S):
137	!!
138	!! REFERENCE(S) :
139	!!
140	!! FLOWCHART :
141	!! \n
142	!_ ================================================================================================================================
143
144	SUBROUTINE vertical_soil_init
145
146	!! 0. Variable and parameter declaration
147	!! 0.1 Local variables
148	INTEGER(i_std), PARAMETER :: nblayermax=500 !! Preset 500 for max. number of soil layers.
149	INTEGER(i_std) :: i, irev, iref, ntmp, ier
150	REAL(r_std) :: zgeoend !! The node position where the geometrical increases of nodes stopped. (m)
151	REAL(r_std) :: dgeoend !! The distance of the node at zgeoend and the node above. (m)
152	REAL(r_std) :: dcst !! Related to refine at bottom. Work in progress...
153	REAL(r_std) :: cstint !! Related to refine at bottom. Work in progress...
154	REAL(r_std) :: hh !! Temporay variable, to calculate the layer thickness for temperature at zmaxh.
155	REAL(r_std) :: ratio !! The ratio of geometric increas.
156	INTEGER(i_std) :: nbrefine !! Related to refine at bottom. Work in progress...
157	INTEGER(i_std) :: nbcst !! Related to refine at bottom. Work in progress...
158	INTEGER(i_std) :: igeo !! Related to refine at bottom. Work in progress...
159	REAL(r_std), DIMENSION(nblayermax+1) :: ztmp !! Depth at the node of each layer (m)
160	REAL(r_std), DIMENSION(nblayermax+1) :: zint !! Depth at the interface of each layer (m)
161	REAL(r_std), DIMENSION(nblayermax+1) :: dtmp !! Distance between the current node and the one above (m)
162	REAL(r_std), DIMENSION(nblayermax+1) :: drefine !! Related to refine at bottom. Work in progress...
163	INTEGER(i_std), DIMENSION(1) :: imin !! Related to refine at bottom. Work in progress...
164
165	!! Variables controling the vertical discretiztaion
166	REAL(r_std) :: depth_topthickness !! Thickness of the top Layer (m)
167	REAL(r_std) :: depth_cstthickness !! Depth at which constant layer thickness start (m)
168	REAL(r_std) :: depth_geom !! Depth at which we resume geometrical increases for temperature (m)
169	REAL(r_std) :: ratio_geom_below !! Ratio of the geometrical series defining the thickness below DEPTH_GEOM.
170	LOGICAL :: refinebottom !! Whether or not the hydrology layers will be refined towards the bottom.
171
172
173	!! 1. Read parameters from run.def file
174
175	!Config Key = DEPTH_MAX_T
176	!Config Desc = Maximum depth of the soil thermodynamics
177	!Config If =
178	!Config Def = 10.0
179	!Config Help = Maximum depth of soil for temperature.
180	!Config Units = m
181	zmaxt = 10.0
182	CALL getin_p("DEPTH_MAX_T",zmaxt)
183
184	!Config Key = DEPTH_MAX_H
185	!Config Desc = Maximum depth of soil moisture
186	!Config If =
187	!Config Def = 2.0 or 4.0 depending on hydrol_cwrr
188	!Config Help = Maximum depth of soil for soil moisture (CWRR).
189	!Config Units = m
190	! Default value for CWRR (only CWRR enters this subroutine)
191	zmaxh=2.0
192	CALL getin_p("DEPTH_MAX_H",zmaxh)
193
194	! Verification
195	IF ( zmaxh > zmaxt) THEN
196	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : zmaxh needs to be smaller than zmaxt.", &
197	"Correction : zmaxh set to zmaxt", " ")
198	ENDIF
199
200	!Config Key = DEPTH_TOPTHICK
201	!Config Desc = Thickness of upper most Layer
202	!Config If =
203	!Config Def = 9.77517107e-04
204	!Config Help = Thickness of top hydrology layer for soil moisture (CWRR).
205	!Config Units = m
206	depth_topthickness = 9.77517107e-04
207	CALL getin_p("DEPTH_TOPTHICK",depth_topthickness)
208
209	!Config Key = DEPTH_CSTTHICK
210	!Config Desc = Depth at which constant layer thickness start
211	!Config If =
212	!Config Def = DEPTH_MAX_H
213	!Config Help = Depth at which constant layer thickness start (smaller than zmaxh/2)
214	!Config Units = m
215	depth_cstthickness=zmaxh
216	CALL getin_p("DEPTH_CSTTHICK",depth_cstthickness)
217
218	IF ( (depth_cstthickness /= zmaxh) .AND. (depth_cstthickness > zmaxh/2.0) ) THEN
219	CALL ipslerr_p(2,'vertical_soil_init',"ERROR : depth_cstthickness needs to be smaller than zmaxh/2.0", &
220	"Correction : Constant thickness disabled", " ")
221	depth_cstthickness = zmaxh
222	ENDIF
223
224	! REFINEBOTTOM option is under work... This option can not be set from the parameter file for the moment.
225	! All the related code for refinebottom except this getin are found in this module but needs to be tested.
226	!!$ !Config Key = REFINEBOTTOM
227	!!$ !Config Desc = Depth at which the hydrology layers will be refined towards the bottom.
228	!!$ !Config If =
229	!!$ !Config Def = .FALSE.
230	!!$ !Config Help = Depth at which the hydrology layers will be refined towards the bottom.
231	!!$ ! This is important when lower boundary conditions is different from a free drainage.
232	!!$ !Config Units = -
233	!!$ !
234	refinebottom = .FALSE.
235	!!$ CALL getin_p("REFINEBOTTOM",refinebottom)
236
237	!Config Key = DEPTH_GEOM
238	!Config Desc = Depth at which we resume geometrical increases for temperature
239	!Config If =
240	!Config Def = DEPTH_MAX_H
241	!Config Help = Depth at which the thickness increases again for temperature.
242	! This depth has to be deeper than the bottom of the hydrology.
243	!Config Units = m
244	depth_geom=zmaxh
245	CALL getin_p("DEPTH_GEOM",depth_geom)
246
247	IF ( depth_geom < zmaxh ) THEN
248	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", &
249	"Correction : setting depth_geom to zmaxh", " ")
250	ENDIF
251
252	!Config Key = RATIO_GEOM_BELOW
253	!Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM
254	!Config If =
255	!Config Def = 2
256	!Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM.
257	! This parameter allows to cover the depth needed for temperature with fewer layers.
258	!Config Units = -
259	ratio_geom_below=2
260	CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below)
261
262	!
263	! Computing the layer depth for soil moisture. This defines the number of layers needed.
264	!
265	ztmp(:) = 0.0
266	dtmp(:) = 0.0
267	DO i=1,nblayermax
268	IF ( ztmp(i) < depth_cstthickness ) THEN
269	ztmp(i+1) = depth_topthickness2.0((2**i)-1)
270	dtmp(i+1) = ztmp(i+1)-ztmp(i)
271	igeo=i+1
272	zgeoend=ztmp(i+1)
273	dgeoend=dtmp(i+1)
274	ELSE
275	ztmp(i+1) = ztmp(i)+dtmp(i)
276	dtmp(i+1) = dtmp(i)
277	ENDIF
278	ENDDO
279	!
280	! refine at bottom. Work in progress...
281	!
282	nbrefine = 1
283	drefine(:) = 0.0
284	IF (refinebottom) THEN
285	!
286	! Compute parameters for the constant increment interval before refining,
287	! If needed !
288	cstint=zmaxh-(2.0*zgeoend)
289	nbcst=MAX(INT(cstint/dgeoend), 0)
290	IF ( nbcst > 0 ) THEN
291	dcst=cstint/nbcst
292	ELSE
293	dcst=dgeoend
294	ENDIF
295	!
296	! If we have to add constant increments
297	!
298	IF ( nbcst > 0 ) THEN
299	! Add linear levels
300	DO i=igeo,igeo+nbcst-1
301	ztmp(i+1) = ztmp(i)+dcst
302	dtmp(i+1) = dcst
303	ENDDO
304	! Refine the levels toward the bottom
305	DO i=igeo+nbcst,2*igeo+nbcst-2
306	irev=(2*igeo+nbcst-1)-i
307	ztmp(i+1) = ztmp(i)+dtmp(irev+1)
308	dtmp(i+1) = ztmp(i+1)-ztmp(i)
309	drefine(nbrefine)=dtmp(i+1)
310	nbrefine=nbrefine+1
311	ENDDO
312	! Without constant increments
313	ELSE
314	imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0)))
315	igeo=imin(1)
316	ztmp(igeo)=zmaxh/2.0
317	dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1)
318	DO i=igeo,2*igeo
319	irev=(2*igeo)-i
320	ztmp(i+1) = ztmp(i)+dtmp(irev)
321	dtmp(i+1) = ztmp(i-1)-ztmp(i)
322	drefine(nbrefine)=dtmp(i+1)
323	nbrefine=nbrefine+1
324	ENDDO
325	ENDIF
326	ENDIF
327
328	! Find the index (nslm) of the node closest to the zmaxh
329	imin=MINLOC(ABS(ztmp(:)-zmaxh))
330	nslm=imin(1)
331	!
332	! ALLOCATE the arrays we need to keep
333	!
334	ALLOCATE(znh(nslm), stat=ier)
335	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znh','','')
336	ALLOCATE(dnh(nslm), stat=ier)
337	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dnh','','')
338	ALLOCATE(dlh(nslm), stat=ier)
339	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlh','','')
340	ALLOCATE(zlh(nslm), stat=ier)
341	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlh','','')
342	!
343	! Extract now the values we need to reach zmaxh
344	!
345	znh(:)=ztmp(1:nslm)
346
347	! Force the last hydrological layer and node to be zmaxh
348	znh(nslm)=zmaxh
349	dnh(:)=dtmp(1:nslm)
350	! Recalculate the distance between the last 2 nodes
351	dnh(nslm) = ztmp(nslm)-ztmp(nslm-1)
352
353	! Calculate the thickness of the layers
354	DO i = 1, nslm-1
355	dlh(i) = (dnh(i) + dnh(i+1)) / 2.0
356	ENDDO
357	dlh(nslm) = dnh(nslm)/2
358
359	WRITE(numout,*) "=========================================================="
360	WRITE(numout,*) "Number of layers in hydrol = ", nslm
361	WRITE(numout,*) "znh = depth of nodes"
362	WRITE(numout,*) znh(:)
363	WRITE(numout,*) "dnh = internode distance"
364	WRITE(numout,*) dnh(:)
365	WRITE(numout,*) "dlh = layer thickness"
366	WRITE(numout,*) dlh(:)
367	WRITE(numout,*) "Total depth in hydrol has been calculated to ", znh(nslm)
368	WRITE(numout,*) "======================================================================"
369
370	!
371	! Extension for the thermodynamics below zmaxh
372	!
373	ntmp = nslm
374	ztmp(:)=0.0
375	! Exception for the first layer. It is at the top in the hydrology and in
376	! the middle for temperature.
377	ztmp(1)=depth_topthickness/2
378	DO i=2,ntmp
379	ztmp(i)=znh(i)
380	ENDDO
381
382	! Exception for the last layer where again temperature needs to be in
383	! the middle of the layer. Also add a layer with the same thickness
384	hh=dnh(ntmp)/2.0
385	ztmp(ntmp)=ztmp(ntmp)-hh/2.0
386	ztmp(ntmp+1)=ztmp(ntmp)+hh*1.5
387	ztmp(ntmp+2)=ztmp(ntmp+1)+hh*2.0
388	ntmp=ntmp+2
389
390	! If we have created a refined region at the bottom of the soil moisture. We need to
391	! unwinde it for temperature below zmaxh
392	IF ( nbrefine > 1 ) THEN
393	DO i=ntmp,ntmp+nbrefine
394	iref=(nbrefine+1)-(i-ntmp)
395	ztmp(i+1)=ztmp(i)+drefine(iref)
396	ENDDO
397	ntmp=ntmp+nbrefine
398	ENDIF
399	!
400	! Resume the geometric increas of thickness
401	DO i=ntmp,nblayermax
402	IF ( ztmp(i) < depth_geom ) THEN
403	ratio = 1.0
404	ELSE
405	ratio = ratio_geom_below
406	ENDIF
407	ztmp(i+1)=ztmp(i)+ratio*(ztmp(i)-ztmp(i-1))
408	ENDDO
409
410	! Compute the depth of the lower interface of each layer.
411	!
412	zint(1) = depth_topthickness
413	DO i=2,nblayermax-1
414	zint(i) = (ztmp(i)+ztmp(i+1))/2.0
415	ENDDO
416	zint(nslm-1) = (znh(nslm-1) + znh(nslm))/2.0
417	zint(nslm) = (znh(nslm))
418	zint(nblayermax) = ztmp(nblayermax)+(ztmp(nblayermax)-ztmp(nblayermax-1))/2.0
419
420	! Determine the total number of layers for thermal (ngrnd).
421	!
422	imin=MINLOC(ABS(zint(:)-zmaxt))
423	ngrnd=imin(1)
424	ALLOCATE(znt(ngrnd), stat=ier)
425	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znt','','')
426	ALLOCATE(dlt(ngrnd), stat=ier)
427	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlt','','')
428	ALLOCATE(zlt(ngrnd), stat=ier)
429	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlt','','')
430
431	! Assign values for znt, zlt, dlt
432	!
433	znt(:)=ztmp(1:ngrnd)
434	zlt(:) = zint(1:ngrnd)
435	dlt(1) = zint(1)
436	DO i=2,ngrnd
437	dlt(i) = zint(i)-zint(i-1)
438	ENDDO
439	! Depth of layers for the hydrology are the same as for thermosoil but only the upper nslm layers are used
440	zlh(:) = zlt(1:nslm)
441
442	! Force the last thermal layer and node to be zmaxt
443	zlt(ngrnd) = zmaxt
444	dlt(ngrnd) = zmaxt-zint(ngrnd-1)
445
446	WRITE(numout,*) "Number of layers in thermosoil = ", ngrnd
447	WRITE(numout,*) "The maximal depth for thermosoil (DEPTH_MAX_T) = ", zlt(ngrnd)
448	WRITE(numout,*) "to be compared with calculated maximal depth for thermosoil (not used) = ", zint(ngrnd)
449	WRITE(numout,*) "Depth of the nodes in thermosoil, znt=",znt
450	WRITE(numout,*) "Thickness of the layers, dlt=", dlt
451	WRITE(numout,*) "Depth of the interface between the layers, zlt=", zlt(1:ngrnd)
452	WRITE(numout,*) "======================================================================"
453
454	END SUBROUTINE vertical_soil_init
455
456	END MODULE vertical_soil

Note: See TracBrowser for help on using the repository browser.

Download in other formats: