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

Last change on this file since 6019 was 5575, checked in by josefine.ghattas, 6 years ago

Correction needed to run on unstructured grid

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