source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/spoc/spoc_regridding/model2.f90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago

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