Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

vertical_soil.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 20.1 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: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/vertical_soil.f90 $
19	!! $Date: 2018-10-03 10:22:28 +0200 (Wed, 03 Oct 2018) $
20	!! $Revision: 5454 $
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) :: 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
208	!Config Help = Maximum depth of soil for soil moisture (CWRR).
209	!Config Units = m
210	zmaxh=2.0
211	CALL getin_p("DEPTH_MAX_H",zmaxh)
212
213	! Verification
214	IF ( zmaxh > zmaxt) THEN
215	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : zmaxh needs to be smaller than zmaxt.", &
216	"Correction : zmaxh set to zmaxt", " ")
217	ENDIF
218
219	!Config Key = DEPTH_TOPTHICK
220	!Config Desc = Thickness of upper most Layer
221	!Config If =
222	!Config Def = 9.77517107e-04
223	!Config Help = Thickness of top hydrology layer for soil moisture (CWRR).
224	!Config Units = m
225	depth_topthickness = 9.77517107e-04
226	CALL getin_p("DEPTH_TOPTHICK",depth_topthickness)
227
228	!Config Key = DEPTH_CSTTHICK
229	!Config Desc = Depth at which constant layer thickness start
230	!Config If =
231	!Config Def = DEPTH_MAX_H
232	!Config Help = Depth at which constant layer thickness start (smaller than zmaxh/2)
233	!Config Units = m
234	depth_cstthickness=zmaxh
235	CALL getin_p("DEPTH_CSTTHICK",depth_cstthickness)
236
237	IF ( (depth_cstthickness /= zmaxh) .AND. (depth_cstthickness > zmaxh/2.0) ) THEN
238	CALL ipslerr_p(2,'vertical_soil_init',"ERROR : depth_cstthickness needs to be smaller than zmaxh/2.0", &
239	"Correction : Constant thickness disabled", " ")
240	depth_cstthickness = zmaxh
241	ENDIF
242
243	!Config Key = REFINEBOTTOM
244	!Config Desc = Depth at which the hydrology layers will be refined towards the bottom.
245	!Config If =
246	!Config Def = .FALSE.
247	!Config Help = Refinebottom is important when lower boundary conditions is different from a free drainage.
248	!Config Units = -
249	refinebottom = .FALSE.
250	CALL getin_p("REFINEBOTTOM",refinebottom)
251
252	!Config Key = DEPTH_GEOM
253	!Config Desc = Depth at which we resume geometrical increases for temperature
254	!Config If =
255	!Config Def = DEPTH_MAX_H
256	!Config Help = Depth at which the thickness increases again for temperature.
257	! This depth has to be deeper than the bottom of the hydrology.
258	!Config Units = m
259	depth_geom=zmaxh
260	CALL getin_p("DEPTH_GEOM",depth_geom)
261
262	IF ( depth_geom < zmaxh ) THEN
263	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", &
264	"Correction : setting depth_geom to zmaxh", " ")
265	ENDIF
266
267	!Config Key = RATIO_GEOM_BELOW
268	!Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM
269	!Config If =
270	!Config Def = 2
271	!Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM.
272	! This parameter allows to cover the depth needed for temperature with fewer layers.
273	!Config Units = -
274	ratio_geom_below=2
275	CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below)
276
277	!
278	! Computing the layer depth for soil moisture. This defines the number of layers needed.
279	!
280	ztmp(:) = 0.0
281	dtmp(:) = 0.0
282	DO i=1,nblayermax
283	IF ( ztmp(i) < depth_cstthickness ) THEN
284	ztmp(i+1) = depth_topthickness2.0((2**i)-1)
285	dtmp(i+1) = ztmp(i+1)-ztmp(i)
286	igeo=i+1
287	zgeoend=ztmp(i+1)
288	dgeoend=dtmp(i+1)
289	ELSE
290	ztmp(i+1) = ztmp(i)+dtmp(i)
291	dtmp(i+1) = dtmp(i)
292	ENDIF
293	ENDDO
294	!
295	! refine at bottom. Work in progress...
296	!
297	nbrefine = 1
298	drefine(:) = 0.0
299	IF (refinebottom) THEN
300	!
301	! Compute parameters for the constant increment interval before refining,
302	! If needed !
303	cstint=zmaxh-(2.0*zgeoend)
304	nbcst=MAX(INT(cstint/dgeoend), 0)
305	IF ( nbcst > 0 ) THEN
306	dcst=cstint/nbcst
307	ELSE
308	dcst=dgeoend
309	ENDIF
310	!
311	! If we have to add constant increments
312	!
313	IF ( nbcst > 0 ) THEN
314	! Add linear levels
315	DO i=igeo,igeo+nbcst-1
316	ztmp(i+1) = ztmp(i)+dcst
317	dtmp(i+1) = dcst
318	ENDDO
319	! Refine the levels toward the bottom
320	DO i=igeo+nbcst,2*igeo+nbcst-2
321	irev=(2*igeo+nbcst-1)-i
322	ztmp(i+1) = ztmp(i)+dtmp(irev+1)
323	dtmp(i+1) = ztmp(i+1)-ztmp(i)
324	drefine(nbrefine)=dtmp(i+1)
325	nbrefine=nbrefine+1
326	ENDDO
327	! Without constant increments
328	ELSE
329	imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0)))
330	igeo=imin(1)
331	ztmp(igeo)=zmaxh/2.0
332	dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1)
333	DO i=igeo,2*igeo
334	irev=(2*igeo)-i
335	ztmp(i+1) = ztmp(i)+dtmp(irev)
336	dtmp(i+1) = ztmp(i-1)-ztmp(i)
337	drefine(nbrefine)=dtmp(i+1)
338	nbrefine=nbrefine+1
339	ENDDO
340	ENDIF
341	nbrefine = nbrefine-1
342	ENDIF
343
344
345
346
347	! Find the index (nslm) of the node closest to the zmaxh
348	imin=MINLOC(ABS(ztmp(:)-zmaxh))
349	nslm=imin(1)
350	!
351	! ALLOCATE the arrays we need to keep
352	!
353	ALLOCATE(znh(nslm), stat=ier)
354	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znh','','')
355	ALLOCATE(dnh(nslm), stat=ier)
356	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dnh','','')
357	ALLOCATE(dlh(nslm), stat=ier)
358	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlh','','')
359	ALLOCATE(zlh(nslm), stat=ier)
360	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlh','','')
361	!
362	! Extract now the values we need to reach zmaxh
363	!
364	znh(:)=ztmp(1:nslm)
365
366	! Force the last hydrological layer and node to be zmaxh
367	znh(nslm)=zmaxh
368	dnh(:)=dtmp(1:nslm)
369	! Recalculate the distance between the last 2 nodes
370	dnh(nslm) = ztmp(nslm)-ztmp(nslm-1)
371
372	! Calculate the thickness of the layers
373	DO i = 1, nslm-1
374	dlh(i) = (dnh(i) + dnh(i+1)) / 2.0
375	ENDDO
376	dlh(nslm) = dnh(nslm)/2
377
378	!
379	! Extension for the thermodynamics below zmaxh
380	!
381	ntmp = nslm
382	ztmp(:)=0.0
383	! Exception for the first layer. It is at the top in the hydrology and in
384	! the middle for temperature.
385	ztmp(1)=depth_topthickness/2
386	DO i=2,ntmp
387	ztmp(i)=znh(i)
388	ENDDO
389
390	! Exception for the last layer where again temperature needs to be in
391	! the middle of the layer. Also add a layer with the same thickness
392	hh=dnh(ntmp)/2.0
393	ztmp(ntmp)=ztmp(ntmp)-hh/2.0
394	ztmp(ntmp+1)=ztmp(ntmp)+hh*1.5
395	ztmp(ntmp+2)=ztmp(ntmp+1)+hh*2.0
396	ntmp=ntmp+2
397
398	! If we have created a refined region at the bottom of the soil moisture. We need to
399	! unwinde it for temperature below zmaxh. Only done for option refinebottom.
400	IF ( nbrefine > 1 ) THEN
401	DO i=ntmp,ntmp+nbrefine-1
402	iref=nbrefine-(i-ntmp)
403	ztmp(i+1)=ztmp(i)+drefine(iref)
404	ENDDO
405	ntmp=ntmp+nbrefine
406	ENDIF
407	!
408	! Resume the geometric increas of thickness
409	DO i=ntmp,nblayermax
410	IF ( ztmp(i) < depth_geom ) THEN
411	ratio = 1.0
412	ELSE
413	ratio = ratio_geom_below
414	ENDIF
415	ztmp(i+1)=ztmp(i)+ratio*(ztmp(i)-ztmp(i-1))
416	ENDDO
417
418	! Compute the depth of the lower interface of each layer.
419	!
420	zint(1) = depth_topthickness
421	DO i=2,nblayermax-1
422	zint(i) = (ztmp(i)+ztmp(i+1))/2.0
423	ENDDO
424	zint(nslm-1) = (znh(nslm-1) + znh(nslm))/2.0
425	zint(nslm) = (znh(nslm))
426	zint(nblayermax) = ztmp(nblayermax)+(ztmp(nblayermax)-ztmp(nblayermax-1))/2.0
427
428	! Determine the total number of layers for thermal (ngrnd).
429	!
430	imin=MINLOC(ABS(zint(:)-zmaxt))
431	ngrnd=imin(1)
432	ALLOCATE(znt(ngrnd), stat=ier)
433	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znt','','')
434	ALLOCATE(dlt(ngrnd), stat=ier)
435	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlt','','')
436	ALLOCATE(zlt(ngrnd), stat=ier)
437	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlt','','')
438
439	! Assign values for znt, zlt, dlt
440	!
441	znt(:)=ztmp(1:ngrnd)
442	zlt(:) = zint(1:ngrnd)
443	dlt(1) = zint(1)
444	DO i=2,ngrnd
445	dlt(i) = zint(i)-zint(i-1)
446	ENDDO
447	! Depth of layers for the hydrology are the same as for thermosoil but only the upper nslm layers are used
448	zlh(:) = zlt(1:nslm)
449
450	! Force the last thermal layer and node to be zmaxt
451	zlt(ngrnd) = zmaxt
452	dlt(ngrnd) = zmaxt-zint(ngrnd-1)
453
454
455
456
457	IF (printlev >=1) THEN
458	WRITE(numout,*) "=== Vertical discretization for hydrology ==================================="
459	WRITE(numout,*) "Number of layers in hydrology: nslm=", nslm
460	WRITE(numout,*) "Total depth in hydrology (DEPTH_MAX_H) (m): znh(nslm)=", znh(nslm)
461	WRITE(numout,*) "Depth of nodes in hydrology (m): znh=", znh
462	WRITE(numout,*) "Depth of lower layer-interface for hydrology (m): zlh=", zlh
463	WRITE(numout,*) "Layer thickness for hydrology (m): dlh=",dlh
464	WRITE(numout,*) "Internode distance for hydrology (m): dnh=",dnh
465	WRITE(numout,*) "=== Vertical discretization for thermosoil =================================="
466	WRITE(numout,*) "Number of layers in thermosoil: ngrnd=", ngrnd
467	WRITE(numout,*) "Total depth for thermosoil (DEPTH_MAX_T) (m): zlt=", zlt(ngrnd)
468	WRITE(numout,*) "to be compared with calculated total depth for thermosoil (not used) zint=", zint(ngrnd)
469	WRITE(numout,*) "Depth of the nodes in thermosoil (m): znt=",znt
470	WRITE(numout,*) "Depth of lower layer-interface for thermosoil (m): zlt=", zlt(1:ngrnd)
471	WRITE(numout,*) "Layer thickness for thermosoil (m): dlt=", dlt
472	WRITE(numout,*) "============================================================================="
473	END IF
474	END SUBROUTINE vertical_soil_init
475
476	END MODULE vertical_soil

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_parameters/vertical_soil.f90 @ 7541

Download in other formats: