Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_driver/testrouting.f90

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

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

Download in other formats: