source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/oasis3-mct/examples/test_interpolation/model1.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: 13.9 KB
Line 
1!------------------------------------------------------------------------
2! Copyright 2018/03, CERFACS, Toulouse, France.
3! All rights reserved. Use is subject to OASIS3 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 grid (P=periodic or R=regional)
30  INTEGER            :: il_overlap_src
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  CHARACTER(len=4)   :: cl_grd_tgt ! name of the target grid
38  NAMELIST /grid_target_characteristics/cl_grd_tgt
39  !
40  CHARACTER(len=14)  :: cl_msk_nam
41  !
42  ! Global grid parameters :
43  INTEGER :: nlon, nlat    ! dimensions in the 2 directions of space
44  INTEGER :: il_size
45  REAL (kind=wp), DIMENSION(:,:), POINTER  :: gg_lon,gg_lat ! lon, lat of the points
46  INTEGER, DIMENSION(:,:), POINTER         :: gg_mask ! mask, 0 == valid point, 1 == masked point
47  !
48  INTEGER :: mype, npes ! rank and  number of pe
49  INTEGER :: localComm  ! local MPI communicator and Initialized
50  INTEGER :: comp_id    ! component identification
51  !
52  INTEGER, DIMENSION(:), ALLOCATABLE :: il_paral ! Decomposition for each proc
53  !
54  INTEGER :: ierror, rank, w_unit
55  INTEGER :: FILE_Debug=2
56  !
57  ! Names of exchanged Fields
58  CHARACTER(len=8), PARAMETER :: var_name = 'FSENDANA' ! 8 characters field sent by model1
59  !
60  ! Used in oasis_def_var and oasis_def_var
61  INTEGER                       :: var_id
62  INTEGER                       :: var_nodims(2) 
63  INTEGER                       :: var_type
64  !
65  !
66  ! Grid parameters definition
67  INTEGER                       :: part_id  ! use to connect the partition to the variables
68  INTEGER                       :: var_sh(4) ! local dimensions of the arrays; 2 x rank (=4)
69  INTEGER :: ibeg, iend, jbeg, jend
70
71  !
72  ! Exchanged local fields arrays
73  REAL (kind=wp),   POINTER     :: field_send(:,:)
74  REAL (kind=wp),   POINTER     :: gradient_i(:,:), gradient_j(:,:), gradient_ij(:,:)
75  REAL (kind=wp),   POINTER     :: grad_lat(:,:), grad_lon(:,:)
76  !
77  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
78  !  INITIALISATION
79  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
80  !
81  CALL oasis_init_comp (comp_id, comp_name, ierror )
82  IF (ierror /= 0) THEN
83      WRITE(0,*) 'oasis_init_comp abort by model1 compid ',comp_id
84      CALL oasis_abort(comp_id,comp_name,'Problem at line 100')
85  ENDIF
86  !
87  ! Unit for output messages : one file for each process
88  CALL MPI_Comm_Rank ( MPI_COMM_WORLD, rank, ierror )
89  IF (ierror /= 0) THEN
90      WRITE(0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
91      CALL oasis_abort(comp_id,comp_name,'Problem at line 107')
92  ENDIF
93  !
94  !
95  !!!!!!!!!!!!!!!!! OASIS_GET_LOCALCOMM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96  !
97  CALL oasis_get_localcomm ( localComm, ierror )
98  IF (ierror /= 0) THEN
99      WRITE (0,*) 'oasis_get_localcomm abort by model1 compid ',comp_id
100      CALL oasis_abort(comp_id,comp_name,'Problem at line 116')
101  ENDIF
102  !
103  ! Get MPI size and rank
104  CALL MPI_Comm_Size ( localComm, npes, ierror )
105  IF (ierror /= 0) THEN
106      WRITE(0,*) 'MPI_comm_size abort by model1 compid ',comp_id
107      CALL oasis_abort(comp_id,comp_name,'Problem at line 123')
108  ENDIF
109  !
110  CALL MPI_Comm_Rank ( localComm, mype, ierror )
111  IF (ierror /= 0) THEN
112      WRITE (0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
113      CALL oasis_abort(comp_id,comp_name,'Problem at line 129')
114  ENDIF
115  !
116  IF ((FILE_Debug == 1) .AND. (mype == 0)) FILE_Debug=2
117  !
118  IF (FILE_Debug <= 1) THEN
119      IF (mype == 0) THEN
120          w_unit = 100 + rank
121          WRITE(chout,'(I3)') w_unit
122          comp_out=comp_name//'.root_'//chout
123          OPEN(w_unit,file=TRIM(comp_out),form='formatted')
124      ELSE
125          w_unit = 15
126          comp_out=comp_name//'.notroot'
127          OPEN(w_unit,file=TRIM(comp_out),form='formatted',position='append')
128      ENDIF
129  ELSE
130      w_unit = 100 + rank
131      WRITE(chout,'(I3)') w_unit
132      comp_out=comp_name//'.out_'//chout
133      OPEN(w_unit,file=TRIM(comp_out),form='formatted')
134  ENDIF
135  !
136  IF (FILE_Debug >= 2) THEN
137      WRITE(w_unit,*) '-----------------------------------------------------------'
138      WRITE(w_unit,*) TRIM(comp_name), ' running with reals compiled as kind ',wp
139      WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' global rank :',rank
140      WRITE(w_unit,*) '----------------------------------------------------------'
141      WRITE(w_unit,*) 'I am the ', TRIM(comp_name), ' ', 'component identifier', comp_id, 'local rank', mype
142      WRITE (w_unit,*) 'Number of processors :',npes
143      WRITE(w_unit,*) '----------------------------------------------------------'
144      CALL FLUSH(w_unit)
145  ENDIF
146  !
147  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
148  !  GRID DEFINITION
149  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
150  !
151  ! Reading global grids.nc and masks.nc netcdf files
152  ! Get arguments giving source grid acronym and field type
153  !
154  OPEN(UNIT=70,FILE='name_grids.dat',FORM='FORMATTED')
155  READ(UNIT=70,NML=grid_source_characteristics)
156  READ(UNIT=70,NML=grid_target_characteristics)
157  CLOSE(70)
158  !
159  IF (FILE_Debug >= 2) THEN
160      WRITE(w_unit,*) 'Source grid name : ',cl_grd_src
161      WRITE(w_unit,*) 'Remapping : ',cl_remap
162      WRITE(w_unit,*) 'Source grid type : ',cl_type_src
163      WRITE(w_unit,*) 'Source grid overlapping pts :',il_overlap_src
164      CALL flush(w_unit)
165  ENDIF
166  !
167  ! Reading dimensions of the global grid
168  CALL read_dimgrid(nlon,nlat,data_gridname,cl_grd_src,w_unit,FILE_Debug)
169  !
170  ! Allocate grid arrays
171  ALLOCATE(gg_lon(nlon,nlat), STAT=ierror )
172  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lon'
173  ALLOCATE(gg_lat(nlon,nlat), STAT=ierror )
174  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gg_lat'
175  ALLOCATE(gg_mask(nlon,nlat), STAT=ierror )
176  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating indice_mask'
177  !
178  ! Read global grid longitudes, latitudes, corners, mask
179  CALL read_grid(nlon,nlat, data_gridname, cl_grd_src, w_unit, FILE_Debug, gg_lon, gg_lat)
180  cl_msk_nam = TRIM(cl_grd_src//'_with_'//cl_grd_tgt)
181  IF (inquire_mask(data_maskname,cl_msk_nam,w_unit, FILE_Debug)) THEN
182     CALL read_mask(nlon,nlat, data_maskname, cl_msk_nam, w_unit, FILE_Debug, gg_mask)
183  ELSE
184     CALL read_mask(nlon,nlat, data_maskname, cl_grd_src, w_unit, FILE_Debug, gg_mask)
185  END IF
186  !
187  IF (FILE_Debug >= 2) THEN
188      WRITE(w_unit,*) 'After grid and mask reading'
189      CALL FLUSH(w_unit)
190  ENDIF
191  !
192  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193  !  PARTITION DEFINITION
194  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
195  !
196#ifdef DECOMP_APPLE
197  il_size = 3
198#elif defined DECOMP_BOX
199  il_size = 5
200#endif
201  ALLOCATE(il_paral(il_size))
202  IF (FILE_Debug >= 2) THEN
203      WRITE(w_unit,*) 'After allocate il_paral, il_size', il_size
204      CALL FLUSH(w_unit)
205  ENDIF
206  !
207  CALL decomp_def (il_paral, il_size, nlon, nlat, mype, npes, w_unit)
208  IF (FILE_Debug >= 2) THEN
209      WRITE(w_unit,*) 'After decomp_def, il_paral = ', il_paral(:)
210      CALL FLUSH(w_unit)
211  ENDIF
212  !
213  CALL oasis_def_partition (part_id, il_paral, ierror)
214  IF (FILE_Debug >= 2) THEN
215      WRITE(w_unit,*) 'After oasis_def_partition'
216      CALL FLUSH(w_unit)
217  ENDIF
218  !
219  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
220  !  COUPLING LOCAL FIELD DECLARATION 
221  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
222  !
223  var_nodims(1) = 2    ! Rank of the field array is 2
224  var_nodims(2) = 1    ! Bundles always 1 for OASIS3
225  var_type = OASIS_Real
226  !
227  var_sh(1) = 1
228  var_sh(2) = il_paral(3)
229  var_sh(3) = 1 
230#ifdef DECOMP_APPLE
231  var_sh(4) = 1
232#elif defined DECOMP_BOX
233  var_sh(4) = il_paral(4)
234#endif
235  !
236  ! Declaration of the field associated with the partition
237  CALL oasis_def_var (var_id, var_name, part_id, &
238     var_nodims, OASIS_Out, var_sh, var_type, ierror)
239  IF (ierror /= 0) THEN
240      WRITE(w_unit,*) 'oasis_def_var abort by model1 compid ',comp_id
241      CALL oasis_abort(comp_id,comp_name,'Problem at line 256')
242  ENDIF
243  IF (FILE_Debug >= 2) THEN
244      WRITE(w_unit,*) 'After def_var'
245      CALL FLUSH(w_unit)
246  ENDIF
247  !
248  DEALLOCATE(il_paral)
249  !
250  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251  !  TERMINATION OF DEFINITION PHASE
252  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
253  !
254  CALL oasis_enddef ( ierror )
255  IF (ierror /= 0) THEN
256      WRITE(w_unit,*) 'oasis_enddef abort by model1 compid ',comp_id
257      CALL oasis_abort(comp_id,comp_name,'Problem at line 272')
258  ENDIF
259  IF (FILE_Debug >= 2) THEN
260      WRITE(w_unit,*) 'After enddef'
261      CALL FLUSH(w_unit)
262  ENDIF
263  !
264  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
265  !  SEND ARRAYS
266  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
267  !
268  ! Allocate the fields send and received by the model1
269  !
270  ALLOCATE(field_send(nlon,nlat), STAT=ierror )
271  IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field1_send'
272  !
273  CALL function_ana(nlon, nlat, gg_lon, gg_lat, field_send)
274  !
275  ibeg=1 ; iend=nlon
276  jbeg=((nlat/npes)*mype)+1 
277  !
278  IF (mype .LT. npes - 1) THEN
279      jend = (nlat/npes)*(mype+1)
280  ELSE
281      jend = nlat 
282  ENDIF
283  !
284  ! Calculate the gradients and send them for bicubic only if grid is not gaussian reduced
285  IF (cl_remap == 'bicu') THEN
286     IF ( trim(cl_type_src) == 'LR') THEN
287        ALLOCATE(gradient_i(nlon,nlat), STAT=ierror )
288        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
289        ALLOCATE(gradient_j(nlon,nlat), STAT=ierror )
290        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
291        ALLOCATE(gradient_ij(nlon,nlat), STAT=ierror )
292        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_ij'
293        call gradient_bicubic(nlon, nlat, field_send, gg_mask, gg_lat, gg_lon, il_overlap_src,  &
294                                  cl_period_src, gradient_i, gradient_j, gradient_ij)
295        IF (FILE_Debug >= 2) THEN
296           WRITE(w_unit,*) 'Bicubic_gradient calculated '
297           CALL FLUSH(w_unit)
298        ENDIF
299        ! For BICUBIC, need to transfer three extra arguments :
300        call oasis_put(var_id, 0, &
301                       RESHAPE(field_send(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
302                       ierror, &
303                       RESHAPE(gradient_i(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
304                       RESHAPE(gradient_j(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
305                       RESHAPE(gradient_ij(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)) )
306     ELSE IF ( trim(cl_type_src) == 'D') THEN
307        call oasis_put(var_id, 0, &
308                       RESHAPE(field_send(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
309                       ierror )
310     ELSE
311        WRITE(w_unit,*) 'Cannot perform bicubic interpolation for type of grid ',cl_type_src
312        CALL oasis_abort(comp_id,comp_name,'Bicubic interpolation impossible for that grid')
313     ENDIF
314
315  ! Calculate the gradients and send them for conserv 2nd only if grid is not gaussian reduced
316  ELSE IF (cl_remap == 'conserv2nd') THEN
317     IF ( trim(cl_type_src) == 'LR') THEN
318        ALLOCATE(grad_lat(nlon,nlat), STAT=ierror )
319        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
320        ALLOCATE(grad_lon(nlon,nlat), STAT=ierror )
321        IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
322        call gradient_conserv(nlon, nlat, field_send, gg_mask, gg_lat, gg_lon, &
323                    & il_overlap_src, cl_period_src, grad_lat, grad_lon)
324        IF (FILE_Debug >= 2) THEN
325           WRITE(w_unit,*) 'Conservative gradient calculated '
326           CALL FLUSH(w_unit)
327        ENDIF
328        ! For CONSERV/SECOND, need to transfer two extra arguments :
329        call oasis_put(var_id, 0, &
330                       RESHAPE(field_send(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
331                       ierror, &
332                       RESHAPE(grad_lat(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
333                       RESHAPE(grad_lon(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)) )
334     ELSE
335        WRITE(w_unit,*) 'Cannot perform second order conserv interpolation for type of grid ',cl_type_src
336        CALL oasis_abort(comp_id,comp_name,'Second order conserv interpolation impossible for that grid')
337     ENDIF
338
339  ELSE
340     call oasis_put(var_id, 0, &
341                    RESHAPE(field_send(ibeg:iend,jbeg:jend),(/var_sh(2),var_sh(4)/)), &
342                    ierror )
343  ENDIF
344  !
345  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
346  !         TERMINATION
347  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
348  IF (FILE_Debug >= 2) THEN
349      WRITE(w_unit,*) 'End of the program, before oasis_terminate'
350      CALL FLUSH(w_unit)
351  ENDIF
352  !
353  CALL oasis_terminate (ierror)
354  IF (ierror /= 0) THEN
355      WRITE(w_unit,*) 'oasis_terminate abort by model1 compid ',comp_id
356      CALL oasis_abort(comp_id,comp_name,'Problem at line 332')
357  ENDIF
358  !
359  !
360END PROGRAM MODEL1
361!
Note: See TracBrowser for help on using the repository browser.