Context Navigation

source: tags/ORCHIDEE_2_0/ORCHIDEE/src_parameters/vertical_soil.f90 @ 6392

Last change on this file since 6392 was 5359, checked in by josefine.ghattas, 6 years ago
Integrated bug correction for refinebottom done in the trunk rev 5358 now also in the tag. See ticket #443
Property svn:keywords set to `Date Revision HeadURL`
File size: 20.0 KB

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

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

Download in other formats: