Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

testrouting.f90 @ 7442

Last change on this file since 7442 was 5080, checked in by chunjing.qiu, 6 years ago
soil freezing, soil moisture, fwet bugs fixed
File size: 19.0 KB

Line
1	! =================================================================================================================================
2	! PROGRAM : testrouting
3	!
4	! CONTACT : jan.polcher _at_ lmd.jussieu.fr
5	!
6	! LICENCE : :)
7	!
8	!>\BRIEF This program tests routing scheme (from routing.f90) which routes the water over the continents
9	!! into the oceans and computes the water stored in floodplains or taken for irrigation.
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! REFERENCE(S) :
16	!!
17	!! SVN :
18	!! $HeadURL : svn://forge.ipsl.jussieu.fr/orchidee/trunk/ORCHIDEE/[somewhere]/testrouting.f90 $
19	!! $Date: 2014-10-27 10:39:00 +0200 (Mon, 27 Oct 2014) $
20	!! $Revision: XXXX $
21	!! \n
22	!_ ================================================================================================================================
23	PROGRAM testrouting
24	!
25	USE ioipsl_para
26	USE pft_parameters
27	USE mod_orchidee_para
28	USE control
29	USE constantes_soil_var
30	USE constantes_var
31	USE constantes
32	USE routing
33	USE timer
34	USE grid
35	!
36	USE getlandseamask
37	!
38	IMPLICIT NONE
39	!
40	INTEGER(i_std) :: nbseg !!
41	!
42	INTEGER(i_std) :: iim !! Size in longitude of coarser grid
43	INTEGER(i_std) :: jjm !! Size in latitude of coarser grid
44	INTEGER(i_std) :: i, j !! Integer variable for loops
45	INTEGER(i_std) :: ibegt, iendt
46	INTEGER(i_std) :: ni !! For checking nbindex
47	REAL(r_std) :: nbyears !! Lenght of simulation in years
48	INTEGER(i_std) :: simlen !! Lenght of simulation: simlen = 36548nbyears
49	REAL(r_std) :: dx, dy !! Lon/Lat resolution of coarser grid
50	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: lon, lat !! Lon/lat of coarser grid
51	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: orog !! New orography after interpolation
52	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: orog_land, orog_loc
53	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: lalo_land
54	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: contfrac_land, contfrac_loc
55	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: contfrac_2d
56	!
57	REAL(r_std), DIMENSION (1) :: lev !! Number of level (requested by restini routine) (unitless)
58	CHARACTER(LEN=80) :: histname !! Name of history file (can not find HISTNAME?)
59	INTEGER(i_std) :: hori_id !! ID of the default horizontal longitude and latitude map.
60	INTEGER(i_std) :: hist_id !! History file identification for ???
61	INTEGER(i_std) :: rest_id !! ID of the restart file
62	!
63	REAL(r_std) :: date0 !! Initial date
64	REAL(r_std) :: date0_rest !! Initial date from restart file
65	REAL(r_std) :: date !! Current date
66	REAL(r_std) :: dt !! Same as dtradia ???
67	REAL(r_std) :: dw !! 86400. ???
68	REAL(r_std) :: one_day_loc
69	!
70	CHARACTER(LEN=40) :: flux_op !! Operations to be performed on fluxes
71	CHARACTER(LEN=40) :: flux_scinsec !! Operation in seconds
72	CHARACTER(LEN=40) :: avescatter, once_wrt !! The various operation to be performed
73	!
74	! Input for routing_main
75	!
76	INTEGER(i_std) :: kjit !! Time step number
77	INTEGER(i_std) :: nbindex !! Number of local continental points
78	REAL(r_std) :: dtradia !! Timestep length
79	!
80	INTEGER(i_std), ALLOCATABLE, DIMENSION (:) :: kindex_g !! Index of land point on 2D map (in local position)
81	INTEGER(i_std), ALLOCATABLE, DIMENSION (:) :: kindex !! index of land point per proc
82	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: runoff !! Grid-point runoff (kg/m^2/dt)
83	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: drainage !! Grid-point drainage (kg/m^2/dt)
84	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: humrel !! Soil moisture stress, root extraction potential (unitless)
85	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: stempdiag !! Diagnostic soil temperature profile
86	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
87	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max !! Max. fraction of vegetation type (LAI -> infty, unitless)
88	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: floodout !! Flow out of floodplains from hydrol
89	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpot !! Potential Transpiration (needed for irrigation)
90	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: precip_rain !! Rainfall (kg/m^2/dt)
91	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: k_litt !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
92	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinf_slope !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
93	INTEGER(i_std) :: hist2_id !! Access to history file 2 (unitless)
94	!
95	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: evapot_corr !! Soil Potential Evaporation
96	!
97	! Output from routing_main
98	!
99	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: returnflow !! The water flow from lakes and swamps which returns to the grid box (kg/m^2/dt)
100	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: irrigation !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
101	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: riverflow !! Outflow of the major rivers, will be located on the continental grid but this should be a coastal point (kg/dt)
102	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: coastalflow !! Outflow on coastal points by small basins, the water which flows in a disperse way into the ocean (kg/dt)
103	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
104	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_res !! Diagnostic of water amount in the floodplains reservoir (kg)
105	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_frac !! Flooded fraction of the grid box (unitless;0-1)
106	!
107	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soil_deficit !! water deficit to reach IRRIG_FULFILL
108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget !! Fraction of vegetation type (unitless, 0-1)
109	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vegstress !! Vegetation moisture stress (only for vegetation growth)
110	!
111	!
112	!_ ================================================================================================================================
113	!
114	!
115	CALL Init_orchidee_para()
116	!
117	CALL getlandseamask_init(iim, jjm, nbindex)
118	ALLOCATE(lon(iim,jjm))
119	ALLOCATE(lat(iim,jjm))
120	ALLOCATE(orog(iim,jjm))
121	ALLOCATE(contfrac_2d(iim,jjm))
122	CALL getlandseamask_read(lon, lat, contfrac_2d, orog)
123	!
124	! ALLOCATE memory needed
125	!
126	ALLOCATE(kindex_g(nbindex))
127	ALLOCATE(lalo_land(nbindex,2))
128	ALLOCATE(contfrac_land(nbindex))
129	ALLOCATE(orog_land(nbindex))
130	!
131	!
132	!
133	ni=0
134	DO j=1,jjm
135	DO i=1,iim
136	IF ( contfrac_2d(i,j) > 0.0 ) THEN
137	ni = ni + 1
138	IF ( ni .GT. nbindex ) THEN
139	WRITE(,) "We are expecting ", nbindex, "point."
140	WRITE(,) "We are at : ", i, j, orog(i,j)
141	STOP 'Too many continental points'
142	ENDIF
143	kindex_g(ni) = (j-1)*iim + i
144	lalo_land(ni,1) = lat(i,j)
145	lalo_land(ni,2) = lon(i,j)
146	contfrac_land(ni) = contfrac_2d(i,j)
147	orog_land(ni) = orog(i,j)
148	ENDIF
149	ENDDO
150	ENDDO
151	!
152	!
153	nbseg = 4
154	!
155	!
156	CALL grid_set_glo(iim, jjm, nbindex)
157	CALL grid_allocate_glo(nbseg)
158	!
159	CALL bcast(nbindex)
160	ALLOCATE(index_g(nbindex))
161	IF ( is_root_prc ) index_g(:) = kindex_g(:)
162	CALL bcast(index_g)
163	!
164	WRITE(,) "GOING INTO Init_orchidee_data_para_driver", nbindex, index_g(1), SIZE(kindex_g)
165	CALL Init_orchidee_data_para_driver(nbindex, index_g)
166	WRITE(,) "OUT OF Init_orchidee_data_para_driver"
167	CALL init_ioipsl_para
168	!
169	WRITE(,) mpi_rank, "DIMENSIONS of grid on processor : iim, jjm, nbindex = ", iim, jjm, nbindex, nbp_loc
170	!
171	CALL grid_init (nbp_loc, nbseg, "RegLonLat", "ForcingGrid")
172	!
173	!==========================================================================
174	!
175	! Transfer the global grid variables to the root proc
176	! _glo -> _g
177	! Variables *_g were allocated with the CALL init_grid
178	!
179	IF ( is_root_prc) THEN
180	!
181	lalo_g(:,:) = lalo_land(:,:)
182	lon_g(:,:) = lon(:,:)
183	lat_g(:,:) = lat(:,:)
184	contfrac_g(:) = contfrac_land(:)
185	!
186	ENDIF
187	!
188	CALL grid_stuff(nbindex, iim, jjm, lon_g, lat_g, index_g)
189	!
190	!
191	! Distribute the grid to all processors
192	!
193	! Redistribute the indeces on all procs (apple distribution of land points)
194	!
195	ALLOCATE(kindex(nbp_loc))
196	ALLOCATE(orog_loc(nbp_loc), contfrac_loc(nbp_loc))
197	CALL bcast(lon_g)
198	CALL bcast(lat_g)
199	CALL scatter(index_g, kindex)
200	CALL scatter(lalo_land, lalo)
201	CALL scatter(orog_land, orog_loc)
202	CALL scatter(contfrac_land, contfrac_loc)
203	!
204	!
205	! Apply the offset needed so that kindex refers to the index of the land point
206	! on the current region, i.e. the local lon lat domain.
207	!
208	kindex(1:nbp_loc)=kindex(1:nbp_loc)-(jj_begin-1)*iim
209	!
210	!
211	!==========================================================================================
212	!
213	! The grid is in place and we can start to prepare the time of integration.
214	!
215	CALL ioconf_calendar("gregorian")
216	!
217	!
218	! Determine initial step or restart one
219	!
220	!Config Key = RESTART_IN
221	!Config Desc = Name of restart file to read at restart
222	!Config If = [-]
223	!Config Def = NONE
224	!Config Help = This function allows to select a restart file with which
225	! the simulation will be initialized.
226	!Config Units = [-]
227	!-
228	!
229	restname_in="NONE"
230	CALL getin('RESTART_IN', restname_in)
231	!
232	!Config Key = RESTART_OUT
233	!Config Desc = Name of restart file to be written at the end of the simulation.
234	!Config If = [-]
235	!Config Def = NONE
236	!Config Help = This function allows to select a restart file which will be written by the
237	! model and which can be used as input for a future restart.
238	!Config Units = [-]
239	!-
240	!
241	restname_out="restart_out.nc"
242	CALL getin('RESTART_OUT', restname_out)
243	!
244	!Config Key = SIMULATION_LEN
245	!Config Desc = Time step length in years for "testrouting"
246	!Config If = [-]
247	!Config Def = 1
248	!Config Help = This is time step length for testrouting
249	!Config Units = [-]
250	nbyears=1.0
251	CALL getin('SIMULATION_LEN', nbyears)
252	!
253	!Config Key = DTRADIA
254	!Config Desc = Time step length for "testrouting"
255	!Config If = [-]
256	!Config Def = 1800.
257	!Config Help = This is time step length for testrouting
258	!Config Units = [-]
259	!-
260	!DTRADIA = 1800.
261	dtradia = 1800.
262	CALL getin('DTRADIA', dtradia)
263	!
264	!- Initial date
265	CALL ymds2ju (2000,1,1,0.0, date0)
266	date0_rest = date0
267	!
268	dt = dtradia
269	dw = 86400.
270	!
271	CALL control_initialize(dt)
272	!
273	CALL ioget_calendar(one_year,one_day_loc)
274	!
275	! We have all we need and we can start to work
276	!
277	!
278	IF (is_root_prc) THEN
279	CALL restini(restname_in, iim, jjm, lon_g, lat_g, 1, lev, &
280	& restname_out, ibegt, date0_rest, dtradia, rest_id, .FALSE.)
281	ELSE
282	rest_id=0
283	ENDIF
284	CALL bcast (ibegt)
285	CALL bcast (date0_rest)
286	CALL bcast (dtradia)
287	!
288	!
289	IF ( INDEX(restname_in, "NONE") > 0 ) THEN
290	kjit = 1
291	ibegt = 1
292	ELSE
293	kjit = ibegt
294	date0 = date0_rest
295	ENDIF
296	WRITE(,) 'Out of restini : kjit=',kjit, " ibegt=", ibegt, " date0=", date0
297	!
298	!- time step length
299	!
300	! Set up the history file
301	!
302	!Config Key = HISTNAME
303	!Config Desc = Name of the history file
304	!Config If = [-]
305	!Config Def = out_testrouting
306	!Config Help = The name of the file which will contain all the diagnostics of routing.
307	!Config Units = [-]
308	!-
309	WRITE(flux_op,'("ave(scatter(X*",F8.1,"))")') one_day_loc/dt
310	WRITE(flux_scinsec,'("ave(scatter(X*",F8.6,"))")') 1.0/dt
311	avescatter = 'ave(scatter(X))'
312	once_wrt = 'once(scatter(X))'
313	!
314	histname="out_testrouting"
315	CALL getin('HISTNAME', histname)
316	!
317	CALL histbeg(histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
318	& kjit-1, date0, dtradia, hori_id, hist_id, domain_id=orch_domain_id)
319	!
320	CALL histdef(hist_id, 'Orog', 'Orography', ' ', &
321	& iim, jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
322	CALL histdef(hist_id, 'Contfrac', 'Fraction of continent', ' ', &
323	& iim, jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
324	CALL histdef(hist_id, 'Areas', 'Mesh areas', 'm2', &
325	& iim,jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
326	!
327	CALL histdef(hist_id, 'riversret', 'Return from endorheic rivers', 'mm/d', &
328	& iim,jjm, hori_id, 1,1,1, -99, 32, flux_op, dt,dw)
329	CALL histdef(hist_id, 'hydrographs', 'Hydrographs of gridbox outflow', 'm^3/s', &
330	& iim,jjm, hori_id, 1,1,1, -99, 32, flux_scinsec, dt,dw)
331	!
332	CALL histdef(hist_id, 'fastr', 'Fast flow reservoir', 'kg/m^2', &
333	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
334	CALL histdef(hist_id, 'slowr', 'Slow flow reservoir', 'kg/m^2', &
335	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
336	CALL histdef(hist_id, 'streamr', 'Stream flow reservoir', 'kg/m^2', &
337	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
338	CALL histdef(hist_id, 'pondr', 'Volume in pond reservoir', 'kg/m^2', &
339	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
340	CALL histdef(hist_id, 'lakevol', 'Volume in lake reservoir', 'kg/m^2', &
341	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
342	!
343	CALL histdef(hist_id, 'basinmap', 'Aproximate map of the river basins', ' ', &
344	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
345	CALL histdef(hist_id, 'nbrivers', 'Number or rivers in the outflow grid box', ' ', &
346	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
347	!
348	CALL histend(hist_id)
349	!
350	! Put a copy of the orography into the restart
351	!
352	!
353	CALL histwrite_p(hist_id, 'Orog', kjit+1, orog_loc, nbp_loc, kindex)
354	CALL histwrite_p(hist_id, 'Contfrac', kjit+1, contfrac_loc, nbp_loc, kindex)
355	CALL histwrite_p(hist_id, 'Areas', kjit+1, area, nbp_loc, kindex)
356	!
357	! Override some settings as testrouting is not as flexible as the rest of ORCHIDEE.
358	!
359	hist2_id=-1
360	almaoutput=.FALSE.
361	!================================================================================
362	!
363	! Set up the routing schemes
364	!
365	!
366	! Allocate all the physical variables
367	!
368	! Input variables
369	ALLOCATE(runoff(nbp_loc), drainage(nbp_loc), humrel(nbp_loc), stempdiag(nbp_loc,nslm), transpot(nbp_loc,nvmc))
370	ALLOCATE(totfrac_nobio(nbp_loc), veget_max(nbp_loc,nvmc), floodout(nbp_loc))
371	ALLOCATE(precip_rain(nbp_loc), k_litt(nbp_loc),reinf_slope(nbp_loc))
372	ALLOCATE(evapot_corr(nbp_loc))
373	! Output variables
374	ALLOCATE(returnflow(nbp_loc), irrigation(nbp_loc), riverflow(nbp_loc), coastalflow(nbp_loc), reinfiltration(nbp_loc))
375	ALLOCATE(flood_frac(nbp_loc), flood_res(nbp_loc))
376	!
377	ALLOCATE(vegstress(nbp_loc, nvmc))
378	ALLOCATE(veget(nbp_loc, nvmc))
379	ALLOCATE(soil_deficit(nbp_loc, nvmc))
380	!
381	! Get some fake value for input arrays
382	!
383	runoff(:) = 1.0
384	drainage(:) = 1.0
385	humrel(:) = 0.75
386	stempdiag(:,:) = 273.5
387	transpot(:,:)=0.0
388	reinfiltration(:)=0.0
389	flood_frac(:)=0.0
390	flood_res(:)=0.0
391	totfrac_nobio = 0.1
392	veget_max = 0.2
393	floodout = 0.0
394	precip_rain = 0.0
395	k_litt = 0.0
396	reinf_slope = 0.1
397	evapot_corr(:) = 10.0
398	!
399	vegstress = 0
400	veget = 0
401	soil_deficit = 0
402	!
403	!
404	!
405	!
406	CALL routing_initialize( kjit, nbp_loc, kindex, &
407	& rest_id, hist_id, hist2_id, lalo, &
408	& neighbours, resolution, contfrac, stempdiag, &
409	& returnflow, reinfiltration, irrigation, riverflow, &
410	& coastalflow, flood_frac, flood_res)
411	!
412	! Do loop over a number of time-steps
413	!
414	simlen = NINT(nbyears365one_day_loc/dtradia)
415	ibegt=kjit
416	iendt=kjit+simlen
417	WRITE(,) "The simulation will go from ", ibegt, " to ", iendt
418	!
419	DO kjit = ibegt,iendt
420	!
421	date = date0 + (kjit-1)*(dtradia/one_day_loc)
422	!
423	IF ( date < date0+1 ) THEN
424	! During one day one kg/m^2d divided up in runoff and drainage
425	runoff(:) = 0.5/48.
426	drainage(:) = 0.5/48.
427	ELSE
428	runoff(:) = 0.0
429	drainage(:) = 0.0
430	ENDIF
431	!
432	vegstress = 0
433	veget = 0
434	evapot_corr = 0
435	soil_deficit = 0
436	!
437	CALL routing_main(kjit, nbp_loc, kindex, &
438	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget, veget_max, soil_deficit, floodout, runoff, &
439	& drainage, transpot, evapot_corr, vegstress, precip_rain, humrel, k_litt, flood_frac, flood_res, &
440	& stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
441	!
442	WRITE(,) "Out of routing at time step = ",kjit,' Seconds since start', (kjit-ibegt)*dtradia
443	!
444	ENDDO
445	!
446	! Shut everything down
447	!
448	CALL routing_finalize(kjit, nbp_loc, rest_id, flood_frac, flood_res)
449	!
450	!
451	CALL histclo
452	IF ( is_root_prc) THEN
453	CALL restclo
454	ENDIF
455	!
456	CALL Finalize_mpi
457	!
458	END PROGRAM testrouting
459	!

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_driver/testrouting.f90 @ 7442

Download in other formats: