Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_oasisdriver/orchoasis_tools.f90 @ 6145

Last change on this file since 6145 was 4503, checked in by jan.polcher, 7 years ago
Some bug correction on the transmission of the rivers to driver2oasis. There is now also protection built in so that it can run without routing.
File size: 33.8 KB

Line
1	MODULE orchoasis_tools
2	!
3	USE defprec
4	USE netcdf
5	!
6	USE ioipsl_para
7	USE mod_orchidee_para
8	!
9	#ifdef OASIS
10	USE mod_oasis
11	!
12	! ORCHIDEE definitions
13	!
14	USE constantes
15	USE constantes_soil
16	USE constantes_var
17	USE constantes_mtc
18	USE pft_parameters
19	!
20	USE control
21	USE routing_reg, ONLY : routing_reg_riverloc, routing_reg_initialize, routing_reg_clear
22	!
23	IMPLICIT NONE
24	!
25	PRIVATE
26	PUBLIC :: orchoasis_time, &
27	& orchoasis_printpoint, orchoasis_defvar, orchoasis_defvecvar, &
28	& orchoasis_putvar, orchoasis_putvecvar, orchoasis_printdate, &
29	& orchoasis_getvar, orchoasis_routinginit, orchoasis_wrtgrid
30	!
31	!
32	! Global variables
33	!
34	!
35	INTEGER(i_std), SAVE :: itau_offset !! This offset is used to phase the
36	!
37	! OASIS PUT VARIABLES
38	!
39	INTEGER(i_std), SAVE :: il_part_id, tair_id, qair_id, zlevtq_id, zlevuv_id
40	INTEGER(i_std), SAVE :: rainf_id, snowf_id, lwdown_id, swnet_id, solarang_id
41	INTEGER(i_std), SAVE :: u_id, v_id, ps_id, cdrag_id
42	!
43	INTEGER(i_std), SAVE :: river_part_id
44	INTEGER(i_std), SAVE :: riverlon_id, riverlat_id, riverflow_id
45	!
46	! OASIS GET VARIABLES
47	!
48	INTEGER(i_std), SAVE :: vevapp_id, fluxsens_id, fluxlat_id, coastal_id, river_id
49	INTEGER(i_std), SAVE :: netco2_id, carblu_id, tsolrad_id, tsolnew_id, qsurf_id
50	INTEGER(i_std), SAVE :: albnir_id, albvis_id, emis_id, z0m_id, z0h_id
51	!
52	LOGICAL, PARAMETER :: check_INPUTS = .FALSE. !! (very) long print of INPUTs in intersurf
53	LOGICAL, SAVE :: check_time = .FALSE.
54	!
55	LOGICAL, SAVE :: landonly = .TRUE.
56	!
57	PUBLIC check_time
58	!
59	CONTAINS
60	!
61	!=============================================================================================
62	!
63	SUBROUTINE orchoasis_time(date_start, dt, nbdt)
64	!
65	!
66	! This subroutine gets the start date of the simulation, the time step and the number
67	! of time steps we need to do until the end of the simulations.
68	!
69	!
70	!
71	REAL(r_std), INTENT(out) :: date_start !! The date at which the simulation starts
72	REAL(r_std), INTENT(out) :: dt !! Time step length in seconds
73	INTEGER(i_std), INTENT(out) :: nbdt !! Number of timesteps to be executed
74	!
75	! Local
76	!
77	CHARACTER(LEN=20) :: str_sdate(2), str_edate(2), tmpstr
78	INTEGER(i_std) :: s_year, s_month, s_day, e_year, e_month, e_day
79	INTEGER(i_std) :: seci, hours, minutes
80	REAL(r_std) :: s_sec, e_sec, dateend, diff_sec, date_end
81	INTEGER(i_std) :: i, ic
82	!
83	!Config Key = START_DATE
84	!Config Desc = Date at which the simulation starts
85	!Config Def = NONE
86	!Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
87	str_sdate = " "
88	CALL getin('START_DATE',str_sdate)
89	!
90	IF ( (INDEX(str_sdate(1),"-") .NE. INDEX(str_sdate(1),"-", .TRUE.)) .AND. &
91	& (INDEX(str_sdate(2),":") .NE. INDEX(str_sdate(2),":", .TRUE.)) ) THEN
92	DO i=1,2
93	tmpstr = str_sdate(1)
94	ic = INDEX(tmpstr,"-")
95	tmpstr(ic:ic) = " "
96	str_sdate(1) = tmpstr
97	tmpstr = str_sdate(2)
98	ic = INDEX(tmpstr,":")
99	tmpstr(ic:ic) = " "
100	str_sdate(2) = tmpstr
101	ENDDO
102	READ (str_sdate(1),*) s_year, s_month, s_day
103	READ (str_sdate(2),*) hours, minutes, seci
104	s_sec = hours3600. + minutes60. + seci
105	ELSE
106	CALL ipslerr(3, "orchoasis_time", "START_DATE incorrectly specified in run.def", str_sdate(1), str_sdate(2))
107	ENDIF
108	CALL ymds2ju (s_year, s_month, s_day, s_sec, date_start)
109	CALL orchoasis_printdate(date_start, "This is after reading the start date")
110	!
111	!Config Key = END_DATE
112	!Config Desc = Date at which the simulation ends
113	!Config Def = NONE
114	!Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
115	str_edate = " "
116	CALL getin('END_DATE',str_edate)
117	!
118	IF ( (INDEX(str_edate(1),"-") .NE. INDEX(str_edate(1),"-", .TRUE.)) .AND. &
119	& (INDEX(str_edate(2),":") .NE. INDEX(str_edate(2),":", .TRUE.)) ) THEN
120	DO i=1,2
121	tmpstr = str_edate(1)
122	ic = INDEX(tmpstr,"-")
123	tmpstr(ic:ic) = " "
124	str_edate(1) = tmpstr
125	tmpstr = str_edate(2)
126	ic = INDEX(tmpstr,":")
127	tmpstr(ic:ic) = " "
128	str_edate(2) = tmpstr
129	ENDDO
130	READ (str_edate(1),*) e_year, e_month, e_day
131	READ (str_edate(2),*) hours, minutes, seci
132	e_sec = hours3600. + minutes60. + seci
133	ELSE
134	CALL ipslerr(3, "orchoasis_time", "END_DATE incorrectly specified in run.def", str_edate(1), str_edate(2))
135	ENDIF
136	CALL ymds2ju (e_year, e_month, e_day, e_sec, date_end)
137	CALL orchoasis_printdate(date_start, "This is after reading the end date")
138	!
139	CALL time_diff (s_year,s_month,s_day,s_sec,e_year,e_month,e_day,e_sec,diff_sec)
140	!
141	!Config Key = DT_SECHIBA
142	!Config Desc = Time step length in seconds for Sechiba component
143	!Config Def = 1800
144	!Config Help =
145	!Config Units = [seconds]
146	dt = 1800
147	CALL getin('DT_SECHIBA', dt)
148	!
149	nbdt = NINT(diff_sec/dt)
150	!
151	END SUBROUTINE orchoasis_time
152	!
153	!=============================================================================================
154	!
155	SUBROUTINE orchoasis_printpoint(julian_day, lon_pt, lat_pt, nbind, lalo, var, message, ktest)
156	!
157	REAL(r_std), INTENT(in) :: julian_day
158	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
159	INTEGER(i_std), INTENT(in) :: nbind
160	REAL(r_std), INTENT(in) :: lalo(:,:)
161	REAL(r_std), INTENT(in) :: var(:)
162	CHARACTER(len=*), INTENT(in) :: message
163	INTEGER(i_std), OPTIONAL, INTENT(out) :: ktest
164	!
165	!
166	!
167	INTEGER(i_std) :: year, month, day, hours, minutes, seci
168	REAL(r_std) :: sec, mindist
169	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist, refdist
170	INTEGER(i_std) :: lon_ind, lat_ind, ind
171	INTEGER(i_std) :: i, imin(1)
172	REAL(r_std), PARAMETER :: mincos = 0.0001
173	REAL(r_std), PARAMETER :: pi = 3.141592653589793238
174	REAL(r_std), PARAMETER :: R_Earth = 6378000.
175	!
176	! Check if there is anything to be done
177	!
178	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
179	IF ( PRESENT(ktest) ) ktest = -1
180	RETURN
181	ENDIF
182	!
183	! Allocate memory
184	!
185	ALLOCATE(dist(nbind))
186	ALLOCATE(refdist(nbind))
187	!
188	! Convert time first
189	!
190	CALL ju2ymds (julian_day, year, month, day, sec)
191	hours = INT(sec/3600)
192	sec = sec - 3600 * hours
193	minutes = INT(sec / 60)
194	sec = sec - 60 * minutes
195	seci = INT(sec)
196	!
197	! Get the location to be analysed
198	!
199	DO i=1,nbind
200	dist(i) = acos( sin(lat_ptpi/180)sin(lalo(i,1)*pi/180) + &
201	& cos(lat_ptpi/180)cos(lalo(i,1)pi/180)cos((lalo(i,2)-lon_pt)pi/180) ) R_Earth
202	ENDDO
203	!
204	! Look for the next grid point closest to the one with the smalest distance.
205	!
206	imin = MINLOC(dist)
207	DO i=1,nbind
208	refdist(i) = acos( sin(lalo(imin(1),1)pi/180)sin(lalo(i,1)*pi/180) + &
209	& cos(lalo(imin(1),1)pi/180)cos(lalo(i,1)pi/180) cos((lalo(i,2)-lalo(imin(1),2))pi/180) ) R_Earth
210	ENDDO
211	refdist(imin(1)) = MAXVAL(refdist)
212	mindist = MINVAL(refdist)
213	!
214	! Are we closer than the closest points ?
215	!
216	IF ( PRESENT(ktest) ) ktest = -1
217	IF ( dist(imin(1)) <= mindist ) THEN
218	!
219	WRITE(numout,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F6.1,',', F6.1,'(i=',I6,') Value = ',F12.4,A38)") &
220	& hours, minutes, seci, lalo(imin(1),2), lalo(imin(1),1), imin(1), var(imin(1)), message
221	!
222	IF ( PRESENT(ktest) ) ktest = imin(1)
223	ENDIF
224	!
225	END SUBROUTINE orchoasis_printpoint
226	!
227	!-------------------------------------------------------------------------------------------------------
228	!-
229	!- orchoasis_wrtgrid
230	!- Subrouting to write the grid to OASIS. Only to be called by root-proc
231	!-------------------------------------------------------------------------------------------------------
232	!
233	SUBROUTINE orchoasis_wrtgrid(oasis_gridname, iim, jjm, nbseg, lon, lat, &
234	& corners, area, mask)
235	!
236	! ARGUMENTS
237	!
238	CHARACTER(LEN=4), INTENT(in) :: oasis_gridname
239	INTEGER(i_std), INTENT(in) :: iim, jjm, nbseg
240	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm)
241	REAL(r_std), INTENT(in) :: corners(iim,jjm,nbseg,2)
242	REAL(r_std), INTENT(in) :: area(iim,jjm)
243	REAL(r_std), INTENT(in) :: mask(iim,jjm)
244	!
245	! LOCAL
246	!
247	CHARACTER(LEN=4) :: oasis_gridname_tmp
248	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: corners_lon, corners_lat
249	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: maskinv
250	INTEGER(i_std) :: flag
251	INTEGER(i_std) :: i, j
252	!
253	IF ( .NOT. is_root_prc) THEN
254	WRITE(numout,*) "orchoasis_wrtgrid has to be called on the roo proc."
255	WRITE(numout,*) "We are here on a non root processor. is_root_prc = ", is_root_prc
256	WRITE(numout,*) "STOP from orchoasis_wrtgrid"
257	STOP
258	ENDIF
259	!
260	CALL oasis_start_grids_writing(flag)
261	!
262	CALL oasis_write_grid(oasis_gridname, iim, jjm, lon, lat)
263	!
264	ALLOCATE(corners_lon(iim, jjm, nbseg), corners_lat(iim, jjm, nbseg))
265	corners_lon(:,:,:) = corners(:,:,:,1)
266	corners_lat(:,:,:) = corners(:,:,:,2)
267	CALL oasis_write_corner(oasis_gridname, iim, jjm, nbseg, corners_lon, corners_lat)
268	!
269	ALLOCATE(maskinv(iim, jjm))
270	DO i=1,iim
271	DO j=1,jjm
272	IF (mask(i,j) == 0) THEN
273	maskinv(i,j) = 1
274	ELSE
275	maskinv(i,j) = 0
276	ENDIF
277	ENDDO
278	ENDDO
279	CALL oasis_write_mask(oasis_gridname, iim, jjm, maskinv)
280	!
281	CALL oasis_write_area(oasis_gridname, iim, jjm, area)
282	!
283	CALL oasis_terminate_grids_writing()
284	!
285	IF ( INDEX(oasis_gridname,"WRF") > 0 ) THEN
286	oasis_gridname_tmp = "twrf"
287	!
288	CALL oasis_start_grids_writing(flag)
289	!
290	CALL oasis_write_grid(oasis_gridname_tmp, iim, jjm, lon, lat)
291	!
292	CALL oasis_write_corner(oasis_gridname_tmp, iim, jjm, nbseg, corners_lon, corners_lat)
293	!
294	CALL oasis_write_mask(oasis_gridname_tmp, iim, jjm, maskinv)
295	!
296	CALL oasis_write_area(oasis_gridname_tmp, iim, jjm, area)
297	!
298	CALL oasis_terminate_grids_writing()
299	ENDIF
300	!
301	DEALLOCATE(corners_lon, corners_lat)
302	DEALLOCATE(maskinv)
303	!
304	END SUBROUTINE orchoasis_wrtgrid
305	!
306	!-------------------------------------------------------------------------------------------------------
307	!-
308	!- orchoasis_defvar
309	!-
310	!-------------------------------------------------------------------------------------------------------
311	SUBROUTINE orchoasis_defvar (mpi_rank, kjpindex)
312	!
313	! ARGUMENTS
314	!
315	INTEGER(i_std), INTENT(in) :: mpi_rank, kjpindex
316	!
317	! LOCAL
318	!
319	INTEGER(i_std) :: ierror, tot_error
320	INTEGER(i_std), DIMENSION(3) :: ig_paral
321	INTEGER(i_std), DIMENSION(2) :: var_nodims, var_shape
322	CHARACTER(LEN=8) :: varname
323	!
324	!
325	!Config Key = COUPLING_LANDONLY
326	!Config Desc = Specifies if the OASIS coupler will provide only data on land or on land and ocean.
327	!Config Def = TRUE
328	!Config Help = If COUPLING_LANDONLY is set to TRUE, then ORCHIDEE expects to receive only land points
329	! from OASIS. On the contrary data on both land and ocean will be received and the arrays
330	! re-indexed so as to keep only land points.
331	landonly = .TRUE.
332	CALL getin('COUPLING_LANDONLY',landonly)
333	!
334	tot_error = 0
335	!
336	IF ( landonly ) THEN
337	ig_paral(1) = 1
338	ig_paral(2) = nbp_mpi_para_begin(mpi_rank)-1
339	ig_paral(3) = kjpindex
340	ELSE
341	ig_paral(1) = 1
342	ig_paral(2) = ij_para_begin(mpi_rank) - 1
343	ig_paral(3) = ij_nb
344	ENDIF
345
346	WRITE(,) mpi_rank, "Start and length =", ig_paral(2:3)
347
348	CALL oasis_def_partition (il_part_id, ig_paral, ierror)
349
350	IF ( landonly ) THEN
351	var_nodims(1) = 1
352	var_nodims(2) = 1
353	var_shape(1) = 1
354	var_shape(1) = kjpindex
355	ELSE
356	var_nodims(1) = 1
357	var_nodims(2) = 1
358	var_shape(1) = 1
359	var_shape(1) = ij_nb
360	ENDIF
361	!
362	! ORCHIDEE's Input variables
363	! ===========================
364	!
365	!
366	! Define levels
367	!
368	varname="HEIGHTTQ"
369	CALL oasis_def_var(zlevtq_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
370	tot_error = tot_error+ierror
371	!
372	varname="HEIGHTUV"
373	CALL oasis_def_var(zlevuv_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
374	tot_error = tot_error+ierror
375	!
376	!
377	varname="TEMPLEV1"
378	CALL oasis_def_var(tair_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
379	tot_error = tot_error+ierror
380	varname="HUMILEV1"
381	CALL oasis_def_var(qair_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
382	tot_error = tot_error+ierror
383	!
384	! Define precipitation variables
385	!
386	varname="RAINFALL"
387	CALL oasis_def_var(rainf_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
388	tot_error = tot_error+ierror
389	varname="SNOWFALL"
390	CALL oasis_def_var(snowf_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
391	tot_error = tot_error+ierror
392	!
393	!
394	! Define precipitation variables
395	!
396	varname="SHORTNET"
397	CALL oasis_def_var(swnet_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
398	tot_error = tot_error+ierror
399	varname="LONWDOWN"
400	CALL oasis_def_var(lwdown_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
401	tot_error = tot_error+ierror
402	varname="SOLARANG"
403	CALL oasis_def_var(solarang_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
404	tot_error = tot_error+ierror
405	!
406	! Define dynamical variables
407	!
408	varname="EASTWIND"
409	CALL oasis_def_var(u_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
410	tot_error = tot_error+ierror
411	varname="NORTWIND"
412	CALL oasis_def_var(v_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
413	tot_error = tot_error+ierror
414	varname="SURFPRES"
415	CALL oasis_def_var(ps_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
416	tot_error = tot_error+ierror
417	varname="MODCDRAG"
418	CALL oasis_def_var(cdrag_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
419	tot_error = tot_error+ierror
420	!
421	! ORCHIDEE's Output variables
422	! ===========================
423	!
424	! Turbulent fluxes
425	!
426	varname="TOTEVAPS"
427	CALL oasis_def_var(vevapp_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
428	tot_error = tot_error+ierror
429	!
430	varname="FLUXSENS"
431	CALL oasis_def_var(fluxsens_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
432	tot_error = tot_error+ierror
433	!
434	varname="FLUXLATE"
435	CALL oasis_def_var(fluxlat_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
436	tot_error = tot_error+ierror
437	!
438	! Discharge to the oceans
439	!
440	varname="COASTFLO"
441	CALL oasis_def_var(coastal_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
442	tot_error = tot_error+ierror
443	!
444	varname="RIVERFLO"
445	CALL oasis_def_var(river_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
446	tot_error = tot_error+ierror
447	!
448	! Carbon fluxes
449	!
450	varname="FLUNECO2"
451	CALL oasis_def_var(netco2_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
452	tot_error = tot_error+ierror
453	!
454	varname="FLULUCO2"
455	CALL oasis_def_var(carblu_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
456	tot_error = tot_error+ierror
457	!
458	! Surface states
459	!
460	varname="TSURFRAD"
461	CALL oasis_def_var(tsolrad_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
462	tot_error = tot_error+ierror
463	!
464	varname="TSURFNEW"
465	CALL oasis_def_var(tsolnew_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
466	tot_error = tot_error+ierror
467	!
468	varname="QSURFNEW"
469	CALL oasis_def_var(qsurf_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
470	tot_error = tot_error+ierror
471	!
472	varname="ALBEDVIS"
473	CALL oasis_def_var(albvis_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
474	tot_error = tot_error+ierror
475	!
476	varname="ALBEDNIR"
477	CALL oasis_def_var(albnir_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
478	tot_error = tot_error+ierror
479	!
480	varname="EMISLONW"
481	CALL oasis_def_var(emis_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
482	tot_error = tot_error+ierror
483	!
484	varname="ROUGHMOM"
485	CALL oasis_def_var(z0m_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
486	tot_error = tot_error+ierror
487	!
488	varname="ROUGHHEA"
489	CALL oasis_def_var(z0h_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
490	tot_error = tot_error+ierror
491	!
492	IF ( tot_error .NE. 0 ) THEN
493	CALL ipslerr(3, "orchoasis_defvar", "Definition of one of the coupling variables failed.", &
494	& "No futher information can be given at this point.", "")
495	ENDIF
496	!
497	END SUBROUTINE orchoasis_defvar
498	!-------------------------------------------------------------------------------------------------------
499	!-
500	!- orchoasis_defvecvar
501	!-
502	!- This routine defines the parition of the vector of river discharge for the ocean model.
503	!-
504	!-------------------------------------------------------------------------------------------------------
505	SUBROUTINE orchoasis_defvecvar (numlargest, nstart, nlen)
506	!
507	! ARGUMENTS
508	!
509	INTEGER, INTENT(IN) :: numlargest, nstart, nlen
510	!
511	! LOCAL
512	!
513	INTEGER(i_std), DIMENSION(3) :: ig_paral
514	CHARACTER(LEN=8) :: varname
515	INTEGER(i_std), DIMENSION(2) :: var_nodims
516	INTEGER(i_std), DIMENSION(2) :: var_shape
517	INTEGER(i_std) :: ier, tot_error
518	!
519	! Define partition for data to be transferred to the ocean.
520	!
521	ig_paral(1) = 1
522	ig_paral(2) = nstart
523	ig_paral(3) = nlen
524	CALL oasis_def_partition (river_part_id, ig_paral, ier)
525	IF ( ier .NE. 0 ) THEN
526	CALL ipslerr(3, "orchoasis_defvar", "Definition of vector partition for rivers failed.", &
527	& "If OASIS has a way to decide the error status it should be put here.", "")
528	ENDIF
529	!
530	var_nodims(1) = 1
531	var_nodims(2) = 1
532	var_shape(1) = 1
533	var_shape(2) = numlargest
534	!
535	!
536	!
537	tot_error = 0
538	varname="LONVECTO"
539	CALL oasis_def_var(riverlon_id, varname, river_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ier)
540	tot_error = tot_error+ier
541	varname="LATVECTO"
542	CALL oasis_def_var(riverlat_id, varname, river_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ier)
543	tot_error = tot_error+ier
544	varname="RIVERVEC"
545	CALL oasis_def_var(riverflow_id, varname, river_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ier)
546	tot_error = tot_error+ier
547	!
548	!
549	IF ( tot_error .NE. 0 ) THEN
550	CALL ipslerr(3, "orchoasis_defvecvar", "End of the definition of coupling variables failed.", &
551	& "If OASIS has a way to decide the error status it should be put here.", "")
552	ENDIF
553	!
554	END SUBROUTINE orchoasis_defvecvar
555
556	!-------------------------------------------------------------------------------------------------------
557	!-
558	!- orchoasis_getvar
559	!-
560	!-------------------------------------------------------------------------------------------------------
561	SUBROUTINE orchoasis_getvar (itau, dt, kjpindex, landindex, zlev_tq, zlev_uv, temp_air, qair, precip_rain, &
562	& precip_snow, swnet, lwdown, sinang, u, v, pb, cdrag)
563	!
564	! ARGUMENTS
565	!
566	INTEGER(i_std), INTENT(in) :: itau, kjpindex
567	INTEGER(i_std), DIMENSION(kjpindex), INTENT(in) :: landindex
568	REAL(r_std), INTENT(in) :: dt
569	!
570	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: zlev_tq, zlev_uv, temp_air, qair, precip_rain
571	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: precip_snow, swnet, lwdown, sinang, u, v, pb
572	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: cdrag
573	!
574	! LOCAL
575	!
576	INTEGER(i_std) :: oasis_info
577	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: buffer
578	!
579	! Dimension the buffer for getting the data from OASIS according to the number of
580	! points we will get : with or without ocean points.
581	!
582	IF ( .NOT. ALLOCATED(buffer) ) THEN
583	IF ( landonly ) THEN
584	ALLOCATE(buffer(kjpindex))
585	ELSE
586	ALLOCATE(buffer(ij_nb))
587	ENDIF
588	ENDIF
589	!
590	! Get first the levels
591	!
592	CALL oasis_get(zlevtq_id, NINT(itau*dt), buffer, oasis_info)
593	IF ( landonly ) THEN
594	zlev_tq = buffer
595	ELSE
596	zlev_tq = buffer(landindex)
597	ENDIF
598	CALL oasis_get(zlevuv_id, NINT(itau*dt), buffer, oasis_info)
599	IF ( landonly ) THEN
600	zlev_uv = buffer
601	ELSE
602	zlev_uv = buffer(landindex)
603	ENDIF
604	!
605	! Get atmospheric state variables
606	!
607	CALL oasis_get(tair_id, NINT(itau*dt), buffer, oasis_info)
608	IF ( landonly ) THEN
609	temp_air = buffer
610	ELSE
611	temp_air = buffer(landindex)
612	ENDIF
613	CALL oasis_get(qair_id, NINT(itau*dt), buffer, oasis_info)
614	IF ( landonly ) THEN
615	qair = buffer
616	ELSE
617	qair = buffer(landindex)
618	ENDIF
619	!
620	! Get precipitation fluxes
621	!
622	CALL oasis_get(rainf_id, NINT(itau*dt), buffer, oasis_info)
623	IF ( landonly ) THEN
624	precip_rain = buffer
625	ELSE
626	precip_rain = buffer(landindex)
627	ENDIF
628	CALL oasis_get(snowf_id, NINT(itau*dt), buffer, oasis_info)
629	IF ( landonly ) THEN
630	precip_snow = buffer
631	ELSE
632	precip_snow = buffer(landindex)
633	ENDIF
634	!
635	! Get Radiation fluxes
636	!
637	CALL oasis_get(swnet_id, NINT(itau*dt), buffer, oasis_info)
638	IF ( landonly ) THEN
639	swnet = buffer
640	ELSE
641	swnet = buffer(landindex)
642	ENDIF
643	CALL oasis_get(lwdown_id, NINT(itau*dt), buffer, oasis_info)
644	IF ( landonly ) THEN
645	lwdown = buffer
646	ELSE
647	lwdown = buffer(landindex)
648	ENDIF
649	CALL oasis_get(solarang_id, NINT(itau*dt), buffer, oasis_info)
650	IF ( landonly ) THEN
651	sinang = buffer
652	ELSE
653	sinang = buffer(landindex)
654	ENDIF
655	!
656	! Get dynamical variables
657	!
658	CALL oasis_get(u_id, NINT(itau*dt), buffer, oasis_info)
659	IF ( landonly ) THEN
660	u = buffer
661	ELSE
662	u = buffer(landindex)
663	ENDIF
664	CALL oasis_get(v_id, NINT(itau*dt), buffer, oasis_info)
665	IF ( landonly ) THEN
666	v = buffer
667	ELSE
668	v = buffer(landindex)
669	ENDIF
670	CALL oasis_get(ps_id, NINT(itau*dt), buffer, oasis_info)
671	IF ( landonly ) THEN
672	pb = buffer
673	ELSE
674	pb = buffer(landindex)
675	ENDIF
676	CALL oasis_get(cdrag_id, NINT(itau*dt), buffer, oasis_info)
677	IF ( landonly ) THEN
678	cdrag = buffer
679	ELSE
680	cdrag = buffer(landindex)
681	ENDIF
682	!
683	END SUBROUTINE orchoasis_getvar
684	!
685	!-------------------------------------------------------------------------------------------------------
686	!-
687	!- orchoasis_putvar
688	!-
689	!-------------------------------------------------------------------------------------------------------
690	SUBROUTINE orchoasis_putvar(itau, dt, kjpindex, landindex, vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
691	& netco2, carblu, tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, z0h)
692	!
693	! ARGUMENTS
694	!
695	INTEGER(i_std), INTENT(in) :: itau, kjpindex
696	INTEGER(i_std), DIMENSION(kjpindex), INTENT(in) :: landindex
697	REAL(r_std), INTENT(in) :: dt
698	!
699	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: vevapp, fluxsens, fluxlat, coastalflow, riverflow
700	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: netco2, carblu, tsol_rad, temp_sol_new, qsurf, emis
701	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: z0m, z0h
702	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: albedo
703	!
704	! LOCAL
705	!
706	INTEGER(i_std) :: oasis_info
707	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: buffer
708	!
709	IF ( .NOT. ALLOCATED(buffer) ) THEN
710	IF ( landonly ) THEN
711	ALLOCATE(buffer(kjpindex))
712	ELSE
713	ALLOCATE(buffer(ij_nb))
714	ENDIF
715	ENDIF
716	!
717	! Set all points to undef. This is the value which will remain over ocean points.
718	!
719	buffer(:) = undef_sechiba
720	!
721	! Turbulent fluxes
722	!
723	IF ( landonly ) THEN
724	buffer(:) = vevapp(:)
725	ELSE
726	buffer(landindex(:)) = vevapp(:)
727	ENDIF
728	CALL oasis_put(vevapp_id, NINT(itau*dt), buffer, oasis_info)
729	IF ( landonly ) THEN
730	buffer(:) = fluxsens(:)
731	ELSE
732	buffer(landindex(:)) = fluxsens(:)
733	ENDIF
734	CALL oasis_put(fluxsens_id, NINT(itau*dt), buffer, oasis_info)
735	IF ( landonly ) THEN
736	buffer(:) = fluxlat(:)
737	ELSE
738	buffer(landindex(:)) = fluxlat(:)
739	ENDIF
740	CALL oasis_put(fluxlat_id, NINT(itau*dt), buffer, oasis_info)
741	!
742	! Water fluxes to the ocean
743	!
744	IF ( landonly ) THEN
745	buffer(:) = coastalflow(:)
746	ELSE
747	buffer(landindex(:)) = coastalflow(:)
748	ENDIF
749	CALL oasis_put(coastal_id, NINT(itau*dt), buffer, oasis_info)
750	IF ( landonly ) THEN
751	buffer(:) = riverflow(:)
752	ELSE
753	buffer(landindex(:)) = riverflow(:)
754	ENDIF
755	CALL oasis_put(river_id, NINT(itau*dt), buffer, oasis_info)
756	!
757	! Carbon
758	!
759	IF ( landonly ) THEN
760	buffer(:) = netco2(:)
761	ELSE
762	buffer(landindex(:)) = netco2(:)
763	ENDIF
764	CALL oasis_put(netco2_id, NINT(itau*dt), buffer, oasis_info)
765	IF ( landonly ) THEN
766	buffer(:) = carblu(:)
767	ELSE
768	buffer(landindex(:)) = carblu(:)
769	ENDIF
770	CALL oasis_put(carblu_id, NINT(itau*dt), buffer, oasis_info)
771	!
772	! Surface conditions
773	!
774	IF ( landonly ) THEN
775	buffer(:) = tsol_rad(:)
776	ELSE
777	buffer(landindex(:)) = tsol_rad(:)
778	ENDIF
779	CALL oasis_put(tsolrad_id, NINT(itau*dt), buffer, oasis_info)
780	IF ( landonly ) THEN
781	buffer(:) = temp_sol_new(:)
782	ELSE
783	buffer(landindex(:)) = temp_sol_new(:)
784	ENDIF
785	CALL oasis_put(tsolnew_id, NINT(itau*dt), buffer, oasis_info)
786	IF ( landonly ) THEN
787	buffer(:) = qsurf(:)
788	ELSE
789	buffer(landindex(:)) = qsurf(:)
790	ENDIF
791	CALL oasis_put(qsurf_id, NINT(itau*dt), buffer, oasis_info)
792	!
793	! Other surface conditions
794	!
795	IF ( landonly ) THEN
796	buffer(:) = albedo(:,ivis)
797	ELSE
798	buffer(landindex(:)) = albedo(:,ivis)
799	ENDIF
800	CALL oasis_put(albvis_id, NINT(itau*dt), buffer, oasis_info)
801	IF ( landonly ) THEN
802	buffer(:) = albedo(:,inir)
803	ELSE
804	buffer(landindex(:)) = albedo(:,inir)
805	ENDIF
806	CALL oasis_put(albnir_id, NINT(itau*dt), buffer, oasis_info)
807	IF ( landonly ) THEN
808	buffer(:) = emis(:)
809	ELSE
810	buffer(landindex(:)) = emis(:)
811	ENDIF
812	CALL oasis_put(emis_id, NINT(itau*dt), buffer, oasis_info)
813	IF ( landonly ) THEN
814	buffer(:) = z0m(:)
815	ELSE
816	buffer(landindex(:)) = z0m(:)
817	ENDIF
818	CALL oasis_put(z0m_id, NINT(itau*dt), buffer, oasis_info)
819	IF ( landonly ) THEN
820	buffer(:) = z0h(:)
821	ELSE
822	buffer(landindex(:)) = z0h(:)
823	ENDIF
824	CALL oasis_put(z0h_id, NINT(itau*dt), buffer, oasis_info)
825	!
826	END SUBROUTINE orchoasis_putvar
827	!-------------------------------------------------------------------------------------------------------
828	!-
829	!- orchoasis_putvecvar
830	!-
831	!-------------------------------------------------------------------------------------------------------
832	SUBROUTINE orchoasis_putvecvar(itau, dt, numl, lonrivers, latrivers, flow)
833	!
834	! ARGUMENTS
835	!
836	INTEGER(i_std), INTENT(in) :: itau, numl
837	REAL(r_std), INTENT(in) :: dt
838	!
839	REAL(r_std), DIMENSION(numl), INTENT(in) :: lonrivers, latrivers, flow
840	!
841	! LOCAL
842	!
843	INTEGER(i_std) :: oasis_info
844	!
845	CALL oasis_put(riverlon_id, NINT(itau*dt), lonrivers, oasis_info)
846	CALL oasis_put(riverlat_id, NINT(itau*dt), latrivers, oasis_info)
847	CALL oasis_put(riverflow_id, NINT(itau*dt), flow, oasis_info)
848	!
849	END SUBROUTINE orchoasis_putvecvar
850	#endif
851	!
852	!=============================================================================================
853	!
854	! A diagnostic routine just to verify that the dates between the various components is the
855	! same.
856	!
857	!=============================================================================================
858	!
859	SUBROUTINE orchoasis_printdate(julian_day, message)
860	!
861	REAL(r_std), INTENT(in) :: julian_day
862	CHARACTER(len=*), INTENT(in) :: message
863	!
864	!
865	!
866	INTEGER(i_std) :: year, month, day, hours, minutes, seci
867	REAL(r_std) :: sec
868	!
869	CALL ju2ymds (julian_day, year, month, day, sec)
870	hours = INT(sec/3600)
871	sec = sec - 3600 * hours
872	minutes = INT(sec / 60)
873	sec = sec - 60 * minutes
874	seci = INT(sec)
875	!
876	WRITE(*,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
877	& year, month, day, hours, minutes, seci, message
878	!
879	END SUBROUTINE orchoasis_printdate
880	!
881	!=============================================================================================
882	!
883	! A routine to initialize the routing scheme before SECHIBA does it so that we can get the
884	! distribution of rivers among the processors and initialize the exchanges with the ocean.
885	!
886	!=============================================================================================
887	!
888	SUBROUTINE orchoasis_routinginit (itausec, nbploc, nbneighb, nbdlt, kindex, &
889	& date0, dt, lalo, &
890	& neighbours, resolution, contfrac, &
891	& numrivers, river_nstart, river_nlen, &
892	& lonrivers, latrivers)
893	!
894	! ARGUMENTS
895	!
896	INTEGER(i_std), INTENT(in) :: itausec, nbploc, nbneighb, nbdlt
897	INTEGER(i_std), INTENT(in) :: kindex(nbploc)
898	REAL(r_std), INTENT(in) :: date0 !! The date at which itau = 0
899	REAL(r_std), INTENT(in) :: dt !! Time step
900	REAL(r_std), INTENT(in) :: lalo(nbploc,2)
901	INTEGER(i_std), INTENT(in) :: neighbours(nbploc,nbneighb) !! Vector of neighbours for each grid point
902	!! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
903	REAL(r_std), INTENT(in) :: resolution(nbploc,2) !! The size of each grid box in X and Y (m)
904	REAL(r_std), INTENT(in) :: contfrac(nbploc) !! Fraction of land in each grid box (unitless;0-1)
905	INTEGER(i_std), INTENT(in) :: numrivers
906	INTEGER(i_std), INTENT(out) :: river_nstart, river_nlen
907	REAL(r_std), INTENT(out) :: lonrivers(numrivers), latrivers(numrivers)
908	!
909	! LOCAL
910	!
911	LOGICAL :: river_routingreg
912	INTEGER(i_std) :: ier
913	INTEGER(i_std) :: restid, histid, hist2id
914	CHARACTER(LEN=80) :: restname_in = 'NONE'
915	CHARACTER(LEN=80) :: restname_out = 'NONE'
916	INTEGER(i_std),PARAMETER :: llm=1
917	REAL(r_std), DIMENSION(llm) :: lev
918	LOGICAL, PARAMETER :: overwrite_time=.TRUE. !! Always override the date from the restart files for SECHIBA and STOMATE.
919	!! The date is taken from the gcm or from the driver restart file.
920	REAL(r_std) :: in_julian
921	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stempdiag
922	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: returnflow
923	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: reinfiltration
924	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: irrigation
925	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: riverflow
926	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: coastalflow
927	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: flood_frac
928	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: flood_res
929	!
930	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
931	!
932	river_routingreg = .FALSE.
933	CALL getin('ROUTING_REGIONAL', river_routingreg)
934	!
935	! Prepare the restart file
936	!
937	CALL getin_p('SECHIBA_restart_in',restname_in)
938	restname_out="TMP_restart.nc"
939	!
940	in_julian = itau2date(itausec, date0, dt)
941	lev(:) = zero
942	IF (is_root_prc) THEN
943	CALL restini( restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
944	& restname_out, itausec, date0, dt, restid, overwrite_time)
945	ELSE
946	restid=0
947	ENDIF
948	histid = -1
949	hist2id = -1
950	!
951	ALLOCATE(stempdiag(nbploc,nbdlt), stat=ier)
952	ALLOCATE(returnflow(nbploc), stat=ier)
953	ALLOCATE(reinfiltration(nbploc), stat=ier)
954	ALLOCATE(irrigation(nbploc), stat=ier)
955	ALLOCATE(riverflow(nbploc), stat=ier)
956	ALLOCATE(coastalflow(nbploc), stat=ier)
957	ALLOCATE(flood_frac(nbploc), stat=ier)
958	ALLOCATE(flood_res(nbploc), stat=ier)
959	!
960	IF ( river_routingreg ) THEN
961
962	stempdiag(:,:) = tp_00
963	returnflow(:) = zero
964	reinfiltration(:) = zero
965	irrigation(:) = zero
966	riverflow(:) = zero
967	coastalflow(:) = zero
968	flood_frac(:) = zero
969	flood_res(:) = zero
970	CALL routing_reg_initialize(itausec, nbploc, kindex, &
971	& restid, histid, hist2id, lalo, &
972	& neighbours, resolution, contfrac, stempdiag, &
973	& returnflow, reinfiltration, irrigation, riverflow, &
974	& coastalflow, flood_frac, flood_res, .TRUE.)
975	!
976	CALL routing_reg_riverloc(river_nstart, river_nlen, lonrivers, latrivers)
977	!
978	ELSE
979	CALL ipslerr(3, "orchoasis_routinginit", "It is not foreseen to call this routing for the global routing", &
980	& "For the moment only the regional routing can be coupled directly to the ocean.", "")
981	ENDIF
982	!
983	DEALLOCATE(stempdiag)
984	DEALLOCATE(returnflow)
985	DEALLOCATE(reinfiltration)
986	DEALLOCATE(irrigation)
987	DEALLOCATE(riverflow)
988	DEALLOCATE(coastalflow)
989	DEALLOCATE(flood_frac)
990	DEALLOCATE(flood_res)
991	!
992	CALL restclo()
993	CALL routing_reg_clear()
994	!
995	END SUBROUTINE orchoasis_routinginit
996	!
997	END MODULE orchoasis_tools

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

Download in other formats: