Context Navigation

intersurf.f90 @ 5143

Last change on this file since 5143 was 4074, checked in by jan.polcher, 7 years ago
Convergence with Trunk version 4061 and in particular introduces the option FROZ_FRAC_CORR to reduce runoff over frozen soils.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 98.0 KB

Line
1	!! This subroutine is the interface between the main program
2	!! (LMDZ or dim2_driver) and SECHIBA.
3	!! - Input fields are gathered to keep just continental points
4	!! - call sechiba_main That's SECHIBA process.
5	!! - Output fields are scattered to complete global fields
6	!!
7	!! @call sechiba_main
8	!! @Version : $Revision$, $Date$
9	!!
10	!! @author Marie-Alice Foujols and Jan Polcher
11	!!
12	!< $HeadURL$
13	!< $Date$
14	!< $Author$
15	!< $Revision$
16	!! IPSL (2006)
17	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
18	!!
19	!f90doc MODULEintersurf
20	MODULE intersurf
21
22	USE IOIPSL
23	USE xios_orchidee
24	USE ioipsl_para
25	USE defprec
26	USE sechiba
27	USE constantes
28	USE constantes_soil
29	USE control
30	USE pft_parameters
31	USE mod_orchidee_para
32	USE solar
33	USE grid
34	USE Orch_Write_field_p
35	USE thermosoilc, ONLY : thermosoilc_levels
36	USE ioipslctrl, ONLY : ioipslctrl_history, ioipslctrl_restini, ok_histsync, max_hist_level, dw
37
38	IMPLICIT NONE
39
40	PRIVATE
41	PUBLIC :: Init_intersurf, intersurf_main, intersurf_main_2d, intersurf_main_gathered, intsurf_time
42	PUBLIC :: intersurf_initialize_2d, intersurf_initialize_gathered
43
44
45	! Interface called from LMDZ
46	INTERFACE intersurf_main
47	MODULE PROCEDURE intersurf_main_gathered
48	END INTERFACE
49
50	!
51	! Global variables
52	!
53	!
54	LOGICAL, SAVE :: lstep_init_intersurf=.TRUE.!! Initialisation has to be done one time
55	!$OMP THREADPRIVATE(lstep_init_intersurf)
56	INTEGER(i_std), SAVE :: printlev_loc !! Write level to this module
57	!$OMP THREADPRIVATE(printlev_loc)
58	INTEGER(i_std), SAVE :: hist_id, rest_id !! IDs for history and restart files
59	!$OMP THREADPRIVATE(hist_id, rest_id)
60	INTEGER(i_std), SAVE :: hist2_id !! ID for the second history files (Hi-frequency ?)
61	!$OMP THREADPRIVATE(hist2_id)
62	INTEGER(i_std), SAVE :: hist_id_stom, hist_id_stom_IPCC, rest_id_stom !! Dito for STOMATE
63	!$OMP THREADPRIVATE(hist_id_stom, hist_id_stom_IPCC, rest_id_stom)
64	!
65	INTEGER(i_std), SAVE :: itau_offset !! This offset is used to phase the
66	! !! calendar of the GCM or the driver.
67	!$OMP THREADPRIVATE(itau_offset)
68	REAL(r_std), SAVE :: date0_shifted
69	!$OMP THREADPRIVATE(date0_shifted)
70	!
71	!! first day of this year
72	REAL(r_std), SAVE :: julian0
73	!$OMP THREADPRIVATE(julian0)
74	!
75	LOGICAL, PARAMETER :: check_INPUTS = .false. !! (very) long print of INPUTs in intersurf
76
77	LOGICAL, SAVE :: ok_q2m_t2m=.TRUE. !! Flag ok if the variables are present in the call in gcm.
78	!$OMP THREADPRIVATE(ok_q2m_t2m)
79	LOGICAL, SAVE :: fatmco2 !! Flag to force the value of atmospheric CO2 for vegetation.
80	!$OMP THREADPRIVATE(fatmco2)
81	REAL(r_std), SAVE :: atmco2 !! atmospheric CO2
82	!$OMP THREADPRIVATE(atmco2)
83	INTEGER(i_std), SAVE :: nb_fields_in=-1 !! Number of fields to give to the GCM
84	!$OMP THREADPRIVATE(nb_fields_in)
85	INTEGER(i_std), SAVE :: nb_fields_out=-1 !! Number of fields to get from the GCM
86	!$OMP THREADPRIVATE(nb_fields_out)
87
88	PUBLIC lstep_init_intersurf
89
90	CONTAINS
91
92	!! =============================================================================================================================
93	!! SUBROUTINE: intersurf_initialize_2d
94	!!
95	!>\BRIEF Initialization and call to sechiba_initialize
96	!!
97	!! DESCRIPTION: Initialization of module variables, read options from parameter file, initialize output files and call to
98	!! sechiba_initialize.
99	!!
100	!! This subroutine is called from dim2_driver before the first call to intersurf_main_2d.
101	!!
102	!! \n
103	!_ ==============================================================================================================================
104
105	SUBROUTINE intersurf_initialize_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
106	lrestart_read, lrestart_write, lon, lat, zcontfrac, zresolution, date0, &
107	zlev, u, v, qair, temp_air, epot_air, ccanopy, &
108	cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
109	precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
110	vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
111	tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m)
112
113	IMPLICIT NONE
114
115	!! 0. Variable and parameter declaration
116	!! 0.1 Input variables
117	INTEGER(i_std),INTENT (in) :: kjit !! Time step number
118	INTEGER(i_std),INTENT (in) :: iim, jjm !! Dimension of input fields
119	INTEGER(i_std),INTENT (in) :: kjpindex !! Number of continental points
120	REAL(r_std),INTENT (in) :: xrdt !! Time step in seconds
121	LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read
122	LOGICAL, INTENT (in) :: lrestart_write !! Logical for _restart_ file to write'
123	REAL(r_std), INTENT (in) :: date0 !! Date at which kjit = 0
124	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points
125	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: u !! Lowest level wind speed
126	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: v !! Lowest level wind speed
127	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: zlev !! Height of first layer
128	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: qair !! Lowest level specific humidity
129	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: precip_rain !! Rain precipitation
130	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: precip_snow !! Snow precipitation
131	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: lwdown !! Down-welling long-wave flux
132	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: swnet !! Net surface short-wave flux
133	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: swdown !! Downwelling surface short-wave flux
134	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: temp_air !! Air temperature in Kelvin
135	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: epot_air !! Air potential energy
136	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: ccanopy !! CO2 concentration in the canopy
137	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution
138	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: peqAcoef !! One for T and another for q
139	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution
140	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: peqBcoef !! One for T and another for q
141	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: pb !! Surface pressure
142	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: lon, lat !! Geographical coordinates
143	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: zcontfrac !! Fraction of continent in the grid
144	REAL(r_std),DIMENSION (iim,jjm,2), INTENT(in) :: zresolution !! resolution in x and y dimensions
145
146	!! 0.2 Output variables
147	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: z0m !! Surface roughness
148	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: coastalflow !! Diffuse flow of water into the ocean (m^3/dt)
149	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: riverflow !! Largest rivers flowing into the ocean (m^3/dt)
150	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: tsol_rad !! Radiative surface temperature
151	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: vevapp !! Total of evaporation
152	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: temp_sol_new !! New soil temperature
153	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: qsurf !! Surface specific humidity
154	REAL(r_std),DIMENSION (iim,jjm,2), INTENT(out) :: albedo !! Albedo
155	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: fluxsens !! Sensible chaleur flux
156	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: fluxlat !! Latent chaleur flux
157	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: emis !! Emissivity
158
159	!! 0.3 Modified variables
160	REAL(r_std),DIMENSION (iim,jjm), INTENT(inout) :: cdrag !! Cdrag
161
162	!! 0.4 Local variables
163	REAL(r_std),DIMENSION (kjpindex) :: zu !! Work array to keep u
164	REAL(r_std),DIMENSION (kjpindex) :: zv !! Work array to keep v
165	REAL(r_std),DIMENSION (kjpindex) :: zzlev !! Work array to keep zlev
166	REAL(r_std),DIMENSION (kjpindex) :: zqair !! Work array to keep qair
167	REAL(r_std),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain
168	REAL(r_std),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow
169	REAL(r_std),DIMENSION (kjpindex) :: zlwdown !! Work array to keep lwdown
170	REAL(r_std),DIMENSION (kjpindex) :: zswnet !! Work array to keep swnet
171	REAL(r_std),DIMENSION (kjpindex) :: zswdown !! Work array to keep swdown
172	REAL(r_std),DIMENSION (kjpindex) :: ztemp_air !! Work array to keep temp_air
173	REAL(r_std),DIMENSION (kjpindex) :: zepot_air !! Work array to keep epot_air
174	REAL(r_std),DIMENSION (kjpindex) :: zccanopy !! Work array to keep ccanopy
175	REAL(r_std),DIMENSION (kjpindex) :: zpetAcoef !! Work array to keep petAcoef
176	REAL(r_std),DIMENSION (kjpindex) :: zpeqAcoef !! Work array to keep peqAcoef
177	REAL(r_std),DIMENSION (kjpindex) :: zpetBcoef !! Work array to keep petBcoef
178	REAL(r_std),DIMENSION (kjpindex) :: zpeqBcoef !! Work array to keep peqVcoef
179	REAL(r_std),DIMENSION (kjpindex) :: zcdrag !! Work array to keep cdrag
180	REAL(r_std),DIMENSION (kjpindex) :: zpb !! Work array to keep pb
181	REAL(r_std),DIMENSION (kjpindex) :: zz0m, zz0h !! Work array to keep zz0m, zz0h
182	REAL(r_std),DIMENSION (kjpindex) :: zcoastal !! Work array to keep coastalflow
183	REAL(r_std),DIMENSION (kjpindex) :: zriver !! Work array to keep riverflow
184	REAL(r_std),DIMENSION (kjpindex) :: znetco2 !! Work array to keep netco2flux
185	REAL(r_std),DIMENSION (kjpindex) :: zcarblu !! Work array to keep fco2_land_use
186	REAL(r_std),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad
187	REAL(r_std),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp
188	REAL(r_std),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new
189	REAL(r_std),DIMENSION (kjpindex) :: zqsurf !! Work array to keep qsurf
190	REAL(r_std),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo
191	REAL(r_std),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens
192	REAL(r_std),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat
193	REAL(r_std),DIMENSION (kjpindex) :: zemis !! Work array to keep emis
194	REAL(r_std),ALLOCATABLE, DIMENSION (:) :: soilth_lev !! Vertical soil axis for thermal scheme (m)
195	INTEGER(i_std) :: i, j, ik
196	INTEGER(i_std) :: ier
197	INTEGER(i_std) :: itau_sechiba
198	REAL(r_std) :: zlev_mean
199	INTEGER :: old_fileout !! old Logical Int for std IO output
200
201
202	CALL ipslnlf_p(new_number=numout,old_number=old_fileout)
203	!
204	! Configuration of SSL specific parameters
205	!
206	CALL control_initialize(xrdt)
207
208	CALL intsurf_time( kjit, date0)
209
210	! Initialize specific write level
211	printlev_loc=get_printlev('instersurf')
212
213	IF ( printlev_loc >=1 ) WRITE(numout,*) 'Initialisation of intersurf_main_2d'
214	OFF_LINE_MODE = .TRUE.
215
216	DO ik=1,kjpindex
217
218	j = ((kindex(ik)-1)/iim) + 1
219	i = (kindex(ik) - (j-1)*iim)
220
221	!- Create the internal coordinate table
222	!-
223	lalo(ik,1) = lat(i,j)
224	lalo(ik,2) = lon(i,j)
225	!
226	!- Store the fraction of the continents only once so that the user
227	!- does not change them afterwards.
228	!-
229	contfrac(ik) = zcontfrac(i,j)
230	ENDDO
231	CALL gather(contfrac,contfrac_g)
232	CALL gather(lalo,lalo_g)
233	CALL gather2D_mpi(lon,lon_g)
234	CALL gather2D_mpi(lat,lat_g)
235	CALL gather2D_mpi(zlev,zlev_g)
236
237	CALL ioipslctrl_restini(kjit, date0, xrdt, rest_id, rest_id_stom, itau_offset, date0_shifted)
238	itau_sechiba = kjit + itau_offset
239
240	!!- Initialize module for output with XIOS
241	!
242	! Get the vertical soil levels for the thermal scheme, to be used in xios_orchidee_init
243	ALLOCATE(soilth_lev(ngrnd), stat=ier)
244	IF (ier /= 0) THEN
245	CALL ipslerr_p(3,'intersurf_main_2d', 'Error in allocation of soilth_lev','','')
246	END IF
247	IF (hydrol_cwrr) THEN
248	soilth_lev(1:ngrnd) = znt(:)
249	ELSE
250	soilth_lev(1:ngrnd) = thermosoilc_levels()
251	END IF
252
253	CALL xios_orchidee_init( MPI_COMM_ORCH, &
254	date0, year, month, day, &
255	lon, lat, soilth_lev )
256
257	!- Initialize IOIPSL sechiba output files
258	CALL ioipslctrl_history(iim, jjm, lon, lat, kindex, kjpindex, itau_sechiba, date0_shifted, xrdt, hist_id, &
259	hist2_id, hist_id_stom, hist_id_stom_IPCC)
260
261	IF ( printlev_loc>=3 ) WRITE(numout,*) 'End of Initialisation of intersurf'
262	!
263	! Shift the time step to phase the two models
264	!
265	itau_sechiba = kjit + itau_offset
266
267	! Update the calendar in xios by sending the new time step
268	! Special case : the model is only in initialization phase and the current itau_sechiba is not a real time step.
269	! Therefor give itau_sechiba+1 to xios to have a correct time axis in output files.
270	CALL xios_orchidee_update_calendar(itau_sechiba+1)
271
272	CALL intsurf_time( itau_sechiba, date0_shifted )
273	!
274	! 1. gather input fields from kindex array
275	! Warning : I'm not sure this interface with one dimension array is the good one
276	!
277	DO ik=1, kjpindex
278
279	j = ((kindex(ik)-1)/iim) + 1
280	i = (kindex(ik) - (j-1)*iim)
281
282	zu(ik) = u(i,j)
283	zv(ik) = v(i,j)
284	zzlev(ik) = zlev(i,j)
285	zqair(ik) = qair(i,j)
286	zprecip_rain(ik) = precip_rain(i,j)*xrdt
287	zprecip_snow(ik) = precip_snow(i,j)*xrdt
288	zlwdown(ik) = lwdown(i,j)
289	zswnet(ik) = swnet(i,j)
290	zswdown(ik) = swdown(i,j)
291	ztemp_air(ik) = temp_air(i,j)
292	zepot_air(ik) = epot_air(i,j)
293	zccanopy(ik) = ccanopy(i,j)
294	zpetAcoef(ik) = petAcoef(i,j)
295	zpeqAcoef(ik) = peqAcoef(i,j)
296	zpetBcoef(ik) = petBcoef(i,j)
297	zpeqBcoef(ik) = peqBcoef(i,j)
298	zcdrag(ik) = cdrag(i,j)
299	zpb(ik) = pb(i,j)
300
301	ENDDO
302	!
303	IF (check_INPUTS) THEN
304	WRITE(numout,*) "Intersurf_main_2D :"
305	WRITE(numout,*) "Time step number = ",kjit
306	WRITE(numout,*) "Dimension of input fields = ",iim, jjm
307	WRITE(numout,*) "Number of continental points = ",kjpindex
308	WRITE(numout,*) "Time step in seconds = ",xrdt
309	WRITE(numout,*) "Logical for _restart_ file to read, write = ",lrestart_read,lrestart_write
310	WRITE(numout,*) "Date at which kjit = 0 = ",date0
311	WRITE(numout,*) "Index for continental points = ",kindex
312	WRITE(numout,*) "Lowest level wind speed North = ",zu
313	WRITE(numout,*) "Lowest level wind speed East = ",zv
314	WRITE(numout,*) "Height of first layer = ",zzlev
315	WRITE(numout,*) "Lowest level specific humidity = ",zqair
316	WRITE(numout,*) "Rain precipitation = ",zprecip_rain
317	WRITE(numout,*) "Snow precipitation = ",zprecip_snow
318	WRITE(numout,*) "Down-welling long-wave flux = ",zlwdown
319	WRITE(numout,*) "Net surface short-wave flux = ",zswnet
320	WRITE(numout,*) "Downwelling surface short-wave flux = ",zswdown
321	WRITE(numout,*) "Air temperature in Kelvin = ",ztemp_air
322	WRITE(numout,*) "Air potential energy = ",zepot_air
323	WRITE(numout,*) "CO2 concentration in the canopy = ",zccanopy
324	WRITE(numout,*) "Coeficients A from the PBL resolution = ",zpetAcoef
325	WRITE(numout,*) "One for T and another for q = ",zpeqAcoef
326	WRITE(numout,*) "Coeficients B from the PBL resolution = ",zpetBcoef
327	WRITE(numout,*) "One for T and another for q = ",zpeqBcoef
328	WRITE(numout,*) "Cdrag = ",zcdrag
329	WRITE(numout,*) "Surface pressure = ",zpb
330	WRITE(numout,*) "Geographical coordinates lon = ", (/ ( lon(ilandindex(ik), jlandindex(ik)), ik=1,kjpindex ) /)
331	WRITE(numout,*) "Geographical coordinates lat = ", (/ ( lat(ilandindex(ik), jlandindex(ik)), ik=1,kjpindex ) /)
332	WRITE(numout,*) "Fraction of continent in the grid = ",contfrac
333	ENDIF
334	!
335	! 2. save the grid
336	!
337	IF ( printlev_loc>=3 ) WRITE(numout,*) 'Save the grid'
338	CALL histwrite_p(hist_id, 'LandPoints', itau_sechiba+1, (/ ( REAL(ik), ik=1,kjpindex ) /), kjpindex, kindex)
339	CALL histwrite_p(hist_id, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
340	IF ( ok_stomate ) THEN
341	CALL histwrite_p(hist_id_stom, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
342	CALL histwrite_p(hist_id_stom_IPCC, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
343	ENDIF
344
345	CALL histwrite_p(hist_id, 'Contfrac', itau_sechiba+1, contfrac, kjpindex, kindex)
346	IF ( is_omp_root .AND. hist_id > 0 ) THEN
347	! Always syncronize output after initialization
348	CALL histsync(hist_id)
349	END IF
350
351	CALL histwrite_p(hist2_id, 'LandPoints', itau_sechiba+1, (/ ( REAL(ik), ik=1,kjpindex ) /), kjpindex, kindex)
352	CALL histwrite_p(hist2_id, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
353	CALL histwrite_p(hist2_id, 'Contfrac', itau_sechiba+1, contfrac, kjpindex, kindex)
354	IF ( is_omp_root .AND. hist2_id > 0 ) THEN
355	! Always syncronize output after initialization
356	CALL histsync(hist2_id)
357	ENDIF
358
359	!
360	! 3. call sechiba for continental points only
361	!
362	IF ( printlev_loc >= 3 ) WRITE(numout,*) 'Calling sechiba_initialize'
363
364	CALL sechiba_initialize( &
365	itau_sechiba, iim*jjm, kjpindex, kindex, &
366	lalo, contfrac, neighbours, resolution, zzlev, &
367	zu, zv, zqair, ztemp_air, ztemp_air, &
368	zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, &
369	zprecip_rain, zprecip_snow, zlwdown, zswnet, zswdown, &
370	zpb, rest_id, hist_id, hist2_id, &
371	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
372	zcoastal, zriver, ztsol_rad, zvevapp, zqsurf, &
373	zz0m, zz0h, zalbedo, zfluxsens, zfluxlat, &
374	zemis, znetco2, zcarblu, ztemp_sol_new, zcdrag)
375
376	IF ( printlev_loc >= 3 ) WRITE(numout,*) 'out of sechiba_initialize'
377	!
378	! 5. scatter output fields
379	!
380	z0m(:,:) = undef_sechiba
381	coastalflow(:,:) = undef_sechiba
382	riverflow(:,:) = undef_sechiba
383	tsol_rad(:,:) = undef_sechiba
384	vevapp(:,:) = undef_sechiba
385	temp_sol_new(:,:) = undef_sechiba
386	qsurf(:,:) = undef_sechiba
387	albedo(:,:,:) = undef_sechiba
388	fluxsens(:,:) = undef_sechiba
389	fluxlat(:,:) = undef_sechiba
390	emis(:,:) = undef_sechiba
391	cdrag(:,:) = undef_sechiba
392
393	DO ik=1, kjpindex
394	j = ((kindex(ik)-1)/iim) + 1
395	i = (kindex(ik) - (j-1)*iim)
396
397	z0m(i,j) = zz0m(ik)
398	coastalflow(i,j) = zcoastal(ik)
399	riverflow(i,j) = zriver(ik)
400	tsol_rad(i,j) = ztsol_rad(ik)
401	vevapp(i,j) = zvevapp(ik)
402	temp_sol_new(i,j) = ztemp_sol_new(ik)
403	qsurf(i,j) = zqsurf(ik)
404	albedo(i,j,1) = zalbedo(ik,1)
405	albedo(i,j,2) = zalbedo(ik,2)
406	fluxsens(i,j) = zfluxsens(ik)
407	fluxlat(i,j) = zfluxlat(ik)
408	emis(i,j) = zemis(ik)
409	cdrag(i,j) = zcdrag(ik)
410
411	ENDDO
412
413	!
414	! 6. Transform the water fluxes into Kg/m^2s and m^3/s
415	!
416	DO ik=1, kjpindex
417
418	j = ((kindex(ik)-1)/iim) + 1
419	i = (kindex(ik) - (j-1)*iim)
420
421	vevapp(i,j) = vevapp(i,j)/xrdt
422	coastalflow(i,j) = coastalflow(i,j)/xrdt
423	riverflow(i,j) = riverflow(i,j)/xrdt
424
425	ENDDO
426
427	IF (is_root_prc) CALL getin_dump
428
429	lstep_init_intersurf = .FALSE.
430	IF (printlev_loc >=3) WRITE (numout,*) ' intersurf_main done '
431	CALL ipslnlf_p(new_number=old_fileout)
432	END SUBROUTINE intersurf_initialize_2d
433
434
435	!! =============================================================================================================================
436	!! SUBROUTINE: intersurf_main_2d
437	!!
438	!>\BRIEF Main subroutine to call ORCHIDEE from dim2_driver using variables on a 2d grid.
439	!!
440	!! DESCRIPTION: This subroutine is the main interface for ORCHIDEE when it is called from the offline driver dim2_driver.
441	!! The variables are all on the 2D grid including ocean points. intersurf_initialize_2d should be called before
442	!! this subroutine is called. This subroutine is called at each time step.
443	!!
444	!! \n
445	!_ ==============================================================================================================================
446
447	SUBROUTINE intersurf_main_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
448	lrestart_read, lrestart_write, lon, lat, zcontfrac, zresolution, date0, &
449	zlev, u, v, qair, temp_air, epot_air, ccanopy, &
450	cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
451	precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
452	vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
453	tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, &
454	coszang)
455
456	IMPLICIT NONE
457
458	!! 0. Variable and parameter declaration
459	!! 0.1 Input variables
460	INTEGER(i_std),INTENT (in) :: kjit !! Time step number
461	INTEGER(i_std),INTENT (in) :: iim, jjm !! Dimension of input fields
462	INTEGER(i_std),INTENT (in) :: kjpindex !! Number of continental points
463	REAL(r_std),INTENT (in) :: xrdt !! Time step in seconds
464	LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read
465	LOGICAL, INTENT (in) :: lrestart_write !! Logical for _restart_ file to write'
466	REAL(r_std), INTENT (in) :: date0 !! Date at which kjit = 0
467	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points
468	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: u !! Lowest level wind speed
469	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: v !! Lowest level wind speed
470	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: zlev !! Height of first layer
471	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: qair !! Lowest level specific humidity
472	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: precip_rain !! Rain precipitation
473	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: precip_snow !! Snow precipitation
474	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: lwdown !! Down-welling long-wave flux
475	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: swnet !! Net surface short-wave flux
476	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: swdown !! Downwelling surface short-wave flux
477	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: coszang !! Cosine of the solar zenith angle (unitless)
478	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: temp_air !! Air temperature in Kelvin
479	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: epot_air !! Air potential energy
480	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: ccanopy !! CO2 concentration in the canopy
481	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution
482	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: peqAcoef !! One for T and another for q
483	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution
484	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: peqBcoef !! One for T and another for q
485	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: pb !! Surface pressure
486	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: lon, lat !! Geographical coordinates
487	REAL(r_std),DIMENSION (iim,jjm), INTENT(in) :: zcontfrac !! Fraction of continent in the grid
488	REAL(r_std),DIMENSION (iim,jjm,2), INTENT(in) :: zresolution !! resolution in x and y dimensions
489
490	!! 0.2 Output variables
491	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: z0m !! Surface roughness for momemtum
492	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: coastalflow !! Diffuse flow of water into the ocean (m^3/dt)
493	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: riverflow !! Largest rivers flowing into the ocean (m^3/dt)
494	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: tsol_rad !! Radiative surface temperature
495	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: vevapp !! Total of evaporation
496	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: temp_sol_new !! New soil temperature
497	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: qsurf !! Surface specific humidity
498	REAL(r_std),DIMENSION (iim,jjm,2), INTENT(out) :: albedo !! Albedo
499	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: fluxsens !! Sensible chaleur flux
500	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: fluxlat !! Latent chaleur flux
501	REAL(r_std),DIMENSION (iim,jjm), INTENT(out) :: emis !! Emissivity
502
503	!! 0.3 Modified variables
504	REAL(r_std),DIMENSION (iim,jjm), INTENT(inout) :: cdrag !! Cdrag
505
506	!! 0.4 Local variables
507	REAL(r_std),DIMENSION (kjpindex) :: zu !! Work array to keep u
508	REAL(r_std),DIMENSION (kjpindex) :: zv !! Work array to keep v
509	REAL(r_std),DIMENSION (kjpindex) :: zzlev !! Work array to keep zlev
510	REAL(r_std),DIMENSION (kjpindex) :: zqair !! Work array to keep qair
511	REAL(r_std),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain
512	REAL(r_std),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow
513	REAL(r_std),DIMENSION (kjpindex) :: zlwdown !! Work array to keep lwdown
514	REAL(r_std),DIMENSION (kjpindex) :: zswnet !! Work array to keep swnet
515	REAL(r_std),DIMENSION (kjpindex) :: zswdown !! Work array to keep swdown
516	REAL(r_std),DIMENSION (kjpindex) :: zcoszang !! Work array to keep coszang
517	REAL(r_std),DIMENSION (kjpindex) :: ztemp_air !! Work array to keep temp_air
518	REAL(r_std),DIMENSION (kjpindex) :: zepot_air !! Work array to keep epot_air
519	REAL(r_std),DIMENSION (kjpindex) :: zccanopy !! Work array to keep ccanopy
520	REAL(r_std),DIMENSION (kjpindex) :: zpetAcoef !! Work array to keep petAcoef
521	REAL(r_std),DIMENSION (kjpindex) :: zpeqAcoef !! Work array to keep peqAcoef
522	REAL(r_std),DIMENSION (kjpindex) :: zpetBcoef !! Work array to keep petBcoef
523	REAL(r_std),DIMENSION (kjpindex) :: zpeqBcoef !! Work array to keep peqVcoef
524	REAL(r_std),DIMENSION (kjpindex) :: zcdrag !! Work array to keep cdrag
525	REAL(r_std),DIMENSION (kjpindex) :: zpb !! Work array to keep pb
526	REAL(r_std),DIMENSION (kjpindex) :: zz0m, zz0h !! Work array to keep zz0m, zz0h
527	REAL(r_std),DIMENSION (kjpindex) :: zcoastal !! Work array to keep coastalflow
528	REAL(r_std),DIMENSION (kjpindex) :: zriver !! Work array to keep riverflow
529	REAL(r_std),DIMENSION (kjpindex) :: znetco2 !! Work array to keep netco2flux
530	REAL(r_std),DIMENSION (kjpindex) :: zcarblu !! Work array to keep fco2_land_use
531	REAL(r_std),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad
532	REAL(r_std),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp
533	REAL(r_std),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new
534	REAL(r_std),DIMENSION (kjpindex) :: zqsurf !! Work array to keep qsurf
535	REAL(r_std),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo
536	REAL(r_std),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens
537	REAL(r_std),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat
538	REAL(r_std),DIMENSION (kjpindex) :: zemis !! Work array to keep emis
539	REAL(r_std),ALLOCATABLE, DIMENSION (:) :: soilth_lev !! Vertical soil axis for thermal scheme (m)
540	INTEGER(i_std) :: i, j, ik
541	INTEGER(i_std) :: ier
542	INTEGER(i_std) :: itau_sechiba
543	REAL(r_std) :: zlev_mean
544	INTEGER :: old_fileout !! old Logical Int for std IO output
545
546	CALL ipslnlf_p(new_number=numout,old_number=old_fileout)
547
548	!
549	! Shift the time step to phase the two models
550	!
551	itau_sechiba = kjit + itau_offset
552	!
553
554	! Update the calendar in xios by sending the new time step
555	CALL xios_orchidee_update_calendar(itau_sechiba)
556
557	CALL intsurf_time( itau_sechiba, date0_shifted )
558	!
559	! 1. gather input fields from kindex array
560	! Warning : I'm not sure this interface with one dimension array is the good one
561	!
562	DO ik=1, kjpindex
563
564	j = ((kindex(ik)-1)/iim) + 1
565	i = (kindex(ik) - (j-1)*iim)
566
567	zu(ik) = u(i,j)
568	zv(ik) = v(i,j)
569	zzlev(ik) = zlev(i,j)
570	zqair(ik) = qair(i,j)
571	zprecip_rain(ik) = precip_rain(i,j)*xrdt
572	zprecip_snow(ik) = precip_snow(i,j)*xrdt
573	zlwdown(ik) = lwdown(i,j)
574	zswnet(ik) = swnet(i,j)
575	zswdown(ik) = swdown(i,j)
576	zcoszang(ik) = coszang(i,j)
577	ztemp_air(ik) = temp_air(i,j)
578	zepot_air(ik) = epot_air(i,j)
579	zccanopy(ik) = ccanopy(i,j)
580	zpetAcoef(ik) = petAcoef(i,j)
581	zpeqAcoef(ik) = peqAcoef(i,j)
582	zpetBcoef(ik) = petBcoef(i,j)
583	zpeqBcoef(ik) = peqBcoef(i,j)
584	zcdrag(ik) = cdrag(i,j)
585	zpb(ik) = pb(i,j)
586
587	ENDDO
588	!
589	IF (check_INPUTS) THEN
590	WRITE(numout,*) "Intersurf_main_2D :"
591	WRITE(numout,*) "Time step number = ",kjit
592	WRITE(numout,*) "Dimension of input fields = ",iim, jjm
593	WRITE(numout,*) "Number of continental points = ",kjpindex
594	WRITE(numout,*) "Time step in seconds = ",xrdt
595	WRITE(numout,*) "Logical for _restart_ file to read, write = ",lrestart_read,lrestart_write
596	WRITE(numout,*) "Date at which kjit = 0 = ",date0
597	WRITE(numout,*) "Index for continental points = ",kindex
598	WRITE(numout,*) "Lowest level wind speed North = ",zu
599	WRITE(numout,*) "Lowest level wind speed East = ",zv
600	WRITE(numout,*) "Height of first layer = ",zzlev
601	WRITE(numout,*) "Lowest level specific humidity = ",zqair
602	WRITE(numout,*) "Rain precipitation = ",zprecip_rain
603	WRITE(numout,*) "Snow precipitation = ",zprecip_snow
604	WRITE(numout,*) "Down-welling long-wave flux = ",zlwdown
605	WRITE(numout,*) "Net surface short-wave flux = ",zswnet
606	WRITE(numout,*) "Downwelling surface short-wave flux = ",zswdown
607	WRITE(numout,*) "Air temperature in Kelvin = ",ztemp_air
608	WRITE(numout,*) "Air potential energy = ",zepot_air
609	WRITE(numout,*) "CO2 concentration in the canopy = ",zccanopy
610	WRITE(numout,*) "Coeficients A from the PBL resolution = ",zpetAcoef
611	WRITE(numout,*) "One for T and another for q = ",zpeqAcoef
612	WRITE(numout,*) "Coeficients B from the PBL resolution = ",zpetBcoef
613	WRITE(numout,*) "One for T and another for q = ",zpeqBcoef
614	WRITE(numout,*) "Cdrag = ",zcdrag
615	WRITE(numout,*) "Surface pressure = ",zpb
616	WRITE(numout,*) "Geographical coordinates lon = ", (/ ( lon(ilandindex(ik), jlandindex(ik)), ik=1,kjpindex ) /)
617	WRITE(numout,*) "Geographical coordinates lat = ", (/ ( lat(ilandindex(ik), jlandindex(ik)), ik=1,kjpindex ) /)
618	WRITE(numout,*) "Fraction of continent in the grid = ",contfrac
619	ENDIF
620
621	!
622	! 3. call sechiba for continental points only
623	!
624	IF ( printlev_loc >= 3 ) WRITE(numout,*) 'Calling sechiba_main'
625
626	CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, &
627	lrestart_read, lrestart_write, &
628	lalo, contfrac, neighbours, resolution, &
629	zzlev, zu, zv, zqair, zqair, ztemp_air, ztemp_air, zepot_air, zccanopy, &
630	zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, &
631	zprecip_rain ,zprecip_snow, zlwdown, zswnet, zswdown, zcoszang, zpb, &
632	zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
633	ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0m, zz0h,&
634	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC )
635
636	IF ( printlev_loc >= 3 ) WRITE(numout,*) 'out of sechiba_main'
637
638	!
639	! 5. scatter output fields
640	!
641	z0m(:,:) = undef_sechiba
642	coastalflow(:,:) = undef_sechiba
643	riverflow(:,:) = undef_sechiba
644	tsol_rad(:,:) = undef_sechiba
645	vevapp(:,:) = undef_sechiba
646	temp_sol_new(:,:) = undef_sechiba
647	qsurf(:,:) = undef_sechiba
648	albedo(:,:,:) = undef_sechiba
649	fluxsens(:,:) = undef_sechiba
650	fluxlat(:,:) = undef_sechiba
651	emis(:,:) = undef_sechiba
652	cdrag(:,:) = undef_sechiba
653	!
654	DO ik=1, kjpindex
655
656
657	j = ((kindex(ik)-1)/iim) + 1
658	i = (kindex(ik) - (j-1)*iim)
659
660	z0m(i,j) = zz0m(ik)
661	coastalflow(i,j) = zcoastal(ik)
662	riverflow(i,j) = zriver(ik)
663	tsol_rad(i,j) = ztsol_rad(ik)
664	vevapp(i,j) = zvevapp(ik)
665	temp_sol_new(i,j) = ztemp_sol_new(ik)
666	qsurf(i,j) = zqsurf(ik)
667	albedo(i,j,1) = zalbedo(ik,1)
668	albedo(i,j,2) = zalbedo(ik,2)
669	fluxsens(i,j) = zfluxsens(ik)
670	fluxlat(i,j) = zfluxlat(ik)
671	emis(i,j) = zemis(ik)
672	cdrag(i,j) = zcdrag(ik)
673
674	ENDDO
675
676	CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /))
677	CALL xios_orchidee_send_field("areas", area)
678	CALL xios_orchidee_send_field("contfrac",contfrac)
679	CALL xios_orchidee_send_field("temp_air",ztemp_air)
680	CALL xios_orchidee_send_field("qair",zqair)
681	CALL xios_orchidee_send_field("swnet",zswnet)
682	CALL xios_orchidee_send_field("swdown",zswdown)
683	CALL xios_orchidee_send_field("pb",zpb)
684
685	IF ( .NOT. almaoutput ) THEN
686	!
687	! scattered during the writing
688	!
689	CALL histwrite_p (hist_id, 'evap', itau_sechiba, zvevapp, kjpindex, kindex)
690	CALL histwrite_p (hist_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
691	CALL histwrite_p (hist_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
692	CALL histwrite_p (hist_id, 'temp_sol', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
693	CALL histwrite_p (hist_id, 'tsol_max', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
694	CALL histwrite_p (hist_id, 'tsol_min', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
695	CALL histwrite_p (hist_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
696	CALL histwrite_p (hist_id, 'fluxlat', itau_sechiba, zfluxlat, kjpindex, kindex)
697	CALL histwrite_p (hist_id, 'swnet', itau_sechiba, zswnet, kjpindex, kindex)
698	CALL histwrite_p (hist_id, 'swdown', itau_sechiba, zswdown, kjpindex, kindex)
699	CALL histwrite_p (hist_id, 'alb_vis', itau_sechiba, zalbedo(:,1), kjpindex, kindex)
700	CALL histwrite_p (hist_id, 'alb_nir', itau_sechiba, zalbedo(:,2), kjpindex, kindex)
701	CALL histwrite_p (hist_id, 'tair', itau_sechiba, ztemp_air, kjpindex, kindex)
702	CALL histwrite_p (hist_id, 'qair', itau_sechiba, zqair, kjpindex, kindex)
703	CALL histwrite_p (hist_id, 'q2m', itau_sechiba, zqair, kjpindex, kindex)
704	CALL histwrite_p (hist_id, 't2m', itau_sechiba, ztemp_air, kjpindex, kindex)
705	CALL histwrite_p (hist2_id, 'evap', itau_sechiba, zvevapp, kjpindex, kindex)
706	CALL histwrite_p (hist2_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
707	CALL histwrite_p (hist2_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
708	CALL histwrite_p (hist2_id, 'temp_sol', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
709	CALL histwrite_p (hist2_id, 'tsol_max', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
710	CALL histwrite_p (hist2_id, 'tsol_min', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
711	CALL histwrite_p (hist2_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
712	CALL histwrite_p (hist2_id, 'fluxlat', itau_sechiba, zfluxlat, kjpindex, kindex)
713	CALL histwrite_p (hist2_id, 'swnet', itau_sechiba, zswnet, kjpindex, kindex)
714	CALL histwrite_p (hist2_id, 'swdown', itau_sechiba, zswdown, kjpindex, kindex)
715	CALL histwrite_p (hist2_id, 'alb_vis', itau_sechiba, zalbedo(:,1), kjpindex, kindex)
716	CALL histwrite_p (hist2_id, 'alb_nir', itau_sechiba, zalbedo(:,2), kjpindex, kindex)
717	CALL histwrite_p (hist2_id, 'tair', itau_sechiba, ztemp_air, kjpindex, kindex)
718	CALL histwrite_p (hist2_id, 'qair', itau_sechiba, zqair, kjpindex, kindex)
719	CALL histwrite_p (hist2_id, 'q2m', itau_sechiba, zqair, kjpindex, kindex)
720	CALL histwrite_p (hist2_id, 't2m', itau_sechiba, ztemp_air, kjpindex, kindex)
721	ELSE
722	CALL histwrite_p (hist_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
723	CALL histwrite_p (hist_id, 'SWnet', itau_sechiba, zswnet, kjpindex, kindex)
724	CALL histwrite_p (hist_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
725	CALL histwrite_p (hist_id, 'Qle', itau_sechiba, zfluxlat, kjpindex, kindex)
726	CALL histwrite_p (hist_id, 'AvgSurfT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
727	CALL histwrite_p (hist_id, 'RadT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
728	CALL histwrite_p (hist_id, 'Tair', itau_sechiba, ztemp_air, kjpindex, kindex)
729	CALL histwrite_p (hist_id, 'Qair', itau_sechiba, zqair, kjpindex, kindex)
730	CALL histwrite_p (hist2_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
731	CALL histwrite_p (hist2_id, 'SWnet', itau_sechiba, zswnet, kjpindex, kindex)
732	CALL histwrite_p (hist2_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
733	CALL histwrite_p (hist2_id, 'Qle', itau_sechiba, zfluxlat, kjpindex, kindex)
734	CALL histwrite_p (hist2_id, 'AvgSurfT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
735	CALL histwrite_p (hist2_id, 'RadT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
736	ENDIF
737	!
738	IF ( is_omp_root ) THEN
739	IF ( (dw .EQ. xrdt) .AND. hist_id > 0 ) THEN
740	! Syncronize output but only if flag ok_histsync is set to true
741	IF (ok_histsync) CALL histsync(hist_id)
742	ENDIF
743	END IF
744
745	!
746	! 6. Transform the water fluxes into Kg/m^2s and m^3/s
747	!
748	DO ik=1, kjpindex
749
750	j = ((kindex(ik)-1)/iim) + 1
751	i = (kindex(ik) - (j-1)*iim)
752
753	vevapp(i,j) = vevapp(i,j)/xrdt
754	coastalflow(i,j) = coastalflow(i,j)/xrdt
755	riverflow(i,j) = riverflow(i,j)/xrdt
756
757	ENDDO
758
759	IF (printlev_loc >=3) WRITE (numout,*) ' intersurf_main done '
760
761	CALL ipslnlf_p(new_number=old_fileout)
762
763	END SUBROUTINE intersurf_main_2d
764
765
766	!! =============================================================================================================================
767	!! SUBROUTINE: init_intersurf
768	!!
769	!>\BRIEF Initialize grid information
770	!!
771	!! DESCRIPTION: This subroutine is called from LMDZ before first call to intersurf_main_gathered or
772	!! intersurf_initialize_gathered
773	!!
774	!! \n
775	!_ ==============================================================================================================================
776
777	SUBROUTINE init_intersurf(nbp_l_lon,nbp_l_lat,kjpindex,kindex,orch_offset,orch_omp_size,orch_omp_rank,COMM)
778
779	USE mod_orchidee_para
780	USE timer
781	IMPLICIT NONE
782
783	INTEGER,INTENT(IN) :: nbp_l_lon
784	INTEGER,INTENT(IN) :: nbp_l_lat
785	INTEGER,INTENT(IN) :: kjpindex
786	INTEGER,INTENT(IN) :: kindex(:)
787	INTEGER,INTENT(IN) :: orch_offset
788	INTEGER,INTENT(IN) :: COMM
789	INTEGER,INTENT(IN) :: orch_omp_size
790	INTEGER,INTENT(IN) :: orch_omp_rank
791
792	INTEGER,DIMENSION(kjpindex) :: kindex_offset
793
794	IF (orch_omp_rank==0) THEN
795	CALL Init_timer
796	CALL start_timer(timer_mpi)
797	CALL grid_set_glo(nbp_l_lon,nbp_l_lat)
798	ENDIF
799	CALL barrier2_omp()
800	CALL init_orchidee_data_para(kjpindex,kindex,orch_offset,orch_omp_size,orch_omp_rank,COMM)
801	CALL Set_stdout_file('out_orchidee')
802
803	IF (is_omp_root) CALL grid_allocate_glo(4)
804	CALL barrier2_omp()
805	CALL init_ioipsl_para
806
807	kindex_offset(:)=kindex(:)+offset
808	CALL gather(kindex_offset,index_g)
809	CALL bcast(index_g)
810
811	WRITE(numout,*) "kjpindex = ",kjpindex
812	WRITE(numout,*) "offset for OMP = ",offset_omp
813	WRITE(numout,*) "Index num local for continental points = ",kindex
814	WRITE(numout,*) "Index num global for continental points = ",kindex_offset
815	IF (is_omp_root) THEN
816	WRITE(numout,*) "ROOT OMP, Index global MPI : ",kindex_mpi(:)
817	ENDIF
818	IF (is_root_prc) THEN
819	WRITE(numout,*) "ROOT global, Index global : ",index_g(:)
820	ENDIF
821
822
823	! Allocation of grid variables
824	CALL grid_init ( kjpindex, 4, "RegLonLat", "2DGrid" )
825
826	END SUBROUTINE init_intersurf
827
828	!! =============================================================================================================================
829	!! SUBROUTINE: intersurf_initialize_gathered
830	!!
831	!>\BRIEF Initialization and call to sechiba_initialize
832	!!
833	!! DESCRIPTION: Initialization of module variables, read options from parameter file, initialize output files and call to
834	!! sechiba_initialize.
835	!!
836	!! This subroutine can be called directly from GCM(LMDZ). If it is not called before the first call to
837	!! intersurf_main_gathered, then it will be done from there. This possibility is done to keep backward
838	!! compatibility with LMDZ.
839	!!
840	!! \n
841	!_ ==============================================================================================================================
842
843	SUBROUTINE intersurf_initialize_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
844	lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
845	zlev, u, v, qair, temp_air, epot_air, &
846	cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
847	precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
848	vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
849	tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, lon_scat_g, lat_scat_g, q2m, t2m, z0h, &
850	field_out_names, fields_out, field_in_names, fields_in)
851
852	USE mod_orchidee_para
853	IMPLICIT NONE
854
855	!! 0. Variable and parameter declaration
856	!! 0.1 Input
857	INTEGER(i_std),INTENT (in) :: kjit !! Time step number
858	INTEGER(i_std),INTENT (in) :: iim_glo !! Dimension of global fields
859	INTEGER(i_std),INTENT (in) :: jjm_glo !! Dimension of global fields
860	INTEGER(i_std),INTENT (in) :: kjpindex !! Number of continental points
861	REAL(r_std),INTENT (in) :: xrdt !! Time step in seconds
862	LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read
863	LOGICAL, INTENT (in) :: lrestart_write !! Logical for _restart_ file to write'
864	REAL(r_std), INTENT (in) :: date0 !! Date at which kjit = 0
865	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points
866	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: u !! Lowest level wind speed
867	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: v !! Lowest level wind speed
868	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: zlev !! Height of first layer
869	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: qair !! Lowest level specific humidity
870	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation
871	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation
872	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: lwdown !! Down-welling long-wave flux
873	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: swnet !! Net surface short-wave flux
874	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: swdown !! Downwelling surface short-wave flux
875	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature in Kelvin
876	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: epot_air !! Air potential energy
877	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution
878	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: peqAcoef !! One for T and another for q
879	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution
880	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: peqBcoef !! One for T and another for q
881	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
882	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: latlon !! Geographical coordinates
883	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: zcontfrac !! Fraction of continent
884	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: zneighbours !! neighbours
885	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: zresolution !! size of the grid box
886	REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN) :: lon_scat_g !! The scattered values for longitude
887	REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN) :: lat_scat_g !! The scattered values for latitudes
888
889	!! 0.2 Output variables
890	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: z0m !! Surface roughness (momentum)
891	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: coastalflow !! Diffuse flow of water into the ocean (m^3/dt)
892	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: riverflow !! Largest rivers flowing into the ocean (m^3/dt)
893	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: tsol_rad !! Radiative surface temperature
894	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: vevapp !! Total of evaporation
895	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: temp_sol_new !! New soil temperature
896	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: qsurf !! Surface specific humidity
897	REAL(r_std),DIMENSION (kjpindex,2), INTENT(out) :: albedo !! Albedo
898	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: fluxsens !! Sensible chaleur flux
899	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: fluxlat !! Latent chaleur flux
900	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: emis !! Emissivity
901
902	!! 0.3 Modified variables
903	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: cdrag !! Cdrag
904
905	!! 0.4 Optional input and output variables
906	REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL :: q2m !! Surface specific humidity
907	REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL :: t2m !! Surface air temperature
908	REAL(r_std),DIMENSION (kjpindex), INTENT(out),OPTIONAL :: z0h !! Surface roughness (heat)
909	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names !! Names for deposit variables to be transported
910	!! from chemistry model by GCM to ORCHIDEE
911	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in !! Fields for deposit variables to be transported
912	!! from chemistry model by GCM to ORCHIDEE
913	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names!! Names for emission variables to be transported
914	!! to chemistry model by GCM from ORCHIDEE
915	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: fields_out !! Fields for emission variables to be transported
916	!! to chemistry model by GCM from ORCHIDEE
917
918
919	!! 0.5 Local variables
920	REAL(r_std),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain
921	REAL(r_std),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow
922	REAL(r_std),DIMENSION (kjpindex) :: zz0m !! Work array to keep zz0m
923	REAL(r_std),DIMENSION (kjpindex) :: zz0h !! Work array to keep zz0h
924	REAL(r_std),DIMENSION (kjpindex) :: zcdrag !! Work array for surface drag
925	REAL(r_std),DIMENSION (kjpindex) :: zcoastal !! Work array to keep coastal flow
926	REAL(r_std),DIMENSION (kjpindex) :: zriver !! Work array to keep river out flow
927	REAL(r_std),DIMENSION (kjpindex) :: znetco2 !! Work array to keep netco2flux
928	REAL(r_std),DIMENSION (kjpindex) :: zcarblu !! Work array to keep fco2_land_use
929	REAL(r_std),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad
930	REAL(r_std),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp
931	REAL(r_std),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new
932	REAL(r_std),DIMENSION (kjpindex) :: zqsurf !! Work array to keep qsurf
933	REAL(r_std),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo
934	REAL(r_std),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens
935	REAL(r_std),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat
936	REAL(r_std),DIMENSION (kjpindex) :: zemis !! Work array to keep emis
937	REAL(r_std),ALLOCATABLE, DIMENSION (:) :: soilth_lev !! Vertical soil axis for thermal scheme (m)
938	REAL(r_std),DIMENSION (kjpindex) :: q2m_loc !! Work array for q2m or qair
939	REAL(r_std),DIMENSION (kjpindex) :: t2m_loc !! Work array for t2m or temp_air
940	REAL(r_std),DIMENSION (:,:),ALLOCATABLE :: lon_scat !! The scattered values for longitude
941	REAL(r_std),DIMENSION (:,:),ALLOCATABLE :: lat_scat !! The scattered values for latitude
942	INTEGER(i_std) :: i, j, ik
943	INTEGER(i_std) :: ier
944	INTEGER(i_std) :: itau_sechiba
945	REAL(r_std) :: mx, zlev_mean
946	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_lon, tmp_lat, tmp_lev
947	INTEGER :: old_fileout !! old Logical Int for std IO output
948	REAL,ALLOCATABLE,DIMENSION(:,:) :: lalo_mpi
949	REAL,ALLOCATABLE,DIMENSION(:) :: zlev_glo1D
950	REAL(r_std),DIMENSION (kjpindex) :: landpoints !! Land point vector
951
952
953	OFF_LINE_MODE = .FALSE.
954	CALL ipslnlf_p(new_number=numout, old_number=old_fileout)
955	!
956	! Configuration of SSL specific parameters
957	!
958	CALL control_initialize(xrdt)
959
960	! Initialize specific write level
961	printlev_loc=get_printlev('instersurf')
962
963	IF ( printlev_loc>=1 ) WRITE(numout,*) 'Entering intersurf_initialize_gathered'
964
965	! Set the variable ok_q2m_t2m=true if q2m and t2m are present in the call from the gcm.
966	! If one of the variables are not present in this subroutine, set ok_q2m_t2m=.FALSE.
967	! Otherwise do not change the current value. Note that the current value for ok_q2m_t2m comes
968	! either from default initialization (true) or from intersurf_main_gathered.
969	IF (.NOT. PRESENT(q2m) .OR. .NOT. PRESENT(t2m)) THEN
970	ok_q2m_t2m=.FALSE.
971	END IF
972
973	IF (ok_q2m_t2m) THEN
974	t2m_loc=t2m
975	q2m_loc=q2m
976	ELSE
977	t2m_loc=temp_air
978	q2m_loc=qair
979	END IF
980
981	CALL intsurf_time( kjit, date0 )
982
983
984	CALL ioget_calendar (one_year, one_day)
985
986	IF (is_omp_root) THEN
987	ALLOCATE(lon_scat(iim_g,jj_nb))
988	ALLOCATE(lat_scat(iim_g,jj_nb))
989	ELSE
990	ALLOCATE(lon_scat(1,1))
991	ALLOCATE(lat_scat(1,1))
992	ENDIF
993
994	CALL init_WriteField_p
995
996	! Create the internal coordinate table
997	!
998	lalo(:,:) = latlon(:,:)
999	CALL gather(lalo,lalo_g)
1000	!
1001	!-
1002	!- Store variable to help describe the grid
1003	!- once the points are gathered.
1004	!-
1005	neighbours(:,:) = zneighbours(:,:)
1006	CALL gather(neighbours,neighbours_g)
1007	!
1008	resolution(:,:) = zresolution(:,:)
1009	CALL gather(resolution,resolution_g)
1010	!
1011	area(:) = resolution(:,1)*resolution(:,2)
1012	CALL gather(area,area_g)
1013	!
1014	!- Store the fraction of the continents only once so that the user
1015	!- does not change them afterwards.
1016	!
1017	contfrac(:) = zcontfrac(:)
1018	CALL gather(contfrac,contfrac_g)
1019	!
1020	!
1021	! Create the internal coordinate table
1022	!
1023	IF ( (.NOT.ALLOCATED(tmp_lon))) THEN
1024	ALLOCATE(tmp_lon(iim_g,jj_nb))
1025	ENDIF
1026	IF ( (.NOT.ALLOCATED(tmp_lat))) THEN
1027	ALLOCATE(tmp_lat(iim_g,jj_nb))
1028	ENDIF
1029	IF ( (.NOT.ALLOCATED(tmp_lev))) THEN
1030	ALLOCATE(tmp_lev(iim_g,jj_nb))
1031	ENDIF
1032	!
1033	! Either we have the scattered coordinates as arguments or
1034	! we have to do the work here.
1035	!
1036	IF (is_omp_root) THEN
1037	lon_scat(:,:)=zero
1038	lat_scat(:,:)=zero
1039	CALL scatter2D_mpi(lon_scat_g,lon_scat)
1040	CALL scatter2D_mpi(lat_scat_g,lat_scat)
1041	lon_scat(:,1)=lon_scat(:,2)
1042	lon_scat(:,jj_nb)=lon_scat(:,2)
1043	lat_scat(:,1)=lat_scat(iim_g,1)
1044	lat_scat(:,jj_nb)=lat_scat(1,jj_nb)
1045
1046	tmp_lon(:,:) = lon_scat(:,:)
1047	tmp_lat(:,:) = lat_scat(:,:)
1048
1049	IF (is_mpi_root) THEN
1050	lon_g(:,:) = lon_scat_g(:,:)
1051	lat_g(:,:) = lat_scat_g(:,:)
1052	ENDIF
1053	ENDIF
1054	!
1055
1056	!ANNE ici possibilite de calculer la grille du modele ? (lignes 965 ??
1057
1058
1059	!ANNE calcul zlev, conserve la version de Martial
1060	IF (is_root_prc) ALLOCATE(zlev_glo1D(nbp_glo))
1061	CALL gather(zlev,zlev_glo1D)
1062
1063	IF (is_root_prc) THEN
1064	DO ik=1, nbp_glo
1065	j = INT( (index_g(ik)-1) / iim_g ) + 1
1066	i = index_g(ik) - (j-1) * iim_g
1067	zlev_g(i,j) = zlev_glo1D(ik)
1068	ENDDO
1069	ENDIF
1070
1071	!Config Key = FORCE_CO2_VEG
1072	!Config Desc = Flag to force the value of atmospheric CO2 for vegetation.
1073	!Config If = Only in coupled mode
1074	!Config Def = FALSE
1075	!Config Help = If this flag is set to true, the ATM_CO2 parameter is used
1076	!Config to prescribe the atmospheric CO2.
1077	!Config This Flag is only use in couple mode.
1078	!Config Units = [FLAG]
1079	fatmco2=.FALSE.
1080	CALL getin_p('FORCE_CO2_VEG',fatmco2)
1081	!
1082	! Next flag is only use in couple mode with a gcm in intersurf.
1083	! In forced mode, it has already been read and set in driver.
1084	IF ( fatmco2 ) THEN
1085	!Config Key = ATM_CO2
1086	!Config IF = FORCE_CO2_VEG (only in coupled mode)
1087	!Config Desc = Value for atm CO2
1088	!Config Def = 350.
1089	!Config Help = Value to prescribe the atm CO2.
1090	!Config For pre-industrial simulations, the value is 286.2 .
1091	!Config 348. for 1990 year.
1092	!Config Units = [ppm]
1093	atmco2=350.
1094	CALL getin_p('ATM_CO2',atmco2)
1095	WRITE(numout,*) 'atmco2 ',atmco2
1096	ENDIF
1097
1098	CALL ioipslctrl_restini(kjit, date0, xrdt, rest_id, rest_id_stom, itau_offset, date0_shifted)
1099	itau_sechiba = kjit + itau_offset
1100
1101	!!- Initialize module for output with XIOS
1102	!
1103	! Get the vertical soil levels for the thermal scheme, to be used in xios_orchidee_init
1104	ALLOCATE(soilth_lev(ngrnd), stat=ier)
1105	IF (ier /= 0) THEN
1106	CALL ipslerr_p(3,'intersurf_main_gathered', 'Error in allocation of soilth_lev','','')
1107	END IF
1108	IF (hydrol_cwrr) THEN
1109	soilth_lev(1:ngrnd) = znt(:)
1110	ELSE
1111	soilth_lev(1:ngrnd) = thermosoilc_levels()
1112	END IF
1113
1114	CALL xios_orchidee_init( MPI_COMM_ORCH, &
1115	date0, year, month, day, &
1116	tmp_lon, tmp_lat, soilth_lev )
1117
1118	!- Initialize IOIPSL sechiba output files
1119	CALL ioipslctrl_history(iim_g, jj_nb, tmp_lon, tmp_lat, kindex, kjpindex, itau_sechiba, &
1120	date0_shifted, xrdt, hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC)
1121
1122	CALL bcast_omp(hist_id)
1123	CALL bcast_omp(hist2_id)
1124	CALL bcast_omp(hist_id_stom)
1125	CALL bcast_omp(hist_id_stom_IPCC)
1126
1127	! Count number of extra output fields to the GCM if it is not already done.
1128	IF (nb_fields_out == -1) THEN
1129	! nb_fields_out is not yet calculated. Do it now.
1130	! This means that the call is done directly from GCM.
1131	IF (PRESENT(field_out_names)) THEN
1132	nb_fields_out=SIZE(field_out_names)
1133	ELSE
1134	nb_fields_out=0
1135	ENDIF
1136	END IF
1137
1138	! Count number of extra input fields to the GCM if it is not already done.
1139	IF (nb_fields_in == -1) THEN
1140	! nb_fields_in is not yet calculated. Do it now.
1141	! This means that the call is done directly from GCM.
1142	IF (PRESENT(field_in_names)) THEN
1143	nb_fields_in=SIZE(field_in_names)
1144	ELSE
1145	nb_fields_in=0
1146	ENDIF
1147	END IF
1148
1149
1150	!
1151	!! Change to be in the orchidee context for XIOS
1152	!
1153	CALL xios_orchidee_change_context("orchidee")
1154
1155	!
1156	! Shift the time step to phase the two models
1157	!
1158	itau_sechiba = kjit + itau_offset
1159
1160	! Update the calendar in xios by sending the new time step
1161	! Special case : the model is only in initialization phase and the current itau_sechiba is not a real time step.
1162	! Therefor give itau_sechiba+1 to xios to have a correct time axis in output files.
1163	CALL xios_orchidee_update_calendar(itau_sechiba+1)
1164
1165	CALL intsurf_time( itau_sechiba, date0_shifted )
1166
1167	!
1168	! 1. Just change the units of some input fields
1169	!
1170	DO ik=1, kjpindex
1171
1172	zprecip_rain(ik) = precip_rain(ik)*xrdt
1173	zprecip_snow(ik) = precip_snow(ik)*xrdt
1174	zcdrag(ik) = cdrag(ik)
1175
1176	ENDDO
1177
1178	! Fields for deposit variables : to be transport from chemistry model by GCM to ORCHIDEE.
1179	! There are currently no fields to be transported into ORCHIDEE in this way
1180	DO i = 1, nb_fields_in
1181	WRITE(numout,*) i," Champ = ",TRIM(field_in_names(i))
1182	SELECT CASE(TRIM(field_in_names(i)))
1183	CASE DEFAULT
1184	CALL ipslerr_p (3,'intsurf_gathered', &
1185	'You ask in GCM an unknown field '//TRIM(field_in_names(i))//&
1186	' to give to ORCHIDEE for this specific version.',&
1187	'This model won''t be able to continue.', &
1188	'(check your tracer parameters in GCM)')
1189	END SELECT
1190	ENDDO
1191
1192	!
1193	! 3. save the grid
1194	!
1195	landpoints(:)=(/ ( REAL(ik), ik=1,kjpindex ) /)
1196	CALL histwrite_p(hist_id, 'LandPoints', itau_sechiba+1, landpoints, kjpindex, kindex)
1197	CALL histwrite_p(hist_id, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
1198	IF ( ok_stomate ) THEN
1199	CALL histwrite_p(hist_id_stom, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
1200	IF ( hist_id_stom_ipcc > 0 ) &
1201	CALL histwrite_p(hist_id_stom_IPCC, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
1202	ENDIF
1203	CALL histwrite_p(hist_id, 'Contfrac', itau_sechiba+1, contfrac, kjpindex, kindex)
1204
1205	! Syncronize output but only if flag ok_histsync is set to true
1206	IF (ok_histsync) THEN
1207	IF (is_omp_root .AND. hist_id > 0) THEN
1208	CALL histsync(hist_id)
1209	END IF
1210	END IF
1211
1212	IF ( hist2_id > 0 ) THEN
1213	CALL histwrite_p(hist2_id, 'LandPoints', itau_sechiba+1, landpoints, kjpindex, kindex)
1214	CALL histwrite_p(hist2_id, 'Areas', itau_sechiba+1, area, kjpindex, kindex)
1215	CALL histwrite_p(hist2_id, 'Contfrac', itau_sechiba+1, contfrac, kjpindex, kindex)
1216
1217	! Syncronize output but only if flag ok_histsync is set to true
1218	IF (ok_histsync .AND. is_omp_root) THEN
1219	CALL histsync(hist2_id)
1220	ENDIF
1221	ENDIF
1222	!
1223
1224	!
1225	! 4. call sechiba for continental points only
1226	!
1227	IF ( printlev_loc>=3 ) WRITE(numout,*) 'Calling sechiba_initialize'
1228
1229	CALL sechiba_initialize( &
1230	itau_sechiba, iim_g*jj_nb, kjpindex, kindex, &
1231	lalo, contfrac, neighbours, resolution, zlev, &
1232	u, v, qair, t2m_loc, temp_air, &
1233	petAcoef, peqAcoef, petBcoef, peqBcoef, &
1234	zprecip_rain, zprecip_snow, lwdown, swnet, swdown, &
1235	pb, rest_id, hist_id, hist2_id, &
1236	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1237	zcoastal, zriver, ztsol_rad, zvevapp, zqsurf, &
1238	zz0m, zz0h, zalbedo, zfluxsens, zfluxlat, zemis, &
1239	znetco2, zcarblu, ztemp_sol_new, zcdrag)
1240
1241	IF ( printlev_loc>=3 ) WRITE(numout,*) 'Out of sechiba_initialize'
1242
1243	!
1244	! 6. scatter output fields
1245	!
1246	DO ik=1, kjpindex
1247	z0m(ik) = zz0m(ik)
1248	coastalflow(ik) = zcoastal(ik)
1249	riverflow(ik) = zriver(ik)
1250	tsol_rad(ik) = ztsol_rad(ik)
1251	vevapp(ik) = zvevapp(ik)
1252	temp_sol_new(ik) = ztemp_sol_new(ik)
1253	qsurf(ik) = zqsurf(ik)
1254	albedo(ik,1) = zalbedo(ik,1)
1255	albedo(ik,2) = zalbedo(ik,2)
1256	fluxsens(ik) = zfluxsens(ik)
1257	fluxlat(ik) = zfluxlat(ik)
1258	emis(ik) = zemis(ik)
1259	cdrag(ik) = zcdrag(ik)
1260	ENDDO
1261	! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1262	IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1263
1264
1265	!
1266	! 7. Transform the water fluxes into Kg/m^2s and m^3/s
1267	!
1268	DO ik=1, kjpindex
1269	vevapp(ik) = vevapp(ik)/xrdt
1270	coastalflow(ik) = coastalflow(ik)/xrdt
1271	riverflow(ik) = riverflow(ik)/xrdt
1272	ENDDO
1273
1274	! Fields for emission variables : to be transport by GCM to chemistry model.
1275	DO i = 1, nb_fields_out
1276	SELECT CASE(TRIM(field_out_names(i)))
1277	CASE("fCO2_land")
1278	fields_out(:,i)=znetco2(:)
1279	CASE("fCO2_land_use")
1280	fields_out(:,i)=zcarblu(:)
1281	CASE DEFAULT
1282	CALL ipslerr_p (3,'intsurf_gathered', &
1283	'You ask from GCM an unknown field '//TRIM(field_out_names(i))//&
1284	' to ORCHIDEE for this specific version.',&
1285	'This model won''t be able to continue.', &
1286	'(check your tracer parameters in GCM)')
1287	END SELECT
1288	END DO
1289
1290	IF(is_root_prc) CALL getin_dump
1291	lstep_init_intersurf = .FALSE.
1292
1293	CALL ipslnlf_p(new_number=old_fileout)
1294	!
1295	!! Change back to be in the LMDZ context for XIOS
1296	!
1297	CALL xios_orchidee_change_context("LMDZ")
1298
1299	END SUBROUTINE intersurf_initialize_gathered
1300
1301
1302	!! =============================================================================================================================
1303	!! SUBROUTINE: intersurf_main_gathered
1304	!!
1305	!>\BRIEF Main subroutine to call ORCHIDEE from the gcm (LMDZ) using variables on a 1D grid with only land points.
1306	!!
1307	!! DESCRIPTION: This subroutine is the main interface for ORCHIDEE when it is called from the gcm (LMDZ).
1308	!! The variables are all gathered before entering this subroutine on the 1D grid with only landpoints.
1309	!!
1310	!! \n
1311	!_ ==============================================================================================================================
1312
1313	SUBROUTINE intersurf_main_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
1314	lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
1315	zlev, u, v, qair, temp_air, epot_air, ccanopy, &
1316	cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1317	precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
1318	vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
1319	tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m,lon_scat_g, lat_scat_g, q2m, t2m, z0h,&
1320	field_out_names, fields_out, field_in_names, fields_in, &
1321	coszang)
1322
1323	USE mod_orchidee_para
1324	IMPLICIT NONE
1325
1326	!! 0. Variable and parameter declaration
1327	!! 0.1 Input variables
1328	INTEGER(i_std),INTENT (in) :: kjit !! Time step number
1329	INTEGER(i_std),INTENT (in) :: iim_glo !! Dimension of global fields
1330	INTEGER(i_std),INTENT (in) :: jjm_glo !! Dimension of global fields
1331	INTEGER(i_std),INTENT (in) :: kjpindex !! Number of continental points
1332	REAL(r_std),INTENT (in) :: xrdt !! Time step in seconds
1333	LOGICAL, INTENT (in) :: lrestart_read !! Logical for _restart_ file to read
1334	LOGICAL, INTENT (in) :: lrestart_write !! Logical for _restart_ file to write'
1335	REAL(r_std), INTENT (in) :: date0 !! Date at which kjit = 0
1336	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: kindex !! Index for continental points
1337	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: u !! Lowest level wind speed
1338	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: v !! Lowest level wind speed
1339	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: zlev !! Height of first layer
1340	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: qair !! Lowest level specific humidity
1341	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation
1342	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation
1343	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: lwdown !! Down-welling long-wave flux
1344	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: swnet !! Net surface short-wave flux
1345	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: swdown !! Downwelling surface short-wave flux
1346	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature in Kelvin
1347	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: epot_air !! Air potential energy
1348	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: ccanopy !! CO2 concentration in the canopy
1349	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: petAcoef !! Coeficients A from the PBL resolution
1350	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: peqAcoef !! One for T and another for q
1351	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: petBcoef !! Coeficients B from the PBL resolution
1352	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: peqBcoef !! One for T and another for q
1353	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
1354	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: latlon !! Geographical coordinates
1355	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: zcontfrac !! Fraction of continent
1356	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: zneighbours !! neighbours
1357	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: zresolution !! size of the grid box
1358	REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN) :: lon_scat_g !! The scattered values for longitude
1359	REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN) :: lat_scat_g !! The scattered values for latitude
1360
1361	!! 0.2 Output variables
1362	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: z0m !! Surface roughness for momentum
1363	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: coastalflow !! Diffuse flow of water into the ocean (m^3/dt)
1364	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: riverflow !! Largest rivers flowing into the ocean (m^3/dt)
1365	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: tsol_rad !! Radiative surface temperature
1366	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: vevapp !! Total of evaporation
1367	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: temp_sol_new !! New soil temperature
1368	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: qsurf !! Surface specific humidity
1369	REAL(r_std),DIMENSION (kjpindex,2), INTENT(out) :: albedo !! Albedo
1370	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: fluxsens !! Sensible chaleur flux
1371	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: fluxlat !! Latent chaleur flux
1372	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: emis !! Emissivity
1373
1374	!! 0.3 Modified variables
1375	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: cdrag !! Cdrag
1376
1377	!! 0.4 Optional input and output variables
1378	REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL:: q2m !! Surface specific humidity
1379	REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL:: t2m !! Surface air temperature
1380	REAL(r_std),DIMENSION (kjpindex), INTENT(out),OPTIONAL:: z0h !! Surface roughness for heat
1381	REAL(r_std), DIMENSION(kjpindex), OPTIONAL, INTENT(in):: coszang !! Cosine of the solar zenith angle (unitless)
1382	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names !! Names for deposit variables to be transported
1383	!! from chemistry model by GCM to ORCHIDEE
1384	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in !! Fields for deposit variables to be transported
1385	!! from chemistry model by GCM to ORCHIDEE
1386	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names !! Names for emission variables to be transported
1387	!! to chemistry model by GCM from ORCHIDEE
1388	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: fields_out !! Fields for emission variables to be transported
1389	!! to chemistry model by GCM from ORCHIDEE
1390
1391	!! 0.5 Local variables
1392	REAL(r_std),DIMENSION (kjpindex) :: zccanopy !! Work array to keep ccanopy
1393	REAL(r_std),DIMENSION (kjpindex) :: zprecip_rain !! Work array to keep precip_rain
1394	REAL(r_std),DIMENSION (kjpindex) :: zprecip_snow !! Work array to keep precip_snow
1395	REAL(r_std),DIMENSION (kjpindex) :: zz0m, zz0h !! Work array to keep zz0m, zz0h
1396	REAL(r_std),DIMENSION (kjpindex) :: zcdrag !! Work array for surface drag
1397	REAL(r_std),DIMENSION (kjpindex) :: zcoastal !! Work array to keep coastal flow
1398	REAL(r_std),DIMENSION (kjpindex) :: zriver !! Work array to keep river out flow
1399	REAL(r_std),DIMENSION (kjpindex) :: znetco2 !! Work array to keep netco2flux
1400	REAL(r_std),DIMENSION (kjpindex) :: zcarblu !! Work array to keep fco2_land_use
1401	REAL(r_std),DIMENSION (kjpindex) :: ztsol_rad !! Work array to keep tsol_rad
1402	REAL(r_std),DIMENSION (kjpindex) :: zvevapp !! Work array to keep vevapp
1403	REAL(r_std),DIMENSION (kjpindex) :: ztemp_sol_new !! Work array to keep temp_sol_new
1404	REAL(r_std),DIMENSION (kjpindex) :: zqsurf !! Work array to keep qsurf
1405	REAL(r_std),DIMENSION (kjpindex,2) :: zalbedo !! Work array to keep albedo
1406	REAL(r_std),DIMENSION (kjpindex) :: zfluxsens !! Work array to keep fluxsens
1407	REAL(r_std),DIMENSION (kjpindex) :: zfluxlat !! Work array to keep fluxlat
1408	REAL(r_std),DIMENSION (kjpindex) :: zemis !! Work array to keep emis
1409	REAL(r_std),DIMENSION (kjpindex) :: zcoszang !! Work array to keep coszang
1410	REAL(r_std),ALLOCATABLE, DIMENSION (:) :: soilth_lev !! Vertical soil axis for thermal scheme (m)
1411	REAL(r_std),DIMENSION (kjpindex) :: q2m_loc !! Work array for q2m or qair
1412	REAL(r_std),DIMENSION (kjpindex) :: t2m_loc !! Work array for t2m or temp_air
1413	REAL(r_std),DIMENSION (:,:),ALLOCATABLE :: lon_scat !! The scattered values for longitude
1414	REAL(r_std),DIMENSION (:,:),ALLOCATABLE :: lat_scat !! The scattered values for latitude
1415	INTEGER(i_std) :: i, j, ik
1416	INTEGER(i_std) :: ier
1417	INTEGER(i_std) :: itau_sechiba
1418	REAL(r_std) :: mx, zlev_mean
1419	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_lon, tmp_lat, tmp_lev
1420	INTEGER :: old_fileout !! old Logical Int for std IO output
1421	REAL,ALLOCATABLE,DIMENSION(:,:) :: lalo_mpi
1422	REAL,ALLOCATABLE,DIMENSION(:) :: zlev_glo1D
1423	REAL(r_std),DIMENSION (kjpindex) :: landpoints !! Local landpoints vector
1424
1425
1426	CALL ipslnlf_p(new_number=numout, old_number=old_fileout)
1427
1428	IF (lstep_init_intersurf) THEN
1429	! Test if q2m and t2m are present
1430	IF (PRESENT(q2m) .AND. PRESENT(t2m)) THEN
1431	ok_q2m_t2m=.TRUE.
1432	ELSE
1433	ok_q2m_t2m=.FALSE.
1434	ENDIF
1435
1436	! Test if field_out_names and field_in_names are present and if so, count
1437	! the number of extra fields to exchange.
1438	IF (PRESENT(field_out_names)) THEN
1439	nb_fields_out=SIZE(field_out_names)
1440	ELSE
1441	nb_fields_out=0
1442	ENDIF
1443
1444	IF (PRESENT(field_in_names)) THEN
1445	nb_fields_in=SIZE(field_in_names)
1446	ELSE
1447	nb_fields_in=0
1448	ENDIF
1449
1450	CALL intersurf_initialize_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
1451	lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
1452	zlev, u, v, qair, temp_air, epot_air, &
1453	cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1454	precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
1455	vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
1456	tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m,lon_scat_g, lat_scat_g, q2m, t2m, zz0h,&
1457	field_out_names, fields_out, field_in_names, fields_in )
1458
1459	! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1460	IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1461	! Return from subroutine intersurf_main_gathered
1462	RETURN
1463	END IF
1464
1465	!
1466	!! Change to be in the orchidee context for XIOS
1467	!
1468	CALL xios_orchidee_change_context("orchidee")
1469
1470	!
1471	! Shift the time step to phase the two models
1472	!
1473	itau_sechiba = kjit + itau_offset
1474
1475	! Update the calendar in xios by sending the new time step
1476	CALL xios_orchidee_update_calendar(itau_sechiba)
1477
1478	CALL intsurf_time( itau_sechiba, date0_shifted )
1479
1480	! Test LMDz inputs
1481	IF (check_INPUTS) THEN
1482	IF ( sec <= 7200 ) THEN
1483	WRITE(numout,*) "sec = ",sec
1484	CALL WriteFieldI_p("Wu",u)
1485	WRITE(numout,*) "Lowest level wind speed :",u
1486	CALL WriteFieldI_p("Wv",v)
1487	WRITE(numout,*) "Lowest level wind speed :",v
1488	CALL WriteFieldI_p("Wzlev",zlev)
1489	WRITE(numout,*) "Height of first layer :",zlev
1490	CALL WriteFieldI_p("Wqair",qair)
1491	WRITE(numout,*) "Lowest level specific humidity :",qair
1492	CALL WriteFieldI_p("Wprecip_rain",precip_rain)
1493	WRITE(numout,*) "Rain precipitation :",precip_rain
1494	CALL WriteFieldI_p("Wprecip_snow",precip_snow)
1495	WRITE(numout,*) "Snow precipitation :",precip_snow
1496	CALL WriteFieldI_p("Wlwdown",lwdown)
1497	WRITE(numout,*) "Down-welling long-wave flux :",lwdown
1498	CALL WriteFieldI_p("Wswnet",swnet)
1499	WRITE(numout,*) "Net surface short-wave flux :",swnet
1500	CALL WriteFieldI_p("Wswdown",swdown)
1501	WRITE(numout,*) "Downwelling surface short-wave flux :",swdown
1502	CALL WriteFieldI_p("Wtemp_air",temp_air)
1503	WRITE(numout,*) "Air temperature in Kelvin :",temp_air
1504	CALL WriteFieldI_p("Wepot_air",epot_air)
1505	WRITE(numout,*) "Air potential energy :",epot_air
1506	CALL WriteFieldI_p("Wccanopy",ccanopy)
1507	WRITE(numout,*) "CO2 concentration in the canopy :",ccanopy
1508	CALL WriteFieldI_p("WpetAcoef",petAcoef)
1509	WRITE(numout,*) "Coeficients A from the PBL resolution :",petAcoef
1510	CALL WriteFieldI_p("WpeqAcoef",peqAcoef)
1511	WRITE(numout,*) "One for T and another for q :",peqAcoef
1512	CALL WriteFieldI_p("WpetBcoef",petBcoef)
1513	WRITE(numout,*) "Coeficients B from the PBL resolution :",petBcoef
1514	CALL WriteFieldI_p("WpeqBcoef",peqBcoef)
1515	WRITE(numout,*) "One for T and another for q :",peqBcoef
1516	CALL WriteFieldI_p("Wcdrag",cdrag)
1517	WRITE(numout,*) "Cdrag :",cdrag
1518	CALL WriteFieldI_p("Wpb",pb)
1519	WRITE(numout,*) "Surface pressure :",pb
1520	CALL WriteFieldI_p("Wzcontfrac",zcontfrac)
1521	WRITE(numout,*) "Fraction of continent :",zcontfrac
1522	IF ( ok_q2m_t2m) THEN
1523	CALL WriteFieldI_p("Wq2m",q2m)
1524	WRITE(numout,*) "Surface specific humidity :",q2m
1525	CALL WriteFieldI_p("Wt2m",t2m)
1526	WRITE(numout,*) "Surface air temperature :",t2m
1527	ENDIF
1528	ENDIF
1529	ENDIF
1530
1531	!
1532	! 1. Just change the units of some input fields
1533	!
1534	DO ik=1, kjpindex
1535
1536	zprecip_rain(ik) = precip_rain(ik)*xrdt
1537	zprecip_snow(ik) = precip_snow(ik)*xrdt
1538	zcdrag(ik) = cdrag(ik)
1539
1540	ENDDO
1541
1542	!>> VOC in coupled mode
1543	IF ( PRESENT(coszang) ) THEN
1544	zcoszang(:) = coszang(:)
1545	ELSE
1546	zcoszang(:) = zero
1547	ENDIF
1548
1549	IF (check_INPUTS) THEN
1550	WRITE(numout,*) "Intersurf_main_gathered :"
1551	WRITE(numout,*) "Time step number = ",kjit
1552	WRITE(numout,*) "Dimension of input fields for local mpi process = ",iim_g, jj_nb
1553	WRITE(numout,*) "Number of continental points = ",kjpindex
1554	WRITE(numout,*) "Time step in seconds = ",xrdt
1555	WRITE(numout,*) "Logical for _restart_ file to read, write = ",lrestart_read,lrestart_write
1556	WRITE(numout,*) "Date at which kjit = 0 = ",date0
1557	WRITE(numout,*) "offset for OMP = ",offset_omp
1558	WRITE(numout,*) "Index for continental points = ",kindex
1559	IF (is_omp_root) THEN
1560	WRITE(numout,*) "ROOT OMP, Index global MPI : ",kindex_mpi(:)
1561	ENDIF
1562	IF (is_root_prc) THEN
1563	WRITE(numout,*) "ROOT global, Index global : ",index_g(:)
1564	ENDIF
1565	WRITE(numout,*) "Lowest level wind speed North = ",u
1566	WRITE(numout,*) "Lowest level wind speed East = ",v
1567	WRITE(numout,*) "Height of first layer = ",zlev
1568	WRITE(numout,*) "Lowest level specific humidity = ",qair
1569	WRITE(numout,*) "Rain precipitation = ",zprecip_rain
1570	WRITE(numout,*) "Snow precipitation = ",zprecip_snow
1571	WRITE(numout,*) "Down-welling long-wave flux = ",lwdown
1572	WRITE(numout,*) "Net surface short-wave flux = ",swnet
1573	WRITE(numout,*) "Downwelling surface short-wave flux = ",swdown
1574	WRITE(numout,*) "Air temperature in Kelvin = ",temp_air
1575	WRITE(numout,*) "Air potential energy = ",epot_air
1576	WRITE(numout,*) "CO2 concentration in the canopy = ",ccanopy
1577	WRITE(numout,*) "Coeficients A from the PBL resolution = ",petAcoef
1578	WRITE(numout,*) "One for T and another for q = ",peqAcoef
1579	WRITE(numout,*) "Coeficients B from the PBL resolution = ",petBcoef
1580	WRITE(numout,*) "One for T and another for q = ",peqBcoef
1581	WRITE(numout,*) "Cdrag = ",zcdrag
1582	WRITE(numout,*) "Surface pressure = ",pb
1583	WRITE(numout,*) "Geographical land coordinates lon = ", lalo(:,2)
1584	WRITE(numout,*) "Geographical land coordinates lat = ", lalo(:,1)
1585	WRITE(numout,*) "Fraction of continent in the grid = ",zcontfrac
1586	ENDIF
1587
1588
1589	! Fields for deposit variables : to be transport from chemistry model by GCM to ORCHIDEE.
1590	DO i = 1, nb_fields_in
1591	WRITE(numout,*) i," Champ = ",TRIM(field_in_names(i))
1592	SELECT CASE(TRIM(field_in_names(i)))
1593	CASE DEFAULT
1594	CALL ipslerr_p (3,'intsurf_gathered', &
1595	'You ask in GCM an unknown field '//TRIM(field_in_names(i))//&
1596	' to give to ORCHIDEE for this specific version.',&
1597	'This model won''t be able to continue.', &
1598	'(check your tracer parameters in GCM)')
1599	END SELECT
1600	ENDDO
1601
1602	!
1603	! 2. modification of co2
1604	!
1605	IF ( fatmco2 ) THEN
1606	zccanopy(:) = atmco2
1607	WRITE (numout,*) 'Modification of the ccanopy value. CO2 = ',atmco2
1608	ELSE
1609	zccanopy(:) = ccanopy(:)
1610	ENDIF
1611
1612	!
1613	! 4. call sechiba for continental points only
1614	!
1615	IF ( printlev_loc>=3 ) WRITE(numout,*) 'Calling sechiba'
1616
1617
1618	IF (ok_q2m_t2m) THEN
1619	t2m_loc=t2m
1620	q2m_loc=q2m
1621	ELSE
1622	t2m_loc=temp_air
1623	q2m_loc=qair
1624	END IF
1625
1626	CALL sechiba_main (itau_sechiba, iim_g*jj_nb, kjpindex, kindex, &
1627	lrestart_read, lrestart_write, &
1628	lalo, contfrac, neighbours, resolution, &
1629	zlev, u, v, qair, q2m_loc, t2m_loc, temp_air, epot_air, zccanopy, &
1630	zcdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1631	zprecip_rain ,zprecip_snow, lwdown, swnet, swdown, zcoszang, pb, &
1632	zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
1633	ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0m, zz0h, &
1634	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC )
1635
1636	IF ( printlev_loc>=3 ) WRITE(numout,*) 'out of SECHIBA'
1637
1638	!
1639	! 6. scatter output fields
1640	!
1641	DO ik=1, kjpindex
1642	z0m(ik) = zz0m(ik)
1643	coastalflow(ik) = zcoastal(ik)
1644	riverflow(ik) = zriver(ik)
1645	tsol_rad(ik) = ztsol_rad(ik)
1646	vevapp(ik) = zvevapp(ik)
1647	temp_sol_new(ik) = ztemp_sol_new(ik)
1648	qsurf(ik) = zqsurf(ik)
1649	albedo(ik,1) = zalbedo(ik,1)
1650	albedo(ik,2) = zalbedo(ik,2)
1651	fluxsens(ik) = zfluxsens(ik)
1652	fluxlat(ik) = zfluxlat(ik)
1653	emis(ik) = zemis(ik)
1654	cdrag(ik) = zcdrag(ik)
1655	ENDDO
1656	! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1657	IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1658
1659	CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /))
1660	CALL xios_orchidee_send_field("areas", area)
1661	CALL xios_orchidee_send_field("contfrac",contfrac)
1662	CALL xios_orchidee_send_field("temp_air",temp_air)
1663	CALL xios_orchidee_send_field("qair",qair)
1664	CALL xios_orchidee_send_field("swnet",swnet)
1665	CALL xios_orchidee_send_field("swdown",swdown)
1666	CALL xios_orchidee_send_field("pb",pb)
1667
1668
1669	IF (ok_q2m_t2m) THEN
1670	CALL xios_orchidee_send_field("t2m",t2m)
1671	CALL xios_orchidee_send_field("q2m",q2m)
1672	ELSE
1673	CALL xios_orchidee_send_field("t2m",temp_air)
1674	CALL xios_orchidee_send_field("q2m",qair)
1675	ENDIF
1676
1677	IF ( .NOT. almaoutput ) THEN
1678	!
1679	! scattered during the writing
1680	!
1681	CALL histwrite_p (hist_id, 'evap', itau_sechiba, zvevapp, kjpindex, kindex)
1682	CALL histwrite_p (hist_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
1683	CALL histwrite_p (hist_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
1684	!
1685	CALL histwrite_p (hist_id, 'temp_sol', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1686
1687	CALL histwrite_p (hist_id, 'tsol_max', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1688	CALL histwrite_p (hist_id, 'tsol_min', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1689	CALL histwrite_p (hist_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
1690	CALL histwrite_p (hist_id, 'fluxlat', itau_sechiba, zfluxlat, kjpindex, kindex)
1691	CALL histwrite_p (hist_id, 'swnet', itau_sechiba, swnet, kjpindex, kindex)
1692	CALL histwrite_p (hist_id, 'swdown', itau_sechiba, swdown, kjpindex, kindex)
1693	CALL histwrite_p (hist_id, 'alb_vis', itau_sechiba, zalbedo(:,1), kjpindex, kindex)
1694	CALL histwrite_p (hist_id, 'alb_nir', itau_sechiba, zalbedo(:,2), kjpindex, kindex)
1695	CALL histwrite_p (hist_id, 'tair', itau_sechiba, temp_air, kjpindex, kindex)
1696	CALL histwrite_p (hist_id, 'qair', itau_sechiba, qair, kjpindex, kindex)
1697	IF (ok_q2m_t2m) THEN
1698	CALL histwrite_p (hist_id, 't2m', itau_sechiba, t2m, kjpindex, kindex)
1699	CALL histwrite_p (hist_id, 'q2m', itau_sechiba, q2m, kjpindex, kindex)
1700	ELSE
1701	CALL histwrite_p (hist_id, 't2m', itau_sechiba, temp_air, kjpindex, kindex)
1702	CALL histwrite_p (hist_id, 'q2m', itau_sechiba, qair, kjpindex, kindex)
1703	ENDIF
1704	!
1705	IF ( hist2_id > 0 ) THEN
1706	CALL histwrite_p (hist2_id, 'evap', itau_sechiba, zvevapp, kjpindex, kindex)
1707	CALL histwrite_p (hist2_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
1708	CALL histwrite_p (hist2_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
1709	!
1710	CALL histwrite_p (hist2_id, 'temp_sol', itau_sechiba, temp_sol_new, kjpindex, kindex)
1711	CALL histwrite_p (hist2_id, 'tsol_max', itau_sechiba, temp_sol_new, kjpindex, kindex)
1712	CALL histwrite_p (hist2_id, 'tsol_min', itau_sechiba, temp_sol_new, kjpindex, kindex)
1713	CALL histwrite_p (hist2_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
1714	CALL histwrite_p (hist2_id, 'fluxlat', itau_sechiba, zfluxlat, kjpindex, kindex)
1715	CALL histwrite_p (hist2_id, 'swnet', itau_sechiba, swnet, kjpindex, kindex)
1716	CALL histwrite_p (hist2_id, 'swdown', itau_sechiba, swdown, kjpindex, kindex)
1717	CALL histwrite_p (hist2_id, 'alb_vis', itau_sechiba, zalbedo(:,1), kjpindex, kindex)
1718	CALL histwrite_p (hist2_id, 'alb_nir', itau_sechiba, zalbedo(:,2), kjpindex, kindex)
1719	CALL histwrite_p (hist2_id, 'tair', itau_sechiba, temp_air, kjpindex, kindex)
1720	CALL histwrite_p (hist2_id, 'qair', itau_sechiba, qair, kjpindex, kindex)
1721	IF (ok_q2m_t2m) THEN
1722	CALL histwrite_p (hist2_id, 't2m', itau_sechiba, t2m, kjpindex, kindex)
1723	CALL histwrite_p (hist2_id, 'q2m', itau_sechiba, q2m, kjpindex, kindex)
1724	ELSE
1725	CALL histwrite_p (hist2_id, 't2m', itau_sechiba, temp_air, kjpindex, kindex)
1726	CALL histwrite_p (hist2_id, 'q2m', itau_sechiba, qair, kjpindex, kindex)
1727	ENDIF
1728	ENDIF
1729	ELSE
1730	CALL histwrite_p (hist_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
1731	CALL histwrite_p (hist_id, 'SWnet', itau_sechiba, swnet, kjpindex, kindex)
1732	CALL histwrite_p (hist_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
1733	CALL histwrite_p (hist_id, 'Qle', itau_sechiba, zfluxlat, kjpindex, kindex)
1734	CALL histwrite_p (hist_id, 'AvgSurfT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1735	CALL histwrite_p (hist_id, 'RadT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1736	!
1737	IF ( hist2_id > 0 ) THEN
1738	CALL histwrite_p (hist2_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
1739	CALL histwrite_p (hist2_id, 'SWnet', itau_sechiba, swnet, kjpindex, kindex)
1740	CALL histwrite_p (hist2_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
1741	CALL histwrite_p (hist2_id, 'Qle', itau_sechiba, zfluxlat, kjpindex, kindex)
1742	CALL histwrite_p (hist2_id, 'AvgSurfT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1743	CALL histwrite_p (hist2_id, 'RadT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1744	ENDIF
1745	ENDIF
1746
1747	! Syncronize output but only if flag ok_histsync is set to true
1748	IF (ok_histsync .AND. is_omp_root) THEN
1749	IF ( (dw .EQ. xrdt) .AND. hist_id > 0 ) THEN
1750	CALL histsync(hist_id)
1751	ENDIF
1752	ENDIF
1753
1754	!
1755	! 7. Transform the water fluxes into Kg/m^2s and m^3/s
1756	!
1757	DO ik=1, kjpindex
1758
1759	vevapp(ik) = vevapp(ik)/xrdt
1760	coastalflow(ik) = coastalflow(ik)/xrdt
1761	riverflow(ik) = riverflow(ik)/xrdt
1762
1763	ENDDO
1764
1765	! Fields for emission variables : to be transport by GCM to chemistry model.
1766	DO i = 1, nb_fields_out
1767	SELECT CASE(TRIM(field_out_names(i)))
1768	CASE("fCO2_land")
1769	fields_out(:,i)=znetco2(:)
1770	CASE("fCO2_land_use")
1771	fields_out(:,i)=zcarblu(:)
1772	CASE DEFAULT
1773	CALL ipslerr_p (3,'intsurf_gathered', &
1774	& 'You ask from GCM an unknown field '//TRIM(field_out_names(i))//&
1775	& ' to ORCHIDEE for this specific version.',&
1776	& 'This model won''t be able to continue.', &
1777	& '(check your tracer parameters in GCM)')
1778	END SELECT
1779	ENDDO
1780
1781	IF (printlev_loc >=3) WRITE (numout,*) ' intersurf_main done '
1782
1783	!Test ORCHIDEE Outputs
1784	IF (check_INPUTS) THEN
1785	IF ( sec <= 7200 ) THEN
1786	CALL WriteFieldI_p("Wcdrag",cdrag)
1787	WRITE(numout,*) "Cdrag :",cdrag
1788	CALL WriteFieldI_p("Wz0m",z0m)
1789	WRITE(numout,*) "Surface roughness :",z0m
1790	CALL WriteFieldI_p("Wcoastalflow",coastalflow)
1791	WRITE(numout,*) "Diffuse flow of water into the ocean (m^3/dt) :",coastalflow
1792	CALL WriteFieldI_p("Wriverflow",riverflow)
1793	WRITE(numout,*) "Largest rivers flowing into the ocean (m^3/dt) :",riverflow
1794	CALL WriteFieldI_p("Wtsol_rad",tsol_rad)
1795	WRITE(numout,*) "Radiative surface temperature :",tsol_rad
1796	CALL WriteFieldI_p("Wvevapp",vevapp)
1797	WRITE(numout,*) "Total of evaporation :",vevapp
1798	CALL WriteFieldI_p("Wtemp_sol_new",temp_sol_new)
1799	WRITE(numout,*) "New soil temperature :",temp_sol_new
1800	CALL WriteFieldI_p("Wqsurf",qsurf)
1801	WRITE(numout,*) "Surface specific humidity :",qsurf
1802	CALL WriteFieldI_p("Walbedo_nir",albedo(:,1))
1803	WRITE(numout,*) "Albedo nir:",albedo(:,1)
1804	CALL WriteFieldI_p("Walbedo_vis",albedo(:,2))
1805	WRITE(numout,*) "Albedo vir :",albedo(:,2)
1806	CALL WriteFieldI_p("Wfluxsens",fluxsens)
1807	WRITE(numout,*) "Sensible chaleur flux :",fluxsens
1808	CALL WriteFieldI_p("Wfluxlat",fluxlat)
1809	WRITE(numout,*) "Latent chaleur flux :",fluxlat
1810	CALL WriteFieldI_p("Wemis",emis)
1811	WRITE(numout,*) "Emissivity :",emis
1812	ENDIF
1813	ENDIF
1814	!
1815	CALL ipslnlf_p(new_number=old_fileout)
1816	!
1817
1818	!
1819	!! Finalize the XIOS orchidee context if it is the last call
1820	!
1821	IF (lrestart_write) THEN
1822	CALL xios_orchidee_context_finalize
1823	END IF
1824	!
1825	!! Change back to be in the LMDZ context for XIOS
1826	!
1827	CALL xios_orchidee_change_context("LMDZ")
1828
1829
1830	END SUBROUTINE intersurf_main_gathered
1831
1832
1833
1834	!! =============================================================================================================================
1835	!! SUBROUTINE: intsurf_time
1836	!!
1837	!>\BRIEF Initalize and update time information
1838	!!
1839	!! DESCRIPTION: Initialize and update time information. This subroutine is called in the initialization phase and at each
1840	!! time step from the different intersurf subroutines.
1841	!!
1842	!! \n
1843	!_ ==============================================================================================================================
1844
1845	SUBROUTINE intsurf_time(istp, date0)
1846
1847	IMPLICIT NONE
1848
1849	INTEGER(i_std), INTENT(in) :: istp !! Time step of the restart file
1850	REAL(r_std), INTENT(in) :: date0 !! The date at which itau = 0
1851
1852	IF (lstep_init_intersurf) THEN
1853	CALL ioget_calendar(calendar_str)
1854	CALL ioget_calendar(one_year, one_day)
1855	CALL tlen2itau('1Y',dt_sechiba,date0,year_length)
1856	IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN
1857	year_spread=un
1858	ELSE
1859	year_spread = one_year/365.2425
1860	ENDIF
1861
1862	IF (printlev_loc >=3) THEN
1863	write(numout,*) "calendar_str =",calendar_str
1864	write(numout,*) "one_year=",one_year,", one_day=",one_day
1865	write(numout,*) "dt_sechiba=",dt_sechiba,", date0=",date0,", year_length=",year_length,", year_spread=",year_spread
1866	ENDIF
1867	ENDIF
1868
1869	IF ( printlev_loc >=4 ) WRITE(numout,*) "---"
1870	! Dans diffuco (ie date0 == date0_shift !!)
1871
1872	IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN
1873	!
1874	! Get Julian date
1875	in_julian = itau2date(istp, date0, dt_sechiba)
1876
1877	! Real date
1878	CALL ju2ymds (in_julian, year, month, day, sec)
1879	!!$ jur=zero
1880	!!$ julian_diff = in_julian
1881	!!$ month_len = ioget_mon_len (year,month)
1882	!!$ IF ( printlev_loc >=4 ) THEN
1883	!!$ write(numout,*) "in_julian, jur, julian_diff=",in_julian, jur, julian_diff
1884	!!$ write(numout,*) 'DATE ymds', year, month,'(',month_len,'d)', day, sec, '-- stp --', istp
1885	!!$ ENDIF
1886
1887	! julian number for january, the first.
1888	CALL ymds2ju (year,1,1,zero, julian0)
1889	julian_diff = in_julian-julian0
1890	! real number of seconds
1891	! sec = (julian_diff-REAL(INT(julian_diff)))*one_day
1892	sec = NINT((julian_diff-REAL(INT(julian_diff)))*one_day)
1893	month_len = ioget_mon_len (year,month)
1894	IF ( printlev_loc >=4 ) THEN
1895	write(numout,*) "2 in_julian, julian0, julian_diff=",in_julian, julian0, julian_diff
1896	write(numout,*) '2 DATE ymds', year, month,'(',month_len,'d)', day, sec, '-- stp --', istp
1897	ENDIF
1898	ELSE
1899	!!$ in_julian = itau2date(istp-1, zero, dt)
1900	!!$ CALL ju2ymds (in_julian, year, month, day, sec)
1901	!!$ jur=zero
1902	!!$ julian_diff = in_julian
1903	!!$ month_len = ioget_mon_len (year,month)
1904	!!$ IF (printlev_loc >=4) THEN
1905	!!$ write(numout,*) "in_julian=",in_julian, jur, julian_diff
1906	!!$ write(numout,*) 'DATE ymds', year, month,'(',month_len,'d)', day, sec, '-- stp --', istp
1907	!!$ ENDIF
1908	!!$
1909	!!$
1910	!!$ CALL ymds2ju (year,1,1,zero, jur)
1911	!!$ julian_diff = in_julian-jur
1912	!!$ CALL ju2ymds (julian_diff, year, month, day, sec)
1913	!!$! sec = (julian_diff-REAL(INT(julian_diff)))*one_day
1914	!!$ sec = NINT((julian_diff-REAL(INT(julian_diff)))*one_day)
1915	!!$ month_len = ioget_mon_len (year,month)
1916	!!$ IF (printlev_loc >=4) THEN
1917	!!$ write(numout,*) "2 in_julian, jur, julian_diff=",in_julian, jur, julian_diff
1918	!!$ write(numout,*) '2 DATE ymds', year, month,'(',month_len,'d)', day, sec, '-- stp --', istp
1919	!!$ ENDIF
1920
1921	!MM
1922	!PB date0 = celui de Soenke (Ã tester avec un autre date0)
1923	! in_julian = itau2date(istp, 153116., dt)
1924	in_julian = itau2date(istp, date0, dt_sechiba)
1925	CALL itau2ymds(istp, dt_sechiba, year, month, day, sec)
1926	CALL ymds2ju (year,1,1,zero, julian0)
1927	julian_diff = in_julian - julian0
1928	month_len = ioget_mon_len (year,month)
1929	IF (printlev_loc >=4) THEN
1930	write(numout,*) "in_julian=",in_julian, julian0, julian_diff
1931	write(numout,*) 'DATE ymds', year, month,'(',month_len,'d)', day, sec, '-- stp --', istp
1932	ENDIF
1933	ENDIF
1934
1935	END SUBROUTINE intsurf_time
1936
1937	END MODULE intersurf

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/intersurf.f90 @ 5143

Download in other formats: