Context Navigation

testrouting.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

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

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_driver/testrouting.f90 @ 6737

Download in other formats: