Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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: 26.0 KB

Line
1	!------------------------------------------------------------------------
2	! Copyright 2021, CERFACS, Toulouse, France.
3	! All rights reserved. Use is subject to OASIS3-MCT license terms.
4	!=============================================================================
5	!
6	PROGRAM model1
7	!
8	USE netcdf
9	USE mod_oasis
10	USE read_all_data
11	USE write_all_fields
12	USE function_ana
13	USE def_parallel_decomposition
14	!
15	IMPLICIT NONE
16	!
17	INCLUDE 'mpif.h'
18	!
19	INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
20	!
21	CHARACTER(len=30), PARAMETER :: data_gridname='grids.nc' ! file with the grids
22	CHARACTER(len=30), PARAMETER :: data_maskname='masks.nc' ! file with the masks
23	CHARACTER(len=30) :: data_filename, field_name
24	!
25	! Component name (6 characters) same as in the namcouple
26	CHARACTER(len=6) :: comp_name = 'model1'
27	CHARACTER(len=128) :: comp_out ! name of the output log file
28	CHARACTER(len=3) :: chout
29	CHARACTER(len=4) :: cl_grd_src, cl_grd_tgt ! name of the source grid
30	CHARACTER(len=11) :: cl_remap ! type of remapping
31	CHARACTER(len=2) :: cl_type_src ! type of the source grid
32	CHARACTER(len=2) :: cl_type_tgt ! type of the target grid
33	CHARACTER(len=8) :: cl_period_src ! periodicity of the source grid (P=periodic or R=regional)
34	INTEGER :: il_overlap_src ! number of overlapping points
35	CHARACTER(len=4) :: cl_library ! remapper (SCRP, ESMF or XIOS)
36	NAMELIST /grid_source_characteristics/cl_grd_src
37	NAMELIST /grid_source_characteristics/cl_remap
38	NAMELIST /grid_source_characteristics/cl_type_src
39	NAMELIST /grid_source_characteristics/cl_period_src
40	NAMELIST /grid_source_characteristics/il_overlap_src
41	NAMELIST /grid_target_characteristics/cl_grd_tgt
42	NAMELIST /grid_target_characteristics/cl_type_tgt
43	NAMELIST /remapper/cl_library
44	!
45	! Grid parameters
46	INTEGER :: il_extentx_s, il_extenty_s, il_offsetx_s, il_offsety_s
47	INTEGER :: il_extentx_t, il_extenty_t, il_offsetx_t, il_offsety_t
48	INTEGER :: il_size_s, il_offset_s
49	INTEGER :: il_size_t, il_offset_t
50	INTEGER :: nlon_s, nlat_s ! dimensions in the 2 space directions
51	INTEGER :: nlon_t, nlat_t ! dimensions in the 2 space directions
52	REAL (kind=wp), ALLOCATABLE :: grid_lon_s(:,:), grid_lat_s(:,:) ! lon, lat of the cell centers
53	REAL (kind=wp), ALLOCATABLE :: grid_lon_t(:,:), grid_lat_t(:,:) ! lon, lat of the cell centers
54	INTEGER, ALLOCATABLE :: grid_msk_s(:,:) ! mask, 0 == valid point, 1 == masked point
55	INTEGER, ALLOCATABLE :: grid_msk_t(:,:), grid_msk_t_global(:,:) ! mask, 0 == valid point, 1 == masked point
56	!
57	INTEGER :: mype, npes ! MPI task rank and number
58	INTEGER :: local_comm ! local MPI communicator
59	INTEGER :: comp_id ! component identification
60	!
61	INTEGER :: il_part_id_s, il_part_id_t
62	INTEGER, DIMENSION(3) :: ig_paral
63	!
64	INTEGER :: ierror, rank, w_unit
65	LOGICAL :: file_debug = .true.
66	!
67	! Names of exchanged Fields
68	CHARACTER(len=8), PARAMETER :: var_name_s = 'FSENDANA' ! 8 characters field sent
69	CHARACTER(len=8), PARAMETER :: var_name_t = 'FRECVANA' ! 8 characters field received
70	!
71	! Used in oasis_def_var and oasis_def_var
72	INTEGER :: var_id_s, var_id_t
73	INTEGER :: var_nodims(2)
74	INTEGER :: var_type
75	!
76	! Grid parameters definition
77	INTEGER :: var_sh(1) ! not used anymore
78	!
79	! Exchanged local fields arrays
80	REAL (kind=wp), ALLOCATABLE :: field_send(:,:)
81	REAL (kind=wp), ALLOCATABLE :: gradient_i(:,:), gradient_j(:,:), gradient_ij(:,:)
82	REAL (kind=wp), ALLOCATABLE :: grad_lat(:,:), grad_lon(:,:)
83	REAL (kind=wp), ALLOCATABLE :: field_recv(:,:), field_recv_global(:,:)
84	REAL (kind=wp), ALLOCATABLE :: grid_lon_global_t(:,:), grid_lat_global_t(:,:)
85	REAL (kind=wp), ALLOCATABLE :: field_ana(:,:), field_error(:,:), field_error_global(:,:)
86	INTEGER, ALLOCATABLE :: mask_error(:,:), mask_error_global(:,:) ! error mask, 0 == masked point, 1 == valid point
87	INTEGER, ALLOCATABLE :: il_size_all(:), il_offset_all(:)
88	!
89	INTEGER :: ic_nmsk, ic_nmskrv
90	REAL (kind=wp) :: min, max ! Minimum and maximum of the error
91	INTEGER, PARAMETER :: echelle = 1 ! To calculate the delta error for plot
92	!
93	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
94	! INITIALISATION
95	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96	!
97	CALL oasis_init_comp (comp_id, comp_name, ierror )
98	IF (ierror /= 0) THEN
99	WRITE(0,*) 'oasis_init_comp abort by model1 compid ',comp_id
100	CALL oasis_abort(comp_id,comp_name,'Problem in oasis_init_comp')
101	ENDIF
102	!
103	! Unit for output messages : one file for each process
104	CALL MPI_Comm_Rank ( MPI_COMM_WORLD, rank, ierror )
105	IF (ierror /= 0) THEN
106	WRITE(0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
107	CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Rank')
108	ENDIF
109	!
110	!!!!!!!!!!!!!!!!! OASIS_GET_LOCALCOMM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111	!
112	CALL oasis_get_localcomm ( local_comm, ierror )
113	IF (ierror /= 0) THEN
114	WRITE (0,*) 'oasis_get_localcomm abort by model1 compid ',comp_id
115	CALL oasis_abort(comp_id,comp_name,'Problem in oasis_get_localcomm')
116	ENDIF
117	!
118	! Get MPI size and rank
119	CALL MPI_Comm_Size ( local_comm, npes, ierror )
120	IF (ierror /= 0) THEN
121	WRITE(0,*) 'MPI_comm_size abort by model1 compid ',comp_id
122	CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Size')
123	ENDIF
124	!
125	CALL MPI_Comm_Rank ( local_comm, mype, ierror )
126	IF (ierror /= 0) THEN
127	WRITE (0,*) 'MPI_Comm_Rank abort by model1 compid ',comp_id
128	CALL oasis_abort(comp_id,comp_name,'Problem in MPI_Comm_Rank')
129	ENDIF
130	!
131	IF (file_debug) THEN
132	w_unit = 100 + rank
133	WRITE(chout,'(I3)') w_unit
134	comp_out=comp_name//'.out_'//chout
135	OPEN(w_unit,file=TRIM(comp_out),form='formatted')
136	ENDIF
137	!
138	IF (file_debug) THEN
139	WRITE(w_unit,*) '-----------------------------------------------------------'
140	WRITE(w_unit,*) TRIM(comp_name), ' running with reals compiled as kind ',wp
141	WRITE (w_unit,*) 'I am component ', TRIM(comp_name), ' global rank :',rank
142	WRITE(w_unit,*) '----------------------------------------------------------'
143	WRITE(w_unit,*) 'I am the ', TRIM(comp_name), ' ', 'component identifier', comp_id, 'local rank', mype
144	WRITE (w_unit,*) 'Number of processors :',npes
145	WRITE(w_unit,*) '----------------------------------------------------------'
146	CALL FLUSH(w_unit)
147	ENDIF
148	!
149	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
150	! GRID ACRONYMS and DIMENSIONS
151	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
152	!
153	! Get arguments giving source grid acronym and characteristics
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	READ(UNIT=70,NML=remapper)
158	CLOSE(70)
159	!
160	IF (file_debug) THEN
161	WRITE(w_unit,*) 'Source and target grid name : ',cl_grd_src, cl_grd_tgt
162	WRITE(w_unit,*) 'Remapping : ',cl_remap
163	WRITE(w_unit,*) 'Source and target grid type : ',cl_type_src, cl_type_tgt
164	WRITE(w_unit,*) 'Source grid overlapping pts :',il_overlap_src
165	WRITE(w_unit,*) 'Remapping library :', cl_library
166	CALL flush(w_unit)
167	ENDIF
168	!
169	! Read dimensions of the source and target global grids
170	CALL read_dimgrid(nlon_s, nlat_s, cl_grd_src, w_unit, file_debug)
171	CALL read_dimgrid(nlon_t, nlat_t, cl_grd_tgt, w_unit, file_debug)
172	!
173	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174	! SOURCE PARTITION DEFINITION
175	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
176	!
177	! Definition of the local partition
178	call def_local_partition(nlon_s, nlat_s, npes, mype, cl_type_src, &
179	il_extentx_s, il_extenty_s, il_size_s, il_offsetx_s, il_offsety_s, il_offset_s)
180	WRITE(w_unit,*) 'Local partition definition'
181	WRITE(w_unit,*) 'il_extentx_s, il_extenty_s, il_size_s, il_offsetx_s, il_offsety_s, il_offset_s = ', &
182	il_extentx_s, il_extenty_s, il_size_s, il_offsetx_s, il_offsety_s, il_offset_s
183	call flush(w_unit)
184	!
185	! APPLE PARTITION
186	ig_paral(1) = 1
187	ig_paral(2) = il_offset_s
188	ig_paral(3) = il_size_s
189	!
190	WRITE(w_unit,*) 'ig_paral = ', ig_paral(:)
191	call flush(w_unit)
192	!
193	CALL oasis_def_partition (il_part_id_s, ig_paral, ierror)
194	!
195	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
196	! TARGET PARTITION DEFINITION
197	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
198	!
199	! Definition of the local partition
200	call def_local_partition(nlon_t, nlat_t, npes, mype, cl_type_tgt, &
201	il_extentx_t, il_extenty_t, il_size_t, il_offsetx_t, il_offsety_t, il_offset_t)
202	WRITE(w_unit,*) 'Local partition definition'
203	WRITE(w_unit,*) 'il_extentx_t, il_extenty_t, il_size_t, il_offsetx_t, il_offsety_t, il_offset_t = ', &
204	il_extentx_t, il_extenty_t, il_size_t, il_offsetx_t, il_offsety_t, il_offset_t
205	call flush(w_unit)
206	!
207	! APPLE PARTITION
208	ig_paral(1) = 1
209	ig_paral(2) = il_offset_t
210	ig_paral(3) = il_size_t
211	!
212	WRITE(w_unit,*) 'ig_paral = ', ig_paral(:)
213	call flush(w_unit)
214	!
215	CALL oasis_def_partition (il_part_id_t, ig_paral, ierror)
216	!
217	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
218	! SOURCE GRID DEFINITION
219	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
220	! Allocation of local grid arrays
221	ALLOCATE(grid_lon_s(il_extentx_s, il_extenty_s), STAT=ierror )
222	ALLOCATE(grid_lat_s(il_extentx_s, il_extenty_s), STAT=ierror )
223	ALLOCATE(grid_msk_s(il_extentx_s, il_extenty_s), STAT=ierror )
224	!
225	! Reading local grid arrays from input file ocean_mesh.nc
226	WRITE(w_unit,*) 'Before read_grid, nlon_s, nlat_s', nlon_s, nlat_s
227	call flush(w_unit)
228	!
229	CALL read_grid(nlon_s, nlat_s, il_offsetx_s+1, il_offsety_s+1, il_extentx_s, il_extenty_s, &
230	cl_grd_src, w_unit, grid_lon_s, grid_lat_s, file_debug)
231	WRITE(w_unit,*) 'After read_grid, nlon_s, nlat_s', nlon_s, nlat_s
232	CALL read_mask(nlon_s, nlat_s, il_offsetx_s+1, il_offsety_s+1, il_extentx_s, il_extenty_s, &
233	cl_grd_src, w_unit, grid_msk_s, file_debug)
234	WRITE(w_unit,*) 'After read_mask, nlon_s, nlat_s', nlon_s, nlat_s
235	call flush(w_unit)
236	!
237	IF (file_debug) THEN
238	WRITE(w_unit,*) 'After grid and mask reading'
239	CALL FLUSH(w_unit)
240	ENDIF
241	!
242	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243	! TARGET GRID DEFINITION
244	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
245	! Allocation of local grid arrays
246	ALLOCATE(grid_lon_t(il_extentx_t, il_extenty_t), STAT=ierror )
247	ALLOCATE(grid_lat_t(il_extentx_t, il_extenty_t), STAT=ierror )
248	ALLOCATE(grid_msk_t(il_extentx_t, il_extenty_t), STAT=ierror )
249	!
250	! Reading local grid arrays from input file ocean_mesh.nc
251	WRITE(w_unit,*) 'Before read_grid, nlon_t, nlat_t', nlon_t, nlat_t
252	call flush(w_unit)
253	!
254	CALL read_grid(nlon_t, nlat_t, il_offsetx_t+1, il_offsety_t+1, il_extentx_t, il_extenty_t, &
255	cl_grd_tgt, w_unit, grid_lon_t, grid_lat_t, file_debug)
256	WRITE(w_unit,*) 'After read_grid, nlon_t, nlat_t', nlon_t, nlat_t
257	CALL read_mask(nlon_t, nlat_t, il_offsetx_t+1, il_offsety_t+1, il_extentx_t, il_extenty_t, &
258	cl_grd_tgt, w_unit, grid_msk_t, file_debug)
259	WRITE(w_unit,*) 'After read_mask, nlon_t, nlat_t', nlon_t, nlat_t
260	call flush(w_unit)
261	!
262	IF (file_debug) THEN
263	WRITE(w_unit,*) 'After grid and mask reading'
264	CALL FLUSH(w_unit)
265	ENDIF
266	!
267	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
268	! COUPLING FIELD DECLARATION
269	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
270	!
271	var_nodims(1) = 2 ! Rank of the field array is 2 (not used anymore)
272	var_nodims(2) = 1 ! Number of bundle fields
273	var_sh(1) = 1 ! (not used anymore)
274	var_type = OASIS_Real
275	!
276	! Declaration of the field associated with the source partition
277	CALL oasis_def_var (var_id_s, var_name_s, il_part_id_s, &
278	var_nodims, OASIS_Out, var_sh, var_type, ierror)
279	!
280	! Declaration of the field associated with the target partition
281	CALL oasis_def_var (var_id_t, var_name_t, il_part_id_t, &
282	var_nodims, OASIS_In, var_sh, var_type, ierror)
283	!
284	IF (ierror /= 0) THEN
285	WRITE(w_unit,*) 'oasis_def_var abort by model1 compid ',comp_id
286	CALL oasis_abort(comp_id,comp_name,'Problem in oasis_def_var')
287	ENDIF
288	IF (file_debug) THEN
289	WRITE(w_unit,*) 'After oasis_def_var'
290	CALL FLUSH(w_unit)
291	ENDIF
292	!
293	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
294	! TERMINATION OF DEFINITION PHASE
295	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
296	!
297	CALL oasis_enddef ( ierror )
298	IF (ierror /= 0) THEN
299	WRITE(w_unit,*) 'oasis_enddef abort by model1 compid ',comp_id
300	CALL oasis_abort(comp_id,comp_name,'Problem in oasis_enddef')
301	ENDIF
302	IF (file_debug) THEN
303	WRITE(w_unit,*) 'After oasis_enddef'
304	CALL FLUSH(w_unit)
305	ENDIF
306	!
307	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
308	! SEND ARRAYS
309	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
310	!
311	! For simplicity, all processes allocate and define the whole global field
312	!
313	! Allocate the field sent by model1
314	ALLOCATE(field_send(il_extentx_s, il_extenty_s), STAT=ierror )
315	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field1_send'
316	!
317	#ifdef FANA1
318	CALL function_ana1(il_extentx_s, il_extenty_s, grid_lon_s, grid_lat_s, field_send)
319	#elif defined FANA2
320	CALL function_ana2(il_extentx_s, il_extenty_s, grid_lon_s, grid_lat_s, field_send)
321	#elif defined FANA3
322	CALL function_vortex(il_extentx_s, il_extenty_s, grid_lon_s, grid_lat_s, field_send)
323	#endif
324	!
325
326	IF (cl_library == 'SCRP') THEN
327	! Special treament for bicubic remapping
328	IF (cl_remap == 'bicu') THEN
329	IF ( trim(cl_type_src) == 'LR') THEN
330	call flush (w_unit)
331	! Calculate the gradients in i, j and ij needed for the bicubic remapping for LR grids
332	! For simplicity, all processes calculate gradients on the whole grid
333	ALLOCATE(gradient_i(il_extentx_s,il_extenty_s), STAT=ierror )
334	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
335	call flush (w_unit)
336	ALLOCATE(gradient_j(il_extentx_s,il_extenty_s), STAT=ierror )
337	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
338	call flush (w_unit)
339	ALLOCATE(gradient_ij(il_extentx_s,il_extenty_s), STAT=ierror )
340	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_ij'
341	call flush (w_unit)
342	WRITE(w_unit,*) 'After allocate gradients'
343	call flush (w_unit)
344	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_ij'
345	call gradient_bicubic(nlon_s, nlat_s, il_offsetx_s+1, il_offsety_s+1, il_extentx_s, il_extenty_s, &
346	cl_grd_src, il_overlap_src, cl_period_src, w_unit, &
347	gradient_i, gradient_j, gradient_ij, file_debug)
348	WRITE(w_unit,*) 'After gradient_bicubic'
349	call flush (w_unit)
350	IF (file_debug) THEN
351	WRITE(w_unit,*) 'Bicubic_gradient calculated '
352	CALL FLUSH(w_unit)
353	ENDIF
354	! Send the local part of the coupling field and gradients
355	call oasis_put(var_id_s, 0, field_send, ierror, &
356	gradient_i, gradient_j, gradient_ij)
357	ELSE IF ( trim(cl_type_src) == 'D') THEN
358	! For Gaussian Reduced grids, a 16-point algorithm is used so gradients
359	! are not needed; send only the local part of the coupling field
360	call oasis_put(var_id_s, 0, field_send, ierror)
361	ELSE
362	! Bicubic remapping is not possible for othe grid types
363	WRITE(w_unit,*) 'Cannot perform bicubic interpolation for type of grid ',cl_type_src
364	CALL oasis_abort(comp_id,comp_name,'Bicubic interpolation impossible for that grid')
365	ENDIF
366	!
367	! Special treament for 2nd order conservative remapping
368	ELSE IF (cl_remap == 'conserv2nd') THEN
369	IF ( trim(cl_type_src) == 'LR') THEN
370	! Calculate the gradients in lat and lon directions needed for 2nd order
371	! conservative remapping for LR grids
372	! For simplicity, all processes calculate gradients on the whole grid
373	ALLOCATE(grad_lat(il_extentx_s,il_extenty_s), STAT=ierror )
374	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_i'
375	ALLOCATE(grad_lon(il_extentx_s,il_extenty_s), STAT=ierror )
376	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating gradient_j'
377	call gradient_conserv(nlon_s, nlat_s, il_offsetx_s+1, il_offsety_s+1, il_extentx_s, il_extenty_s, &
378	cl_grd_src, il_overlap_src, cl_period_src, w_unit, &
379	grad_lon, grad_lat, file_debug)
380	IF (file_debug) THEN
381	WRITE(w_unit,*) 'Conservative gradient calculated '
382	CALL FLUSH(w_unit)
383	ENDIF
384	!
385	! Send the local part of the coupling field and gradients
386	call oasis_put(var_id_s, 0, field_send, ierror, &
387	grad_lat, grad_lon)
388	ELSE
389	! 2nd order conservative not implemented for grids other than LR
390	WRITE(w_unit,*) 'Cannot perform second order conserv interpolation for type of grid ',cl_type_src
391	CALL oasis_abort(comp_id,comp_name,'Second order conserv interpolation impossible for that grid')
392	ENDIF
393	ELSE
394	! Standard oasis_put for other types of remappings for SCROP
395	CALL oasis_put(var_id_s, 0, field_send, ierror)
396	ENDIF
397	ELSE
398	CALL oasis_put(var_id_s, 0, field_send, ierror)
399	ENDIF
400	IF (file_debug) THEN
401	WRITE(w_unit,*) 'After oasis_put'
402	CALL FLUSH(w_unit)
403	ENDIF
404	!
405	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
406	! RECEIVE ARRAYS AND CALCULATE THE ERROR
407	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
408	!
409	! Allocate the field received by the model
410	!
411	ALLOCATE(field_recv(il_extentx_t, il_extenty_t), STAT=ierror )
412	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_recv'
413	ALLOCATE(field_ana(il_extentx_t, il_extenty_t),STAT=ierror )
414	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_ana'
415	ALLOCATE(field_error(il_extentx_t, il_extenty_t),STAT=ierror )
416	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_error'
417	ALLOCATE(mask_error(il_extentx_t, il_extenty_t),STAT=ierror )
418	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating mask_error'
419	!
420	#ifdef FANA1
421	CALL function_ana1(il_extentx_t, il_extenty_t, grid_lon_t, grid_lat_t, field_ana)
422	#elif defined FANA2
423	CALL function_ana2(il_extentx_t, il_extenty_t, grid_lon_t, grid_lat_t, field_ana)
424	#elif defined FANA3
425	CALL function_vortex(il_extentx_t, il_extenty_t, grid_lon_t, grid_lat_t, field_ana)
426	#endif
427	!
428	! Get the field FRECVANA
429	field_recv=-1.0
430	CALL oasis_get(var_id_t, 0, field_recv, ierror)
431	!
432	IF (ierror .NE. OASIS_Ok .AND. ierror .LT. OASIS_Recvd) THEN
433	WRITE (w_unit,*) 'oasis_get abort by model1 compid ',comp_id
434	CALL oasis_abort(comp_id,comp_name,'Problem in oasis_get')
435	ENDIF
436	IF (file_debug) THEN
437	WRITE(w_unit,'(A39,1X,F6.4,1X,F6.4)') 'MINVAL(field_recv),MAXVAL(field_recv)=',MINVAL(field_recv),MAXVAL(field_recv)
438	ENDIF
439	!
440	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
441	! Modify the field so to have:
442	! - on masked points: field value of 10000, error of 0.
443	! - on non-masked points that did not receive any interpolated value: field value of 1.e20, error of -1.e20
444	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
445	!
446	mask_error=1
447	WHERE (grid_msk_t == 1) ! masked points
448	field_error=0.
449	field_recv=10000.
450	mask_error=0
451	ELSEWHERE ! non-masked points
452	WHERE (field_recv /= 0.) ! non-masked points that received an interpolated value (works only if function is not 0 anywhere)
453	field_error = ABS(((field_ana - field_recv)/field_recv))*100
454	ELSEWHERE ! non-masked points that did not receive an interpolated value
455	field_error=-1.e20
456	field_recv=1.e20
457	mask_error=0
458	END WHERE
459	END WHERE
460	!
461	IF (file_debug) THEN
462	WRITE (w_unit,*) 'After calculating the interpolation error'
463	CALL FLUSH(w_unit)
464	ENDIF
465	!
466	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
467	! Gather field received and error on master process to calculate min and max and write them to files
468	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
469	!
470	IF (file_debug) THEN
471	WRITE (w_unit,*) 'Before gathering the field received and the error '
472	CALL FLUSH(w_unit)
473	ENDIF
474
475	! Gather the global field and global error on the master process
476	IF (mype == 0) THEN
477	ALLOCATE(il_size_all(npes))
478	ALLOCATE(il_offset_all(npes))
479	ALLOCATE(field_recv_global(nlon_t, nlat_t),STAT=ierror )
480	ALLOCATE(field_error_global(nlon_t, nlat_t),STAT=ierror )
481	IF ( ierror /= 0 ) WRITE(w_unit,*) 'Error allocating field_recv_global'
482	ALLOCATE(mask_error_global(nlon_t, nlat_t),STAT=ierror )
483	ALLOCATE(grid_msk_t_global(nlon_t, nlat_t),STAT=ierror )
484	ALLOCATE(grid_lon_global_t(nlon_t, nlat_t),STAT=ierror )
485	ALLOCATE(grid_lat_global_t(nlon_t, nlat_t),STAT=ierror )
486	field_recv_global(:,:) = 0.0
487	grid_lon_global_t(:,:) = 0.0
488	grid_lat_global_t(:,:) = 0.0
489	ENDIF
490	CALL MPI_GATHER(il_size_t, 1, MPI_INTEGER, il_size_all, 1, MPI_INTEGER, 0, local_comm, ierror)
491	WRITE (w_unit,*) 'After MPI_GATHER il_size_all target'
492	CALL FLUSH(w_unit)
493	CALL MPI_GATHER(il_offset_t, 1, MPI_INTEGER, il_offset_all, 1, MPI_INTEGER, 0, local_comm, ierror)
494	WRITE (w_unit,*) 'After MPI_GATHER il_offset_all target'
495	CALL FLUSH(w_unit)
496	CALL MPI_Gatherv(field_recv, il_size_t, MPI_DOUBLE_PRECISION, field_recv_global, il_size_all , il_offset_all, MPI_DOUBLE_PRECISION, 0, local_comm, ierror)
497	WRITE (w_unit,*) 'After MPI_GATHERv field_recv_global'
498	CALL FLUSH(w_unit)
499	CALL MPI_Gatherv(field_error, il_size_t, MPI_DOUBLE_PRECISION, field_error_global, il_size_all , il_offset_all, MPI_DOUBLE_PRECISION, 0, local_comm, ierror)
500	WRITE (w_unit,*) 'After MPI_GATHERv field_error_global'
501	CALL FLUSH(w_unit)
502	CALL MPI_Gatherv(mask_error, il_size_t, MPI_INTEGER, mask_error_global, il_size_all , il_offset_all, MPI_INTEGER, 0, local_comm, ierror)
503	WRITE (w_unit,*) 'After MPI_GATHERv mask_error_global'
504	CALL FLUSH(w_unit)
505	CALL MPI_Gatherv(grid_msk_t, il_size_t, MPI_INTEGER, grid_msk_t_global, il_size_all , il_offset_all, MPI_INTEGER, 0, local_comm, ierror)
506	WRITE (w_unit,*) 'After MPI_GATHERv mask_error_global'
507	CALL FLUSH(w_unit)
508	CALL MPI_Gatherv(grid_lon_t, il_size_t, MPI_DOUBLE_PRECISION, grid_lon_global_t, il_size_all , il_offset_all, MPI_DOUBLE_PRECISION, 0, local_comm, ierror)
509	WRITE (w_unit,*) 'After MPI_GATHERv grid_lon_global_t'
510	CALL FLUSH(w_unit)
511	CALL MPI_Gatherv(grid_lat_t, il_size_t, MPI_DOUBLE_PRECISION, grid_lat_global_t, il_size_all , il_offset_all, MPI_DOUBLE_PRECISION, 0, local_comm, ierror)
512	WRITE (w_unit,*) 'After MPI_GATHERv grid_lat_global_t'
513	CALL FLUSH(w_unit)
514	!
515	IF (mype == 0) THEN
516	data_filename='FRECVANA.nc'
517	field_name='FRECVANA'
518	!
519	CALL write_field(nlon_t, nlat_t, data_filename, field_name, w_unit, file_debug, grid_lon_global_t, grid_lat_global_t, field_recv_global)
520	IF (file_debug) THEN
521	WRITE (w_unit,*) 'After writing the field received'
522	CALL FLUSH(w_unit)
523	ENDIF
524	data_filename='REGRID_ERROR.nc'
525	field_name='ERROR'
526	!
527	CALL write_field(nlon_t, nlat_t, data_filename, field_name, w_unit, file_debug, grid_lon_global_t, grid_lat_global_t, field_error_global)
528	IF (file_debug) THEN
529	WRITE (w_unit,*) 'After writing the field received'
530	CALL FLUSH(w_unit)
531	ENDIF
532	!
533	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
534	! Calculate error min and max on non-masked points that received an interpolated value
535	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
536	!
537	min=MINVAL(field_error_global, MASK=mask_error_global>0)
538	IF (file_debug) THEN
539	WRITE(w_unit,*) 'Min (%) and its location in the error field : ',min
540	WRITE(w_unit,*) MINLOC(field_error_global)
541	CALL FLUSH(w_unit)
542	ENDIF
543	!
544	max=MAXVAL(field_error_global, MASK=mask_error_global>0)
545	IF (file_debug) THEN
546	WRITE(w_unit,'(A47,1X,F6.4)')'Max (%) and its location in the error field : ',max
547	WRITE(w_unit,*) MAXLOC(field_error_global)
548	CALL FLUSH(w_unit)
549	ENDIF
550	!
551	IF (file_debug) THEN
552	ic_nmsk=nlon_t*nlat_t-SUM(grid_msk_t_global)
553	WRITE(w_unit,*) 'Number of non-masked points :',ic_nmsk
554	ic_nmskrv=SUM(mask_error_global)
555	WRITE(w_unit,*) 'Number of non-masked points that received a value :',ic_nmskrv
556	WRITE(w_unit,'(A60,1X,F6.4)') 'Error mean on non masked points that received a value (%): ', SUM(ABS(field_error_global), MASK=mask_error_global>0)/ic_nmskrv
557	WRITE(w_unit,*) 'Delta error (%/echelle) :',(max - min)/echelle
558	WRITE(w_unit,*) 'End calculation of stat on the error'
559	CALL FLUSH(w_unit)
560	ENDIF
561	ELSE
562	WRITE(w_unit,*) 'Before oasis_terminate'
563	CALL FLUSH(w_unit)
564	ENDIF
565	!
566	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
567	! TERMINATION
568	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
569	IF (file_debug) THEN
570	WRITE(w_unit,*) 'End of the program, before oasis_terminate'
571	CALL FLUSH(w_unit)
572	ENDIF
573	!
574	CALL oasis_terminate (ierror)
575	IF (ierror /= 0) THEN
576	WRITE(w_unit,*) 'oasis_terminate abort by model1 compid ',comp_id
577	ENDIF
578	!
579	!
580	END PROGRAM MODEL1
581	!

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/regrid_environment/src/model1.F90 @ 6331

Download in other formats: