Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_parameters/vertical_soil.f90 @ 8398

Last change on this file since 8398 was 4282, checked in by josefine.ghattas, 7 years ago
Added module comments following coding guidelines.
Property svn:keywords set to `Date Revision HeadURL`
File size: 20.3 KB

Line
1	! ================================================================================================================================
2	! MODULE : vertical_soil
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.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 = 10.0
199	!Config Help = Maximum depth of soil for temperature.
200	!Config Units = m
201	zmaxt = 10.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	! REFINEBOTTOM option is under work... This option can not be set from the parameter file for the moment.
245	! All the related code for refinebottom except this getin are found in this module but needs to be tested.
246	!!$ !Config Key = REFINEBOTTOM
247	!!$ !Config Desc = Depth at which the hydrology layers will be refined towards the bottom.
248	!!$ !Config If =
249	!!$ !Config Def = .FALSE.
250	!!$ !Config Help = Depth at which the hydrology layers will be refined towards the bottom.
251	!!$ ! This is important when lower boundary conditions is different from a free drainage.
252	!!$ !Config Units = -
253	!!$ !
254	refinebottom = .FALSE.
255	!!$ CALL getin_p("REFINEBOTTOM",refinebottom)
256
257	!Config Key = DEPTH_GEOM
258	!Config Desc = Depth at which we resume geometrical increases for temperature
259	!Config If =
260	!Config Def = DEPTH_MAX_H
261	!Config Help = Depth at which the thickness increases again for temperature.
262	! This depth has to be deeper than the bottom of the hydrology.
263	!Config Units = m
264	depth_geom=zmaxh
265	CALL getin_p("DEPTH_GEOM",depth_geom)
266
267	IF ( depth_geom < zmaxh ) THEN
268	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", &
269	"Correction : setting depth_geom to zmaxh", " ")
270	ENDIF
271
272	!Config Key = RATIO_GEOM_BELOW
273	!Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM
274	!Config If =
275	!Config Def = 2
276	!Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM.
277	! This parameter allows to cover the depth needed for temperature with fewer layers.
278	!Config Units = -
279	ratio_geom_below=2
280	CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below)
281
282	!
283	! Computing the layer depth for soil moisture. This defines the number of layers needed.
284	!
285	ztmp(:) = 0.0
286	dtmp(:) = 0.0
287	DO i=1,nblayermax
288	IF ( ztmp(i) < depth_cstthickness ) THEN
289	ztmp(i+1) = depth_topthickness2.0((2**i)-1)
290	dtmp(i+1) = ztmp(i+1)-ztmp(i)
291	igeo=i+1
292	zgeoend=ztmp(i+1)
293	dgeoend=dtmp(i+1)
294	ELSE
295	ztmp(i+1) = ztmp(i)+dtmp(i)
296	dtmp(i+1) = dtmp(i)
297	ENDIF
298	ENDDO
299	!
300	! refine at bottom. Work in progress...
301	!
302	nbrefine = 1
303	drefine(:) = 0.0
304	IF (refinebottom) THEN
305	!
306	! Compute parameters for the constant increment interval before refining,
307	! If needed !
308	cstint=zmaxh-(2.0*zgeoend)
309	nbcst=MAX(INT(cstint/dgeoend), 0)
310	IF ( nbcst > 0 ) THEN
311	dcst=cstint/nbcst
312	ELSE
313	dcst=dgeoend
314	ENDIF
315	!
316	! If we have to add constant increments
317	!
318	IF ( nbcst > 0 ) THEN
319	! Add linear levels
320	DO i=igeo,igeo+nbcst-1
321	ztmp(i+1) = ztmp(i)+dcst
322	dtmp(i+1) = dcst
323	ENDDO
324	! Refine the levels toward the bottom
325	DO i=igeo+nbcst,2*igeo+nbcst-2
326	irev=(2*igeo+nbcst-1)-i
327	ztmp(i+1) = ztmp(i)+dtmp(irev+1)
328	dtmp(i+1) = ztmp(i+1)-ztmp(i)
329	drefine(nbrefine)=dtmp(i+1)
330	nbrefine=nbrefine+1
331	ENDDO
332	! Without constant increments
333	ELSE
334	imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0)))
335	igeo=imin(1)
336	ztmp(igeo)=zmaxh/2.0
337	dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1)
338	DO i=igeo,2*igeo
339	irev=(2*igeo)-i
340	ztmp(i+1) = ztmp(i)+dtmp(irev)
341	dtmp(i+1) = ztmp(i-1)-ztmp(i)
342	drefine(nbrefine)=dtmp(i+1)
343	nbrefine=nbrefine+1
344	ENDDO
345	ENDIF
346	ENDIF
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
401	IF ( nbrefine > 1 ) THEN
402	DO i=ntmp,ntmp+nbrefine
403	iref=(nbrefine+1)-(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 >=2) 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: