source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/oasis3-mct/examples/test_interpolation/model2.f90 @ 5725

Last change on this file since 5725 was 5725, checked in by aclsce, 3 years ago

Added new oasis3-MCT version to be used to handle ensembles simulations with XIOS.
File size: 14.3 KB
Line 
1!------------------------------------------------------------------------
2! Copyright 2018/03, CERFACS, Toulouse, France.
3! All rights reserved. Use is subject to OASIS3 license terms.
4!=============================================================================
5!
6PROGRAM model2
7  !
8  USE netcdf
9  USE mod_oasis
10  USE read_all_data
11  USE write_all_fields
12  !
13  IMPLICIT NONE
14  !
15  INCLUDE 'mpif.h'
16  !
17  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
18  !
19  CHARACTER(len=30), PARAMETER   :: data_gridname='grids.nc' ! file with the grids
20  CHARACTER(len=30), PARAMETER   :: data_maskname='masks.nc' ! file with the masks
21  CHARACTER(len=30)              :: data_filename, field_name
22  !
23  ! Component name (6 characters) same as in the namcouple
24  CHARACTER(len=6)   :: comp_name = 'model2'
25  CHARACTER(len=128) :: comp_out ! name of the output log file
26  CHARACTER(len=3)   :: chout
27  CHARACTER(len=4)   :: cl_grd_src     ! name of the source grid
28  CHARACTER(len=11)  :: cl_remap       ! type of remapping
29  CHARACTER(len=2)   :: cl_type_src    ! type of the source grid
30  CHARACTER(len=8)   :: cl_period_src  ! Periodicity of grid (P=periodic or R=regional)
31  INTEGER            :: il_overlap_src
32  NAMELIST /grid_source_characteristics/cl_grd_src
33  NAMELIST /grid_source_characteristics/cl_remap
34  NAMELIST /grid_source_characteristics/cl_type_src
35  NAMELIST /grid_source_characteristics/cl_period_src
36  NAMELIST /grid_source_characteristics/il_overlap_src
37  !
38  CHARACTER(len=4)   :: cl_grd_tgt ! name of the target grid
39  NAMELIST /grid_target_characteristics/cl_grd_tgt
40  !
41  CHARACTER(len=14)  :: cl_msk_nam
42  !
43  ! Global grid parameters :
44  INTEGER :: nlon, nlat    ! dimensions in the 2 directions of space
45  INTEGER :: il_size
46  INTEGER, PARAMETER :: echelle=1            ! To calculate th delta error for plot
47  REAL (kind=wp), DIMENSION(:,:), POINTER    :: gg_lon,gg_lat
48  INTEGER, DIMENSION(:,:), POINTER           :: gg_mask ! mask, 0 == valid point, 1 == masked point
49  !
50  INTEGER :: mype, npes ! rank and number of pe
51  INTEGER :: localComm  ! local MPI communicator and Initialized
52  INTEGER :: comp_id    ! component identification
53  !
54  INTEGER, DIMENSION(:), ALLOCATABLE :: il_paral ! Decomposition for each proc
55  !
56  INTEGER :: ierror, rank, w_unit
57  INTEGER :: ic_nmsk, ic_nmskrv
58  INTEGER :: FILE_Debug=2
59  !
60  ! Names of exchanged Fields
61  CHARACTER(len=8), PARAMETER :: var_name = 'FRECVANA' ! 8 characters field received by the atmospheremodel2 from model1
62  !
63  ! Used in oasis_def_var and oasis_def_var
64  INTEGER                       :: var_id 
65  INTEGER                       :: var_nodims(2) 
66  INTEGER                       :: var_type
67  !
68  REAL (kind=wp), PARAMETER     :: field_ini = -1. ! initialisation of received fields
69  !
70  ! Grid parameter definition
71  INTEGER                       :: part_id  ! use to connect the partition to the variables
72  INTEGER                       :: var_sh(4) ! local dimensions of the arrays; 2 x rank (=4)
73  INTEGER :: ibeg, iend, jbeg, jend
74  !
75  ! Local fields arrays used in routines oasis_put and oasis_get
76  REAL (kind=wp), POINTER       :: field_recv1d(:,:), field_recv(:,:), field_ana(:,:), gg_error(:,:)
77  INTEGER, POINTER              :: mask_error(:,:) ! error mask, 0 == masked point, 1 == valid point
78  !
79  ! Min and Max of the error of interpolation
80  REAL (kind=wp)             :: min,max
81  !
82  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
83  !   INITIALISATION
84  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85  !
86  CALL mpi_init(ierror)
87  !
88  CALL oasis_init_comp (comp_id, comp_name, ierror )
89  IF (ierror /= 0) THEN
90      WRITE(0,*) 'oasis_init_comp abort by model2 compid ',comp_id
91      CALL oasis_abort(comp_id,comp_name,'Problem at line 109')
92  ENDIF
93  !
94  ! Unit for output messages : one file for each process
95  CALL MPI_Comm_Rank ( MPI_COMM_WORLD, rank, ierror )
96  IF (ierror /= 0) THEN
97      WRITE(0,*) 'MPI_Comm_Rank abort by model2 compid ',comp_id
98      CALL oasis_abort(comp_id,comp_name,'Problem at line 116')
99  ENDIF
100  !
101  !
102  !!!!!!!!!!!!!!!!! OASIS_GET_LOCALCOMM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
103  !
104  CALL oasis_get_localcomm ( localComm, ierror )
105  IF (ierror /= 0) THEN
106      WRITE (0,*) 'oasis_get_localcomm abort by model2 compid ',comp_id
107      CALL oasis_abort(comp_id,comp_name,'Problem at line 125')
108  ENDIF
109  !
110  ! Get MPI size and rank
111  CALL MPI_Comm_Size ( localComm, npes, ierror )
112  IF (ierror /= 0) THEN
113      WRITE(0,*) 'MPI_comm_size abort by model2 compid ',comp_id
114      CALL oasis_abort(comp_id,comp_name,'Problem at line 132')
115  ENDIF
116  !
117  CALL MPI_Comm_Rank ( localComm, mype, ierror )
118  IF (ierror /= 0) THEN
119      WRITE (0,*) 'MPI_Comm_Rank abort by model2 compid ',comp_id
120      CALL oasis_abort(comp_id,comp_name,'Problem at line 138')
121  ENDIF
122  !
123  IF ((FILE_Debug == 1) .AND. (mype == 0)) FILE_Debug=2
124  !
125  IF (FILE_Debug <= 1) THEN
126      IF (mype == 0) THEN
127          w_unit = 100 + rank
128          WRITE(chout,'(I3)') w_unit
129          comp_out=comp_name//'.root_'//chout
130          OPEN(w_unit,file=TRIM(comp_out),form='formatted')
131      ELSE
132          w_unit = 15
133          comp_out=comp_name//'.notroot'
134          OPEN(w_unit,file=TRIM(comp_out),form='formatted',position='append')
135      ENDIF
136  ELSE
137      w_unit = 100 + rank
138      WRITE(chout,'(I3)') w_unit
139      comp_out=comp_name//'.out_'//chout
140      OPEN(w_unit,file=TRIM(comp_out),form='formatted')
141  ENDIF
142  !
143  IF (FILE_Debug >= 2) THEN
144      WRITE (w_unit,*) '-----------------------------------------------------------'
145      WRITE (w_unit,*) TRIM(comp_name), ' running with reals compiled as kind ',wp
146      WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' global rank :',rank
147      WRITE (w_unit,*) '----------------------------------------------------------'
148      WRITE(w_unit,*) 'I am the ', TRIM(comp_name), ' ', 'component identifier', comp_id, 'local rank', mype
149      WRITE (w_unit,*) 'Number of processors :',npes
150      WRITE (w_unit,*) '----------------------------------------------------------'
151      CALL FLUSH(w_unit)
152  ENDIF
153  !
154  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155  !  GRID DEFINITION
156  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
157  !
158  ! Reading global grids.nc and masks.nc netcdf files
159  ! Get arguments giving source grid acronym and field type
160  !
161  OPEN(UNIT=70,FILE='name_grids.dat',FORM='FORMATTED')
162  READ(UNIT=70,NML=grid_source_characteristics)
163  READ(UNIT=70,NML=grid_target_characteristics)
164  CLOSE(70)
165  !
166  IF (FILE_Debug >= 2) THEN
167      WRITE(w_unit,*) 'Target grid name : ',cl_grd_tgt
168      CALL flush(w_unit)
169  ENDIF
170  !
171  ! Reading dimensions of the grid
172  CALL read_dimgrid(nlon,nlat,data_gridname,cl_grd_tgt,w_unit,FILE_Debug)
173  !
174  ! Allocate grid arrays
175  ALLOCATE(gg_lon(nlon,nlat), STAT=ierror )
176  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lon'
177  ALLOCATE(gg_lat(nlon,nlat), STAT=ierror )
178  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lat'
179  ALLOCATE(gg_mask(nlon,nlat), STAT=ierror )
180  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating indice_mask'
181  ALLOCATE(gg_error(nlon,nlat),STAT=ierror )
182  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_error'
183  ALLOCATE(mask_error(nlon,nlat),STAT=ierror )
184  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating mask_error'
185  !
186  ! Read global grid longitudes, latitudes and mask
187  !
188  CALL read_grid(nlon,nlat, data_gridname, cl_grd_tgt, w_unit, FILE_Debug, gg_lon,gg_lat)
189  cl_msk_nam = TRIM(cl_grd_tgt//'_with_'//cl_grd_src)
190  IF (inquire_mask(data_maskname,cl_msk_nam,w_unit, FILE_Debug)) THEN
191     CALL read_mask(nlon,nlat, data_maskname, cl_msk_nam, w_unit, FILE_Debug, gg_mask)
192  ELSE
193     CALL read_mask(nlon,nlat, data_maskname, cl_grd_tgt, w_unit, FILE_Debug, gg_mask)
194  END IF
195  !
196  IF (FILE_Debug >= 2) THEN
197      WRITE(w_unit,*) 'After grid and mask reading'
198      CALL FLUSH(w_unit)
199  ENDIF
200  !
201  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
202  !  PARTITION DEFINITION
203  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
204  !
205  il_size = 3
206  ALLOCATE(il_paral(il_size))
207  IF (FILE_Debug >= 2) THEN
208      WRITE(w_unit,*) 'After allocate il_paral, il_size', il_size
209      CALL FLUSH(w_unit)
210  ENDIF
211  !
212  il_paral(1)=0
213  il_paral(2)=0
214  IF (mype == 0) THEN
215     il_paral(3)=nlon*nlat
216  ELSE
217     il_paral(3)=0
218  END IF
219  !       
220  IF (FILE_Debug >= 2) THEN
221      WRITE(w_unit,*) 'il_paral = ', il_paral(:)
222      CALL FLUSH(w_unit)
223  ENDIF
224  !       
225  CALL oasis_def_partition (part_id, il_paral, ierror)
226  IF (FILE_Debug >= 2) THEN
227      WRITE(w_unit,*) 'After oasis_def_partition'
228      CALL FLUSH(w_unit)
229  ENDIF
230  !
231  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
232  ! COUPLING LOCAL FIELD DECLARATION
233  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
234  !
235  var_nodims(1) = 2    ! Rank of the field array is 2
236  var_nodims(2) = 1    ! Bundles always 1 for OASIS3
237  var_type = OASIS_Real
238  !
239  var_sh(1) = 1
240  var_sh(2) = il_paral(3)
241  var_sh(3) = 1 
242  var_sh(4) = 1
243  !
244  CALL oasis_def_var (var_id,var_name, part_id, &
245     var_nodims, OASIS_In, var_sh, var_type, ierror)
246  IF (ierror /= 0) THEN
247      WRITE (w_unit,*) 'oasis_def_var abort by model2 compid ',comp_id
248      CALL oasis_abort(comp_id,comp_name,'Problem at line 266')
249  ENDIF
250  IF (FILE_Debug >= 2) THEN
251      WRITE(w_unit,*) 'After def_var'
252      CALL FLUSH(w_unit)
253  ENDIF
254  !
255  DEALLOCATE(il_paral)
256  !
257  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
258  !         TERMINATION OF DEFINITION PHASE
259  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
260  !
261  CALL oasis_enddef ( ierror )
262  IF (ierror /= 0) THEN
263      WRITE (w_unit,*) 'oasis_enddef abort by model2 compid ',comp_id
264      CALL oasis_abort(comp_id,comp_name,'Problem at line 281')
265  ENDIF
266  IF (FILE_Debug >= 2) THEN
267      WRITE(w_unit,*) 'After enddef'
268      CALL FLUSH(w_unit)
269  ENDIF
270  !
271IF (mype == 0) THEN
272  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273  ! RECEIVE ARRAYS AND CALCULATE THE ERROR
274  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
275  !
276  ! Allocate the fields send and received by the model
277  !
278  ALLOCATE(field_recv1d(var_sh(2), var_sh(4)), STAT=ierror )
279  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_recv1d'
280  ALLOCATE(field_recv(nlon,nlat), STAT=ierror )
281  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_recv'
282  ALLOCATE(field_ana(nlon,nlat),STAT=ierror )
283  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_ana'
284  !
285  CALL function_ana(nlon, nlat, gg_lon, gg_lat, field_ana)
286  !
287  ! Get the field FRECVANA
288  field_recv1d=field_ini
289  CALL oasis_get(var_id, 0, field_recv1d, ierror)
290  field_recv = RESHAPE(field_recv1d,(/nlon,nlat/))
291  DEALLOCATE(field_recv1d)
292  !
293  IF (ierror .NE. OASIS_Ok .AND. ierror .LT. OASIS_Recvd) THEN
294      WRITE (w_unit,*) 'oasis_get abort by model2 compid ',comp_id
295      CALL oasis_abort(comp_id,comp_name,'Problem at line 327')
296  ENDIF
297  IF (FILE_Debug >= 2) THEN
298      WRITE(w_unit,*) 'MINVAL(field_recv),MAXVAL(field_recv)=',MINVAL(field_recv),MAXVAL(field_recv)
299  ENDIF
300  !
301  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
302  ! Modify the field so to have:
303  !   - on masked points: value of 10000, error of -10000
304  !   - on non-masked points that did not receive any interpolated value: value of 1.e20, error of -1.e20
305  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306  !
307  gg_error=-10000
308  mask_error=1
309  WHERE (gg_mask == 1)     ! masked points
310        field_recv=10000.
311        mask_error=0
312  ELSEWHERE                     ! non-masked points
313        WHERE (field_recv /= 0.)   ! non-masked points that received an interpolated value
314           gg_error = ABS(((field_ana - field_recv)/field_recv))*100
315        ELSEWHERE   ! non-masked points that did not receive an interpolated value
316           gg_error=-1.e20
317           field_recv=1.e20
318           mask_error=0
319        END WHERE
320  END WHERE   
321  !
322  IF (FILE_Debug >= 2) THEN
323      WRITE (w_unit,*) 'After calculating the interpolation error'
324      CALL FLUSH(w_unit)
325  ENDIF
326  !
327  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
328  !  Write the error and the field in a NetCDF file by proc0
329  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
330  !
331  data_filename='error_interp.nc'
332  field_name='error_interp'
333  CALL write_field(nlon, nlat, data_filename, field_name, w_unit, FILE_Debug, gg_lon, gg_lat, gg_error)
334  !
335  data_filename='FRECVANA.nc'
336  field_name='FRECVANA' 
337  CALL write_field(nlon, nlat, data_filename, field_name, w_unit, FILE_Debug, gg_lon, gg_lat, field_recv)
338  !
339  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
340  ! Calculate error min and max on non-masked points that received an interpolated value
341  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
342  !
343  min=MINVAL(gg_error, MASK=mask_error>0)
344  IF (FILE_Debug >= 2) THEN
345     WRITE(w_unit,*) 'Min (%) and its location in the error field : ',min
346     WRITE(w_unit,*) MINLOC(gg_error)
347     CALL FLUSH(w_unit)
348  ENDIF
349  !
350  max=MAXVAL(gg_error, MASK=mask_error>0)
351  IF (FILE_Debug >= 2) THEN
352      WRITE(w_unit,*)'Max (%) and its location in the error field : ',max
353      WRITE(w_unit,*) MAXLOC(gg_error)
354      CALL FLUSH(w_unit)
355  ENDIF
356  !
357  IF (FILE_Debug >= 2) THEN
358      ic_nmsk=nlon*nlat-SUM(gg_mask)
359      WRITE(w_unit,*) 'Number of non-masked points :',ic_nmsk
360      ic_nmskrv=SUM(mask_error)
361      WRITE(w_unit,*) 'Number of non-masked points that received a value :',ic_nmskrv
362      WRITE(w_unit,*) 'Error mean on non masked points that received a value (%): ', SUM(ABS(gg_error), MASK=mask_error>0)/ic_nmskrv
363      WRITE(w_unit,*) 'Delta error (%/echelle) :',(max - min)/echelle
364      WRITE(w_unit,*) 'End calculation of stat on the error'
365      CALL FLUSH(w_unit)
366  ENDIF
367
368  CALL SYSTEM("ln -s "//TRIM(comp_out)//" model2.out")
369ENDIF ! of proc0
370  !
371  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
372  !         TERMINATION
373  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
374  IF (FILE_Debug >= 2) THEN
375      WRITE (w_unit,*) 'End of the program, before oasis_terminate'
376      CALL FLUSH(w_unit)
377  ENDIF
378  !
379  CALL oasis_terminate (ierror)
380  IF (ierror /= 0) THEN
381      WRITE (w_unit,*) 'oasis_terminate abort by model2 compid ',comp_id
382      CALL oasis_abort(comp_id,comp_name,'Problem at line 477')
383  ENDIF
384  !
385  CALL mpi_finalize(ierror)
386  !
387END PROGRAM MODEL2
388!
Note: See TracBrowser for help on using the repository browser.