source: branches/publications/ORCHIDEE_GLUC_r6545/src_oasisdriver/driver2oasis.f90 @ 6737

Last change on this file since 6737 was 4359, checked in by albert.jornet, 7 years ago

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