Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/testrouting.f90 @ 7262

Last change on this file since 7262 was 5559, checked in by josefine.ghattas, 6 years ago
Enhancement on grid module: merge variables grid_type and GridType? meaning the same. Only grid_type is used now, this one is choosen because it is the same as used in LMDZ.
File size: 18.4 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	!
109	!_ ================================================================================================================================
110	!
111	!
112	CALL Init_orchidee_para()
113	!
114	CALL getlandseamask_init(iim, jjm, nbindex)
115	ALLOCATE(lon(iim,jjm))
116	ALLOCATE(lat(iim,jjm))
117	ALLOCATE(orog(iim,jjm))
118	ALLOCATE(contfrac_2d(iim,jjm))
119	CALL getlandseamask_read(lon, lat, contfrac_2d, orog)
120	!
121	! ALLOCATE memory needed
122	!
123	ALLOCATE(kindex_g(nbindex))
124	ALLOCATE(lalo_land(nbindex,2))
125	ALLOCATE(contfrac_land(nbindex))
126	ALLOCATE(orog_land(nbindex))
127	!
128	!
129	!
130	ni=0
131	DO j=1,jjm
132	DO i=1,iim
133	IF ( contfrac_2d(i,j) > 0.0 ) THEN
134	ni = ni + 1
135	IF ( ni .GT. nbindex ) THEN
136	WRITE(,) "We are expecting ", nbindex, "point."
137	WRITE(,) "We are at : ", i, j, orog(i,j)
138	STOP 'Too many continental points'
139	ENDIF
140	kindex_g(ni) = (j-1)*iim + i
141	lalo_land(ni,1) = lat(i,j)
142	lalo_land(ni,2) = lon(i,j)
143	contfrac_land(ni) = contfrac_2d(i,j)
144	orog_land(ni) = orog(i,j)
145	ENDIF
146	ENDDO
147	ENDDO
148	!
149	!
150	nbseg = 4
151	!
152	!
153	CALL grid_set_glo(iim, jjm, nbindex)
154	CALL grid_allocate_glo(nbseg)
155	!
156	CALL bcast(nbindex)
157	ALLOCATE(index_g(nbindex))
158	IF ( is_root_prc ) index_g(:) = kindex_g(:)
159	CALL bcast(index_g)
160	!
161	WRITE(,) "GOING INTO Init_orchidee_data_para_driver", nbindex, index_g(1), SIZE(kindex_g)
162	CALL Init_orchidee_data_para_driver(nbindex, index_g)
163	WRITE(,) "OUT OF Init_orchidee_data_para_driver"
164	CALL init_ioipsl_para
165	!
166	WRITE(,) mpi_rank, "DIMENSIONS of grid on processor : iim, jjm, nbindex = ", iim, jjm, nbindex, nbp_loc
167	!
168	CALL grid_init (nbp_loc, nbseg, regular_lonlat, "ForcingGrid")
169	!
170	!==========================================================================
171	!
172	! Transfer the global grid variables to the root proc
173	! _glo -> _g
174	! Variables *_g were allocated with the CALL init_grid
175	!
176	IF ( is_root_prc) THEN
177	!
178	lalo_g(:,:) = lalo_land(:,:)
179	lon_g(:,:) = lon(:,:)
180	lat_g(:,:) = lat(:,:)
181	contfrac_g(:) = contfrac_land(:)
182	!
183	ENDIF
184	!
185	CALL grid_stuff(nbindex, iim, jjm, lon_g, lat_g, index_g)
186	!
187	!
188	! Distribute the grid to all processors
189	!
190	! Redistribute the indeces on all procs (apple distribution of land points)
191	!
192	ALLOCATE(kindex(nbp_loc))
193	ALLOCATE(orog_loc(nbp_loc), contfrac_loc(nbp_loc))
194	CALL bcast(lon_g)
195	CALL bcast(lat_g)
196	CALL scatter(index_g, kindex)
197	CALL scatter(lalo_land, lalo)
198	CALL scatter(orog_land, orog_loc)
199	CALL scatter(contfrac_land, contfrac_loc)
200	!
201	!
202	! Apply the offset needed so that kindex refers to the index of the land point
203	! on the current region, i.e. the local lon lat domain.
204	!
205	kindex(1:nbp_loc)=kindex(1:nbp_loc)-(jj_begin-1)*iim
206	!
207	!
208	!==========================================================================================
209	!
210	! The grid is in place and we can start to prepare the time of integration.
211	!
212	CALL ioconf_calendar("gregorian")
213	!
214	!
215	! Determine initial step or restart one
216	!
217	!Config Key = RESTART_IN
218	!Config Desc = Name of restart file to read at restart
219	!Config If = [-]
220	!Config Def = NONE
221	!Config Help = This function allows to select a restart file with which
222	! the simulation will be initialized.
223	!Config Units = [-]
224	!-
225	!
226	restname_in="NONE"
227	CALL getin('RESTART_IN', restname_in)
228	!
229	!Config Key = RESTART_OUT
230	!Config Desc = Name of restart file to be written at the end of the simulation.
231	!Config If = [-]
232	!Config Def = NONE
233	!Config Help = This function allows to select a restart file which will be written by the
234	! model and which can be used as input for a future restart.
235	!Config Units = [-]
236	!-
237	!
238	restname_out="restart_out.nc"
239	CALL getin('RESTART_OUT', restname_out)
240	!
241	!Config Key = SIMULATION_LEN
242	!Config Desc = Time step length in years for "testrouting"
243	!Config If = [-]
244	!Config Def = 1
245	!Config Help = This is time step length for testrouting
246	!Config Units = [-]
247	nbyears=1.0
248	CALL getin('SIMULATION_LEN', nbyears)
249	!
250	!Config Key = DTRADIA
251	!Config Desc = Time step length for "testrouting"
252	!Config If = [-]
253	!Config Def = 1800.
254	!Config Help = This is time step length for testrouting
255	!Config Units = [-]
256	!-
257	!DTRADIA = 1800.
258	dtradia = 1800.
259	CALL getin('DTRADIA', dtradia)
260	!
261	!- Initial date
262	CALL ymds2ju (2000,1,1,0.0, date0)
263	date0_rest = date0
264	!
265	dt = dtradia
266	dw = 86400.
267	!
268	CALL control_initialize
269	!
270	CALL ioget_calendar(one_year,one_day_loc)
271	!
272	! We have all we need and we can start to work
273	!
274	!
275	IF (is_root_prc) THEN
276	CALL restini(restname_in, iim, jjm, lon_g, lat_g, 1, lev, &
277	& restname_out, ibegt, date0_rest, dtradia, rest_id, .FALSE.)
278	ELSE
279	rest_id=0
280	ENDIF
281	CALL bcast (ibegt)
282	CALL bcast (date0_rest)
283	CALL bcast (dtradia)
284	!
285	!
286	IF ( INDEX(restname_in, "NONE") > 0 ) THEN
287	kjit = 1
288	ibegt = 1
289	ELSE
290	kjit = ibegt
291	date0 = date0_rest
292	ENDIF
293	WRITE(,) 'Out of restini : kjit=',kjit, " ibegt=", ibegt, " date0=", date0
294	!
295	!- time step length
296	!
297	! Set up the history file
298	!
299	!Config Key = HISTNAME
300	!Config Desc = Name of the history file
301	!Config If = [-]
302	!Config Def = out_testrouting
303	!Config Help = The name of the file which will contain all the diagnostics of routing.
304	!Config Units = [-]
305	!-
306	WRITE(flux_op,'("ave(scatter(X*",F8.1,"))")') one_day_loc/dt
307	WRITE(flux_scinsec,'("ave(scatter(X*",F8.6,"))")') 1.0/dt
308	avescatter = 'ave(scatter(X))'
309	once_wrt = 'once(scatter(X))'
310	!
311	histname="out_testrouting"
312	CALL getin('HISTNAME', histname)
313	!
314	CALL histbeg(histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
315	& kjit-1, date0, dtradia, hori_id, hist_id, domain_id=orch_domain_id)
316	!
317	CALL histdef(hist_id, 'Orog', 'Orography', ' ', &
318	& iim, jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
319	CALL histdef(hist_id, 'Contfrac', 'Fraction of continent', ' ', &
320	& iim, jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
321	CALL histdef(hist_id, 'Areas', 'Mesh areas', 'm2', &
322	& iim,jjm, hori_id, 1,1,1, -99, 32, once_wrt, dt, dw)
323	!
324	CALL histdef(hist_id, 'riversret', 'Return from endorheic rivers', 'mm/d', &
325	& iim,jjm, hori_id, 1,1,1, -99, 32, flux_op, dt,dw)
326	CALL histdef(hist_id, 'hydrographs', 'Hydrographs of gridbox outflow', 'm^3/s', &
327	& iim,jjm, hori_id, 1,1,1, -99, 32, flux_scinsec, dt,dw)
328	!
329	CALL histdef(hist_id, 'fastr', 'Fast flow reservoir', 'kg/m^2', &
330	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
331	CALL histdef(hist_id, 'slowr', 'Slow flow reservoir', 'kg/m^2', &
332	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
333	CALL histdef(hist_id, 'streamr', 'Stream flow reservoir', 'kg/m^2', &
334	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
335	CALL histdef(hist_id, 'pondr', 'Volume in pond reservoir', 'kg/m^2', &
336	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
337	CALL histdef(hist_id, 'lakevol', 'Volume in lake reservoir', 'kg/m^2', &
338	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
339	!
340	CALL histdef(hist_id, 'basinmap', 'Aproximate map of the river basins', ' ', &
341	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
342	CALL histdef(hist_id, 'nbrivers', 'Number or rivers in the outflow grid box', ' ', &
343	& iim,jjm, hori_id, 1,1,1, -99, 32, avescatter, dt,dw)
344	!
345	CALL histend(hist_id)
346	!
347	! Put a copy of the orography into the restart
348	!
349	!
350	CALL histwrite_p(hist_id, 'Orog', kjit+1, orog_loc, nbp_loc, kindex)
351	CALL histwrite_p(hist_id, 'Contfrac', kjit+1, contfrac_loc, nbp_loc, kindex)
352	CALL histwrite_p(hist_id, 'Areas', kjit+1, area, nbp_loc, kindex)
353	!
354	! Override some settings as testrouting is not as flexible as the rest of ORCHIDEE.
355	!
356	hist2_id=-1
357	almaoutput=.FALSE.
358	!================================================================================
359	!
360	! Set up the routing schemes
361	!
362	!
363	! Allocate all the physical variables
364	!
365	! Input variables
366	ALLOCATE(runoff(nbp_loc), drainage(nbp_loc), humrel(nbp_loc), stempdiag(nbp_loc,nslm), transpot(nbp_loc,nvmc))
367	ALLOCATE(totfrac_nobio(nbp_loc), veget_max(nbp_loc,nvmc), floodout(nbp_loc))
368	ALLOCATE(precip_rain(nbp_loc), k_litt(nbp_loc),reinf_slope(nbp_loc))
369	ALLOCATE(evapot_corr(nbp_loc))
370	! Output variables
371	ALLOCATE(returnflow(nbp_loc), irrigation(nbp_loc), riverflow(nbp_loc), coastalflow(nbp_loc), reinfiltration(nbp_loc))
372	ALLOCATE(flood_frac(nbp_loc), flood_res(nbp_loc))
373	!
374	! Get some fake value for input arrays
375	!
376	runoff(:) = 1.0
377	drainage(:) = 1.0
378	humrel(:) = 0.75
379	stempdiag(:,:) = 273.5
380	transpot(:,:)=0.0
381	reinfiltration(:)=0.0
382	flood_frac(:)=0.0
383	flood_res(:)=0.0
384	totfrac_nobio = 0.1
385	veget_max = 0.2
386	floodout = 0.0
387	precip_rain = 0.0
388	k_litt = 0.0
389	reinf_slope = 0.1
390	evapot_corr(:) = 10.0
391	!
392	!
393	!
394	CALL routing_initialize( kjit, nbp_loc, kindex, &
395	& rest_id, hist_id, hist2_id, lalo, &
396	& neighbours, resolution, contfrac, stempdiag, &
397	& returnflow, reinfiltration, irrigation, riverflow, &
398	& coastalflow, flood_frac, flood_res)
399	!
400	! Do loop over a number of time-steps
401	!
402	simlen = NINT(nbyears365one_day_loc/dtradia)
403	ibegt=kjit
404	iendt=kjit+simlen
405	WRITE(,) "The simulation will go from ", ibegt, " to ", iendt
406	!
407	DO kjit = ibegt,iendt
408	!
409	date = date0 + (kjit-1)*(dtradia/one_day_loc)
410	!
411	IF ( date < date0+1 ) THEN
412	! During one day one kg/m^2d divided up in runoff and drainage
413	runoff(:) = 0.5/48.
414	drainage(:) = 0.5/48.
415	ELSE
416	runoff(:) = 0.0
417	drainage(:) = 0.0
418	ENDIF
419	!
420	CALL routing_main(kjit, nbp_loc, kindex, &
421	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
422	& drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
423	& stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
424	!
425	WRITE(,) "Out of routing at time step = ",kjit,' Seconds since start', (kjit-ibegt)*dtradia
426	!
427	ENDDO
428	!
429	! Shut everything down
430	!
431	CALL routing_finalize(kjit, nbp_loc, rest_id, flood_frac, flood_res)
432	!
433	!
434	CALL histclo
435	IF ( is_root_prc) THEN
436	CALL restclo
437	ENDIF
438	!
439	CALL Finalize_mpi
440	!
441	END PROGRAM testrouting
442	!

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

Download in other formats: