Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_parameters/vertical_soil.f90 @ 5491

Last change on this file since 5491 was 3564, checked in by albert.jornet, 8 years ago

Merge: from [3313:3545/trunk/ORCHIDEE]
Clean: output subroutine variables compile warning messages are solved.

Done in branches/ORCHIDEE-MICT/ORCHIDEE_MICT_TRUNK

Property svn:keywords set to Date Revision HeadURL

File size: 19.5 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 ^
104	!! \| \|
105	!! \| X znt(1) \| depth_topthickness, dlt(1)
106	!! \|dnh(2) \|
107	!! ----------------- \| ------------------- 2nd zlt(1) v ^
108	!! \| \|
109	!! X znh(2) v ^ X znt(2) \|dlt(2)
110	!! \| \|
111	!! ----------------- \|dnh(3) ------------------- 3rd zlt(2) 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 = 38
179	!Config Help = Maximum depth of soil for temperature.
180	!Config Units = m
181	!For carbon permafrost, to get ndeep=32 layers, we set DEPTH_TMAX=38. and RATIO_GEOM_BELOW=1.05
182	zmaxt = 38.0
183	CALL getin_p("DEPTH_MAX_T",zmaxt)
184
185	!Config Key = DEPTH_MAX_H
186	!Config Desc = Maximum depth of soil moisture
187	!Config If =
188	!Config Def = 2.0 or 4.0 depending on hydrol_cwrr
189	!Config Help = Maximum depth of soil for soil moisture (CWRR).
190	!Config Units = m
191	! Default value for CWRR (only CWRR enters this subroutineDEPTH_TMAX=38.)
192	zmaxh=2.0
193	CALL getin_p("DEPTH_MAX_H",zmaxh)
194
195	! Verification
196	IF ( zmaxh > zmaxt) THEN
197	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : zmaxh needs to be smaller than zmaxt.", &
198	"Correction : zmaxh set to zmaxt", " ")
199	ENDIF
200
201	!Config Key = DEPTH_TOPTHICK
202	!Config Desc = Thickness of upper most Layer
203	!Config If =
204	!Config Def = 9.77517107e-04
205	!Config Help = Thickness of top hydrology layer for soil moisture (CWRR).
206	!Config Units = m
207	depth_topthickness = 9.77517107e-04
208	CALL getin_p("DEPTH_TOPTHICK",depth_topthickness)
209
210	!Config Key = DEPTH_CSTTHICK
211	!Config Desc = Depth at which constant layer thickness start
212	!Config If =
213	!Config Def = DEPTH_MAX_H
214	!Config Help = Depth at which constant layer thickness start (smaller than zmaxh/2)
215	!Config Units = m
216	depth_cstthickness=zmaxh
217	CALL getin_p("DEPTH_CSTTHICK",depth_cstthickness)
218
219	IF ( (depth_cstthickness /= zmaxh) .AND. (depth_cstthickness > zmaxh/2.0) ) THEN
220	CALL ipslerr_p(2,'vertical_soil_init',"ERROR : depth_cstthickness needs to be smaller than zmaxh/2.0", &
221	"Correction : Constant thickness disabled", " ")
222	depth_cstthickness = zmaxh
223	ENDIF
224
225	! REFINEBOTTOM option is under work... This option can not be set from the parameter file for the moment.
226	! All the related code for refinebottom except this getin are found in this module but needs to be tested.
227	!!$ !Config Key = REFINEBOTTOM
228	!!$ !Config Desc = Depth at which the hydrology layers will be refined towards the bottom.
229	!!$ !Config If =
230	!!$ !Config Def = .FALSE.
231	!!$ !Config Help = Depth at which the hydrology layers will be refined towards the bottom.
232	!!$ ! This is important when lower boundary conditions is different from a free drainage.
233	!!$ !Config Units = -
234	!!$ !
235	refinebottom = .FALSE.
236	!!$ CALL getin_p("REFINEBOTTOM",refinebottom)
237
238	!Config Key = DEPTH_GEOM
239	!Config Desc = Depth at which we resume geometrical increases for temperature
240	!Config If =
241	!Config Def = DEPTH_MAX_H
242	!Config Help = Depth at which the thickness increases again for temperature.
243	! This depth has to be deeper than the bottom of the hydrology.
244	!Config Units = m
245	depth_geom=zmaxh
246	CALL getin_p("DEPTH_GEOM",depth_geom)
247
248	IF ( depth_geom < zmaxh ) THEN
249	CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", &
250	"Correction : setting depth_geom to zmaxh", " ")
251	ENDIF
252
253	!Config Key = RATIO_GEOM_BELOW
254	!Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM
255	!Config If =
256	!Config Def = 2
257	!Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM.
258	! This parameter allows to cover the depth needed for temperature with fewer layers.
259	!Config Units = -
260	! ratio_geom_below=2
261	!For carbon permafrost, to get ndeep=32 layers, we set DEPTH_TMAX=38. and RATIO_GEOM_BELOW=1.05
262	ratio_geom_below=1.05
263	CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below)
264
265	!
266	! Computing the layer depth for soil moisture. This defines the number of layers needed.
267	!
268	ztmp(:) = 0.0
269	dtmp(:) = 0.0
270	DO i=1,nblayermax
271	IF ( ztmp(i) < depth_cstthickness ) THEN
272	ztmp(i+1) = depth_topthickness2.0((2**i)-1)
273	dtmp(i+1) = ztmp(i+1)-ztmp(i)
274	igeo=i+1
275	zgeoend=ztmp(i+1)
276	dgeoend=dtmp(i+1)
277	ELSE
278	ztmp(i+1) = ztmp(i)+dtmp(i)
279	dtmp(i+1) = dtmp(i)
280	ENDIF
281	ENDDO
282	!
283	! refine at bottom. Work in progress...
284	!
285	nbrefine = 1
286	drefine(:) = 0.0
287	IF (refinebottom) THEN
288	!
289	! Compute parameters for the constant increment interval before refining,
290	! If needed !
291	cstint=zmaxh-(2.0*zgeoend)
292	nbcst=MAX(INT(cstint/dgeoend), 0)
293	IF ( nbcst > 0 ) THEN
294	dcst=cstint/nbcst
295	ELSE
296	dcst=dgeoend
297	ENDIF
298	!
299	! If we have to add constant increments
300	!
301	IF ( nbcst > 0 ) THEN
302	! Add linear levels
303	DO i=igeo,igeo+nbcst-1
304	ztmp(i+1) = ztmp(i)+dcst
305	dtmp(i+1) = dcst
306	ENDDO
307	! Refine the levels toward the bottom
308	DO i=igeo+nbcst,2*igeo+nbcst-2
309	irev=(2*igeo+nbcst-1)-i
310	ztmp(i+1) = ztmp(i)+dtmp(irev+1)
311	dtmp(i+1) = ztmp(i+1)-ztmp(i)
312	drefine(nbrefine)=dtmp(i+1)
313	nbrefine=nbrefine+1
314	ENDDO
315	! Without constant increments
316	ELSE
317	imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0)))
318	igeo=imin(1)
319	ztmp(igeo)=zmaxh/2.0
320	dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1)
321	DO i=igeo,2*igeo
322	irev=(2*igeo)-i
323	ztmp(i+1) = ztmp(i)+dtmp(irev)
324	dtmp(i+1) = ztmp(i-1)-ztmp(i)
325	drefine(nbrefine)=dtmp(i+1)
326	nbrefine=nbrefine+1
327	ENDDO
328	ENDIF
329	ENDIF
330
331	! Find the index (nslm) of the node closest to the zmaxh
332	imin=MINLOC(ABS(ztmp(:)-zmaxh))
333	nslm=imin(1)
334	!
335	! ALLOCATE the arrays we need to keep
336	!
337	ALLOCATE(znh(nslm), stat=ier)
338	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znh','','')
339	ALLOCATE(dnh(nslm), stat=ier)
340	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dnh','','')
341	ALLOCATE(dlh(nslm), stat=ier)
342	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlh','','')
343	ALLOCATE(zlh(nslm), stat=ier)
344	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlh','','')
345	!
346	! Extract now the values we need to reach zmaxh
347	!
348	znh(:)=ztmp(1:nslm)
349
350	! Force the last hydrological layer and node to be zmaxh
351	znh(nslm)=zmaxh
352	dnh(:)=dtmp(1:nslm)
353	! Recalculate the distance between the last 2 nodes
354	dnh(nslm) = ztmp(nslm)-ztmp(nslm-1)
355
356	! Calculate the thickness of the layers
357	DO i = 1, nslm-1
358	dlh(i) = (dnh(i) + dnh(i+1)) / 2.0
359	ENDDO
360	dlh(nslm) = dnh(nslm)/2
361
362	WRITE(numout,*) "=========================================================="
363	WRITE(numout,*) "Number of layers in hydrol = ", nslm
364	WRITE(numout,*) "znh = depth of nodes"
365	WRITE(numout,*) znh(:)
366	WRITE(numout,*) "dnh = internode distance"
367	WRITE(numout,*) dnh(:)
368	WRITE(numout,*) "dlh = layer thickness"
369	WRITE(numout,*) dlh(:)
370	WRITE(numout,*) "Total depth in hydrol has been calculated to ", znh(nslm)
371	WRITE(numout,*) "======================================================================"
372
373	!
374	! Extension for the thermodynamics below zmaxh
375	!
376	ntmp = nslm
377	ztmp(:)=0.0
378	! Exception for the first layer. It is at the top in the hydrology and in
379	! the middle for temperature.
380	ztmp(1)=depth_topthickness/2
381	DO i=2,ntmp
382	ztmp(i)=znh(i)
383	ENDDO
384
385	! Exception for the last layer where again temperature needs to be in
386	! the middle of the layer. Also add a layer with the same thickness
387	hh=dnh(ntmp)/2.0
388	ztmp(ntmp)=ztmp(ntmp)-hh/2.0
389	ztmp(ntmp+1)=ztmp(ntmp)+hh*1.5
390	ztmp(ntmp+2)=ztmp(ntmp+1)+hh*2.0
391	ntmp=ntmp+2
392
393	! If we have created a refined region at the bottom of the soil moisture. We need to
394	! unwinde it for temperature below zmaxh
395	IF ( nbrefine > 1 ) THEN
396	DO i=ntmp,ntmp+nbrefine
397	iref=(nbrefine+1)-(i-ntmp)
398	ztmp(i+1)=ztmp(i)+drefine(iref)
399	ENDDO
400	ntmp=ntmp+nbrefine
401	ENDIF
402	!
403	! Resume the geometric increas of thickness
404	DO i=ntmp,nblayermax
405	IF ( ztmp(i) < depth_geom ) THEN
406	ratio = 1.0
407	ELSE
408	ratio = ratio_geom_below
409	ENDIF
410	ztmp(i+1)=ztmp(i)+ratio*(ztmp(i)-ztmp(i-1))
411	ENDDO
412
413	! Compute the depth of the lower interface of each layer.
414	!
415	zint(1) = depth_topthickness
416	DO i=2,nblayermax-1
417	zint(i) = (ztmp(i)+ztmp(i+1))/2.0
418	ENDDO
419	zint(nslm-1) = (znh(nslm-1) + znh(nslm))/2.0
420	zint(nslm) = (znh(nslm))
421	zint(nblayermax) = ztmp(nblayermax)+(ztmp(nblayermax)-ztmp(nblayermax-1))/2.0
422
423	! Determine the total number of layers for thermal (ngrnd).
424	!
425	imin=MINLOC(ABS(zint(:)-zmaxt))
426	ngrnd=imin(1)
427	ALLOCATE(znt(ngrnd), stat=ier)
428	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znt','','')
429	ALLOCATE(dlt(ngrnd), stat=ier)
430	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlt','','')
431	ALLOCATE(zlt(ngrnd), stat=ier)
432	IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlt','','')
433
434	! Assign values for znt, zlt, dlt
435	!
436	znt(:)=ztmp(1:ngrnd)
437	zlt(:) = zint(1:ngrnd)
438	dlt(1) = zint(1)
439	DO i=2,ngrnd
440	dlt(i) = zint(i)-zint(i-1)
441	ENDDO
442	! Depth of layers for the hydrology are the same as for thermosoil but only the upper nslm layers are used
443	zlh(:) = zlt(1:nslm)
444
445	! Force the last thermal layer and node to be zmaxt
446	zlt(ngrnd) = zmaxt
447	dlt(ngrnd) = zmaxt-zint(ngrnd-1)
448
449	WRITE(numout,*) "Number of layers in thermosoil = ", ngrnd
450	WRITE(numout,*) "The maximal depth for thermosoil (DEPTH_MAX_T) = ", zlt(ngrnd)
451	WRITE(numout,*) "to be compared with calculated maximal depth for thermosoil (not used) = ", zint(ngrnd)
452	WRITE(numout,*) "Depth of the nodes in thermosoil, znt=",znt
453	WRITE(numout,*) "Thickness of the layers, dlt=", dlt
454	WRITE(numout,*) "Depth of the interface between the layers, zlt=", zlt(1:ngrnd)
455	WRITE(numout,*) "======================================================================"
456
457	END SUBROUTINE vertical_soil_init
458
459	END MODULE vertical_soil

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

Download in other formats: