source: branches/publications/ORCHIDEE_gmd-2018-261/src_oasisdriver/driver2oasis.f90 @ 7367

Last change on this file since 7367 was 3447, checked in by josefine.ghattas, 9 years ago

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