source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/spoc/spoc_regridding/model1.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.5 KB
Line 
1!------------------------------------------------------------------------
2! Copyright 2021, CERFACS, Toulouse, France.
3! All rights reserved. Use is subject to OASIS3-MCT license terms.
4!=============================================================================
5!
6PROGRAM model1
7  !
8  USE netcdf
9  USE mod_oasis
10  USE read_all_data
11  !
12  IMPLICIT NONE
13  !
14  INCLUDE 'mpif.h'
15  !
16  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
17  !
18  CHARACTER(len=30), PARAMETER   :: data_gridname='grids.nc' ! file with the grids
19  CHARACTER(len=30), PARAMETER   :: data_maskname='masks.nc' ! file with the masks
20  CHARACTER(len=30)              :: data_filename, field_name
21  !
22  ! Component name (6 characters) same as in the namcouple
23  CHARACTER(len=6)   :: comp_name = 'model1'
24  CHARACTER(len=128) :: comp_out       ! name of the output log file
25  CHARACTER(len=3)   :: chout
26  CHARACTER(len=4)   :: cl_grd_src     ! name of the source grid
27  CHARACTER(len=11)  :: cl_remap       ! type of remapping
28  CHARACTER(len=2)   :: cl_type_src    ! type of the source grid
29  CHARACTER(len=8)   :: cl_period_src  ! periodicity of the source grid (P=periodic or R=regional)
30  INTEGER            :: il_overlap_src ! number of overlapping points
31  NAMELIST /grid_source_characteristics/cl_grd_src
32  NAMELIST /grid_source_characteristics/cl_remap
33  NAMELIST /grid_source_characteristics/cl_type_src
34  NAMELIST /grid_source_characteristics/cl_period_src
35  NAMELIST /grid_source_characteristics/il_overlap_src
36  !
37  ! Global grid parameters :
38  INTEGER :: nlon, nlat    ! dimensions in the 2 space directions
39  INTEGER :: il_size       ! dimension of the partition description array
40  REAL (kind=wp), DIMENSION(:,:), POINTER  :: gg_lon,gg_lat ! grid point longitudes and latitudes
41  INTEGER, DIMENSION(:,:), POINTER         :: gg_mask       ! grid point mask
42  !
43  INTEGER :: mype, npes ! MPI task rank and number
44  INTEGER :: local_comm  ! local MPI communicator
45  INTEGER :: comp_id    ! component identification
46  !
47  INTEGER, DIMENSION(:), ALLOCATABLE :: il_paral ! Decomposition for each proc
48  !
49  INTEGER :: ierror, rank, w_unit
50  LOGICAL :: file_debug = .true.
51  !
52  ! Names of exchanged Fields
53  CHARACTER(len=8), PARAMETER :: var_name = 'FSENDANA' ! 8 characters field sent by model1
54  !
55  ! Used in oasis_def_var and oasis_def_var
56  INTEGER                       :: var_id
57  INTEGER                       :: var_nodims(2) 
58  INTEGER                       :: var_type
59  !
60  ! Grid parameters definition
61  INTEGER                       :: part_id   ! use to connect the partition to the variables
62  INTEGER                       :: var_sh(1) ! not used anymore
63  INTEGER :: ibeg, iend, jbeg, jend
64  !
65  ! Exchanged local fields arrays
66  REAL (kind=wp),   POINTER     :: field_send(:,:)
67  REAL (kind=wp),   POINTER     :: gradient_i(:,:), gradient_j(:,:), gradient_ij(:,:)
68  REAL (kind=wp),   POINTER     :: grad_lat(:,:), grad_lon(:,:)
69  !
70  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
71  !  INITIALISATION
72  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73  !
74  CALL oasis_init_comp (comp_id, comp_name, ierror )
75  IF (ierror /= 0) THEN
76      WRITE(0,*) 'oasis_init_comp abort by model1 compid ',comp_id
77      CALL oasis_abort(comp_id,comp_name,'Problem in oasis_init_comp')
78  ENDIF
79  !
80  ! Unit for output messages : one file for each process
81  CALL MPI_Comm_Rank ( MPI_COMM_WORLD, rank, ierror )
82  IF (ierror /= 0) THEN
83      WRITE(0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
84      CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Rank')
85  ENDIF
86  !
87  !!!!!!!!!!!!!!!!! OASIS_GET_LOCALCOMM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88  !
89  CALL oasis_get_localcomm ( local_comm, ierror )
90  IF (ierror /= 0) THEN
91      WRITE (0,*) 'oasis_get_localcomm abort by model1 compid ',comp_id
92      CALL oasis_abort(comp_id,comp_name,'Problem in oasis_get_localcomm')
93  ENDIF
94  !
95  ! Get MPI size and rank
96  CALL MPI_Comm_Size ( local_comm, npes, ierror )
97  IF (ierror /= 0) THEN
98      WRITE(0,*) 'MPI_comm_size abort by model1 compid ',comp_id
99      CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Size')
100  ENDIF
101  !
102  CALL MPI_Comm_Rank ( local_comm, mype, ierror )
103  IF (ierror /= 0) THEN
104      WRITE (0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
105      CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Rank')
106  ENDIF
107  !
108  IF (file_debug) THEN
109      w_unit = 100 + rank
110      WRITE(chout,'(I3)') w_unit
111      comp_out=comp_name//'.out_'//chout
112      OPEN(w_unit,file=TRIM(comp_out),form='formatted')
113  ENDIF
114  !
115  IF (file_debug) THEN
116      WRITE(w_unit,*) '-----------------------------------------------------------'
117      WRITE(w_unit,*) TRIM(comp_name), ' running with reals compiled as kind ',wp
118      WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' global rank :',rank
119      WRITE(w_unit,*) '----------------------------------------------------------'
120      WRITE(w_unit,*) 'I am the ', TRIM(comp_name), ' ', 'component identifier', comp_id, 'local rank', mype
121      WRITE (w_unit,*) 'Number of processors :',npes
122      WRITE(w_unit,*) '----------------------------------------------------------'
123      CALL FLUSH(w_unit)
124  ENDIF
125  !
126  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
127  !  GRID DEFINITION
128  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
129  !
130  ! For simplicity, all processes read the whole global grid and mask files grids.nc and masks.nc
131  !
132  ! Get arguments giving source grid acronym and characteristics
133  OPEN(UNIT=70,FILE='name_grids.dat',FORM='FORMATTED')
134  READ(UNIT=70,NML=grid_source_characteristics)
135  CLOSE(70)
136  !
137  IF (file_debug) THEN
138      WRITE(w_unit,*) 'Source grid name : ',cl_grd_src
139      WRITE(w_unit,*) 'Remapping : ',cl_remap
140      WRITE(w_unit,*) 'Source grid type : ',cl_type_src
141      WRITE(w_unit,*) 'Source grid overlapping pts :',il_overlap_src
142      CALL flush(w_unit)
143  ENDIF
144  !
145  ! Read dimensions of the global grid
146  CALL read_dimgrid(nlon,nlat,data_gridname,cl_grd_src,w_unit,file_debug)
147  !
148  ! Allocate grid arrays
149  ALLOCATE(gg_lon(nlon,nlat), STAT=ierror )
150  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lon'
151  ALLOCATE(gg_lat(nlon,nlat), STAT=ierror )
152  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lat'
153  ALLOCATE(gg_mask(nlon,nlat), STAT=ierror )
154  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating indice_mask'
155  !
156  ! Read global grid longitudes, latitudes and mask
157  CALL read_grid(nlon,nlat, data_gridname, cl_grd_src, w_unit, file_debug, gg_lon, gg_lat)
158  CALL read_mask(nlon,nlat, data_maskname, cl_grd_src, w_unit, file_debug, gg_mask)
159  !
160  IF (file_debug) THEN
161      WRITE(w_unit,*) 'After grid and mask reading'
162      CALL FLUSH(w_unit)
163  ENDIF
164  !
165  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
166  !  PARTITION DEFINITION (each process treats a local part of the coupling field)
167  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
168  !
169#ifdef DECOMP_APPLE
170  il_size = 3
171#elif defined DECOMP_BOX
172  il_size = 5
173#endif
174  ALLOCATE(il_paral(il_size))
175  IF (file_debug) THEN
176      WRITE(w_unit,*) 'After allocate il_paral, il_size', il_size
177      CALL FLUSH(w_unit)
178  ENDIF
179  !
180  CALL decomp_def (il_paral, il_size, nlon, nlat, mype, npes, w_unit)
181  IF (file_debug) THEN
182      WRITE(w_unit,*) 'After decomp_def, il_paral = ', il_paral(:)
183      CALL FLUSH(w_unit)
184  ENDIF
185  !
186  CALL oasis_def_partition (part_id, il_paral, ierror)
187  IF (file_debug) THEN
188      WRITE(w_unit,*) 'After oasis_def_partition'
189      CALL FLUSH(w_unit)
190  ENDIF
191  !
192  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193  !  COUPLING FIELD DECLARATION 
194  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
195  !
196  var_nodims(1) = 2    ! Rank of the field array is 2 (not used anymore)
197  var_nodims(2) = 1    ! Number of bundle fields
198  var_sh(1) = 1        ! (not used anymore)
199  var_type = OASIS_Real
200  !
201  ! Declaration of the field associated with the partition
202  CALL oasis_def_var (var_id, var_name, part_id, &
203     var_nodims, OASIS_Out, var_sh, var_type, ierror)
204  IF (ierror /= 0) THEN
205      WRITE(w_unit,*) 'oasis_def_var abort by model1 compid ',comp_id
206      CALL oasis_abort(comp_id,comp_name,'Problem in oasis_def_var')
207  ENDIF
208  IF (file_debug) THEN
209      WRITE(w_unit,*) 'After oasis_def_var'
210      CALL FLUSH(w_unit)
211  ENDIF
212  !
213  DEALLOCATE(il_paral)
214  !
215  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
216  !  TERMINATION OF DEFINITION PHASE
217  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
218  !
219  CALL oasis_enddef ( ierror )
220  IF (ierror /= 0) THEN
221      WRITE(w_unit,*) 'oasis_enddef abort by model1 compid ',comp_id
222      CALL oasis_abort(comp_id,comp_name,'Problem in oasis_enddef')
223  ENDIF
224  IF (file_debug) THEN
225      WRITE(w_unit,*) 'After oasis_enddef'
226      CALL FLUSH(w_unit)
227  ENDIF
228  !
229  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
230  !  SEND ARRAYS
231  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
232  !
233  ! For simplicity, all processes allocate and define the whole global field
234  !
235  ! Allocate the field sent by model1
236  ALLOCATE(field_send(nlon,nlat), STAT=ierror )
237  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field1_send'
238  !
239  CALL function_ana(nlon, nlat, gg_lon, gg_lat, field_send)
240  !
241  ! Define indices corresponding to the local part of the coupling field
242  ibeg=1 ; iend=nlon
243  jbeg=((nlat/npes)*mype)+1 
244  !
245  IF (mype .LT. npes - 1) THEN
246      jend = (nlat/npes)*(mype+1)
247  ELSE
248      jend = nlat 
249  ENDIF
250  WRITE(w_unit,*) 'ibeg, iend, jbeg, jend', ibeg, iend, jbeg, jend
251  !
252  ! Special treament for bicubic remapping
253  IF (cl_remap == 'bicu') THEN
254     IF ( trim(cl_type_src) == 'LR') THEN
255        ! Calculate the gradients in i, j and ij needed for the bicubic remapping for LR grids
256        ! For simplicity, all processes calculate gradients on the whole grid
257        ALLOCATE(gradient_i(nlon,nlat), STAT=ierror )
258        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
259        ALLOCATE(gradient_j(nlon,nlat), STAT=ierror )
260        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
261        ALLOCATE(gradient_ij(nlon,nlat), STAT=ierror )
262        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_ij'
263        call gradient_bicubic(nlon, nlat, field_send, gg_mask, gg_lat, gg_lon, il_overlap_src,  &
264                                  cl_period_src, gradient_i, gradient_j, gradient_ij)
265        IF (file_debug) THEN
266           WRITE(w_unit,*) 'Bicubic_gradient calculated '
267           CALL FLUSH(w_unit)
268        ENDIF
269        ! Send the local part of the coupling field and gradients
270        call oasis_put(var_id, 0, &
271                       field_send(ibeg:iend,jbeg:jend), &
272                       ierror, &
273                       gradient_i(ibeg:iend,jbeg:jend), &
274                       gradient_j(ibeg:iend,jbeg:jend), &
275                       gradient_ij(ibeg:iend,jbeg:jend) )
276     ELSE IF ( trim(cl_type_src) == 'D') THEN
277        ! For Gaussian Reduced grids, a 16-point algorithm is used so gradients
278        ! are not needed; send only the local part of the coupling field
279        call oasis_put(var_id, 0, &
280                       field_send(ibeg:iend,jbeg:jend), &
281                       ierror )
282     ELSE
283        ! Bicubic remapping is not possible for othe grid types
284        WRITE(w_unit,*) 'Cannot perform bicubic interpolation for type of grid ',cl_type_src
285        CALL oasis_abort(comp_id,comp_name,'Bicubic interpolation impossible for that grid')
286     ENDIF
287  !
288  ! Special treament for 2nd order conservative remapping
289  ELSE IF (cl_remap == 'conserv2nd') THEN
290     IF ( trim(cl_type_src) == 'LR') THEN
291        ! Calculate the gradients in lat and lon directions needed for 2nd order
292        ! conservative remapping for LR grids
293        ! For simplicity, all processes calculate gradients on the whole grid
294        ALLOCATE(grad_lat(nlon,nlat), STAT=ierror )
295        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
296        ALLOCATE(grad_lon(nlon,nlat), STAT=ierror )
297        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
298        call gradient_conserv(nlon, nlat, field_send, gg_mask, gg_lat, gg_lon, &
299                    & il_overlap_src, cl_period_src, grad_lat, grad_lon)
300        IF (file_debug) THEN
301           WRITE(w_unit,*) 'Conservative gradient calculated '
302           CALL FLUSH(w_unit)
303        ENDIF
304        ! Send the local part of the coupling field and gradients
305        call oasis_put(var_id, 0, &
306                       field_send(ibeg:iend,jbeg:jend), &
307                       ierror, &
308                       grad_lat(ibeg:iend,jbeg:jend), &
309                       grad_lon(ibeg:iend,jbeg:jend) )
310     ELSE
311        ! 2nd order conservative not implemented for grids other than LR
312        WRITE(w_unit,*) 'Cannot perform second order conserv interpolation for type of grid ',cl_type_src
313        CALL oasis_abort(comp_id,comp_name,'Second order conserv interpolation impossible for that grid')
314     ENDIF
315  ! Standard oasis_put for other types of remappings
316  ELSE
317     call oasis_put(var_id, 0, &
318                    field_send(ibeg:iend,jbeg:jend),&
319                    ierror )
320  ENDIF
321  !
322  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
323  !         TERMINATION
324  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
325  IF (file_debug) THEN
326      WRITE(w_unit,*) 'End of the program, before oasis_terminate'
327      CALL FLUSH(w_unit)
328  ENDIF
329  !
330  CALL oasis_terminate (ierror)
331  IF (ierror /= 0) THEN
332      WRITE(w_unit,*) 'oasis_terminate abort by model1 compid ',comp_id
333  ENDIF
334  !
335  !
336END PROGRAM MODEL1
337!
Note: See TracBrowser for help on using the repository browser.