source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/intersurf.f90 @ 8543

Last change on this file since 8543 was 8032, checked in by josefine.ghattas, 13 months ago

As done in the trunk [8031] : Corrections to run in debug mode for unsctructured grid. Done by Y. Meurdesoif.
See ticket #923

  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 88.9 KB
Line 
1! ================================================================================================================================
2!  MODULE       : intersurf
3!
4!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
5!
6!  LICENCE      : IPSL (2006)
7!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF   Subroutines for interfacing between the driver dim2_driver and sechiba and between LMDZ and sechiba
10!!
11!!\n DESCRIPTION:    This module contains subroutines for the interfacing between ORCHIDEE and LMDZ and between the driver
12!!                   dim2_driver and the reste of ORCHIDEE. Follwoing subroutines exists:
13!!
14!!              - intersurf_initialize_2d       : Intitialization and call to sechiba_initialize, called from dim2_driver
15!!                                                before the first call to intersurf_main_2d.
16!!              - intersurf_main_2d             : Main subroutine will call sechiba_main, called from dim2_driver at each time-step
17!!              - init_intersurf                : Initialize grid information when coupling with LMDZ. This subroutine is called
18!!                                                from LMDZ before call to intersurf_initialize_gathered
19!!              - intersurf_initialize_gathered : Intitialization and call to sechiba_initialize, called from LMDZ
20!!                                                before the first call to intersurf_main_2d.
21!!              - intersurf_main_gathered       : Main subroutine will call sechiba_main, called from LMDZ at each time-step
22!!              - intsurf_time                  : Initialize and update time information, called in the initialization phase
23!!                                                and at each time step from the different intersurf subroutines.
24!!
25!! RECENT CHANGE(S):
26!!
27!! REFERENCE(S) : None
28!!
29!! SVN          :
30!! $HeadURL$
31!! $Date$
32!! $Revision$
33!! \n
34!_ ================================================================================================================================
35MODULE intersurf
36
37  USE IOIPSL
38  USE xios_orchidee
39  USE ioipsl_para 
40  USE defprec
41  USE sechiba
42  USE constantes
43  USE constantes_soil
44  USE control
45  USE pft_parameters
46  USE mod_orchidee_para
47  USE solar
48  USE grid
49  USE time, ONLY : one_year, one_day, dt_sechiba, time_initialize, time_nextstep, julian_diff
50  USE time, ONLY : year_end, month_end, day_end, sec_end
51  USE ioipslctrl, ONLY : ioipslctrl_history, ioipslctrl_restini, ioipslctrl_restclo, ok_histsync, max_hist_level, dw
52
53  IMPLICIT NONE
54
55  PRIVATE
56  PUBLIC :: Init_intersurf, intersurf_main, intersurf_main_2d, intersurf_main_gathered
57  PUBLIC :: intersurf_initialize_2d, intersurf_initialize_gathered, intersurf_clear
58
59
60  ! Interface called from LMDZ in older verisons
61  INTERFACE intersurf_main
62    MODULE PROCEDURE intersurf_main_gathered
63  END INTERFACE
64
65  LOGICAL, SAVE                                      :: lstep_init_intersurf=.TRUE.!! Initialisation has to be done one time
66!$OMP THREADPRIVATE(lstep_init_intersurf)
67  INTEGER(i_std), SAVE                               :: printlev_loc            !! Write level to this module
68!$OMP THREADPRIVATE(printlev_loc)
69  INTEGER(i_std), SAVE                               :: hist_id, rest_id        !! IDs for history and restart files
70!$OMP THREADPRIVATE(hist_id, rest_id)
71  INTEGER(i_std), SAVE                               :: hist2_id                !! ID for the second history files (Hi-frequency ?)
72!$OMP THREADPRIVATE(hist2_id)
73  INTEGER(i_std), SAVE                               :: hist_id_stom, hist_id_stom_IPCC, rest_id_stom !! Dito for STOMATE
74!$OMP THREADPRIVATE(hist_id_stom, hist_id_stom_IPCC, rest_id_stom)
75  INTEGER(i_std), SAVE                               :: itau_offset  !! This offset is used to phase the
76                                                                     !! calendar of the GCM or the driver.
77!$OMP THREADPRIVATE(itau_offset)
78  REAL(r_std), SAVE                                  :: date0_shifted
79!$OMP THREADPRIVATE(date0_shifted)
80  REAL(r_std), SAVE :: julian0                       !! first day of this year
81!$OMP THREADPRIVATE(julian0)
82  LOGICAL, SAVE                                      :: fatmco2           !! Flag to force the value of atmospheric CO2 for vegetation.
83!$OMP THREADPRIVATE(fatmco2)
84  REAL(r_std), SAVE                                  :: atmco2            !! atmospheric CO2
85!$OMP THREADPRIVATE(atmco2)
86  REAL(r_std), SAVE                                  :: coeff_rel         !! Coefficient for time filter on riverflow and coastalflow
87!$OMP THREADPRIVATE(coeff_rel)
88  REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE       :: riverflow_cpl0    !! Value from previous time step for riverflow
89!$OMP THREADPRIVATE(riverflow_cpl0)
90  REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE       :: coastalflow_cpl0  !! Value from previous time step for coastalflow
91!$OMP THREADPRIVATE(coastalflow_cpl0)
92  INTEGER(i_std), SAVE                               :: nbcf_from_LMDZ    !! Number of optional coupled fields(fields_out) received from LMDZ, for ESM configuration
93!$OMP THREADPRIVATE(nbcf_from_LMDZ)
94  INTEGER(i_std), SAVE                               :: nbcf_into_LMDZ    !! Number of optional coupled fields(fields_in) sent to LMDZ, for ESM configuration
95!$OMP THREADPRIVATE(nbcf_into_LMDZ)
96  CHARACTER(LEN=30), ALLOCATABLE, DIMENSION(:), SAVE  :: field_out_names_loc !! Names of optional fields(fields_out) received from LMDZ, for ESM configuration
97!$OMP THREADPRIVATE(field_out_names_loc)
98  CHARACTER(LEN=30), ALLOCATABLE, DIMENSION(:), SAVE  :: field_in_names_loc  !! Number of optional fields(fields_in) sent to LMDZ, for ESM configuration
99!$OMP THREADPRIVATE(field_in_names_loc)
100
101  PUBLIC lstep_init_intersurf
102 
103CONTAINS
104
105!!  =============================================================================================================================
106!! SUBROUTINE:    intersurf_initialize_2d
107!!
108!>\BRIEF          Initialization and call to sechiba_initialize
109!!
110!! DESCRIPTION:   Initialization of module variables, read options from parameter file, initialize output files and call to
111!!                sechiba_initialize.
112!!
113!!                This subroutine is called from dim2_driver before the first call to intersurf_main_2d.
114!!
115!! \n
116!_ ==============================================================================================================================
117
118  SUBROUTINE intersurf_initialize_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
119       lrestart_read, lrestart_write, lon, lat, zcontfrac, zresolution, date0, &
120       zlev, u, v, qair, temp_air, epot_air, ccanopy, &
121       cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
122       precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
123       vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
124       tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m)
125
126    IMPLICIT NONE
127
128    !! 0. Variable and parameter declaration
129    !! 0.1 Input variables
130    INTEGER(i_std),INTENT (in)                            :: kjit            !! Time step number
131    INTEGER(i_std),INTENT (in)                            :: iim, jjm        !! Dimension of input fields
132    INTEGER(i_std),INTENT (in)                            :: kjpindex        !! Number of continental points
133    REAL(r_std),INTENT (in)                               :: xrdt            !! Time step in seconds
134    LOGICAL, INTENT (in)                                 :: lrestart_read    !! Logical for _restart_ file to read
135    LOGICAL, INTENT (in)                                 :: lrestart_write   !! Logical for _restart_ file to write'
136    REAL(r_std), INTENT (in)                              :: date0           !! Date at which kjit = 0
137    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)      :: kindex          !! Index for continental points
138    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: u             !! Lowest level wind speed
139    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: v             !! Lowest level wind speed
140    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: zlev          !! Height of first layer
141    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: qair          !! Lowest level specific humidity
142    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: precip_rain   !! Rain precipitation
143    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: precip_snow   !! Snow precipitation
144    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: lwdown        !! Down-welling long-wave flux
145    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: swnet         !! Net surface short-wave flux
146    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: swdown        !! Downwelling surface short-wave flux
147    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: temp_air      !! Air temperature in Kelvin
148    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: epot_air      !! Air potential energy
149    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: ccanopy       !! CO2 concentration in the canopy
150    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: petAcoef      !! Coeficients A from the PBL resolution
151    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: peqAcoef      !! One for T and another for q
152    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: petBcoef      !! Coeficients B from the PBL resolution
153    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: peqBcoef      !! One for T and another for q
154    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: pb            !! Surface pressure
155    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: lon, lat      !! Geographical coordinates
156    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: zcontfrac     !! Fraction of continent in the grid
157    REAL(r_std),DIMENSION (iim,jjm,2), INTENT(in)           :: zresolution   !! resolution in x and y dimensions
158
159    !! 0.2 Output variables
160    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: z0m            !! Surface roughness
161    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: coastalflow   !! Diffuse flow of water into the ocean (m^3/s)
162    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: riverflow     !! Largest rivers flowing into the ocean (m^3/s)
163    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: tsol_rad      !! Radiative surface temperature
164    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: vevapp        !! Total of evaporation
165    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: temp_sol_new  !! New soil temperature
166    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: qsurf         !! Surface specific humidity
167    REAL(r_std),DIMENSION (iim,jjm,2), INTENT(out)          :: albedo        !! Albedo
168    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: fluxsens      !! Sensible chaleur flux
169    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: fluxlat       !! Latent chaleur flux
170    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: emis          !! Emissivity
171
172    !! 0.3 Modified variables
173    REAL(r_std),DIMENSION (iim,jjm), INTENT(inout)          :: cdrag         !! Cdrag
174
175    !! 0.4 Local variables
176    REAL(r_std),DIMENSION (kjpindex)                      :: zu            !! Work array to keep u
177    REAL(r_std),DIMENSION (kjpindex)                      :: zv            !! Work array to keep v
178    REAL(r_std),DIMENSION (kjpindex)                      :: zzlev         !! Work array to keep zlev
179    REAL(r_std),DIMENSION (kjpindex)                      :: zqair         !! Work array to keep qair
180    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
181    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
182    REAL(r_std),DIMENSION (kjpindex)                      :: zlwdown       !! Work array to keep lwdown
183    REAL(r_std),DIMENSION (kjpindex)                      :: zswnet        !! Work array to keep swnet
184    REAL(r_std),DIMENSION (kjpindex)                      :: zswdown       !! Work array to keep swdown
185    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_air     !! Work array to keep temp_air
186    REAL(r_std),DIMENSION (kjpindex)                      :: zepot_air     !! Work array to keep epot_air
187    REAL(r_std),DIMENSION (kjpindex)                      :: zccanopy      !! Work array to keep ccanopy
188    REAL(r_std),DIMENSION (kjpindex)                      :: zpetAcoef     !! Work array to keep petAcoef
189    REAL(r_std),DIMENSION (kjpindex)                      :: zpeqAcoef     !! Work array to keep peqAcoef
190    REAL(r_std),DIMENSION (kjpindex)                      :: zpetBcoef     !! Work array to keep petBcoef
191    REAL(r_std),DIMENSION (kjpindex)                      :: zpeqBcoef     !! Work array to keep peqVcoef
192    REAL(r_std),DIMENSION (kjpindex)                      :: zcdrag        !! Work array to keep cdrag
193    REAL(r_std),DIMENSION (kjpindex)                      :: zpb           !! Work array to keep pb
194    REAL(r_std),DIMENSION (kjpindex)                      :: zz0m, zz0h    !! Work array to keep zz0m, zz0h
195    REAL(r_std),DIMENSION (kjpindex)                      :: zcoastal      !! Work array to keep coastalflow (m^3/dt)
196    REAL(r_std),DIMENSION (kjpindex)                      :: zriver        !! Work array to keep riverflow (m^3/dt)
197    REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
198    REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
199    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
200    REAL(r_std),DIMENSION (kjpindex)                      :: zqsurf        !! Work array to keep qsurf
201    REAL(r_std),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
202    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
203    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
204    REAL(r_std),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis
205    INTEGER(i_std)                                       :: i, j, ik
206    INTEGER(i_std)                                       :: ier
207    INTEGER(i_std)                                       :: itau_sechiba
208    INTEGER                                              :: old_fileout   !! old Logical Int for std IO output
209
210
211    IF (printlev >= 2) WRITE(numout,*) 'Start intersurf_initialize_2d'
212    CALL ipslnlf_p(new_number=numout,old_number=old_fileout)
213
214    ! Initialize variables in module time
215    CALL time_initialize(kjit, date0, xrdt, "END")
216
217    !  Configuration of SSL specific parameters
218    CALL control_initialize
219
220    ! Initialize specific write level
221    printlev_loc=get_printlev('intersurf')
222   
223    OFF_LINE_MODE = .TRUE. 
224   
225    DO ik=1,kjpindex
226       
227       j = ((kindex(ik)-1)/iim) + 1
228       i = (kindex(ik) - (j-1)*iim)
229       
230       !- Create the internal coordinate table
231       !-
232       lalo(ik,1) = lat(i,j)
233       lalo(ik,2) = lon(i,j)
234       !
235       !- Store the fraction of the continents only once so that the user
236       !- does not change them afterwards.
237       !-
238       contfrac(ik) = zcontfrac(i,j)
239    ENDDO
240    CALL gather(contfrac,contfrac_g)
241    CALL gather(lalo,lalo_g)
242    CALL gather2D_mpi(lon,lon_g)
243    CALL gather2D_mpi(lat,lat_g)
244
245    CALL bcast(lalo_g)
246    CALL bcast(contfrac_g)
247   
248    CALL ioipslctrl_restini(kjit, date0, xrdt, rest_id, rest_id_stom, itau_offset, date0_shifted)
249    itau_sechiba = kjit + itau_offset
250   
251    !!- Initialize module for output with XIOS
252    !
253    CALL xios_orchidee_init( MPI_COMM_ORCH,                        &
254         date0,   year_end,  month_end,     day_end,  julian_diff, &
255         lon,     lat )
256
257    CALL sechiba_xios_initialize
258   
259    CALL xios_orchidee_close_definition
260    IF (printlev >= 2) WRITE(numout,*) 'After xios_orchidee_close_definition'
261
262    !- Initialize IOIPSL sechiba output files
263    CALL ioipslctrl_history(iim, jjm, lon, lat,  kindex, kjpindex, itau_sechiba, date0_shifted, xrdt, hist_id, &
264         hist2_id, hist_id_stom, hist_id_stom_IPCC)
265 
266    !
267    !  Shift the time step to phase the two models
268    !
269    itau_sechiba = kjit + itau_offset
270   
271    ! Update the calendar in xios by sending the new time step
272    ! Special case : the model is only in initialization phase and the current itau_sechiba is not a real time step.
273    ! Therefor give itau_sechiba+1 to xios to have a correct time axis in output files.
274    CALL xios_orchidee_update_calendar(itau_sechiba+1)
275
276    !
277    ! 1. gather input fields from kindex array
278    !    Warning : I'm not sure this interface with one dimension array is the good one
279    !
280    DO ik=1, kjpindex
281     
282       j = ((kindex(ik)-1)/iim) + 1
283       i = (kindex(ik) - (j-1)*iim)
284       
285       zu(ik)           = u(i,j)
286       zv(ik)           = v(i,j)
287       zzlev(ik)        = zlev(i,j)
288       zqair(ik)        = qair(i,j)
289       zprecip_rain(ik) = precip_rain(i,j)*xrdt
290       zprecip_snow(ik) = precip_snow(i,j)*xrdt
291       zlwdown(ik)      = lwdown(i,j)
292       zswnet(ik)       = swnet(i,j)
293       zswdown(ik)      = swdown(i,j)
294       ztemp_air(ik)    = temp_air(i,j)
295       zepot_air(ik)    = epot_air(i,j)
296       zccanopy(ik)     = ccanopy(i,j)
297       zpetAcoef(ik)    = petAcoef(i,j)
298       zpeqAcoef(ik)    = peqAcoef(i,j)
299       zpetBcoef(ik)    = petBcoef(i,j)
300       zpeqBcoef(ik)    = peqBcoef(i,j)
301       zcdrag(ik)       = cdrag(i,j)
302       zpb(ik)          = pb(i,j)
303       
304    ENDDO
305
306    !
307    ! 2. save the grid
308    !
309    CALL histwrite_p(hist_id, 'LandPoints',  itau_sechiba+1, (/ ( REAL(ik), ik=1,kjpindex ) /), kjpindex, kindex)
310    CALL histwrite_p(hist_id, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
311    IF ( ok_stomate ) THEN
312       CALL histwrite_p(hist_id_stom, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
313       CALL histwrite_p(hist_id_stom_IPCC, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
314    ENDIF
315   
316    CALL histwrite_p(hist_id, 'Contfrac',  itau_sechiba+1, contfrac, kjpindex, kindex)
317    IF ( is_omp_root .AND. hist_id > 0 ) THEN
318       ! Always syncronize output after initialization
319       CALL histsync(hist_id)
320    END IF
321   
322    CALL histwrite_p(hist2_id, 'LandPoints',  itau_sechiba+1, (/ ( REAL(ik), ik=1,kjpindex ) /), kjpindex, kindex)
323    CALL histwrite_p(hist2_id, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
324    CALL histwrite_p(hist2_id, 'Contfrac',  itau_sechiba+1, contfrac, kjpindex, kindex)
325    IF ( is_omp_root .AND. hist2_id > 0 ) THEN
326       ! Always syncronize output after initialization
327       CALL histsync(hist2_id)
328    ENDIF
329   
330    !
331    ! 3. call sechiba for continental points only
332    !
333    IF (printlev_loc >= 3) WRITE(numout,*) 'Before call to sechiba_initialize'
334   
335    CALL sechiba_initialize( &
336         itau_sechiba, iim*jjm,      kjpindex,      kindex,                     &
337         lalo,         contfrac,     neighbours,    resolution,  zzlev,         &
338         zu,           zv,           zqair,         ztemp_air,                  &
339         zpetAcoef,    zpeqAcoef,    zpetBcoef,     zpeqBcoef,                  &
340         zprecip_rain, zprecip_snow, zlwdown,       zswnet,      zswdown,       &
341         zpb,          rest_id,      hist_id,       hist2_id,                   &
342         rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                         &
343         zcoastal,     zriver,       ztsol_rad,     zvevapp,     zqsurf,        &
344         zz0m,         zz0h,         zalbedo,      zfluxsens,     zfluxlat,     &
345         zemis,        ztemp_sol_new,zcdrag)
346
347   
348    IF (printlev_loc >= 3) WRITE(numout,*) 'After call to sechiba_initialize'
349    !
350    ! 5. scatter output fields
351    !
352    z0m(:,:)           = undef_sechiba
353    coastalflow(:,:)  = undef_sechiba
354    riverflow(:,:)    = undef_sechiba
355    tsol_rad(:,:)     = undef_sechiba
356    vevapp(:,:)       = undef_sechiba
357    temp_sol_new(:,:) = undef_sechiba 
358    qsurf(:,:)        = undef_sechiba 
359    albedo(:,:,:)     = undef_sechiba
360    fluxsens(:,:)     = undef_sechiba
361    fluxlat(:,:)      = undef_sechiba
362    emis(:,:)         = undef_sechiba 
363    cdrag(:,:)        = undef_sechiba 
364   
365    DO ik=1, kjpindex
366       j = ((kindex(ik)-1)/iim) + 1
367       i = (kindex(ik) - (j-1)*iim)
368       
369       z0m(i,j)           = zz0m(ik)
370       coastalflow(i,j)  = zcoastal(ik)
371       riverflow(i,j)    = zriver(ik)
372       tsol_rad(i,j)     = ztsol_rad(ik)
373       vevapp(i,j)       = zvevapp(ik)
374       temp_sol_new(i,j) = ztemp_sol_new(ik)
375       qsurf(i,j)        = zqsurf(ik)
376       albedo(i,j,1)     = zalbedo(ik,1)
377       albedo(i,j,2)     = zalbedo(ik,2)
378       fluxsens(i,j)     = zfluxsens(ik)
379       fluxlat(i,j)      = zfluxlat(ik)
380       emis(i,j)         = zemis(ik)
381       cdrag(i,j)        = zcdrag(ik)
382    ENDDO
383
384    !
385    ! 6. Transform the water fluxes into Kg/m^2s and m^3/s
386    !
387    DO ik=1, kjpindex
388   
389       j = ((kindex(ik)-1)/iim) + 1
390       i = (kindex(ik) - (j-1)*iim)
391
392       vevapp(i,j) = vevapp(i,j)/xrdt
393       coastalflow(i,j) = coastalflow(i,j)/xrdt
394       riverflow(i,j) = riverflow(i,j)/xrdt
395
396    ENDDO
397
398    IF (is_root_prc) CALL getin_dump
399
400    lstep_init_intersurf = .FALSE.
401    CALL ipslnlf_p(new_number=old_fileout)
402    IF (printlev_loc >= 1) WRITE (numout,*) 'End intersurf_initialize_2d'
403
404  END SUBROUTINE intersurf_initialize_2d
405
406
407!!  =============================================================================================================================
408!! SUBROUTINE:    intersurf_main_2d
409!!
410!>\BRIEF          Main subroutine to call ORCHIDEE from dim2_driver using variables on a 2d grid.
411!!
412!! DESCRIPTION:   This subroutine is the main interface for ORCHIDEE when it is called from the offline driver dim2_driver.
413!!                The variables are all on the 2D grid including ocean points. intersurf_initialize_2d should be called before
414!!                this subroutine is called. This subroutine is called at each time step.
415!!
416!! \n
417!_ ==============================================================================================================================
418
419  SUBROUTINE intersurf_main_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
420       lrestart_read, lrestart_write, lon, lat, zcontfrac, zresolution, date0, &
421       zlev, u, v, qair, temp_air, epot_air, ccanopy, &
422       cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
423       precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
424       vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
425       tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, &
426       coszang)
427
428    IMPLICIT NONE
429
430    !! 0. Variable and parameter declaration
431    !! 0.1 Input variables
432    INTEGER(i_std),INTENT (in)                              :: kjit            !! Time step number
433    INTEGER(i_std),INTENT (in)                              :: iim, jjm        !! Dimension of input fields
434    INTEGER(i_std),INTENT (in)                              :: kjpindex        !! Number of continental points
435    REAL(r_std),INTENT (in)                                 :: xrdt            !! Time step in seconds
436    LOGICAL, INTENT (in)                                    :: lrestart_read   !! Logical for _restart_ file to read
437    LOGICAL, INTENT (in)                                    :: lrestart_write  !! Logical for _restart_ file to write'
438    REAL(r_std), INTENT (in)                                :: date0           !! Date at which kjit = 0
439    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)        :: kindex          !! Index for continental points
440    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: u             !! Lowest level wind speed
441    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: v             !! Lowest level wind speed
442    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: zlev          !! Height of first layer
443    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: qair          !! Lowest level specific humidity
444    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: precip_rain   !! Rain precipitation
445    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: precip_snow   !! Snow precipitation
446    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: lwdown        !! Down-welling long-wave flux
447    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: swnet         !! Net surface short-wave flux
448    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: swdown        !! Downwelling surface short-wave flux
449    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: coszang       !! Cosine of the solar zenith angle (unitless)
450    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: temp_air      !! Air temperature in Kelvin
451    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: epot_air      !! Air potential energy
452    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: ccanopy       !! CO2 concentration in the canopy
453    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: petAcoef      !! Coeficients A from the PBL resolution
454    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: peqAcoef      !! One for T and another for q
455    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: petBcoef      !! Coeficients B from the PBL resolution
456    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: peqBcoef      !! One for T and another for q
457    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: pb            !! Surface pressure
458    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: lon, lat      !! Geographical coordinates
459    REAL(r_std),DIMENSION (iim,jjm), INTENT(in)             :: zcontfrac     !! Fraction of continent in the grid
460    REAL(r_std),DIMENSION (iim,jjm,2), INTENT(in)           :: zresolution   !! resolution in x and y dimensions
461
462    !! 0.2 Output variables
463    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: z0m            !! Surface roughness for momemtum
464    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: coastalflow   !! Diffuse flow of water into the ocean (m^3/s)
465    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: riverflow     !! Largest rivers flowing into the ocean (m^3/s)
466    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: tsol_rad      !! Radiative surface temperature
467    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: vevapp        !! Total of evaporation
468    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: temp_sol_new  !! New soil temperature
469    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: qsurf         !! Surface specific humidity
470    REAL(r_std),DIMENSION (iim,jjm,2), INTENT(out)          :: albedo        !! Albedo
471    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: fluxsens      !! Sensible chaleur flux
472    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: fluxlat       !! Latent chaleur flux
473    REAL(r_std),DIMENSION (iim,jjm), INTENT(out)            :: emis          !! Emissivity
474
475    !! 0.3 Modified variables
476    REAL(r_std),DIMENSION (iim,jjm), INTENT(inout)          :: cdrag         !! Cdrag
477
478    !! 0.4 Local variables
479    REAL(r_std),DIMENSION (kjpindex)                      :: zu            !! Work array to keep u
480    REAL(r_std),DIMENSION (kjpindex)                      :: zv            !! Work array to keep v
481    REAL(r_std),DIMENSION (kjpindex)                      :: zzlev         !! Work array to keep zlev
482    REAL(r_std),DIMENSION (kjpindex)                      :: zqair         !! Work array to keep qair
483    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
484    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
485    REAL(r_std),DIMENSION (kjpindex)                      :: zlwdown       !! Work array to keep lwdown
486    REAL(r_std),DIMENSION (kjpindex)                      :: zswnet        !! Work array to keep swnet
487    REAL(r_std),DIMENSION (kjpindex)                      :: zswdown       !! Work array to keep swdown
488    REAL(r_std),DIMENSION (kjpindex)                      :: zcoszang      !! Work array to keep coszang
489    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_air     !! Work array to keep temp_air
490    REAL(r_std),DIMENSION (kjpindex)                      :: zepot_air     !! Work array to keep epot_air
491    REAL(r_std),DIMENSION (kjpindex)                      :: zccanopy      !! Work array to keep ccanopy
492    REAL(r_std),DIMENSION (kjpindex)                      :: zpetAcoef     !! Work array to keep petAcoef
493    REAL(r_std),DIMENSION (kjpindex)                      :: zpeqAcoef     !! Work array to keep peqAcoef
494    REAL(r_std),DIMENSION (kjpindex)                      :: zpetBcoef     !! Work array to keep petBcoef
495    REAL(r_std),DIMENSION (kjpindex)                      :: zpeqBcoef     !! Work array to keep peqVcoef
496    REAL(r_std),DIMENSION (kjpindex)                      :: zcdrag        !! Work array to keep cdrag
497    REAL(r_std),DIMENSION (kjpindex)                      :: zpb           !! Work array to keep pb
498    REAL(r_std),DIMENSION (kjpindex)                      :: zz0m, zz0h    !! Work array to keep zz0m, zz0h
499    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zveget        !! Work array to keep veget
500    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zlai          !! Work array to keep lai
501    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zheight       !! Work array to keep height
502    REAL(r_std),DIMENSION (kjpindex)                      :: zcoastal      !! Work array to keep coastalflow (m^3/dt)
503    REAL(r_std),DIMENSION (kjpindex)                      :: zriver        !! Work array to keep riverflow (m^3/dt)
504    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux (gC/m2/one_day)
505    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use (gC/m2/one_day)
506    REAL(r_std),DIMENSION (kjpindex)                      :: zcarbwh       !! Work array to keep fco2_wh (gC/m2/one_day)
507    REAL(r_std),DIMENSION (kjpindex)                      :: zcarbha       !! Work array to keep fco2_ha (gC/m2/one_day)
508    REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
509    REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
510    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
511    REAL(r_std),DIMENSION (kjpindex)                      :: zqsurf        !! Work array to keep qsurf
512    REAL(r_std),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
513    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
514    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
515    REAL(r_std),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis
516    INTEGER(i_std)                                        :: i, j, ik
517    INTEGER(i_std)                                        :: ier
518    INTEGER(i_std)                                        :: itau_sechiba
519    INTEGER                                               :: old_fileout   !! old Logical Int for std IO output
520
521    IF (printlev_loc >= 3) WRITE(numout,*) 'Start intersurf_main_2d'
522    CALL ipslnlf_p(new_number=numout,old_number=old_fileout)
523
524    !
525    !  Shift the time step to phase the two models
526    !
527    itau_sechiba = kjit + itau_offset
528    !
529
530    ! Update the calendar in xios by sending the new time step
531    CALL xios_orchidee_update_calendar(itau_sechiba)
532
533    ! Update the calendar and all time variables in module time
534    CALL time_nextstep(itau_sechiba)
535    !
536    ! 1. gather input fields from kindex array
537    !    Warning : I'm not sure this interface with one dimension array is the good one
538    !
539    DO ik=1, kjpindex
540     
541       j = ((kindex(ik)-1)/iim) + 1
542       i = (kindex(ik) - (j-1)*iim)
543       
544       zu(ik)           = u(i,j)
545       zv(ik)           = v(i,j)
546       zzlev(ik)        = zlev(i,j)
547       zqair(ik)        = qair(i,j)
548       zprecip_rain(ik) = precip_rain(i,j)*xrdt
549       zprecip_snow(ik) = precip_snow(i,j)*xrdt
550       zlwdown(ik)      = lwdown(i,j)
551       zswnet(ik)       = swnet(i,j)
552       zswdown(ik)      = swdown(i,j)
553       zcoszang(ik)     = coszang(i,j)
554       ztemp_air(ik)    = temp_air(i,j)
555       zepot_air(ik)    = epot_air(i,j)
556       zccanopy(ik)     = ccanopy(i,j)
557       zpetAcoef(ik)    = petAcoef(i,j)
558       zpeqAcoef(ik)    = peqAcoef(i,j)
559       zpetBcoef(ik)    = petBcoef(i,j)
560       zpeqBcoef(ik)    = peqBcoef(i,j)
561       zcdrag(ik)       = cdrag(i,j)
562       zpb(ik)          = pb(i,j)
563       
564    ENDDO
565
566    !
567    ! 3. call sechiba for continental points only
568    !
569    IF (printlev_loc >= 3) WRITE(numout,*) 'Before call to sechiba_main from intersurf_main_2d'
570
571    CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, &
572         lrestart_read, lrestart_write, &
573         lalo, contfrac, neighbours, resolution, &
574         zzlev, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
575         zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, &
576         zprecip_rain ,zprecip_snow,  zlwdown, zswnet, zswdown, zcoszang, zpb, &
577         zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
578         znetco2, zcarblu, zcarbwh, zcarbha, &
579         ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0m, zz0h,&
580         zveget, zlai, zheight, &
581         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC ) 
582   
583    IF (printlev_loc >= 3) WRITE(numout,*) 'After call to sechiba_main'
584
585    !
586    ! 5. scatter output fields
587    !
588    z0m(:,:)          = undef_sechiba
589    coastalflow(:,:)  = undef_sechiba
590    riverflow(:,:)    = undef_sechiba
591    tsol_rad(:,:)     = undef_sechiba
592    vevapp(:,:)       = undef_sechiba
593    temp_sol_new(:,:) = undef_sechiba 
594    qsurf(:,:)        = undef_sechiba 
595    albedo(:,:,:)     = undef_sechiba
596    fluxsens(:,:)     = undef_sechiba
597    fluxlat(:,:)      = undef_sechiba
598    emis(:,:)         = undef_sechiba 
599    cdrag(:,:)        = undef_sechiba 
600    !
601    DO ik=1, kjpindex
602     
603   
604       j = ((kindex(ik)-1)/iim) + 1
605       i = (kindex(ik) - (j-1)*iim)
606
607       z0m(i,j)           = zz0m(ik)
608       coastalflow(i,j)  = zcoastal(ik)
609       riverflow(i,j)    = zriver(ik)
610       tsol_rad(i,j)     = ztsol_rad(ik)
611       vevapp(i,j)       = zvevapp(ik)
612       temp_sol_new(i,j) = ztemp_sol_new(ik)
613       qsurf(i,j)        = zqsurf(ik)
614       albedo(i,j,1)     = zalbedo(ik,1)
615       albedo(i,j,2)     = zalbedo(ik,2)
616       fluxsens(i,j)     = zfluxsens(ik)
617       fluxlat(i,j)      = zfluxlat(ik)
618       emis(i,j)         = zemis(ik)
619       cdrag(i,j)        = zcdrag(ik)
620
621    ENDDO
622
623    CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /))
624    CALL xios_orchidee_send_field("areas", area)
625    CALL xios_orchidee_send_field("contfrac",contfrac)
626    CALL xios_orchidee_send_field("temp_air",ztemp_air)
627    CALL xios_orchidee_send_field("qair",zqair)
628    CALL xios_orchidee_send_field("swnet",zswnet)
629    CALL xios_orchidee_send_field("swdown",zswdown)
630    CALL xios_orchidee_send_field("pb",zpb)
631    CALL xios_orchidee_send_field("znetco2",znetco2/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
632    CALL xios_orchidee_send_field("zcarblu",zcarblu/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
633    CALL xios_orchidee_send_field("zcarbwh",zcarbwh/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
634    CALL xios_orchidee_send_field("zcarbha",zcarbha/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
635
636   
637    IF ( .NOT. almaoutput ) THEN
638       !
639       !  scattered during the writing
640       !
641       CALL histwrite_p (hist_id, 'evap',     itau_sechiba, zvevapp, kjpindex, kindex)
642       CALL histwrite_p (hist_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
643       CALL histwrite_p (hist_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
644       CALL histwrite_p (hist_id, 'temp_sol', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
645       CALL histwrite_p (hist_id, 'tsol_max', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
646       CALL histwrite_p (hist_id, 'tsol_min', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
647       CALL histwrite_p (hist_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
648       CALL histwrite_p (hist_id, 'fluxlat',  itau_sechiba, zfluxlat, kjpindex, kindex)
649       CALL histwrite_p (hist_id, 'swnet',    itau_sechiba, zswnet, kjpindex, kindex)
650       CALL histwrite_p (hist_id, 'swdown',   itau_sechiba, zswdown, kjpindex, kindex)
651       CALL histwrite_p (hist_id, 'alb_vis',  itau_sechiba, zalbedo(:,1), kjpindex, kindex)
652       CALL histwrite_p (hist_id, 'alb_nir',  itau_sechiba, zalbedo(:,2), kjpindex, kindex)
653       CALL histwrite_p (hist_id, 'tair',     itau_sechiba, ztemp_air, kjpindex, kindex)
654       CALL histwrite_p (hist_id, 'qair',     itau_sechiba, zqair, kjpindex, kindex)
655       CALL histwrite_p (hist2_id, 'evap',     itau_sechiba, zvevapp, kjpindex, kindex)
656       CALL histwrite_p (hist2_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
657       CALL histwrite_p (hist2_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
658       CALL histwrite_p (hist2_id, 'temp_sol', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
659       CALL histwrite_p (hist2_id, 'tsol_max', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
660       CALL histwrite_p (hist2_id, 'tsol_min', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
661       CALL histwrite_p (hist2_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
662       CALL histwrite_p (hist2_id, 'fluxlat',  itau_sechiba, zfluxlat, kjpindex, kindex)
663       CALL histwrite_p (hist2_id, 'swnet',    itau_sechiba, zswnet, kjpindex, kindex)
664       CALL histwrite_p (hist2_id, 'swdown',   itau_sechiba, zswdown, kjpindex, kindex)
665       CALL histwrite_p (hist2_id, 'alb_vis',  itau_sechiba, zalbedo(:,1), kjpindex, kindex)
666       CALL histwrite_p (hist2_id, 'alb_nir',  itau_sechiba, zalbedo(:,2), kjpindex, kindex)
667       CALL histwrite_p (hist2_id, 'tair',     itau_sechiba, ztemp_air, kjpindex, kindex)
668       CALL histwrite_p (hist2_id, 'qair',     itau_sechiba, zqair, kjpindex, kindex)
669    ELSE
670       CALL histwrite_p (hist_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
671       CALL histwrite_p (hist_id, 'SWnet',    itau_sechiba, zswnet, kjpindex, kindex)
672       CALL histwrite_p (hist_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
673       CALL histwrite_p (hist_id, 'Qle',  itau_sechiba, zfluxlat, kjpindex, kindex)
674       CALL histwrite_p (hist_id, 'AvgSurfT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
675       CALL histwrite_p (hist_id, 'RadT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
676       CALL histwrite_p (hist_id, 'Tair', itau_sechiba, ztemp_air, kjpindex, kindex)
677       CALL histwrite_p (hist_id, 'Qair', itau_sechiba, zqair, kjpindex, kindex)
678       CALL histwrite_p (hist2_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
679       CALL histwrite_p (hist2_id, 'SWnet',    itau_sechiba, zswnet, kjpindex, kindex)
680       CALL histwrite_p (hist2_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
681       CALL histwrite_p (hist2_id, 'Qle',  itau_sechiba, zfluxlat, kjpindex, kindex)
682       CALL histwrite_p (hist2_id, 'AvgSurfT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
683       CALL histwrite_p (hist2_id, 'RadT', itau_sechiba, ztemp_sol_NEW, kjpindex, kindex)
684    ENDIF
685    !
686    IF ( is_omp_root ) THEN
687       IF ( (dw .EQ. xrdt) .AND. hist_id > 0 ) THEN
688          ! Syncronize output but only if flag ok_histsync is set to true
689          IF (ok_histsync) CALL histsync(hist_id)
690       ENDIF
691    END IF
692
693    !
694    ! 6. Transform the water fluxes into Kg/m^2s and m^3/s
695    !
696    DO ik=1, kjpindex
697       
698       j = ((kindex(ik)-1)/iim) + 1
699       i = (kindex(ik) - (j-1)*iim)
700       
701       vevapp(i,j) = vevapp(i,j)/xrdt
702       coastalflow(i,j) = coastalflow(i,j)/xrdt
703       riverflow(i,j) = riverflow(i,j)/xrdt
704       
705    ENDDO
706   
707    CALL ipslnlf_p(new_number=old_fileout)
708    IF (printlev_loc >= 3) WRITE (numout,*) 'End intersurf_main_2d'
709
710  END SUBROUTINE intersurf_main_2d
711
712
713!!  =============================================================================================================================
714!! SUBROUTINE:    init_intersurf
715!!
716!>\BRIEF          Initialize grid information
717!!
718!! DESCRIPTION:   This subroutine is called from LMDZ before first call to intersurf_main_gathered or
719!!                intersurf_initialize_gathered
720!!
721!! \n
722!_ ==============================================================================================================================
723
724  SUBROUTINE init_intersurf(nbp_l_lon,nbp_l_lat,kjpindex,kindex,orch_offset,orch_omp_size,orch_omp_rank,COMM, grid)
725
726    USE mod_orchidee_para
727    USE timer
728    IMPLICIT NONE
729
730    INTEGER,INTENT(IN)  :: nbp_l_lon
731    INTEGER,INTENT(IN)  :: nbp_l_lat
732    INTEGER,INTENT(IN)  :: kjpindex
733    INTEGER,INTENT(IN)  :: kindex(:)
734    INTEGER,INTENT(IN)  :: orch_offset
735    INTEGER,INTENT(IN)  :: COMM
736    INTEGER,INTENT(IN)  :: orch_omp_size
737    INTEGER,INTENT(IN)  :: orch_omp_rank
738    INTEGER(i_std), INTENT(in), OPTIONAL :: grid          !! grid type : regular_lonlat or unstructured (dynamico)
739
740    INTEGER,DIMENSION(kjpindex)  :: kindex_offset
741
742
743    IF (printlev >= 1) WRITE(*,*) 'Start ORCHIDEE'
744
745    IF (orch_omp_rank==0) THEN
746      CALL Init_timer
747      CALL start_timer(timer_mpi)
748      CALL grid_set_glo(nbp_l_lon,nbp_l_lat)
749    ENDIF
750    CALL barrier2_omp()   
751    CALL init_orchidee_data_para(kjpindex,kindex,orch_offset,orch_omp_size,orch_omp_rank,COMM)
752
753    ! Initialize specific write level and open text output file
754    printlev_loc=get_printlev('intersurf')
755
756    IF (printlev >= 1) WRITE(numout,*) 'Start ORCHIDEE intitalization phase'
757   
758    IF (PRESENT(grid)) THEN
759       IF (grid==unstructured) THEN
760          CALL grid_allocate_glo(6)
761       ELSE
762          CALL grid_allocate_glo(4)
763       END IF
764    ELSE
765       CALL grid_allocate_glo(4)
766    END IF
767   
768    CALL barrier2_omp()
769    CALL init_ioipsl_para
770         
771    kindex_offset(:)=kindex(:)+offset
772    CALL gather(kindex_offset,index_g)
773    CALL bcast(index_g) 
774
775    IF (printlev >= 2) THEN
776       WRITE(numout,*) "kjpindex = ",kjpindex
777       WRITE(numout,*) "offset for OMP = ",offset_omp
778       WRITE(numout,*) "Index num local for continental points = ",kindex
779       WRITE(numout,*) "Index num global for continental points = ",kindex_offset
780       IF (is_omp_root) THEN
781          WRITE(numout,*) "ROOT OMP, Index global MPI : ",kindex_mpi(:)
782       ENDIF
783       IF (is_root_prc) THEN
784          WRITE(numout,*) "ROOT global, Index global : ",index_g(:)
785       ENDIF
786       WRITE(numout,*) "NbNeighb = ", NbNeighb
787    END IF
788
789
790  END SUBROUTINE init_intersurf
791
792!!  =============================================================================================================================
793!! SUBROUTINE:    intersurf_initialize_gathered
794!!
795!>\BRIEF          Initialization and call to sechiba_initialize
796!!
797!! DESCRIPTION:   Initialization of module variables, read options from parameter file, initialize output files and call to
798!!                sechiba_initialize.
799!!
800!!                This subroutine can be called directly from GCM(LMDZ). If it is not called before the first call to
801!!                intersurf_main_gathered, then it will be done from there. This possibility is done to keep backward
802!!                compatibility with LMDZ.
803!!
804!! \n
805!_ ==============================================================================================================================
806
807  SUBROUTINE intersurf_initialize_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
808       lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
809       zlev,  u, v, qair, temp_air, epot_air, &
810       cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
811       precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
812       vevapp, fluxsens, fluxlat, coastalflow_cpl, riverflow_cpl, &
813       tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, lon_scat_g, lat_scat_g, &
814       q2m, t2m, z0h, nvm_out, &
815       field_out_names, field_in_names, &
816       grid, bounds_latlon, cell_area, ind_cell_glo)
817
818    USE mod_orchidee_para
819    IMPLICIT NONE
820
821    !! 0. Variable and parameter declaration
822    !! 0.1 Input
823    INTEGER(i_std),INTENT (in)                             :: kjit           !! Time step number
824    INTEGER(i_std),INTENT (in)                             :: iim_glo        !! Dimension of global fields
825    INTEGER(i_std),INTENT (in)                             :: jjm_glo        !! Dimension of global fields
826    INTEGER(i_std),INTENT (in)                             :: kjpindex       !! Number of continental points
827    REAL(r_std),INTENT (in)                                :: xrdt           !! Time step in seconds
828    LOGICAL, INTENT (in)                                   :: lrestart_read  !! Logical for _restart_ file to read
829    LOGICAL, INTENT (in)                                   :: lrestart_write !! Logical for _restart_ file to write'
830    REAL(r_std), INTENT (in)                               :: date0          !! Date at which kjit = 0
831    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)       :: kindex         !! Index for continental points
832    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: u              !! Lowest level wind speed
833    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: v              !! Lowest level wind speed
834    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: zlev           !! Height of first layer
835    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: qair           !! Lowest level specific humidity
836    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: precip_rain    !! Rain precipitation
837    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: precip_snow    !! Snow precipitation
838    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: lwdown         !! Down-welling long-wave flux
839    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: swnet          !! Net surface short-wave flux
840    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: swdown         !! Downwelling surface short-wave flux
841    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: temp_air       !! Air temperature in Kelvin
842    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: epot_air       !! Air potential energy
843    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: petAcoef       !! Coeficients A from the PBL resolution
844    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: peqAcoef       !! One for T and another for q
845    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: petBcoef       !! Coeficients B from the PBL resolution
846    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: peqBcoef       !! One for T and another for q
847    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: pb             !! Surface pressure
848    REAL(r_std),DIMENSION (:,:), INTENT(in)                :: latlon         !! Geographical coordinates, this variable has different dimension when using grid unstructured
849    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: zcontfrac      !! Fraction of continent
850    INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: zneighbours    !! neighbours
851    REAL(r_std),DIMENSION (kjpindex,2), INTENT(in)         :: zresolution    !! size of the grid box
852    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN)    :: lon_scat_g     !! Longitudes on the global 2D grid including ocean
853    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN)    :: lat_scat_g     !! Latitudes on the global 2D grid including ocean
854
855    !! 0.2 Output variables
856    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: z0m            !! Surface roughness (momentum)
857    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: coastalflow_cpl!! Diffuse flow of water into the ocean (m^3/s)
858    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: riverflow_cpl  !! Largest rivers flowing into the ocean (m^3/s)
859    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: tsol_rad       !! Radiative surface temperature
860    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: vevapp         !! Total of evaporation
861    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: temp_sol_new   !! New soil temperature
862    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: qsurf          !! Surface specific humidity
863    REAL(r_std),DIMENSION (kjpindex,2), INTENT(out)        :: albedo         !! Albedo
864    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: fluxsens       !! Sensible chaleur flux
865    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: fluxlat        !! Latent chaleur flux
866    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: emis           !! Emissivity
867
868    !! 0.3 Modified variables
869    REAL(r_std),DIMENSION (kjpindex), INTENT(inout)        :: cdrag          !! Cdrag
870
871    !! 0.4 Optional input and output variables
872    REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL :: q2m            !! Surface specific humidity, obsolete: this variable is not used any more, kept only for retrocomp with LMDZ
873    REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL :: t2m            !! Surface air temperature, obsolete: this variable is not used any more, kept only for retrocomp with LMDZ
874    REAL(r_std),DIMENSION (kjpindex), INTENT(out),OPTIONAL :: z0h            !! Surface roughness (heat)
875    INTEGER(i_std), INTENT(out), OPTIONAL                  :: nvm_out        !! Number of vegetation types, PFTs
876    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN)    :: field_out_names!! Name of optional fields sent to LMDZ (fields_out), for ESM configuration
877    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN)    :: field_in_names !! Name of optional fields received from LMDZ (fields_in), for ESM configuration
878    INTEGER(i_std)                  , INTENT(in), OPTIONAL :: grid           !! grid type : regular_lonlat or unstructured (dynamico)
879    REAL(r_std),DIMENSION (:,:,:)   , INTENT(in), OPTIONAL :: bounds_latlon  !! boundaries Geographical coordinates
880    REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL :: cell_area      !! cell area
881    INTEGER(i_std),DIMENSION (:)    , INTENT(in), OPTIONAL :: ind_cell_glo   !! cell order (for XIOS)
882
883
884    !! 0.5 Local variables
885    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
886    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
887    REAL(r_std),DIMENSION (kjpindex)                      :: zz0m          !! Work array to keep zz0m
888    REAL(r_std),DIMENSION (kjpindex)                      :: zz0h          !! Work array to keep zz0h
889    REAL(r_std),DIMENSION (kjpindex)                      :: zcdrag        !! Work array for surface drag
890    REAL(r_std),DIMENSION (kjpindex)                      :: zcoastal      !! Work array to keep coastal flow
891    REAL(r_std),DIMENSION (kjpindex)                      :: zriver        !! Work array to keep river out flow
892    REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
893    REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
894    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
895    REAL(r_std),DIMENSION (kjpindex)                      :: zqsurf        !! Work array to keep qsurf
896    REAL(r_std),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
897    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
898    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
899    REAL(r_std),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis
900    INTEGER(i_std)                                        :: i, j, ik
901    INTEGER(i_std)                                        :: ier
902    INTEGER(i_std)                                        :: itau_sechiba
903    INTEGER(i_std)                                        :: grid_loc      !! Local variable for grid type
904    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)              :: tmp_lon       !! Longitudes for local MPI process. Only available on master OMP.
905    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)              :: tmp_lat       !! Latitudes for local MPI process. Only available on master OMP.
906    INTEGER                                               :: old_fileout   !! old Logical Int for std IO output
907    REAL,ALLOCATABLE,DIMENSION(:,:)                       :: lalo_mpi
908    REAL(r_std),DIMENSION (kjpindex)                      :: landpoints    !! Land point vector
909    REAL(r_std)                                           :: tau_outflow   !! Number of days for temoral filter on river- and coastalflow
910
911 
912    ! Initialize specific write level
913    printlev_loc=get_printlev('intersurf')
914
915    IF (printlev_loc >= 1) WRITE(numout,*) 'Entering intersurf_initialize_gathered'
916    IF (printlev_loc >= 2) WRITE(numout,*) 'using printlev_loc for intersurf:', printlev_loc
917
918    OFF_LINE_MODE = .FALSE. 
919    CALL ipslnlf_p(new_number=numout, old_number=old_fileout)
920
921    ! Initialize variables in module time
922    CALL time_initialize( kjit, date0, xrdt, "END" )
923
924    !  Configuration of SSL specific parameters
925    CALL control_initialize
926
927   
928    ! Check if variable grid is in the argument list
929    IF (PRESENT(grid)) THEN
930      grid_loc=grid
931    ELSE
932      grid_loc=regular_lonlat
933    END IF   
934    IF (printlev >= 2) WRITE(numout,*) 'intersurf_initialize_gathered: grid=', grid_loc
935
936    IF (grid_loc==unstructured) THEN
937      ! First check if all variables needed are present in the call from LMDZ
938      IF (.NOT. PRESENT(bounds_latlon)) CALL ipslerr_p(3,'intersurf_initialize_gathered', &
939        'Error in call intersurf_initialize_gathered from LMDZ','bounds_latlon was not present','')
940
941      IF (.NOT. PRESENT(cell_area)) CALL ipslerr_p(3,'intersurf_initialize_gathered', &
942        'Error in call intersurf_initialize_gathered from LMDZ','cell_area was not present','')
943
944      IF (.NOT. PRESENT(ind_cell_glo)) CALL ipslerr_p(3,'intersurf_initialize_gathered', &
945        'Error in call intersurf_initialize_gathered from LMDZ','ind_cell_glo was not present','')
946
947      ! Allocation of grid variables for unstructured grid
948      IF (printlev_loc >= 3) WRITE(numout,*) 'intersurf_initialize_gathered: call grid_init_unstructured for unstructured grid'
949      CALL grid_init_unstructured(kjpindex,             ij_omp_nb,            latlon(:,2),  latlon(:,1), &
950                                  bounds_latlon(:,:,2), bounds_latlon(:,:,1), cell_area,    ind_cell_glo)
951    ELSE
952      ! Allocation of grid variables for rectilinear lon_lat grid
953      IF (printlev_loc >= 3) WRITE(numout,*) 'intersurf_initialize_gathered: call grid_init for rectilinear lon_lat grid'
954      CALL grid_init ( kjpindex, 4, regular_lonlat, "2DGrid" )
955    ENDIF
956
957
958    !
959    !  Create the internal coordinate table
960    !
961    IF (grid_loc==unstructured) THEN
962      DO ik=1, kjpindex
963        lalo(ik,:)=latlon(kindex(ik),:)
964      ENDDO   
965    ELSE
966      lalo(:,:) = latlon(:,:)
967    ENDIF
968   
969    CALL gather(lalo,lalo_g)
970    CALL bcast(lalo_g)
971    !
972    !-
973    !- Store variable to help describe the grid
974    !- once the points are gathered.
975    !-
976    neighbours(:,:) = zneighbours(:,:)
977    CALL gather(neighbours,neighbours_g)
978    CALL bcast(neighbours_g)
979    !
980    resolution(:,:) = zresolution(:,:)
981    CALL gather(resolution,resolution_g)
982    CALL bcast(resolution_g)
983    !
984    IF (grid_type==regular_lonlat) area(:) = resolution(:,1)*resolution(:,2)
985    CALL gather(area,area_g)
986    CALL bcast(area_g)
987    !
988    !- Store the fraction of the continents only once so that the user
989    !- does not change them afterwards.
990    !
991    contfrac(:) = zcontfrac(:)
992    CALL gather(contfrac,contfrac_g)
993    CALL bcast(contfrac_g)
994    !
995    !
996    !  Create the internal coordinate table
997    !
998    IF ( (.NOT.ALLOCATED(tmp_lon))) THEN
999       ALLOCATE(tmp_lon(iim_g,jj_nb))
1000    ENDIF
1001    IF ( (.NOT.ALLOCATED(tmp_lat))) THEN
1002       ALLOCATE(tmp_lat(iim_g,jj_nb))
1003    ENDIF
1004
1005
1006    IF (is_omp_root) THEN
1007       ! Extract from global variables the longitudes and latitudes
1008       ! for the local MPI process using the indices for the latitude bands(jj_begin, jj_end).
1009       tmp_lon(:,:)=lon_scat_g(:,jj_begin:jj_end)
1010       tmp_lat(:,:)=lat_scat_g(:,jj_begin:jj_end)
1011       
1012       ! Save the global variables only on mpi root, to be used in other modules
1013       IF (is_mpi_root) THEN
1014          lon_g(:,:) = lon_scat_g(:,:)
1015          lat_g(:,:) = lat_scat_g(:,:)
1016       ENDIF
1017    ENDIF
1018   
1019    !Config Key  = FORCE_CO2_VEG
1020    !Config Desc = Flag to force the value of atmospheric CO2 for vegetation.
1021    !Config If   = Only in coupled mode
1022    !Config Def  = FALSE
1023    !Config Help = If this flag is set to true, the ATM_CO2 parameter is used
1024    !Config        to prescribe the atmospheric CO2.
1025    !Config        This Flag is only use in couple mode.
1026    !Config Units = [FLAG]
1027    fatmco2=.FALSE.
1028    CALL getin_p('FORCE_CO2_VEG',fatmco2)
1029    !
1030    ! Next flag is only use in couple mode with a gcm in intersurf.
1031    ! In forced mode, it has already been read and set in driver.
1032    IF ( fatmco2 ) THEN
1033       atmco2=350.
1034       CALL getin_p('ATM_CO2',atmco2)
1035       WRITE(numout,*) 'atmco2 ',atmco2
1036    ENDIF
1037   
1038    CALL ioipslctrl_restini(kjit, date0, xrdt, rest_id, rest_id_stom, itau_offset, date0_shifted)
1039    itau_sechiba = kjit + itau_offset
1040   
1041    !!- Initialize module for output with XIOS
1042    !
1043    CALL xios_orchidee_init( MPI_COMM_ORCH,                   &
1044         date0,    year_end, month_end,     day_end, julian_diff, &
1045         tmp_lon,  tmp_lat )
1046
1047    IF (printlev >= 3) WRITE(numout,*) 'intersurf_initialize_gathered: After xios_orchidee_init'
1048    CALL sechiba_xios_initialize
1049    IF (printlev >= 3) WRITE(numout,*) 'intersurf_initialize_gathered: After sechiba_xios_initialize'
1050
1051    CALL xios_orchidee_close_definition
1052    IF (printlev >= 2) WRITE(numout,*) 'intersurf_initialize_gathered: After xios_orchidee_close_definition'
1053
1054    !- Initialize IOIPSL sechiba output files
1055    CALL ioipslctrl_history(iim_g, jj_nb, tmp_lon, tmp_lat,  kindex, kjpindex, itau_sechiba, &
1056         date0_shifted, xrdt, hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC)
1057   
1058    CALL bcast_omp(hist_id)
1059    CALL bcast_omp(hist2_id)
1060    CALL bcast_omp(hist_id_stom)
1061    CALL bcast_omp(hist_id_stom_IPCC)
1062
1063    !
1064    !! Change to be in the orchidee context for XIOS
1065    !
1066    CALL xios_orchidee_change_context("orchidee")
1067
1068    !
1069    !  Shift the time step to phase the two models
1070    !
1071    itau_sechiba = kjit + itau_offset
1072   
1073    ! Update the calendar in xios by sending the new time step
1074    ! Special case : the model is only in initialization phase and the current itau_sechiba is not a real time step.
1075    ! Therefor give itau_sechiba+1 to xios to have a correct time axis in output files.
1076    CALL xios_orchidee_update_calendar(itau_sechiba+1)
1077   
1078    !
1079    ! 1. Just change the units of some input fields
1080    !
1081    DO ik=1, kjpindex
1082       
1083       zprecip_rain(ik) = precip_rain(ik)*xrdt
1084       zprecip_snow(ik) = precip_snow(ik)*xrdt
1085       zcdrag(ik)       = cdrag(ik)
1086    ENDDO
1087
1088     
1089    !
1090    ! 3. save the grid
1091    !
1092    landpoints(:)=(/ ( REAL(ik), ik=1,kjpindex ) /)
1093    CALL histwrite_p(hist_id, 'LandPoints',  itau_sechiba+1, landpoints, kjpindex, kindex)
1094    CALL histwrite_p(hist_id, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
1095    IF ( ok_stomate ) THEN
1096       CALL histwrite_p(hist_id_stom, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
1097       IF ( hist_id_stom_ipcc > 0 ) &
1098            CALL histwrite_p(hist_id_stom_IPCC, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
1099    ENDIF
1100    CALL histwrite_p(hist_id, 'Contfrac',  itau_sechiba+1, contfrac, kjpindex, kindex)
1101   
1102    ! Syncronize output but only if flag ok_histsync is set to true       
1103    IF (ok_histsync) THEN
1104       IF (is_omp_root .AND. hist_id > 0) THEN
1105          CALL histsync(hist_id)
1106       END IF
1107    END IF
1108   
1109    IF ( hist2_id > 0 ) THEN
1110       CALL histwrite_p(hist2_id, 'LandPoints',  itau_sechiba+1, landpoints, kjpindex, kindex)
1111       CALL histwrite_p(hist2_id, 'Areas',  itau_sechiba+1, area, kjpindex, kindex)
1112       CALL histwrite_p(hist2_id, 'Contfrac',  itau_sechiba+1, contfrac, kjpindex, kindex)
1113       
1114       ! Syncronize output but only if flag ok_histsync is set to true
1115       IF (ok_histsync .AND. is_omp_root) THEN
1116          CALL histsync(hist2_id)
1117       ENDIF
1118    ENDIF
1119    !
1120   
1121    !
1122    ! 4. call sechiba for continental points only
1123    !
1124    IF ( printlev_loc>=3 ) WRITE(numout,*) 'Before call to sechiba_initialize'
1125
1126    CALL sechiba_initialize( &
1127         itau_sechiba, iim_g*jj_nb,  kjpindex,      kindex,                     &
1128         lalo,         contfrac,     neighbours,    resolution,  zlev,          &
1129         u,            v,            qair,          temp_air,                   &
1130         petAcoef,     peqAcoef,     petBcoef,      peqBcoef,                   &
1131         zprecip_rain, zprecip_snow, lwdown,        swnet,       swdown,        &
1132         pb,           rest_id,      hist_id,       hist2_id,                   &
1133         rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                         &
1134         zcoastal,     zriver,       ztsol_rad,     zvevapp,     zqsurf,        &
1135         zz0m,          zz0h,        zalbedo,      zfluxsens,     zfluxlat,    zemis,         &
1136         ztemp_sol_new, zcdrag)
1137
1138   
1139    IF ( printlev_loc>=3 ) WRITE(numout,*) 'After call to sechiba_initialize'
1140
1141    !
1142    ! 6. scatter output fields
1143    !
1144    DO ik=1, kjpindex
1145       z0m(ik)          = zz0m(ik)
1146       tsol_rad(ik)     = ztsol_rad(ik)
1147       vevapp(ik)       = zvevapp(ik)/xrdt ! Transform into kg/m^2/s
1148       temp_sol_new(ik) = ztemp_sol_new(ik)
1149       qsurf(ik)        = zqsurf(ik)
1150       albedo(ik,1)     = zalbedo(ik,1)
1151       albedo(ik,2)     = zalbedo(ik,2)
1152       fluxsens(ik)     = zfluxsens(ik)
1153       fluxlat(ik)      = zfluxlat(ik)
1154       emis(ik)         = zemis(ik)
1155       cdrag(ik)        = zcdrag(ik)
1156    ENDDO
1157
1158    ! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1159    IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1160
1161
1162    !
1163    ! 7. Count optional fields for the coupling with LMDZ, used in the ESM configuration
1164    !
1165    ! Treat fields sent to LMDZ
1166    IF (PRESENT(field_out_names)) THEN
1167       nbcf_from_LMDZ=SIZE(field_out_names)
1168       ALLOCATE(field_out_names_loc(nbcf_from_LMDZ), stat=ier)
1169       IF (ier /= 0) CALL ipslerr_p(3,'intersurf_initialized_gathered', 'Error in allocation of field_out_names_loc','','')
1170       field_out_names_loc=field_out_names
1171    ELSE
1172       nbcf_from_LMDZ=0
1173    ENDIF
1174
1175    ! Treat fields received from LMDZ
1176    IF (PRESENT(field_in_names)) THEN
1177       nbcf_into_LMDZ=SIZE(field_in_names)
1178       ALLOCATE(field_in_names_loc(nbcf_into_LMDZ), stat=ier)
1179       IF (ier /= 0) CALL ipslerr_p(3,'intersurf_initialized_gathered', 'Error in allocation of field_in_names_loc','','')
1180       field_in_names_loc=field_in_names
1181
1182       ! Coherence test: veget_update>0 needed for fCO2_fLuc and fCO2_nbp
1183       DO i = 1, nbcf_into_LMDZ
1184          IF ((field_in_names_loc(i) == 'fCO2_fLuc') .OR. (field_in_names_loc(i) == 'fCO2_nbp')) THEN
1185             IF (veget_update .LE. 0) THEN
1186                CALL ipslerr_p(3,'intersurf_initialize_gathered', &
1187                     'Error in set up. VEGET_UPDATE must be > 0 when using the field: '//TRIM(field_in_names_loc(i)),'','')
1188             END IF
1189          END IF
1190       END DO
1191    ELSE
1192       nbcf_into_LMDZ=0
1193    ENDIF
1194    IF ( printlev_loc>=2 ) THEN
1195       WRITE(numout,*) 'intersurf_initialized_gathered --- nbcf_from_LMDZ, nbcf_into_LMDZ = ',nbcf_from_LMDZ, nbcf_into_LMDZ
1196       IF (nbcf_into_LMDZ > 0) THEN
1197          WRITE(numout,*) 'intersurf_initialized_gathered --- field_in_names_loc = ', field_in_names_loc
1198       END IF
1199       IF (nbcf_from_LMDZ > 0) THEN
1200          WRITE(numout,*) 'intersurf_initialized_gathered --- field_out_names_loc = ', field_out_names_loc
1201       END IF
1202    END IF
1203
1204    !
1205    ! 8. Distribution of coastal- and riverflow on several time steps
1206    !
1207    !Config Key  = TAU_OUTFLOW
1208    !Config Desc = Number of days over which the coastal- and riverflow will be distributed
1209    !Config If   = Only in coupled mode
1210    !Config Def  = 0
1211    !Config Help = The default value 0 makes the distribution instanteneous
1212    !Config Units = [days]
1213    tau_outflow = 0
1214    CALL getin_p('TAU_OUTFLOW',tau_outflow)
1215    IF (tau_outflow <=xrdt/one_day) THEN
1216       coeff_rel = 1.0
1217    ELSE
1218       coeff_rel = (1.0 - exp(-xrdt/(tau_outflow*one_day)))
1219    END IF
1220    IF (printlev_loc >=2)  WRITE(numout,*),'tau_outflow, coeff_rel = ', tau_outflow, coeff_rel
1221
1222    ! Allocate and read riverflow_cpl0 from restart file. Initialize to 0 if the variable was not found in the restart file.
1223    ALLOCATE(riverflow_cpl0(kjpindex), stat=ier) 
1224    IF (ier /= 0) CALL ipslerr_p(3,'intersurf_initialize_gathered', 'Error in allocation of riverflow_cpl0','','')
1225    CALL restget_p (rest_id, 'riverflow_cpl0', nbp_glo, 1, 1, kjit, .TRUE., riverflow_cpl0, "gather", nbp_glo, index_g)
1226    IF ( ALL(riverflow_cpl0(:) == val_exp) ) riverflow_cpl0(:)=0
1227
1228    ! Allocate and read coastalflow_cpl0 from restart file. Initialize to 0 if the variable was not found in the restart file.
1229    ALLOCATE(coastalflow_cpl0(kjpindex), stat=ier) 
1230    IF (ier /= 0) CALL ipslerr_p(3,'intersurf_initialize_gathered', 'Error in allocation of coastalflow_cpl0','','')
1231    CALL restget_p (rest_id, 'coastalflow_cpl0', nbp_glo, 1, 1, kjit, .TRUE., coastalflow_cpl0, "gather", nbp_glo, index_g)
1232    IF ( ALL(coastalflow_cpl0(:) == val_exp) ) coastalflow_cpl0(:)=0
1233
1234    ! Do not applay the filter now in initialization phase.
1235    ! These variables will never be used anyway in the initialization phase.
1236    ! Transform into m^3/s
1237    riverflow_cpl = zriver/xrdt
1238    coastalflow_cpl = zcoastal/xrdt
1239   
1240    !
1241    ! 9.
1242    !
1243    ! Copy the number of vegetation types to local output variable
1244    IF (PRESENT(nvm_out)) nvm_out=nvm
1245
1246    IF(is_root_prc) CALL getin_dump
1247    lstep_init_intersurf = .FALSE.
1248   
1249    CALL ipslnlf_p(new_number=old_fileout)
1250    !
1251    !! Change back to be in the LMDZ context for XIOS
1252    !
1253    CALL xios_orchidee_change_context("LMDZ")
1254
1255    IF (printlev_loc >= 2) WRITE (numout,*) 'End intersurf_initialize_gathered'
1256    IF (printlev_loc >= 1) WRITE (numout,*) 'Initialization phase for ORCHIDEE is finished.'
1257
1258  END SUBROUTINE intersurf_initialize_gathered
1259
1260
1261!!  =============================================================================================================================
1262!! SUBROUTINE:    intersurf_main_gathered
1263!!
1264!>\BRIEF          Main subroutine to call ORCHIDEE from the gcm (LMDZ) using variables on a 1D grid with only land points.
1265!!
1266!! DESCRIPTION:   This subroutine is the main interface for ORCHIDEE when it is called from the gcm (LMDZ).
1267!!                The variables are all gathered before entering this subroutine on the 1D grid with only landpoints.
1268!!
1269!! \n
1270!_ ==============================================================================================================================
1271
1272  SUBROUTINE intersurf_main_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
1273       lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
1274       zlev,  u, v, qair, temp_air, epot_air, ccanopy, &
1275       cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1276       precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
1277       vevapp, fluxsens, fluxlat, coastalflow_cpl, riverflow_cpl, &
1278       tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m,lon_scat_g, lat_scat_g, q2m, t2m, z0h, &
1279       veget, lai, height, &
1280       fields_out, fields_in,  &
1281       coszang)
1282
1283    USE mod_orchidee_para
1284    IMPLICIT NONE
1285
1286    !! 0. Variable and parameter declaration
1287    !! 0.1 Input variables
1288    INTEGER(i_std),INTENT (in)                            :: kjit            !! Time step number
1289    INTEGER(i_std),INTENT (in)                            :: iim_glo         !! Dimension of global fields
1290    INTEGER(i_std),INTENT (in)                            :: jjm_glo         !! Dimension of global fields
1291    INTEGER(i_std),INTENT (in)                            :: kjpindex        !! Number of continental points
1292    REAL(r_std),INTENT (in)                               :: xrdt            !! Time step in seconds
1293    LOGICAL, INTENT (in)                                  :: lrestart_read   !! Logical for _restart_ file to read
1294    LOGICAL, INTENT (in)                                  :: lrestart_write  !! Logical for _restart_ file to write'
1295    REAL(r_std), INTENT (in)                              :: date0           !! Date at which kjit = 0
1296    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)      :: kindex          !! Index for continental points
1297    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: u               !! Lowest level wind speed
1298    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: v               !! Lowest level wind speed
1299    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: zlev            !! Height of first layer
1300    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: qair            !! Lowest level specific humidity
1301    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: precip_rain     !! Rain precipitation
1302    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: precip_snow     !! Snow precipitation
1303    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: lwdown          !! Down-welling long-wave flux
1304    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: swnet           !! Net surface short-wave flux
1305    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: swdown          !! Downwelling surface short-wave flux
1306    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: temp_air        !! Air temperature in Kelvin
1307    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: epot_air        !! Air potential energy
1308    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: ccanopy         !! CO2 concentration in the canopy
1309    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: petAcoef        !! Coeficients A from the PBL resolution
1310    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: peqAcoef        !! One for T and another for q
1311    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: petBcoef        !! Coeficients B from the PBL resolution
1312    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: peqBcoef        !! One for T and another for q
1313    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: pb              !! Surface pressure
1314    REAL(r_std),DIMENSION (:,:)     , INTENT(in)          :: latlon          !! Geographical coordinates
1315    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: zcontfrac       !! Fraction of continent
1316    INTEGER(i_std),DIMENSION (:,:), INTENT(in)            :: zneighbours     !! neighbours
1317    REAL(r_std),DIMENSION (kjpindex,2), INTENT(in)        :: zresolution     !! size of the grid box
1318    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN)   :: lon_scat_g      !! The scattered values for longitude
1319    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN)   :: lat_scat_g      !! The scattered values for latitude
1320
1321    !! 0.2 Output variables
1322    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: z0m             !! Surface roughness for momentum
1323    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: coastalflow_cpl !! Diffuse flow of water into the ocean, time filtered (m^3/s)
1324    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: riverflow_cpl   !! Largest rivers flowing into the ocean, time filtered (m^3/s)
1325    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: tsol_rad        !! Radiative surface temperature
1326    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: vevapp          !! Total of evaporation
1327    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: temp_sol_new    !! New soil temperature
1328    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: qsurf           !! Surface specific humidity
1329    REAL(r_std),DIMENSION (kjpindex,2), INTENT(out)       :: albedo          !! Albedo
1330    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: fluxsens        !! Sensible chaleur flux
1331    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: fluxlat         !! Latent chaleur flux
1332    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: emis            !! Emissivity
1333
1334    !! 0.3 Modified variables
1335    REAL(r_std),DIMENSION (kjpindex), INTENT(inout)       :: cdrag           !! Cdrag
1336
1337    !! 0.4 Optional input and output variables
1338    REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL:: q2m             !! Surface specific humidity, obsolete: this variable is not used any more, kept only for retrocomp with LMDZ
1339    REAL(r_std),DIMENSION (kjpindex), INTENT(in), OPTIONAL:: t2m             !! Surface air temperature, obsolete: this variable is not used any more, kept only for retrocomp with LMDZ
1340    REAL(r_std),DIMENSION (kjpindex), INTENT(out),OPTIONAL:: z0h             !! Surface roughness for heat
1341    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out), OPTIONAL :: veget     !! Fraction of vegetation type (unitless, 0-1)
1342    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out), OPTIONAL :: lai       !! Leaf area index (m^2 m^{-2}
1343    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out), OPTIONAL :: height    !! Vegetation Height (m)
1344    REAL(r_std), DIMENSION(kjpindex), OPTIONAL, INTENT(in):: coszang         !! Cosine of the solar zenith angle (unitless)
1345    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(OUT)     :: fields_in       !! Optional fields sent to LMDZ, for ESM configuration
1346    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN)      :: fields_out      !! Optional fields received from LMDZ, for ESM configuration
1347
1348
1349    !! 0.5 Local variables
1350    REAL(r_std),DIMENSION (kjpindex)                      :: zccanopy      !! Work array to keep ccanopy
1351    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
1352    REAL(r_std),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
1353    REAL(r_std),DIMENSION (kjpindex)                      :: zz0m, zz0h    !! Work array to keep zz0m, zz0h
1354    REAL(r_std),DIMENSION (kjpindex)                      :: zcdrag        !! Work array for surface drag
1355    REAL(r_std),DIMENSION (kjpindex)                      :: zcoastal      !! Work array to keep coastal flow
1356    REAL(r_std),DIMENSION (kjpindex)                      :: zriver        !! Work array to keep river out flow
1357    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux (gC/m2/one_day)
1358    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use (gC/m2/one_day)
1359    REAL(r_std),DIMENSION (kjpindex)                      :: zcarbwh       !! Work array to keep fco2_wh (gC/m2/one_day)
1360    REAL(r_std),DIMENSION (kjpindex)                      :: zcarbha       !! Work array to keep fco2_ha (gC/m2/one_day)
1361    REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
1362    REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
1363    REAL(r_std),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
1364    REAL(r_std),DIMENSION (kjpindex)                      :: zqsurf        !! Work array to keep qsurf
1365    REAL(r_std),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
1366    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
1367    REAL(r_std),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
1368    REAL(r_std),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis
1369    REAL(r_std),DIMENSION (kjpindex)                      :: zcoszang      !! Work array to keep coszang
1370    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zveget        !! Work array to keep veget
1371    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zlai          !! Work array to keep lai
1372    REAL(r_std),DIMENSION (kjpindex,nvm)                  :: zheight       !! Work array to keep height
1373    INTEGER(i_std)                                        :: i, j, ik
1374    INTEGER(i_std)                                        :: ier
1375    INTEGER(i_std)                                        :: itau_sechiba
1376    INTEGER                                               :: old_fileout   !! old Logical Int for std IO output
1377    REAL,ALLOCATABLE,DIMENSION(:,:)                       :: lalo_mpi
1378    REAL(r_std),DIMENSION (kjpindex)                      :: landpoints    !! Local landpoints vector
1379
1380    IF (printlev_loc >= 3) WRITE(numout,*) 'Start intersurf_main_gathered'
1381    CALL ipslnlf_p(new_number=numout, old_number=old_fileout)
1382   
1383    IF (lstep_init_intersurf) THEN
1384       ! intersurf_initialize_gathered was not yet called. An old version of
1385       ! LMDZ is used. Now prepare and call intersurf_initialize_gathered
1386
1387       CALL intersurf_initialize_gathered (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
1388            lrestart_read, lrestart_write, latlon, zcontfrac, zneighbours, zresolution, date0, &
1389            zlev,  u, v, qair, temp_air, epot_air, &
1390            cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1391            precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
1392            vevapp, fluxsens, fluxlat, coastalflow_cpl, riverflow_cpl, &
1393            tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m,lon_scat_g, lat_scat_g, &
1394            zz0h)
1395
1396       ! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1397       IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1398       ! Return from subroutine intersurf_main_gathered
1399       RETURN
1400    END IF
1401
1402    !
1403    !! Change to be in the orchidee context for XIOS
1404    !
1405    CALL xios_orchidee_change_context("orchidee")
1406   
1407    !
1408    !  Shift the time step to phase the two models
1409    !
1410    itau_sechiba = kjit + itau_offset
1411   
1412    ! Update the calendar in xios by sending the new time step
1413    CALL xios_orchidee_update_calendar(itau_sechiba)
1414   
1415    ! Update the calendar and all time variables in module time
1416    CALL time_nextstep(itau_sechiba)
1417   
1418    !
1419    ! 1. Just change the units of some input fields
1420    !
1421    DO ik=1, kjpindex
1422       
1423       zprecip_rain(ik) = precip_rain(ik)*xrdt
1424       zprecip_snow(ik) = precip_snow(ik)*xrdt
1425       zcdrag(ik)       = cdrag(ik)
1426       
1427    ENDDO
1428
1429    !>> VOC in coupled mode
1430    IF ( PRESENT(coszang) )  THEN
1431       zcoszang(:) = coszang(:)
1432    ELSE
1433       zcoszang(:) = zero
1434    ENDIF
1435 
1436    ! Optional fields received from LMDZ, for ESM configuration
1437    ! No fields are currently received from LMDZ
1438    DO i = 1, nbcf_from_LMDZ 
1439       SELECT CASE(TRIM(field_out_names_loc(i)))
1440       ! CASE("...")
1441          ! Things to do for this case...
1442        CASE("atmco2")
1443          IF (printlev_loc > 4) WRITE(numout,*) &
1444               'Retreive the atmco2 field, in this case equal to fields_out(:,i)=ccanopy(:)'
1445       CASE DEFAULT 
1446          CALL ipslerr_p (3,'intersurf_main_gathered', &
1447               'The selected field received from LMDZ '//TRIM(field_out_names_loc(i))//' is not implemented.','','')
1448       END SELECT
1449    ENDDO
1450
1451    !
1452    ! 2. modification of co2
1453    !
1454    IF ( fatmco2 ) THEN
1455       zccanopy(:) = atmco2
1456       WRITE (numout,*) 'Modification of the ccanopy value. CO2 = ',atmco2
1457    ELSE
1458       zccanopy(:) = ccanopy(:)
1459    ENDIF
1460
1461    !
1462    ! 4. call sechiba for continental points only
1463    !
1464    IF ( printlev_loc >= 3 ) WRITE(numout,*) 'Before call to sechiba_main from intersurf_main_gathered'
1465   
1466
1467    CALL sechiba_main (itau_sechiba, iim_g*jj_nb, kjpindex, kindex, &
1468         lrestart_read, lrestart_write, &
1469         lalo, contfrac, neighbours, resolution, &
1470         zlev, u, v, qair, temp_air, epot_air, zccanopy, &
1471         zcdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1472         zprecip_rain ,zprecip_snow,  lwdown, swnet, swdown, zcoszang, pb, &
1473         zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
1474         znetco2, zcarblu, zcarbwh, zcarbha, &
1475         ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0m, zz0h, &
1476         zveget, zlai, zheight, &
1477         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC ) 
1478
1479    IF ( printlev_loc>=3 ) WRITE(numout,*) 'After call to sechiba_main'
1480
1481    !
1482    ! 6. scatter output fields
1483    !
1484    DO ik=1, kjpindex
1485       z0m(ik)          = zz0m(ik)
1486       tsol_rad(ik)     = ztsol_rad(ik)
1487       temp_sol_new(ik) = ztemp_sol_new(ik)
1488       qsurf(ik)        = zqsurf(ik)
1489       albedo(ik,1)     = zalbedo(ik,1)
1490       albedo(ik,2)     = zalbedo(ik,2)
1491       fluxsens(ik)     = zfluxsens(ik)
1492       fluxlat(ik)      = zfluxlat(ik)
1493       emis(ik)         = zemis(ik)
1494       cdrag(ik)        = zcdrag(ik)
1495       ! Transform the water fluxes into Kg/m^2/s
1496       vevapp(ik)       = zvevapp(ik)/xrdt
1497    ENDDO
1498
1499    ! Copy variables only if the optional variables are present in the call to the subroutine
1500    IF (PRESENT(veget))  veget(:,:)  = zveget(:,:) 
1501    IF (PRESENT(lai))    lai(:,:)    = zlai(:,:) 
1502    IF (PRESENT(height)) height(:,:) = zheight(:,:) 
1503
1504
1505    !
1506    ! 7. Set variables to be sent to LMDZ, for ESM configuration
1507    !    field_in_names_loc should be consitent with the names used in coupling_fields.def read by LMDZ
1508    !
1509    IF (nbcf_into_LMDZ.GT.0) THEN
1510       DO ik=1, kjpindex
1511          DO i = 1, nbcf_into_LMDZ
1512             SELECT CASE(TRIM(field_in_names_loc(i)))
1513             CASE("fCO2_nep")
1514                ! Transform from gC/m2/one_day into kgC/m2/s
1515                fields_in(ik,i) = znetco2(ik)/1.e3/one_day
1516             CASE("fCO2_fLuc")
1517                ! Transform from gC/m2/one_day into kgC/m2/s
1518                fields_in(ik,i) = zcarblu(ik)/1.e3/one_day 
1519             CASE("fCO2_fWoodharvest")
1520                ! Transform from gC/m2/one_day into kgC/m2/s
1521                fields_in(ik,i) = zcarbwh(ik)/1.e3/one_day
1522             CASE("fCO2_fHarvest")
1523                ! Transform from gC/m2/one_day into kgC/m2/s
1524                fields_in(ik,i) = zcarbha(ik)/1.e3/one_day
1525             CASE("fCO2_nbp")
1526                fields_in(ik,i) = znetco2(ik)/1.e3/one_day + zcarblu(ik)/1.e3/one_day + &
1527                       zcarbwh(ik)/1.e3/one_day + zcarbha(ik)/1.e3/one_day
1528             CASE DEFAULT
1529                CALL ipslerr_p (3,'intersurf_main_gathered', &
1530                     'The selected field to be sent to LMDZ '//TRIM(field_in_names_loc(i))//' is not yet implemented.','','')
1531             END SELECT
1532          END  DO
1533       ENDDO
1534    END IF
1535   
1536
1537    ! Applay time filter to distribut the river- and coastalflow over a longer time period.
1538    ! When coeff_rel=1(default case when tau_outflow=0), the distribution is instanteneous.
1539    ! Use TAU_OUTFLOW in run.def to set the number of days of distribution.
1540    ! The water fluxes zriver and zcoastal coming from sechiba are at the same time transfromed
1541    ! from m^3/dt_sechiba into m^3/s by dividing with the sechiba time step (xrdt).
1542    riverflow_cpl = coeff_rel*zriver/xrdt + (1.-coeff_rel)*riverflow_cpl0
1543    riverflow_cpl0 = riverflow_cpl
1544
1545    coastalflow_cpl = coeff_rel*zcoastal/xrdt + (1.-coeff_rel)*coastalflow_cpl0
1546    coastalflow_cpl0 = coastalflow_cpl 
1547
1548   
1549    ! z0h is a optional output variable. Check first if it is present in the call from LMDZ.
1550    IF ( PRESENT(z0h) ) z0h(:) = zz0h(:)
1551       
1552    CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /))
1553    CALL xios_orchidee_send_field("areas", area)
1554    CALL xios_orchidee_send_field("contfrac",contfrac)
1555    CALL xios_orchidee_send_field("temp_air",temp_air)
1556    CALL xios_orchidee_send_field("qair",qair)
1557    CALL xios_orchidee_send_field("swnet",swnet)
1558    CALL xios_orchidee_send_field("swdown",swdown)
1559    CALL xios_orchidee_send_field("pb",pb)
1560    CALL xios_orchidee_send_field("riverflow_cpl",riverflow_cpl)
1561    CALL xios_orchidee_send_field("coastalflow_cpl",coastalflow_cpl)
1562    CALL xios_orchidee_send_field("znetco2",znetco2/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
1563    CALL xios_orchidee_send_field("zcarblu",zcarblu/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
1564    CALL xios_orchidee_send_field("zcarbwh",zcarbwh/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
1565    CALL xios_orchidee_send_field("zcarbha",zcarbha/1.e3/one_day)    ! gC/m2/one_day -> kgC/m2/s
1566
1567   
1568    IF ( .NOT. almaoutput ) THEN
1569       !
1570       !  scattered during the writing
1571       !           
1572       CALL histwrite_p (hist_id, 'evap',     itau_sechiba, zvevapp, kjpindex, kindex)
1573       CALL histwrite_p (hist_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
1574       CALL histwrite_p (hist_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
1575       CALL histwrite_p (hist_id, 'temp_sol', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1576       CALL histwrite_p (hist_id, 'tsol_max', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1577       CALL histwrite_p (hist_id, 'tsol_min', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1578       CALL histwrite_p (hist_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
1579       CALL histwrite_p (hist_id, 'fluxlat',  itau_sechiba, zfluxlat,  kjpindex, kindex)
1580       CALL histwrite_p (hist_id, 'swnet',    itau_sechiba, swnet,    kjpindex, kindex)
1581       CALL histwrite_p (hist_id, 'swdown',   itau_sechiba, swdown,   kjpindex, kindex)
1582       CALL histwrite_p (hist_id, 'alb_vis',  itau_sechiba, zalbedo(:,1), kjpindex, kindex)
1583       CALL histwrite_p (hist_id, 'alb_nir',  itau_sechiba, zalbedo(:,2), kjpindex, kindex)
1584       CALL histwrite_p (hist_id, 'tair',     itau_sechiba, temp_air, kjpindex, kindex)
1585       CALL histwrite_p (hist_id, 'qair',     itau_sechiba, qair, kjpindex, kindex)
1586       
1587       IF ( hist2_id > 0 ) THEN
1588          CALL histwrite_p (hist2_id, 'evap',     itau_sechiba, zvevapp, kjpindex, kindex)
1589          CALL histwrite_p (hist2_id, 'coastalflow',itau_sechiba, zcoastal, kjpindex, kindex)
1590          CALL histwrite_p (hist2_id, 'riverflow',itau_sechiba, zriver, kjpindex, kindex)
1591          CALL histwrite_p (hist2_id, 'temp_sol', itau_sechiba, temp_sol_new, kjpindex, kindex)
1592          CALL histwrite_p (hist2_id, 'tsol_max', itau_sechiba, temp_sol_new, kjpindex, kindex)
1593          CALL histwrite_p (hist2_id, 'tsol_min', itau_sechiba, temp_sol_new, kjpindex, kindex)
1594          CALL histwrite_p (hist2_id, 'fluxsens', itau_sechiba, zfluxsens, kjpindex, kindex)
1595          CALL histwrite_p (hist2_id, 'fluxlat',  itau_sechiba, zfluxlat,  kjpindex, kindex)
1596          CALL histwrite_p (hist2_id, 'swnet',    itau_sechiba, swnet,    kjpindex, kindex)
1597          CALL histwrite_p (hist2_id, 'swdown',   itau_sechiba, swdown,   kjpindex, kindex)
1598          CALL histwrite_p (hist2_id, 'alb_vis',  itau_sechiba, zalbedo(:,1), kjpindex, kindex)
1599          CALL histwrite_p (hist2_id, 'alb_nir',  itau_sechiba, zalbedo(:,2), kjpindex, kindex)
1600          CALL histwrite_p (hist2_id, 'tair',     itau_sechiba, temp_air, kjpindex, kindex)
1601          CALL histwrite_p (hist2_id, 'qair',     itau_sechiba, qair, kjpindex, kindex)
1602       ENDIF
1603    ELSE
1604       CALL histwrite_p (hist_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
1605       CALL histwrite_p (hist_id, 'SWnet',    itau_sechiba, swnet, kjpindex, kindex)
1606       CALL histwrite_p (hist_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
1607       CALL histwrite_p (hist_id, 'Qle',  itau_sechiba, zfluxlat, kjpindex, kindex)
1608       CALL histwrite_p (hist_id, 'AvgSurfT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1609       CALL histwrite_p (hist_id, 'RadT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1610       IF ( hist2_id > 0 ) THEN
1611          CALL histwrite_p (hist2_id, 'Evap', itau_sechiba, zvevapp, kjpindex, kindex)
1612          CALL histwrite_p (hist2_id, 'SWnet',    itau_sechiba, swnet, kjpindex, kindex)
1613          CALL histwrite_p (hist2_id, 'Qh', itau_sechiba, zfluxsens, kjpindex, kindex)
1614          CALL histwrite_p (hist2_id, 'Qle',  itau_sechiba, zfluxlat, kjpindex, kindex)
1615          CALL histwrite_p (hist2_id, 'AvgSurfT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1616          CALL histwrite_p (hist2_id, 'RadT', itau_sechiba, ztemp_sol_new, kjpindex, kindex)
1617       ENDIF
1618    ENDIF
1619   
1620    ! Syncronize output but only if flag ok_histsync is set to true
1621    IF (ok_histsync .AND. is_omp_root) THEN
1622       IF ( (dw .EQ. xrdt) .AND. hist_id > 0 ) THEN
1623          CALL histsync(hist_id)
1624       ENDIF
1625    ENDIF
1626   
1627 
1628
1629    ! Write variables to restart file
1630    IF (lrestart_write) THEN
1631       CALL restput_p (rest_id, 'riverflow_cpl0', nbp_glo, 1, 1, kjit, riverflow_cpl0, 'scatter',  nbp_glo, index_g)
1632       CALL restput_p (rest_id, 'coastalflow_cpl0', nbp_glo, 1, 1, kjit, coastalflow_cpl0, 'scatter',  nbp_glo, index_g)
1633    END IF
1634
1635    !
1636    CALL ipslnlf_p(new_number=old_fileout)
1637    !       
1638
1639    !
1640    !! Finalize the XIOS orchidee context if it is the last call
1641    !
1642    IF (lrestart_write) THEN
1643       CALL ioipslctrl_restclo
1644       CALL xios_orchidee_context_finalize
1645    END IF
1646    !
1647    !! Change back to be in the LMDZ context for XIOS
1648    !
1649    CALL xios_orchidee_change_context("LMDZ")
1650
1651    IF (printlev_loc >= 3) WRITE (numout,*) 'End intersurf_main_gathered'
1652
1653  END SUBROUTINE intersurf_main_gathered
1654
1655!! ================================================================================================================================
1656!! SUBROUTINE   : intersurf_clear
1657!!
1658!>\BRIEF         Clear intersurf module and underlaying modules
1659!!
1660!! DESCRIPTION  :  Deallocate memory and reset initialization variables to there original values.
1661!!                 Call the clearing for sechiba module.
1662!!
1663!_ ================================================================================================================================
1664  SUBROUTINE intersurf_clear
1665    CALL sechiba_clear
1666  END SUBROUTINE intersurf_clear
1667
1668END MODULE intersurf
Note: See TracBrowser for help on using the repository browser.