source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/examples/test_interpolation/model2.F90 @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 5 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 16.6 KB
Line 
1!------------------------------------------------------------------------
2! Copyright 2010, CERFACS, Toulouse, France.
3! All rights reserved. Use is subject to OASIS3 license terms.
4!=============================================================================
5!
6PROGRAM model2
7  !
8  ! Use for netCDF library
9  USE netcdf
10  ! Use for OASIS communication library
11  USE mod_oasis
12  !
13  ! Use module to read the data
14  USE read_all_data
15  !
16  ! Use module to write the data
17  USE write_all_fields
18  !
19  IMPLICIT NONE
20  !
21  INCLUDE 'mpif.h'
22  !
23  ! By default OASIS3-MCT exchanges data in double precision,
24  ! but conversion to or from single precision data is supported in the interface
25#ifdef NO_USE_DOUBLE_PRECISION
26  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(6,37)   ! real
27#elif defined USE_DOUBLE_PRECISION
28  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
29#endif
30  !
31  CHARACTER(len=30), PARAMETER   :: data_gridname='grids.nc' ! file with the grids
32  CHARACTER(len=30), PARAMETER   :: data_maskname='masks.nc' ! file with the masks
33  CHARACTER(len=30)              :: data_filename, field_name
34  !
35  ! Component name (6 characters) same as in the namcouple
36  CHARACTER(len=6)   :: comp_name = 'model2'
37  CHARACTER(len=128) :: comp_out ! name of the output log file
38  CHARACTER(len=3)   :: chout
39  CHARACTER(len=4)   :: cl_grd_src ! name of the source grid
40  CHARACTER(len=4)   :: cl_grd_tgt ! name of the target grid
41  NAMELIST /grid_source_characteristics/cl_grd_src
42  NAMELIST /grid_target_characteristics/cl_grd_tgt
43  !
44  ! Global grid parameters :
45  INTEGER :: nlon, nlat, ntot    ! dimensions in the 2 directions of space + total size
46  INTEGER :: il_size
47  INTEGER :: nc             ! number of corners in the (i,j) plan
48  !
49  REAL (kind=wp), DIMENSION(:,:), POINTER    :: gg_lon,gg_lat
50  INTEGER, DIMENSION(:,:), POINTER           :: gg_mask ! mask, 0 == valid point, 1 == masked point
51  !
52  INTEGER :: mype, npes ! rank and number of pe
53  INTEGER :: localComm  ! local MPI communicator and Initialized
54  INTEGER :: comp_id    ! component identification
55  !
56  INTEGER, DIMENSION(:), ALLOCATABLE :: il_paral ! Decomposition for each proc
57  !
58  INTEGER :: ierror, rank, w_unit
59  INTEGER :: i, j, ij, ic_msk
60  INTEGER :: FILE_Debug=1
61  !
62  ! Names of exchanged Fields
63  CHARACTER(len=8), PARAMETER :: var_name = 'FRECVANA' ! 8 characters field received by the atmosphere from the ocean
64  !
65  ! Used in oasis_def_var and oasis_def_var
66  INTEGER                       :: var_id 
67  INTEGER                       :: var_nodims(2) 
68  INTEGER                       :: var_type
69  !
70  REAL (kind=wp), PARAMETER     :: field_ini = -1. ! initialisation of received fields
71  INTEGER,        PARAMETER     :: err_msk = -10000
72  !
73  INTEGER               ::  ib
74  INTEGER, PARAMETER    ::  il_nb_time_steps = 1 ! number of time steps
75  INTEGER, PARAMETER    ::  delta_t = 3600      ! time step
76  !
77  INTEGER                       :: il_flag          ! Flag for grid writing
78  !
79  INTEGER                       :: itap_sec ! Time
80  !
81  ! Grid parameter definition
82  INTEGER                       :: part_id  ! use to connect the partition to the variables
83                                            ! in oasis_def_var
84  INTEGER                       :: var_sh(4) ! local dimensions of the arrays to the pe
85                                             ! 2 x field rank (= 4 because fields are of rank = 2)
86  INTEGER :: indi_beg, indi_end, indj_beg, indj_end
87  !
88  ! Local fields arrays used in routines oasis_put and oasis_get
89  REAL (kind=wp), POINTER       :: field_recv(:,:), field_ana(:,:), error(:,:)
90  !
91  ! Global array to plot the error by proc 0 and calculate min and max
92  REAL (kind=wp), POINTER       :: global_error(:)
93  INTEGER, POINTER              :: global_shapes(:,:)
94  INTEGER, POINTER              :: displs(:),rcounts(:)
95  !
96  ! Min and Max of the error of interpolation
97  REAL (kind=wp)             :: min,max
98  !
99  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
100  !   INITIALISATION
101  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102  !
103  CALL mpi_init(ierror)
104  !!!!!!!!!!!!!!!!! OASIS_INIT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105  !
106  CALL oasis_init_comp (comp_id, comp_name, ierror )
107  IF (ierror /= 0) THEN
108      WRITE(0,*) 'oasis_init_comp abort by model2 compid ',comp_id
109      CALL oasis_abort(comp_id,comp_name,'Problem at line 109')
110  ENDIF
111  !
112  ! Unit for output messages : one file for each process
113  CALL MPI_Comm_Rank ( MPI_COMM_WORLD, rank, ierror )
114  IF (ierror /= 0) THEN
115      WRITE(0,*) 'MPI_Comm_Rank abort by model2 compid ',comp_id
116      CALL oasis_abort(comp_id,comp_name,'Problem at line 116')
117  ENDIF
118  !
119  !
120  !!!!!!!!!!!!!!!!! OASIS_GET_LOCALCOMM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121  !
122  CALL oasis_get_localcomm ( localComm, ierror )
123  IF (ierror /= 0) THEN
124      WRITE (0,*) 'oasis_get_localcomm abort by model2 compid ',comp_id
125      CALL oasis_abort(comp_id,comp_name,'Problem at line 125')
126  ENDIF
127  !
128  ! Get MPI size and rank
129  CALL MPI_Comm_Size ( localComm, npes, ierror )
130  IF (ierror /= 0) THEN
131      WRITE(0,*) 'MPI_comm_size abort by model2 compid ',comp_id
132      CALL oasis_abort(comp_id,comp_name,'Problem at line 132')
133  ENDIF
134  !
135  CALL MPI_Comm_Rank ( localComm, mype, ierror )
136  IF (ierror /= 0) THEN
137      WRITE (0,*) 'MPI_Comm_Rank abort by model2 compid ',comp_id
138      CALL oasis_abort(comp_id,comp_name,'Problem at line 138')
139  ENDIF
140  !
141  !
142  IF ((FILE_Debug == 1) .AND. (mype == 0)) FILE_Debug=2
143  !
144  IF (FILE_Debug <= 1) THEN
145      IF (mype == 0) THEN
146          w_unit = 100 + rank
147          WRITE(chout,'(I3)') w_unit
148          comp_out=comp_name//'.root_'//chout
149          OPEN(w_unit,file=TRIM(comp_out),form='formatted')
150      ELSE
151          w_unit = 15
152          comp_out=comp_name//'.notroot'
153          OPEN(w_unit,file=TRIM(comp_out),form='formatted',position='append')
154      ENDIF
155  ELSE
156      w_unit = 100 + rank
157      WRITE(chout,'(I3)') w_unit
158      comp_out=comp_name//'.out_'//chout
159      OPEN(w_unit,file=TRIM(comp_out),form='formatted')
160  ENDIF
161  !
162  IF (FILE_Debug >= 2) THEN
163      OPEN(w_unit,file=TRIM(comp_out),form='formatted')
164      WRITE (w_unit,*) '-----------------------------------------------------------'
165      WRITE (w_unit,*) TRIM(comp_name), ' Running with reals compiled as kind =',wp
166      WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' rank :',rank
167      WRITE (w_unit,*) '----------------------------------------------------------'
168      WRITE(w_unit,*) 'I am the', TRIM(comp_name), ' ', 'comp', comp_id, 'local rank', mype
169      WRITE (w_unit,*) 'Number of porcessors :',npes
170      CALL FLUSH(w_unit)
171  ENDIF
172  !
173  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174  !  GRID DEFINITION
175  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
176  !
177  ! Reading global grids.nc and masks.nc netcdf files
178  ! Get arguments giving source grid acronym and field type
179  !
180  OPEN(UNIT=70,FILE='name_grids.dat',FORM='FORMATTED')
181  READ(UNIT=70,NML=grid_source_characteristics)
182  READ(UNIT=70,NML=grid_target_characteristics)
183  CLOSE(70)
184  !
185  IF (FILE_Debug >= 2) THEN
186      WRITE(w_unit,*) 'Target grid name :',cl_grd_tgt
187      CALL flush(w_unit)
188  ENDIF
189  !
190  ! Reading dimensions of the grid
191  CALL read_dimgrid(nlon,nlat,data_gridname,cl_grd_tgt,w_unit,FILE_Debug)
192  nc=4
193  !
194  ! Allocation
195  ALLOCATE(gg_lon(nlon,nlat), STAT=ierror )
196  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lon'
197  ALLOCATE(gg_lat(nlon,nlat), STAT=ierror )
198  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lat'
199  ALLOCATE(gg_mask(nlon,nlat), STAT=ierror )
200  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating indice_mask'
201  !
202  ! Read global grid longitudes, latitudes and mask
203  !
204  CALL read_grid(nlon,nlat, data_gridname, cl_grd_tgt, w_unit, FILE_Debug, &
205                 gg_lon,gg_lat)
206  CALL read_mask(nlon,nlat, data_maskname, cl_grd_tgt, w_unit, FILE_Debug, &
207                 gg_mask)
208  !
209  IF (FILE_Debug >= 2) THEN
210      WRITE(w_unit,*) 'After grids writing'
211      CALL FLUSH(w_unit)
212  ENDIF
213  !
214  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215  !  PARTITION DEFINITION
216  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
217  !
218  ! Definition of the partition of the grid
219  ntot=nlon*nlat
220#ifdef DECOMP_APPLE
221  il_size = 3
222#elif defined DECOMP_BOX
223  il_size = 5
224#endif
225  ALLOCATE(il_paral(il_size))
226  IF (FILE_Debug >= 2) THEN
227      WRITE(w_unit,*) 'After allocate il_paral, il_size', il_size
228      CALL FLUSH(w_unit)
229  ENDIF
230  !
231  CALL decomp_def (part_id,il_paral,il_size,nlon,nlat,mype,npes,w_unit)
232  IF (FILE_Debug >= 2) THEN
233      WRITE(w_unit,*) 'After decomp_def, il_paral = ', il_paral(:)
234      CALL FLUSH(w_unit)
235  ENDIF
236  !
237  CALL oasis_def_partition (part_id, il_paral, ierror)
238  !
239  !
240  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
241  ! DEFINITION OF THE LOCAL FIELDS 
242  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243  !
244  !!!!!!!!!!!!!!!!!! OASIS_DEF_VAR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
245  !
246  ! Define transient variables
247  !
248  var_nodims(1) = 2    ! Rank of the field array is 2
249  var_nodims(2) = 1    ! Bundles always 1 for OASIS3
250  var_type = OASIS_Real
251  !
252  var_sh(1) = 1
253  var_sh(2) = il_paral(3)
254  var_sh(3) = 1 
255#ifdef DECOMP_APPLE
256  var_sh(4) = 1
257#elif defined DECOMP_BOX
258  var_sh(4) = il_paral(4)
259#endif
260  !
261  ! Declaration of the field associated with the partition of the grid
262  CALL oasis_def_var (var_id,var_name, part_id, &
263     var_nodims, OASIS_In, var_sh, var_type, ierror)
264  IF (ierror /= 0) THEN
265      WRITE (w_unit,*) 'oasis_def_var abort by model2 compid ',comp_id
266      CALL oasis_abort(comp_id,comp_name,'Problem at line 266')
267  ENDIF
268  !
269  !
270  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
271  !         TERMINATION OF DEFINITION PHASE
272  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273  !  All processes involved in the coupling must call oasis_enddef;
274  !  here all processes are involved in coupling
275  !
276  !!!!!!!!!!!!!!!!!! OASIS_ENDDEF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
277  !
278  CALL oasis_enddef ( ierror )
279  IF (ierror /= 0) THEN
280      WRITE (w_unit,*) 'oasis_enddef abort by model2 compid ',comp_id
281      CALL oasis_abort(comp_id,comp_name,'Problem at line 281')
282  ENDIF
283  !
284  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
285  ! RECEIVE ARRAYS AND CALCULATE THE ERROR
286  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
287  !
288  ! Allocate the fields send and received by the model
289  !
290  ALLOCATE(field_recv(var_sh(2), var_sh(4)), STAT=ierror )
291  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_recv'
292  ALLOCATE(field_ana(var_sh(2), var_sh(4)),STAT=ierror )
293  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_ana'
294  ALLOCATE(error(var_sh(2), var_sh(4)),STAT=ierror )
295  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating error'
296  !
297  !!!!!!!!!!!!!!!!!!!!!!!!!!!!OASIS_PUT/OASIS_GET !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
298  !
299  indi_beg=1 ; indi_end=nlon
300  indj_beg=((nlat/npes)*mype)+1 
301  !
302  IF (mype .LT. npes - 1) THEN
303      indj_end = (nlat/npes)*(mype+1)
304  ELSE
305      indj_end = nlat 
306  ENDIF
307  !
308  ! Data exchange in time loop
309  !
310  ib=1
311  itap_sec = delta_t * (ib-1) ! Time
312  !
313  CALL function_ana(var_sh(2),&
314                    var_sh(4), &
315                    RESHAPE(gg_lon(indi_beg:indi_end,indj_beg:indj_end),(/ var_sh(2), var_sh(4) /)), &
316                    RESHAPE(gg_lat(indi_beg:indi_end,indj_beg:indj_end),(/ var_sh(2), var_sh(4) /)), &
317                    field_ana,ib)
318  !
319  ! Get the field FRECVANA
320  field_recv=field_ini
321  CALL oasis_get(var_id,itap_sec, field_recv, ierror)
322  IF (FILE_Debug >= 2) THEN
323      WRITE(w_unit,*) 'tcx recvf1 ',itap_sec,MINVAL(field_recv),MAXVAL(field_recv)
324  ENDIF
325  IF ( ierror .NE. OASIS_Ok .AND. ierror .LT. OASIS_Recvd) THEN
326      WRITE (w_unit,*) 'oasis_get abort by model2 compid ',comp_id
327      CALL oasis_abort(comp_id,comp_name,'Problem at line 327')
328  ENDIF
329  !
330  !
331  IF (FILE_Debug >= 2) THEN
332      WRITE (w_unit,*) 'Calculate the error of interpolation'
333      CALL FLUSH(w_unit)
334  ENDIF
335  !
336  error=err_msk
337  !
338  WHERE (RESHAPE(gg_mask(indi_beg:indi_end,indj_beg:indj_end),(/ var_sh(2), var_sh(4) /)) == 0)
339      error = ABS(((field_ana - field_recv)/field_recv))*100
340  END WHERE
341  !
342  !
343  !!!!!!!!!!!!!!!!!!!!!!! Write the error and the field in a NetCDF file by proc 0 !!!!!!!!!!!!!!!!!!!!
344  ! field_recv is 0 on masked points => put error = err_msk on these points
345  !
346  IF (mype == 0) THEN
347      !
348      ALLOCATE(global_error(ntot),STAT=ierror )
349      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating global_error'
350      ALLOCATE(global_shapes(4,npes),STAT=ierror )
351      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating global_shapes'
352      ALLOCATE(displs(npes),STAT=ierror )
353      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating displs'
354      ALLOCATE(rcounts(npes),STAT=ierror )
355      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating rcounts'
356      ! 
357      global_error=0.
358      global_shapes=0
359      displs=0
360      rcounts=0
361      !
362  ENDIF
363  !     
364  !
365  CALL MPI_GATHER(var_sh,4,MPI_INTEGER,global_shapes,4,MPI_INTEGER,0,localComm,ierror)
366  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error collecting shapes'
367  IF (mype == 0) THEN
368      IF (FILE_Debug >= 2) THEN
369          WRITE (w_unit,*) 'Shapes of all the processes :',global_shapes
370          CALL FLUSH(w_unit)
371      ENDIF
372  ENDIF
373      !
374  IF (mype == 0) THEN
375      displs(1)=0
376      DO i=2,npes
377        displs(i) = displs(i-1) + global_shapes(2,i-1) * global_shapes(4,i-1)
378      ENDDO
379      DO i=1,npes
380        rcounts(i) = global_shapes(2,i) * global_shapes(4,i)
381      ENDDO
382      !
383      IF (FILE_Debug >= 2) THEN
384          WRITE(w_unit,*) 'displs :', displs
385          WRITE(w_unit,*) 'rcounts :', rcounts
386          CALL FLUSH(w_unit)
387      ENDIF
388  ENDIF
389      !
390#ifdef NO_USE_DOUBLE_PRECISION
391      CALL MPI_gatherv(error, var_sh(2)*var_sh(4),MPI_REAL,&
392                       global_error,rcounts,displs,MPI_REAL,0,localComm,ierror)
393      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error collecting errors'
394      !
395#elif defined USE_DOUBLE_PRECISION
396      CALL MPI_gatherv(error, var_sh(2)*var_sh(4),MPI_REAL8,&
397                       global_error,rcounts,displs,MPI_REAL8,0,localComm,ierror)
398      IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error collecting errors'
399      !
400#endif
401      !
402  IF (mype == 0) THEN
403      !
404      data_filename='error.nc'
405      field_name='error'
406      CALL write_field(nlon,nlat, &
407                       data_filename, field_name,  &
408                       w_unit, FILE_Debug, &
409                       gg_lon, gg_lat, global_error)
410  ENDIF
411  !
412  !!!!!!!!!!!!!!!!!!!! Min and Max of the error on non masked points calculated by proc 0 !!!!!!!!!!!!!!!!
413  !
414  IF (mype == 0) THEN
415      ic_msk=0
416      DO i=1,nlon
417        DO j=1,nlat
418          ij=i+(j-1)*(nlon)
419          IF ( gg_mask(i,j) == 1 ) THEN
420              global_error(ij) = -err_msk
421              ic_msk = ic_msk + 1
422          ENDIF
423        ENDDO
424      ENDDO
425      min=MINVAL(REAL(global_error))
426      IF (FILE_Debug >= 2) THEN
427          WRITE(w_unit,*) 'Min and its location in the error field : ',min
428          WRITE(w_unit,*) MINLOC(REAL(global_error))
429          CALL FLUSH(w_unit)
430      ENDIF
431      DO i=1,nlon
432        DO j=1,nlat
433          ij=i+(j-1)*(nlon)
434          IF ( gg_mask(i,j) == 1 ) THEN
435              global_error(ij) = err_msk
436          ENDIF
437        ENDDO
438      ENDDO
439      max=MAXVAL(REAL(global_error))
440      IF (FILE_Debug >= 2) THEN
441          WRITE(w_unit,*)'Max and its location in the error field : ',max
442          WRITE(w_unit,*) MAXLOC(REAL(global_error))
443          CALL FLUSH(w_unit)
444      ENDIF
445      DO i=1,nlon
446        DO j=1,nlat
447          ij=i+(j-1)*(nlon)
448          IF ( gg_mask(i,j) == 1 ) THEN
449              global_error(ij) = 0.
450          ENDIF
451        ENDDO
452      ENDDO
453      IF (FILE_Debug >= 2) THEN
454          WRITE(w_unit,*) 'Number of masked points :',ic_msk
455          WRITE(w_unit,*) 'Error mean on non masked points: ', &
456             SUM(ABS(global_error))/((ntot)-ic_msk)
457          WRITE(w_unit,*) 'End calculation of stat on the error'
458          CALL FLUSH(w_unit)
459      ENDIF
460  ENDIF
461  !
462  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
463  !         TERMINATION
464  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
465  IF (FILE_Debug >= 2) THEN
466      WRITE (w_unit,*) 'End of the program, before oasis_terminate'
467      CALL FLUSH(w_unit)
468  ENDIF
469  !
470  !!!!!!!!!!!!!!!!!! OASIS_ENDDEF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
471  !
472  ! Collective call to terminate the coupling exchanges
473  !
474  CALL oasis_terminate (ierror)
475  IF (ierror /= 0) THEN
476      WRITE (w_unit,*) 'oasis_terminate abort by model2 compid ',comp_id
477      CALL oasis_abort(comp_id,comp_name,'Problem at line 477')
478  ENDIF
479  !
480  CALL mpi_finalize(ierror)
481  !
482END PROGRAM MODEL2
483!
Note: See TracBrowser for help on using the repository browser.