source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/intersurf.f90 @ 7541

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