source: branches/ORCHIDEE_2_2/ORCHIDEE/src_oasisdriver/driver2oasis.f90

Last change on this file was 8377, checked in by jan.polcher, 6 months ago

The modifications performed on the trunk for the coupling to WRF and the interpolation by XIOS on curviliean grids has been backported.
File size: 35.9 KB
Line 
1PROGRAM driver2oasis
2  !---------------------------------------------------------------------
3  !-
4  !- Reads the forcing file in the ALMA format and feeds the data to the
5  !- OASIS coupler.
6  !-
7  !---------------------------------------------------------------------
8  USE defprec
9  !
10  USE netcdf
11  !
12  USE constantes
13  !
14  USE ioipsl_para
15  USE mod_orchidee_para
16  !
17  USE grid
18  USE time
19  USE timer
20  USE constantes
21  USE constantes_soil
22  USE forcing_tools
23  USE globgrd
24  !
25  USE mod_oasis
26  USE timer
27  USE xios
28  !-
29  IMPLICIT NONE
30  !-
31  INCLUDE 'mpif.h'
32  !-
33  CHARACTER(LEN=80) :: gridfilename
34  CHARACTER(LEN=80), DIMENSION(100) :: forfilename
35  CHARACTER(LEN=6)   :: comp_name = 'driver'
36  !
37  !
38  CHARACTER(LEN=8)  :: model_guess
39  INTEGER(i_std)    :: iim_glo, jjm_glo, file_id
40  !-
41  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon_glo, lat_glo, area_glo
42  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_glo
43  INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: maskinv_glo
44  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: corners_glo
45  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: corners_lon, corners_lat
46  INTEGER(i_std) :: nbindex_g, kjpindex
47  INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: kindex_g
48  REAL(r_std), DIMENSION(2) :: zoom_lon, zoom_lat
49  CHARACTER(LEN=20) :: calendar, calXIOS
50  !-
51  !- Variables local to each processors.
52  !-
53  INTEGER(i_std) :: i, j, ik, nbdt, first_point
54  INTEGER(i_std) :: nb_forcefile
55  INTEGER(i_std) :: itau, itau_offset, itau_sechiba
56  REAL(r_std)    :: date_start, dt
57  REAL(r_std)    :: timestep_interval(2), timestep_int_next(2), julian
58  INTEGER(i_std) :: rest_id, rest_id_stom
59  INTEGER(i_std) :: hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC
60  INTEGER(i_std)                 :: yeardr, monthdr, daydr    !! Time origin date information
61  REAL(r_std)                    :: secdr                     !! Time origin date information
62!-
63!- input fields
64!-
65  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: u             !! Lowest level wind speed
66  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: v             !! Lowest level wind speed
67  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: pb            !! Surface pressure
68  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: zlev_tq       !! Height of layer for T and q
69  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: zlev_uv       !! Height of layer for u and v
70  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: temp_air      !! Air temperature in Kelvin
71  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: qair          !! Lowest level specific humidity
72  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: ccanopy       !! CO2 concentration in the canopy
73  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: cdrag         !! Cdrag
74  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: precip_rain   !! Rain precipitation
75  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: precip_snow   !! Snow precipitation
76  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: lwdown        !! Down-welling long-wave flux
77  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: swdown        !! Downwelling surface short-wave flux
78  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: swnet         !! Net surface short-wave flux
79  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: solarang      !! Cosine of solar zenith angle
80  !-
81  !- output fields
82  !-
83  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0m           !! Momentum Surface roughness
84  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0h           !! Temperature Surface roughness
85  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: coastalflow   !! Diffuse flow of water into the ocean (m^3/dt)
86  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: riverflow     !! Largest rivers flowing into the ocean (m^3/dt)
87  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: tsol_rad      !! Radiative surface temperature
88  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: vevapp        !! Total of evaporation
89  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: temp_sol_new  !! New soil temperature
90  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: qsurf         !! Surface specific humidity
91  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: albedo        !! Albedo
92  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: fluxsens      !! Sensible chaleur flux
93  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: fluxlat       !! Latent chaleur flux
94  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: emis          !! Emissivity
95  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: netco2        !! netco2flux
96  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: carblu        !! fco2_land_use
97  !-
98  !- Declarations for OASIS
99  !-
100  INTEGER(i_std) :: glo_rank, glo_size ! rank and  number of pe
101  INTEGER(i_std) :: loc_rank, loc_size ! rank and  number of pe
102  INTEGER(i_std) :: LOCAL_OASIS_COMM, TMP_COMM ! local MPI communicator and Initialized
103  INTEGER(i_std) :: comp_id    ! component identification
104  INTEGER(i_std) :: ierror, flag, oasis_info
105  CHARACTER(LEN=4) :: drv_gridname
106  CHARACTER(LEN=8) :: varname
107  INTEGER(i_std), DIMENSION(3) :: ig_paral
108  ! OASIS Output
109  INTEGER(i_std) :: il_part_id, tair_id, qair_id, zlevtq_id, zlevuv_id
110  INTEGER(i_std) :: rainf_id, snowf_id, swnet_id, lwdown_id, solarang_id
111  INTEGER(i_std) :: u_id, v_id, ps_id, cdrag_id
112  ! OASIS Input
113  INTEGER(i_std) :: vevapp_id, fluxsens_id, fluxlat_id, coastal_id, river_id
114  INTEGER(i_std) :: netco2_id, carblu_id, tsolrad_id, tsolnew_id, qsurf_id
115  INTEGER(i_std) :: albnir_id, albvis_id, emis_id, z0m_id, z0h_id
116  !
117  INTEGER(i_std), DIMENSION(2) :: var_nodims, var_shape
118  INTEGER(i_std) :: nbarg, iret, helpmsg = 0
119  CHARACTER(LEN=10) :: arg
120  LOGICAL :: initmode = .FALSE.
121  !-
122  !- XIOS
123  !-
124  LOGICAL                        :: xios_driver_ok = .FALSE.
125  TYPE(xios_context)             :: ctx_hdl               !! Handel for driver2oasis
126  CHARACTER(len=*),PARAMETER     :: id="driver"           !! Id for initialization of driver2oasis in XIOS
127  TYPE(xios_duration)            :: dtime_xios
128  TYPE(xios_date)                :: start_date
129  TYPE(xios_date)                :: time_origin
130  TYPE(xios_fieldgroup)          :: fieldgroup_handle
131  TYPE(xios_field)               :: field_handle
132  TYPE(xios_file)                :: file_handle
133  !-
134  ! -
135  INTEGER(i_std) :: freq_diag=10, debug_lev=1
136  INTEGER(i_std) :: w_unit=737
137  !
138  ! Timer variables
139  !
140  LOGICAL, PARAMETER :: timemeasure=.TRUE.
141  REAL(r_std) :: waitput_cputime=0.0, waitget_cputime=0.0, preparation_cputime=0.0
142  REAL(r_std) :: waitput_walltime=0.0, waitget_walltime=0.0, preparation_walltime=0.0
143  !
144  ! Print point
145  !
146!!  REAL(r_std), DIMENSION(2) :: testpt=(/44.8,-25.3/)
147!!  REAL(r_std), DIMENSION(2) :: testpt=(/44.8,-18.3/)
148!!  REAL(r_std), DIMENSION(2) :: testpt=(/-60.25,-5.25/)
149!!  REAL(r_std), DIMENSION(2) :: testpt=(/-5.25,41.25/)
150  REAL(r_std), DIMENSION(2) :: testpt=(/9999.99,9999.99/)
151!!  REAL(r_std), DIMENSION(2) :: testpt=(/46.7,10.3/)
152!!  REAL(r_std), DIMENSION(2) :: testpt=(/0.25,49.25/)
153  !
154  INTEGER iargc, getarg
155  EXTERNAL iargc, getarg 
156  !-
157  !---------------------------------------------------------------------------------------
158  !-
159  !- The code has 2 execution mode :
160  !-          1) initialisation of the grid file based on the forcing file
161  !-          2) Reading the forcing file and sending the data out with OASIS
162  !-
163  !-
164  nbarg = iargc()
165  IF ( nbarg > 1 ) THEN
166     helpmsg = 1
167  ELSE IF ( nbarg == 1 ) THEN
168     iret = getarg(1,arg)
169     SELECT CASE(arg)
170        !
171     CASE('-h')
172        helpmsg = 1
173     CASE('-init')
174        initmode = .TRUE.
175     CASE DEFAULT
176        helpmsg = 1
177        !
178     END SELECT
179  ELSE
180     initmode = .FALSE.
181  ENDIF
182  !
183  ! Does the user just want help ?
184  !
185  IF ( helpmsg > 0 ) THEN
186     WRITE(*,*) "USAGE : driver2oasis [-init] " 
187     WRITE(*,*) "             The program will read the forcing file provided by variable" 
188     WRITE(*,*) "             FORCING_FILE in the run.def file and do one of 2 things :"
189     WRITE(*,*) "        "
190     WRITE(*,*) "      -init  a grid description file will be generated for the specified "
191     WRITE(*,*) "             forcing file. This grid description file will be written into"
192     WRITE(*,*) "             the file name provided by the variable GRID_FILE "
193     WRITE(*,*) "             of the run.def. "
194     WRITE(*,*) "       "
195     WRITE(*,*) "             If no arguments are provided driver2oasis will initiate and "
196     WRITE(*,*) "             send the forcing data out via OASIS."
197     STOP "HELP from driver2oasis"
198  ENDIF
199  !-
200  !- Open output file for driver
201  !-
202  OPEN(UNIT=w_unit, FILE="out_driver_mono", FORM="formatted")
203  !---------------------------------------------------------------------------------------
204  !-
205  !- Get the general information we need
206  !-
207  !---------------------------------------------------------------------------------------
208  !-
209!!  CALL getin_name("run.def")
210  !
211  !Config Key   = FORCING_FILE
212  !Config Desc  = Name of file containing the forcing data
213  !Config If    = [-]
214  !Config Def   = forcing_file.nc
215  !Config Help  = This is the name of the file which should be opened
216  !Config         for reading the forcing data of the dim0 model.
217  !Config         The format of the file has to be netCDF and COADS
218  !Config         compliant.
219  !Config Units = [FILE]
220  !-
221  forfilename(:)=" "
222  forfilename(1)='forcing_file.nc'
223  CALL getin('FORCING_FILE', forfilename)
224  !
225  !Config Key   = GRID_FILE
226  !Config Desc  = Name of file containing the forcing data
227  !Config If    = [-]
228  !Config Def   = grid_file.nc
229  !Config Help  = This is the name of the file from which we will read
230  !Config         or write into it the description of the grid from
231  !Config         the forcing file.
232  !Config         compliant.
233  !Config Units = [FILE]
234  !-
235  gridfilename='grid_file.nc'
236  CALL getin('GRID_FILE', gridfilename)
237  !-
238  !- Define the zoom
239  !-
240  zoom_lon=(/-180,180/)
241  zoom_lat=(/-90,90/)
242  !
243  !Config Key   = LIMIT_WEST
244  !Config Desc  = Western limit of region
245  !Config If    = [-]
246  !Config Def   = -180.
247  !Config Help  = Western limit of the region we are
248  !Config         interested in. Between -180 and +180 degrees
249  !Config         The model will use the smalest regions from
250  !Config         region specified here and the one of the forcing file.
251  !Config Units = [Degrees]
252  !-
253  CALL getin_p('LIMIT_WEST',zoom_lon(1))
254  !-
255  !Config Key   = LIMIT_EAST
256  !Config Desc  = Eastern limit of region
257  !Config If    = [-]
258  !Config Def   = 180.
259  !Config Help  = Eastern limit of the region we are
260  !Config         interested in. Between -180 and +180 degrees
261  !Config         The model will use the smalest regions from
262  !Config         region specified here and the one of the forcing file.
263  !Config Units = [Degrees]
264  !-
265  CALL getin_p('LIMIT_EAST',zoom_lon(2))
266  !-
267  !Config Key   = LIMIT_NORTH
268  !Config Desc  = Northern limit of region
269  !Config If    = [-]
270  !Config Def   = 90.
271  !Config Help  = Northern limit of the region we are
272  !Config         interested in. Between +90 and -90 degrees
273  !Config         The model will use the smalest regions from
274  !Config         region specified here and the one of the forcing file.
275  !Config Units = [Degrees]
276  !-
277  CALL getin_p('LIMIT_NORTH',zoom_lat(2))
278  !-
279  !Config Key   = LIMIT_SOUTH
280  !Config Desc  = Southern limit of region
281  !Config If    = [-]
282  !Config Def   = -90.
283  !Config Help  = Southern limit of the region we are
284  !Config         interested in. Between 90 and -90 degrees
285  !Config         The model will use the smalest regions from
286  !Config         region specified here and the one of the forcing file.
287  !Config Units = [Degrees]
288  !-
289  CALL getin_p('LIMIT_SOUTH',zoom_lat(1))
290  IF ( (zoom_lon(1)+180 < EPSILON(zoom_lon(1))) .AND. (zoom_lon(2)-180 < EPSILON(zoom_lon(2))) .AND.&
291       &(zoom_lat(1)+90 < EPSILON(zoom_lat(1))) .AND. (zoom_lat(2)-90 < EPSILON(zoom_lat(2))) ) THEN
292     !
293     !Config Key   = WEST_EAST
294     !Config Desc  = Longitude interval to use from the forcing data
295     !Config If    = [-]
296     !Config Def   = -180, 180
297     !Config Help  = This function allows to zoom into the forcing data
298     !Config Units = [degrees east]
299     !-
300     CALL getin('WEST_EAST', zoom_lon)
301     !
302     !Config Key   = SOUTH_NORTH
303     !Config Desc  = Latitude interval to use from the forcing data
304     !Config If    = [-]
305     !Config Def   = -90, 90
306     !Config Help  = This function allows to zoom into the forcing data
307     !Config Units = [degrees north]
308     !-
309     CALL getin('SOUTH_NORTH', zoom_lat)
310  ENDIF
311  !-
312  debug_lev=1
313  CALL getin('BAVARD', debug_lev)
314  IF ( debug_lev .EQ. 4 ) THEN
315     freq_diag=1
316  ELSE IF ( debug_lev .EQ. 3 ) THEN
317     freq_diag=10
318  ELSE IF ( debug_lev .EQ. 2 ) THEN
319     freq_diag=24
320  ELSE
321     freq_diag=96
322  ENDIF
323  !
324  xios_driver_ok=.TRUE.
325  WRITE(w_unit,*) 'Using XIOS :', xios_driver_ok
326  !-
327  !---------------------------------------------------------------------------------------
328  !-
329  !- We go into the mode which initialises the grid of the forcing file and writes it
330  !- for future usage by the driver and ORCHIDEE.
331  !-
332  !---------------------------------------------------------------------------------------
333  !-
334  IF ( initmode ) THEN
335
336     WRITE(w_unit,*) 'Forcing files : ', forfilename(:)
337     WRITE(w_unit,*) 'Grid files : ', gridfilename
338
339     nb_forcefile = 0
340     DO ik=1,100
341        IF ( INDEX(forfilename(ik), '.nc') > 0 ) nb_forcefile = nb_forcefile+1
342     ENDDO
343     !
344     ! This mode of driver2oasis is monoproc and thus we need to force is_root_prc
345     !
346     is_root_prc = .TRUE.
347     !
348     CALL globgrd_getdomsz("NONE", iim_glo, jjm_glo, nbindex_g, model_guess, file_id, forfilename, zoom_lon, zoom_lat)
349     !
350     CALL forcing_getglogrid(nb_forcefile, forfilename, iim_glo, jjm_glo, nbindex_g, .TRUE.)
351
352     CALL forcing_zoomgrid(zoom_lon, zoom_lat, forfilename(1), model_guess, .TRUE.)
353 
354     CALL globgrd_writegrid(gridfilename)
355
356     STOP "Grid file sucessfully generated"
357  ENDIF
358  !-
359  !---------------------------------------------------------------------------------------
360  !-
361  !- This mode will read the grid file and the forcing file and produce the
362  !- data to be handed over to OASIS.
363  !-
364  !---------------------------------------------------------------------------------------
365  !---------------------------------------------------------------------------------------
366  !-
367  !- Define MPI communicator and set-up OASIS
368  !-
369  CALL oasis_init_comp(comp_id, comp_name, ierror)
370  CALL oasis_get_localcomm(TMP_COMM, ierror)
371  !
372  IF ( xios_driver_ok ) THEN
373     CALL xios_initialize(id,local_comm=TMP_COMM, return_comm=LOCAL_OASIS_COMM)
374  ELSE
375     LOCAL_OASIS_COMM = TMP_COMM
376  ENDIF
377  !
378  CALL Init_orchidee_para(LOCAL_OASIS_COMM, .FALSE.)
379  !-
380  !-
381  !---------------------------------------------------------------------------------------
382  !-
383  !- Get the grid associated to the forcing file ... it should be generated by
384  !- this code with a special option before the OASIS coupled case is run.
385  !-
386  !---------------------------------------------------------------------------------------
387  !-
388  IF ( timemeasure ) THEN
389     CALL init_timer
390     CALL start_timer(timer_global)
391  ENDIF
392  !
393  CALL globgrd_getdomsz(gridfilename, iim_glo, jjm_glo, nbindex_g, model_guess, file_id)
394  !-
395  !- Allocation of memory
396  !- variables over the entire grid (thus in x,y)
397  ALLOCATE(lon_glo(iim_glo, jjm_glo))
398  ALLOCATE(lat_glo(iim_glo, jjm_glo))
399  ALLOCATE(mask_glo(iim_glo, jjm_glo))
400  ALLOCATE(area_glo(iim_glo, jjm_glo))
401  ALLOCATE(corners_glo(iim_glo, jjm_glo, 4, 2))
402  !
403  ! Gathered variables
404  ALLOCATE(kindex_g(nbindex_g))
405  ALLOCATE(contfrac(nbindex_g))
406  !-
407  !-
408  !-
409  CALL globgrd_getgrid(file_id, iim_glo, jjm_glo, nbindex_g, model_guess, &
410       &               lon_glo, lat_glo, mask_glo, area_glo, corners_glo,&
411       &               kindex_g, contfrac, calendar)
412  !-
413  !- Set the calendar and get some information
414  !-
415  CALL ioconf_calendar(calendar)
416  CALL ioget_calendar(one_year, one_day)
417  !-
418  !- lalo needs to be created before going into the parallel region
419  !-
420  ALLOCATE(lalo(nbindex_g,2))
421  DO ik=1,nbindex_g
422     !
423     j = ((kindex_g(ik)-1)/iim_glo)+1
424     i = (kindex_g(ik)-(j-1)*iim_glo)
425     !
426     IF ( i > iim_glo .OR. j > jjm_glo ) THEN
427        WRITE(w_unit,*) "Error in the indexing (ik, kindex_g, i, j) : ", ik, kindex_g(ik), i, j
428        STOP "ERROR in driver2oasis"
429     ENDIF
430     !
431     lalo(ik,1) = lat_glo(i,j)
432     lalo(ik,2) = lon_glo(i,j)
433     !
434  ENDDO
435  !
436  WRITE(w_unit,*) "Rank", mpi_rank, " Before opening forcingfile. All land points : ",  nbindex_g
437  WRITE(w_unit,*) "Rank", mpi_rank, " from ", iim_glo, " point in Lon. and ", jjm_glo, "in Lat."
438  !
439  ! Set-up the paralelisation so that all gather and scatters work properly on this monoproc task.
440  !
441  CALL grid_set_glo(iim_glo, jjm_glo, nbindex_g)
442  CALL grid_allocate_glo(4)
443  CALL bcast(nbindex_g)
444  CALL bcast(kindex_g)
445  !
446  WRITE(numout,*) "Rank", mpi_rank, "Into Init_orchidee_data_para_driver with ", nbindex_g,index_g(1)
447  !
448  CALL Init_orchidee_data_para_driver(nbindex_g,kindex_g)
449  CALL init_ioipsl_para 
450  !-
451  !---------------------------------------------------------------------------------------
452  !-
453  !- Open the forcing file and get the time information.
454  !-
455  !---------------------------------------------------------------------------------------
456  !-
457  ! Add w_unit as last arguments in order to get some print_out from the forcing_open routine.
458  !
459  CALL forcing_open(forfilename, iim_glo,  jjm_glo, lon_glo, lat_glo, nbindex_g, zoom_lon, zoom_lat, &
460       &            kindex_g, nbindex_g, w_unit, model_guess)
461  CALL forcing_integration_time(date_start, dt, nbdt)
462  !
463  !
464  ALLOCATE(zlev_tq(nbindex_g), zlev_uv(nbindex_g))
465  ALLOCATE(u(nbindex_g), v(nbindex_g), pb(nbindex_g))
466  ALLOCATE(temp_air(nbindex_g))
467  ALLOCATE(qair(nbindex_g))
468  ALLOCATE(ccanopy(nbindex_g))
469  ALLOCATE(cdrag(nbindex_g))
470  ALLOCATE(precip_rain(nbindex_g))
471  ALLOCATE(precip_snow(nbindex_g))
472  ALLOCATE(swdown(nbindex_g))
473  ALLOCATE(swnet(nbindex_g))
474  ALLOCATE(lwdown(nbindex_g))
475  ALLOCATE(solarang(nbindex_g))
476  ALLOCATE(vevapp(nbindex_g))
477  ALLOCATE(fluxsens(nbindex_g))
478  ALLOCATE(fluxlat(nbindex_g))
479  ALLOCATE(coastalflow(nbindex_g))
480  ALLOCATE(riverflow(nbindex_g))
481  ALLOCATE(netco2(nbindex_g))
482  ALLOCATE(carblu(nbindex_g))
483  ALLOCATE(tsol_rad(nbindex_g))
484  ALLOCATE(temp_sol_new(nbindex_g))
485  ALLOCATE(qsurf(nbindex_g))
486  ALLOCATE(albedo(nbindex_g,2))
487  ALLOCATE(emis(nbindex_g))
488  ALLOCATE(z0m(nbindex_g))
489  ALLOCATE(z0h(nbindex_g))
490  !
491  !-
492  !---------------------------------------------------------------------------------------
493  !-
494  !- OASIS Diagnostics
495  !-
496  !---------------------------------------------------------------------------------------
497  !-
498  ! Unit for output messages : one file for each process
499  !-
500  CALL MPI_Comm_Size (LOCAL_OASIS_COMM, loc_size, ierror )
501  CALL MPI_Comm_Size (MPI_COMM_WORLD, glo_size, ierror )
502  IF (ierror /= 0) THEN
503     WRITE(w_unit,*) 'MPI_comm_size abort by model1 compid ',comp_id
504  ENDIF
505  CALL MPI_Comm_Rank (LOCAL_OASIS_COMM, loc_rank, ierror )
506  CALL MPI_Comm_Rank (MPI_COMM_WORLD, glo_rank, ierror )
507  IF (ierror /= 0) THEN
508     WRITE(0,*) 'MPI_Comm_Rank abort by orchidee comp_id ',comp_id
509  ENDIF
510  !
511  WRITE (w_unit,*) '-----------------------------------------------------------'
512  WRITE (w_unit,*) TRIM(comp_name), ' Running with reals compiled as kind =',r_std
513  WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' Local rank :',loc_rank, &
514       &           " Global rank : ", glo_rank
515  WRITE (w_unit,*) '----------------------------------------------------------'
516  !
517  !
518  WRITE (w_unit,*) 'DD I am the ', TRIM(comp_name), 'local rank', loc_rank, " with ", iim_glo*jjm_glo, "points."
519  WRITE (w_unit,*) 'DD Local number of processors :', loc_size, " Global value :", glo_size
520  WRITE (w_unit,*) 'DD Local MPI communicator is :', LOCAL_OASIS_COMM,  MPI_COMM_WORLD
521  !-
522  !-
523  !---------------------------------------------------------------------------------------
524  !-
525  !- Send the grid to OASIS ... only on the root processor
526  !-
527  !---------------------------------------------------------------------------------------
528  !-
529  IF (loc_rank == 0) THEN
530     !
531     ! TOCOMPLETE - Put here OASIS grid, corner, areas and mask writing calls !
532     !
533     drv_gridname = "DRIV"
534     !
535     CALL oasis_start_grids_writing(flag)
536     !
537     CALL oasis_write_grid(drv_gridname, iim_glo, jjm_glo, lon_glo, lat_glo)
538     !
539     ALLOCATE(corners_lon(iim_glo, jjm_glo, 4), corners_lat(iim_glo, jjm_glo, 4))
540     corners_lon(:,:,:) = corners_glo(:,:,:,1)
541     corners_lat(:,:,:) = corners_glo(:,:,:,2)
542     CALL oasis_write_corner(drv_gridname, iim_glo, jjm_glo, 4, corners_lon, corners_lat)
543     DEALLOCATE(corners_lon, corners_lat)
544     !
545     ALLOCATE(maskinv_glo(iim_glo, jjm_glo))
546     DO i=1,iim_glo
547        DO j=1,jjm_glo
548           IF (mask_glo(i,j) == 0) THEN
549              maskinv_glo(i,j) = 1
550           ELSE
551              maskinv_glo(i,j) = 0
552           ENDIF
553        ENDDO
554     ENDDO
555     CALL oasis_write_mask(drv_gridname, iim_glo, jjm_glo, maskinv_glo)
556     !
557     CALL oasis_write_area(drv_gridname, iim_glo, jjm_glo, area_glo)
558     !
559     CALL oasis_terminate_grids_writing()
560  ENDIF
561  !
562  !-
563  !---------------------------------------------------------------------------------------
564  !-
565  !- Declare the variables
566  !-
567  !---------------------------------------------------------------------------------------
568  !-
569  ig_paral(1) = 0
570  ig_paral(2) = 0
571  ig_paral(3) = nbindex_g
572
573  CALL oasis_def_partition (il_part_id, ig_paral, ierror)
574
575  var_nodims(1) = 1
576  var_nodims(2) = 1
577  var_shape(1) = 1
578  var_shape(1) = nbindex_g
579  !
580  ! Variables SENT to ORCHIDEE
581  ! ==========================
582  !
583  ! Define levels
584  !
585  varname="ZLTQDRIV"
586  CALL oasis_def_var(zlevtq_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
587  varname="ZLUVDRIV"
588  CALL oasis_def_var(zlevuv_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
589  !
590  ! Define scalar atmospheric variables
591  !
592  varname="TAIRDRIV"
593  CALL oasis_def_var(tair_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
594  varname="QAIRDRIV"
595  CALL oasis_def_var(qair_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
596  !
597  ! Define precipitation variables
598  !
599  varname="RAINDRIV"
600  CALL oasis_def_var(rainf_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
601  varname="SNOWDRIV"
602  CALL oasis_def_var(snowf_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
603  !
604  ! Define radiation variables
605  !
606  varname="SWNEDRIV"
607  CALL oasis_def_var(swnet_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
608  varname="LWDODRIV"
609  CALL oasis_def_var(lwdown_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
610  varname="SOLADRIV"
611  CALL oasis_def_var(solarang_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
612  !
613  ! Finaly pressure and wind
614  !
615  varname="UWINDRIV"
616  CALL oasis_def_var(u_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
617  varname="VWINDRIV"
618  CALL oasis_def_var(v_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
619  varname="PRESDRIV"
620  CALL oasis_def_var(ps_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
621  varname="DRAGDRIV"
622  CALL oasis_def_var(cdrag_id, varname, il_part_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror)
623  !
624  !
625  ! Variables RECEIVED from ORCHIDEE
626  ! ================================
627  !
628  !
629  ! Turbulent fluxes
630  !
631  varname="EVAPDRIV"
632  CALL oasis_def_var(vevapp_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
633  !
634  varname="SENSDRIV"
635  CALL oasis_def_var(fluxsens_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
636  !
637  varname="LATEDRIV"
638  CALL oasis_def_var(fluxlat_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
639  !
640  ! Discharge to the oceans
641  !
642  varname="COASDRIV"
643  CALL oasis_def_var(coastal_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
644  !
645  varname="RIVEDRIV"
646  CALL oasis_def_var(river_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
647  !
648  ! Carbon fluxes
649  !
650  varname="NECODRIV"
651  CALL oasis_def_var(netco2_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
652  !
653  varname="LUCODRIV"
654  CALL oasis_def_var(carblu_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
655  !
656  ! Surface states
657  !
658  varname="TRADDRIV"
659  CALL oasis_def_var(tsolrad_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
660  !
661  varname="TNEWDRIV"
662  CALL oasis_def_var(tsolnew_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
663  !
664  varname="QSURDRIV"
665  CALL oasis_def_var(qsurf_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
666  !
667  varname="ANIRDRIV"
668  CALL oasis_def_var(albnir_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
669  !
670  varname="AVISDRIV"
671  CALL oasis_def_var(albvis_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
672  !
673  varname="EMISDRIV"
674  CALL oasis_def_var(emis_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
675  !
676  varname="Z0MODRIV"
677  CALL oasis_def_var(z0m_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
678  !
679  varname="Z0HEDRIV"
680  CALL oasis_def_var(z0h_id, varname, il_part_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
681  !
682  CALL oasis_enddef (ierror)
683  !-
684  WRITE(w_unit,*) 'OASIS Vars declared'
685  !
686  !---------------------------------------------------------------------------------------
687  !-
688  !- Set-up XIOS
689  !-
690  !---------------------------------------------------------------------------------------
691  !
692  IF ( xios_driver_ok ) THEN
693     CALL xios_context_initialize("driver2oasis", LOCAL_OASIS_COMM)
694     CALL xios_get_handle("driver2oasis",ctx_hdl)
695     CALL xios_set_current_context(ctx_hdl)
696     !
697     CALL ju2ymds(date_start, yeardr, monthdr, daydr, secdr)
698     !
699     dtime_xios%second=dt
700     IF (INDEX(calendar, 'gregor') > 0 .OR. INDEX(calendar, 'Gregor') > 0 .OR. &
701          &    calendar == 'standard') THEN
702        calXIOS = 'gregorian'
703     ELSE IF (calendar == 'noleap') THEN
704        calXIOS = 'noleap'
705     ELSE IF (calendar == '360d') THEN
706        calXIOS = 'd360'
707     END IF
708     
709     CALL xios_define_calendar(type=calXIOS, start_date=xios_date(yeardr,monthdr,daydr,0,0,0), &
710          time_origin=xios_date(yeardr,monthdr,daydr,0,0,0), timestep=dtime_xios)
711     CALL xios_set_domain_attr("domain_landpoints", ni_glo=iim_glo, nj_glo=jjm_glo, type="curvilinear")
712     CALL xios_set_domain_attr("domain_landpoints",lonvalue_2d=lon_glo,latvalue_2d=lat_glo)
713     CALL xios_close_context_definition()
714  ENDIF
715  WRITE(w_unit,*) 'XIOS Closed Context'
716  !
717  !---------------------------------------------------------------------------------------
718  !-
719  !- Get a first set of forcing data
720  !-
721  !---------------------------------------------------------------------------------------
722  !-
723  timestep_interval(1) = date_start
724  timestep_interval(2) = date_start + (dt/one_day)
725  CALL forcing_getvalues(timestep_interval, dt, zlev_tq, zlev_uv, temp_air, qair, &
726       &                 precip_rain, precip_snow, swdown, lwdown, solarang, u, v, pb)
727  ! An atmsopheric model will provide this information. If zero ORCHIDEE will compute it.
728  cdrag(:) = 0.0
729  ! Transform the swdown into a swnet
730  albedo(:,:)=0.13
731  swnet(:) = (1.-(albedo(:,1)+albedo(:,2))/2.)*swdown(:)
732  !-
733  !---------------------------------------------------------------------------------------
734  !-
735  !- Go into the time loop
736  !-
737  !---------------------------------------------------------------------------------------
738  !-
739  IF ( timemeasure ) THEN
740     WRITE(w_unit,*) '------> CPU Time for start-up of driver : ',Get_cpu_Time(timer_global)
741     WRITE(w_unit,*) '------> Real Time for start-up of driver : ',Get_real_Time(timer_global)
742     CALL stop_timer(timer_global)
743     CALL start_timer(timer_global)
744  ENDIF
745  !-
746  DO itau = 0,nbdt-1
747     !
748     timestep_interval(1) = date_start + itau*(dt/one_day)
749     timestep_interval(2) = date_start + (itau+1)*(dt/one_day)
750     julian = date_start + (itau+0.5)*(dt/one_day)
751     !
752     CALL xios_set_current_context(ctx_hdl)
753     CALL xios_update_calendar(itau+1)
754     !
755     ! Write some diagnostics to look at the shape of the maps
756     !
757     IF ( itau == 0 ) THEN
758       CALL globgrd_writevar(iim_glo, jjm_glo, lon_glo, lat_glo, &
759          &                  nbindex_g, lalo, temp_air, "TAIR", "Tair_Driver_0.nc")
760     ENDIF
761     !
762     ! Put first the height of the atmospheric levels
763     !
764     CALL forcing_printpoint(julian, testpt(1), testpt(2), zlev_tq, "TQ height")
765     CALL oasis_put(zlevtq_id, NINT(itau*dt), zlev_tq, oasis_info)
766     CALL oasis_put(zlevuv_id, NINT(itau*dt), zlev_uv, oasis_info)
767     !
768     !
769     CALL forcing_printpoint(julian, testpt(1), testpt(2), temp_air, "Air temperature")
770     CALL forcing_printpoint(julian, testpt(1), testpt(2), qair, "Air humidity")
771     !
772     CALL oasis_put(tair_id, NINT(itau*dt), temp_air, oasis_info)
773     CALL oasis_put(qair_id, NINT(itau*dt), qair, oasis_info)
774     !
775     ! Precipitation fluxes
776     !
777     CALL forcing_printpoint(julian, testpt(1), testpt(2), precip_rain*one_day, "Rainfall")
778     CALL forcing_printpoint(julian, testpt(1), testpt(2), precip_snow*one_day, "Snowfall")
779     CALL oasis_put(rainf_id, NINT(itau*dt), precip_rain, oasis_info)
780     CALL oasis_put(snowf_id, NINT(itau*dt), precip_snow, oasis_info)
781     !
782     ! Radiation fluxes
783     !
784     CALL forcing_printpoint(julian, testpt(1), testpt(2), swnet, "Net solar")
785     CALL oasis_put(swnet_id, NINT(itau*dt), swnet, oasis_info)
786     CALL oasis_put(lwdown_id, NINT(itau*dt), lwdown, oasis_info)
787     CALL forcing_printpoint(julian, testpt(1), testpt(2), lwdown, "Downward Longwave") 
788     CALL oasis_put(solarang_id, NINT(itau*dt), solarang, oasis_info)
789     !
790     ! Dynamical fields
791     !
792     CALL forcing_printpoint(julian, testpt(1), testpt(2), u, "East-ward wind")
793     CALL forcing_printpoint(julian, testpt(1), testpt(2), pb, "Surface Pressure")
794     CALL oasis_put(u_id, NINT(itau*dt), u, oasis_info)
795     CALL oasis_put(v_id, NINT(itau*dt), v, oasis_info)
796     CALL oasis_put(ps_id, NINT(itau*dt), pb, oasis_info)
797     CALL oasis_put(cdrag_id, NINT(itau*dt), cdrag, oasis_info)
798     !
799     IF ( timemeasure ) THEN
800        waitput_cputime = waitput_cputime + Get_cpu_Time(timer_global)
801        waitput_walltime = waitput_walltime + Get_real_Time(timer_global)
802        CALL stop_timer(timer_global)
803        CALL start_timer(timer_global)
804     ENDIF
805     !
806     ! Get the forcing data for the next time step before we go into the blocking gets
807     !
808     IF ( itau < (nbdt-1) ) THEN
809        timestep_int_next(1) = date_start + (itau+1)*(dt/one_day)
810        timestep_int_next(2) = date_start + (itau+2)*(dt/one_day)
811        CALL forcing_getvalues(timestep_int_next, dt, zlev_tq, zlev_uv, temp_air, qair, &
812             &                 precip_rain, precip_snow, swdown, lwdown, solarang, u, v, pb)
813        ! An atmsopheric model will provide this information. If zero ORCHIDEE will compute it.
814        cdrag(:) = 0.0
815        ! Compute swnet with the albedo computed in the previous call to ORCHIDEE
816        swnet(:) = (1.-(albedo(:,1)+albedo(:,2))/2.)*swdown(:)
817        !
818     ENDIF
819     !
820     ! Timing of computations for next forcing data
821     !
822     IF ( timemeasure ) THEN
823        preparation_cputime = preparation_cputime + Get_cpu_Time(timer_global)
824        preparation_walltime = preparation_walltime + Get_real_Time(timer_global)
825        CALL stop_timer(timer_global)
826        CALL start_timer(timer_global)
827     ENDIF
828     !
829     ! Go into the blocking gets from OASIS
830     !
831     CALL oasis_get(vevapp_id, NINT((itau)*dt), vevapp, oasis_info)
832     CALL oasis_get(fluxsens_id, NINT((itau)*dt), fluxsens, oasis_info)
833     CALL oasis_get(fluxlat_id, NINT((itau)*dt), fluxlat, oasis_info)
834     !
835     CALL oasis_get(coastal_id, NINT((itau)*dt), coastalflow, oasis_info)
836     CALL oasis_get(river_id, NINT((itau)*dt), riverflow, oasis_info)
837     !
838     CALL oasis_get(netco2_id, NINT((itau)*dt), netco2, oasis_info)
839     CALL oasis_get(carblu_id, NINT((itau)*dt), carblu, oasis_info)
840     !
841     CALL oasis_get(tsolrad_id, NINT((itau)*dt), tsol_rad, oasis_info)
842     CALL oasis_get(tsolnew_id, NINT((itau)*dt), temp_sol_new, oasis_info)
843     CALL oasis_get(qsurf_id, NINT((itau)*dt), qsurf, oasis_info)
844     !
845     CALL oasis_get(albvis_id, NINT((itau)*dt), albedo(:,1), oasis_info)
846     CALL oasis_get(albnir_id, NINT((itau)*dt), albedo(:,2), oasis_info)
847     CALL oasis_get(emis_id, NINT((itau)*dt), emis, oasis_info)
848     CALL oasis_get(z0m_id, NINT((itau)*dt), z0m, oasis_info)
849     CALL oasis_get(z0h_id, NINT((itau)*dt), z0h, oasis_info)
850     CALL forcing_printpoint(julian, testpt(1), testpt(2), vevapp, "Recived evap")
851     !
852     ! Timing of waits
853     !
854     IF ( timemeasure ) THEN
855        waitget_cputime = waitget_cputime + Get_cpu_Time(timer_global)
856        waitget_walltime = waitget_walltime + Get_real_Time(timer_global)
857        CALL stop_timer(timer_global)
858        CALL start_timer(timer_global)
859     ENDIF
860     !
861     !---------------------------------------------------------------------------------------
862     ! Some dianostics
863     !---------------------------------------------------------------------------------------
864     !
865     IF ( MOD(itau, freq_diag) == 0 ) THEN
866        !
867        CALL forcing_printdate(timestep_interval(1), "Diagnostics START", w_unit)
868        !
869        WRITE(w_unit,*) MINVAL(temp_sol_new), " << TEMP_SOL_NEW << ", MAXVAL(temp_sol_new)
870        WRITE(w_unit,*) MINVAL(vevapp), " << VEVAPP << ", MAXVAL(vevapp)
871        WRITE(w_unit,*) MINVAL(fluxsens), " << FLUXSENS << ", MAXVAL(fluxsens)
872        WRITE(w_unit,*) MINVAL(fluxlat), " << FLUXLAT << ", MAXVAL(fluxlat)
873        WRITE(w_unit,*) MINVAL(coastalflow), " <<  COASTALFLOW << ", MAXVAL(coastalflow)
874        WRITE(w_unit,*) MINVAL(riverflow), " << RIVERFLOW << ", MAXVAL(riverflow)
875        WRITE(w_unit,*) MINVAL(netco2), " << NETCO2 << ", MAXVAL(netco2)
876        WRITE(w_unit,*) MINVAL(carblu), " << CARBLU << ", MAXVAL(carblu)
877        WRITE(w_unit,*) MINVAL(tsol_rad), " << TSOL_RAD, << ", MAXVAL(tsol_rad)
878        WRITE(w_unit,*) MINVAL(temp_sol_new), " << TEMP_SOL_NEW << ", MAXVAL(temp_sol_new)
879        WRITE(w_unit,*) MINVAL(qsurf), " << QSURF << ", MAXVAL(qsurf)
880        WRITE(w_unit,*) MINVAL(albedo(:,1)), " << ALBEDO VIS << ", MAXVAL(albedo(:,1))
881        WRITE(w_unit,*) MINVAL(albedo(:,2)), " << ALBEDO NIR << ", MAXVAL(albedo(:,2))
882        WRITE(w_unit,*) MINVAL(emis), " << EMIS << ", MAXVAL(emis)
883        WRITE(w_unit,*) MINVAL(z0m), " << Z0M << ", MAXVAL(z0m)
884        !
885        CALL forcing_printdate(timestep_interval(2), "Diagnostics END", w_unit)
886        !
887     ENDIF
888     !
889  ENDDO
890  !-
891  !-
892  !-
893  IF ( timemeasure ) THEN
894     WRITE(w_unit,*) '------> Total CPU Time waiting for put to ORCH : ',waitput_cputime
895     WRITE(w_unit,*) '------> Total Real Time waiting for put to ORCH : ',waitput_walltime
896     WRITE(w_unit,*) '------> Total CPU Time for preparing forcing : ', preparation_cputime
897     WRITE(w_unit,*) '------> Total Real Time for preparing forcing : ', preparation_walltime
898     WRITE(w_unit,*) '------> Total CPU Time waiting for get from ORCH : ',waitget_cputime
899     WRITE(w_unit,*) '------> Total Real Time waiting for get from ORCH : ',waitget_walltime
900     CALL stop_timer(timer_global)
901  ENDIF
902  !-
903  !---------------------------------------------------------------------------------------
904  !-
905  !- Close OASIS and MPI
906  !-
907  !---------------------------------------------------------------------------------------
908  !-
909  CALL forcing_close()
910  WRITE(w_unit,*) "Closed forcing"
911  !-
912  IF ( xios_driver_ok ) THEN
913     CALL xios_context_finalize()
914     CALL xios_finalize()
915     WRITE(w_unit,*) "XIOS finalized"
916     CALL oasis_terminate(ierror)
917  ELSE
918     CALL oasis_terminate(ierror)
919  ENDIF
920  WRITE(w_unit,*) "OASIS closed"
921  !-
922  CLOSE(UNIT=w_unit)
923  !-
924END PROGRAM driver2oasis
Note: See TracBrowser for help on using the repository browser.