source: branches/publications/ORCHIDEE_gmd-2018-57/src_oasisdriver/driver2oasis.f90 @ 6145

Last change on this file since 6145 was 4513, checked in by jan.polcher, 7 years ago

This commit corrects an error in the output of the largest rivers together with the location of the estuaries.

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