New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

grid.f90 in branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src/grid.f90 @ 4213

Last change on this file since 4213 was 4213, checked in by cbricaud, 10 years ago
first draft of the CONFIGURATION MANAGER demonstrator
File size: 94.9 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! MODULE: grid
6	!
7	! DESCRIPTION:
8	!> @brief grid manager <br/>
9	!>
10	!> @details
11	!>
12	!> @author
13	!> J.Paul
14	! REVISION HISTORY:
15	!> @date Nov, 2013 - Initial Version
16	!
17	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
18	!> @todo
19	!----------------------------------------------------------------------
20	MODULE grid
21	USE netcdf
22	USE kind ! F90 kind parameter
23	USE fct ! basic usefull function
24	USE global ! global parameter
25	USE phycst ! physical constant
26	USE logger ! log file manager
27	USE file ! file manager
28	USE var ! variable manager
29	USE dim ! dimension manager
30	USE dom ! domain manager
31	USE iom ! I/O manager
32	USE mpp ! MPP manager
33	USE iom_mpp ! MPP I/O manager
34	IMPLICIT NONE
35	PRIVATE
36	! NOTE_avoid_public_variables_if_possible
37
38	! type and variable
39
40	! function and subroutine
41	PUBLIC :: grid_check_dom !< check domain validity
42	PUBLIC :: grid_get_coarse_index !< get closest coarse grid indices of fine grid domain.
43	PUBLIC :: grid_is_global !< check if grid is global or not
44	PUBLIC :: grid_get_closest !< return closest coarse grid point from another point
45	PUBLIC :: grid_distance !< compute grid distance to a point
46	PUBLIC :: grid_get_fine_offset !< get fine grid offset
47	PUBLIC :: grid_check_coincidence !< check fine and coarse grid coincidence
48	PUBLIC :: grid_get_perio !< return NEMO periodicity index
49	PUBLIC :: grid_get_pivot !< return NEMO pivot point index
50	PUBLIC :: grid_add_ghost !< add ghost cell at boundaries.
51	PUBLIC :: grid_del_ghost !< delete ghost cell at boundaries.
52	PUBLIC :: grid_get_ghost !< return ghost cell factor
53	PUBLIC :: grid_split_domain !<
54	PUBLIC :: grid_fill_small_dom !<
55
56	PRIVATE :: grid_get_coarse_index_ff
57	PRIVATE :: grid_get_coarse_index_cf
58	PRIVATE :: grid_get_coarse_index_fc
59	PRIVATE :: grid_get_coarse_index_cc
60	PRIVATE :: grid__get_ghost_f
61	PRIVATE :: grid__get_ghost_ll
62	PRIVATE :: grid__check_corner
63
64	INTERFACE grid_get_ghost
65	MODULE PROCEDURE grid__get_ghost_ll
66	MODULE PROCEDURE grid__get_ghost_f
67	END INTERFACE grid_get_ghost
68
69	INTERFACE grid_get_coarse_index
70	MODULE PROCEDURE grid_get_coarse_index_ff
71	MODULE PROCEDURE grid_get_coarse_index_cf
72	MODULE PROCEDURE grid_get_coarse_index_fc
73	MODULE PROCEDURE grid_get_coarse_index_cc
74	END INTERFACE grid_get_coarse_index
75
76	CONTAINS
77	!-------------------------------------------------------------------
78	!> @brief
79	!> This funtion return NEMO pivot point index of the input variable.
80	!> - F-point : 0
81	!> - T-point : 1
82	!>
83	!> @warning
84	!> - variable must be nav_lon or nav_lat
85	!> - do not work with ORCA2 grid (T-point)
86	!>
87	!> @author J.Paul
88	!> - Nov, 2013- Subroutine written
89	!
90	!> @todo
91	!> - improve check between T or F pivot.
92	!
93	!> @param[in] td_file : file structure
94	!> @param[in] cd_varname : variable name
95	!> @return NEMO pivot point index
96	!-------------------------------------------------------------------
97	!> @code
98	INTEGER(i4) FUNCTION grid_get_pivot(td_file)
99	IMPLICIT NONE
100	! Argument
101	TYPE(TFILE), INTENT(IN) :: td_file
102
103	! local variable
104	TYPE(TVAR) :: tl_var
105	INTEGER(i4) :: il_varid
106	INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
107
108	! loop indices
109	INTEGER(i4) :: ji
110
111	INTEGER(i4) :: it1
112	INTEGER(i4) :: it2
113	INTEGER(i4) :: jt1
114	INTEGER(i4) :: jt2
115
116	INTEGER(i4) :: if1
117	INTEGER(i4) :: if2
118	INTEGER(i4) :: jf1
119	INTEGER(i4) :: jf2
120	!----------------------------------------------------------------
121	! initialise
122	grid_get_pivot=-1
123
124	! look for suitable variable
125	il_varid=0
126	DO ji=1,td_file%i_nvar
127	IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
128	SELECT CASE(TRIM(fct_lower(td_file%t_var(ji)%c_stdname)) )
129	CASE('longitude','latitude')
130	CASE DEFAULT
131	il_varid=ji
132	EXIT
133	END SELECT
134	ENDDO
135
136	IF( il_varid/=0 )THEN
137	IF( ASSOCIATED(td_file%t_var(il_varid)%d_value) )THEN
138	CALL logger_debug("GRID GET PIVOT: ASSOCIATED")
139	tl_var=td_file%t_var(il_varid)
140	ELSE
141	! read variable
142	il_dim(:)=td_file%t_var(il_varid)%t_dim(:)%i_len
143
144	CALL logger_debug("GRID GET PIVOT: read variable")
145	tl_var=iom_read_var(td_file, td_file%t_var(il_varid)%c_name, &
146	& id_start=(/1,il_dim(2)-3,1,1/), &
147	& id_count=(/3,4,1,1/) )
148	ENDIF
149
150	CALL logger_debug("GRID GET PIVOT: use variable "//TRIM(tl_var%c_name))
151
152	IF( ASSOCIATED(tl_var%d_value) )THEN
153
154	CALL logger_debug("GRID GET PIVOT: point "//TRIM(tl_var%c_point))
155	! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
156	it1=1 ; jt1=4
157	it2=3 ; jt2=2
158
159	! F-point pivot !case of ORCA05 grid
160	if1=1 ; jf1=4
161	if2=2 ; jf2=3
162
163	SELECT CASE(TRIM(tl_var%c_point))
164	CASE('T')
165	IF( ABS(tl_var%d_value(it1,jt1,1,1)) == &
166	& ABS(tl_var%d_value(it2,jt2,1,1)) )THEN
167	CALL logger_info("GRID GET PIVOT: T-pivot")
168	grid_get_pivot=1
169	ELSEIF( ABS(tl_var%d_value(if1,jf1,1,1)) == &
170	& ABS(tl_var%d_value(if2,jf2,1,1)) )THEN
171	CALL logger_info("GRID GET PIVOT: F-pivot")
172	grid_get_pivot=0
173	ELSE
174	CALL logger_error("GRID GET PIVOT: something wrong when "//&
175	& "computing pivot point")
176	ENDIF
177	CASE('U')
178	IF( ABS(tl_var%d_value(it1 ,jt1,1,1)) == &
179	& ABS(tl_var%d_value(it2-1,jt2,1,1)) )THEN
180	CALL logger_info("GRID GET PIVOT: T-pivot")
181	grid_get_pivot=1
182	ELSEIF( ABS(tl_var%d_value(if1 ,jf1,1,1)) == &
183	& ABS(tl_var%d_value(if2-1,jf2,1,1)) )THEN
184	CALL logger_info("GRID GET PIVOT: F-pivot")
185	grid_get_pivot=0
186	ELSE
187	CALL logger_error("GRID GET PIVOT: something wrong when "//&
188	& "computing pivot point")
189	ENDIF
190	CASE('V')
191	IF( ABS(tl_var%d_value(it1,jt1 ,1,1)) == &
192	& ABS(tl_var%d_value(it2,jt2-1,1,1)) )THEN
193	CALL logger_info("GRID GET PIVOT: T-pivot")
194	grid_get_pivot=1
195	ELSEIF( ABS(tl_var%d_value(if1,jf1 ,1,1)) == &
196	& ABS(tl_var%d_value(if2,jf2-1,1,1)) )THEN
197	CALL logger_info("GRID GET PIVOT: F-pivot")
198	grid_get_pivot=0
199	ELSE
200	CALL logger_error("GRID GET PIVOT: something wrong when "//&
201	& "computing pivot point")
202	ENDIF
203	CASE('F')
204	IF( ABS(tl_var%d_value(it1 ,jt1 ,1,1)) == &
205	& ABS(tl_var%d_value(it2-1,jt2-1,1,1)) )THEN
206	CALL logger_info("GRID GET PIVOT: T-pivot")
207	grid_get_pivot=1
208	ELSEIF( ABS(tl_var%d_value(if1 ,jf1 ,1,1)) == &
209	& ABS(tl_var%d_value(if2-1,jf2-1,1,1)) )THEN
210	CALL logger_info("GRID GET PIVOT: F-pivot")
211	grid_get_pivot=0
212	ELSE
213	CALL logger_error("GRID GET PIVOT: something wrong when "//&
214	& "computing pivot point")
215	ENDIF
216	END SELECT
217	ELSE
218	CALL logger_error("GRID GET PIVOT: can't compute pivot point. "//&
219	& "no value associated to variable "//TRIM(tl_var%c_name) )
220	ENDIF
221
222	ELSE
223	CALL logger_error("GRID GET PIVOT: no suitable variable to compute "//&
224	& "pivot point in file "//TRIM(td_file%c_name))
225	ENDIF
226
227	END FUNCTION grid_get_pivot
228	!> @endcode
229	!-------------------------------------------------------------------
230	!> @brief
231	!> This funtion return NEMO periodicity index of the input file.
232	!> The variable used must be on T point.
233	!>
234	!> @note the NEMO periodicity index can't be compute from coordinates file,
235	!> neither with mpp files.
236	!>
237	!> 0: closed boundaries
238	!> 1: cyclic east-west boundary
239	!> 2: symmetric boundary condition across the equator
240	!> 3: North fold boundary (with a F-point pivot)
241	!> 4: North fold boundary (with a F-point pivot) and cyclic east-west boundary
242	!> 5: North fold boundary (with a T-point pivot)
243	!> 6: North fold boundary (with a T-point pivot) and cyclic east-west boundary
244	!>
245	!> @author J.Paul
246	!> - Nov, 2013- Subroutine written
247	!
248	!> @todo
249	!> - improve check between T or F pivot.
250	!> - manage mpp case (read only border files)
251	!
252	!> @param[in] td_file : file structure
253	!> @param[in] id_pivot : pivot point
254	!> @return NEMO periodicity index
255	!-------------------------------------------------------------------
256	!> @code
257	INTEGER(i4) FUNCTION grid_get_perio(td_file, id_pivot)
258	IMPLICIT NONE
259
260	! Argument
261	TYPE(TFILE), INTENT(IN) :: td_file
262	INTEGER(i4), INTENT(IN) :: id_pivot
263
264	! local variable
265	TYPE(TVAR) :: tl_var
266	INTEGER(i4) :: il_varid
267	INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
268
269	! loop indices
270	INTEGER(i4) :: ji
271	!----------------------------------------------------------------
272
273	! initialise
274	grid_get_perio=-1
275
276	IF( id_pivot < 0 .OR. id_pivot > 1 )THEN
277	CALL logger_error("GRID GET PERIO: invalid pivot point index. "//&
278	& "you should use grid_get_pivot to compute it")
279	ENDIF
280
281	! look for suitable variable
282	il_varid=0
283	DO ji=1,td_file%i_nvar
284	IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
285	SELECT CASE(TRIM(fct_lower(td_file%t_var(ji)%c_stdname)) )
286	CASE('longitude','latitude')
287	CASE DEFAULT
288	il_varid=ji
289	EXIT
290	END SELECT
291	ENDDO
292
293	IF( il_varid==0 )THEN
294
295	CALL logger_error("GRID GET PERIO: no suitable variable to compute "//&
296	& " periodicity in file "//TRIM(td_file%c_name))
297	ELSE
298	il_dim(:)=td_file%t_var(il_varid)%t_dim(:)%i_len
299
300	IF( ASSOCIATED(td_file%t_var(il_varid)%d_value) )THEN
301	tl_var=td_file%t_var(il_varid)
302	ELSE
303	! read variable
304	tl_var=iom_read_var(td_file, td_file%t_var(il_varid)%c_name, &
305	& id_start=(/1,1,1,1/), &
306	& id_count=(/il_dim(1),il_dim(2),1,1/) )
307	ENDIF
308
309	IF(ALL(tl_var%d_value( 1 , : ,1,1)/=tl_var%d_fill).AND.&
310	& ALL(tl_var%d_value(il_dim(1), : ,1,1)/=tl_var%d_fill).AND.&
311	& ALL(tl_var%d_value( : , 1 ,1,1)/=tl_var%d_fill).AND.&
312	& ALL(tl_var%d_value( : ,il_dim(2),1,1)/=tl_var%d_fill))THEN
313	! no boundary closed
314	CALL logger_warn("GRID GET PERIO: can't determined periodicity. "//&
315	& "there is no boundary closed for variable "//&
316	& TRIM(tl_var%c_name)//" in file "//&
317	& TRIM(td_file%c_name) )
318	ELSE
319	! check periodicity
320	IF(ANY(tl_var%d_value( 1 ,:,1,1)/=tl_var%d_fill).OR.&
321	& ANY(tl_var%d_value(il_dim(1),:,1,1)/=tl_var%d_fill))THEN
322	! East-West cyclic (1,4,6)
323
324	IF( ANY(tl_var%d_value(:, 1, 1, 1) /= tl_var%d_fill) )THEN
325	! South boundary not closed
326
327	CALL logger_error("GRID GET PERIO: should have been an "//&
328	& "impossible case")
329	CALL logger_debug("GRID GET PERIO: East_West cyclic")
330	CALL logger_debug("GRID GET PERIO: South boundary not closed")
331
332	ELSE
333	! South boundary closed (1,4,6)
334	CALL logger_info("GRID GET PERIO: South boundary closed")
335
336	IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill) )THEN
337	! North boundary not closed (4,6)
338	CALL logger_info("GRID GET PERIO: North boundary not closed")
339	! check pivot
340	SELECT CASE(id_pivot)
341	CASE(0)
342	! F pivot
343	grid_get_perio=4
344	CASE(1)
345	! T pivot
346	grid_get_perio=6
347	CASE DEFAULT
348	CALL logger_error("GRID GET PERIO: invalid pivot ")
349	END SELECT
350	ELSE
351	! North boundary closed
352	CALL logger_info("GRID GET PERIO: North boundary closed")
353	grid_get_perio=1 ! North and South boundaries closed
354	ENDIF
355
356	ENDIF
357
358	ELSE
359	! East-West boundaries closed (0,2,3,5)
360	CALL logger_info("GRID GET PERIO: East West boundaries closed")
361
362	IF( ANY(tl_var%d_value(:, 1, 1, 1) /= tl_var%d_fill) )THEN
363	! South boundary not closed (2)
364	CALL logger_info("GRID GET PERIO: South boundary not closed")
365
366	IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill))THEN
367	! North boundary not closed
368	CALL logger_error("GRID GET PERIO: should have been "//&
369	& "an impossible case")
370	CALL logger_debug("GRID GET PERIO: East West boundaries "//&
371	& "closed")
372	CALL logger_debug("GRID GET PERIO: South boundary not closed")
373	CALL logger_debug("GRID GET PERIO: North boundary not closed")
374	ELSE
375	! North boundary closed
376	grid_get_perio=2 ! East-West and North boundaries closed
377	ENDIF
378
379	ELSE
380	! South boundary closed (0,3,5)
381	CALL logger_info("GRID GET PERIO: South boundary closed")
382
383	IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill))THEN
384	! North boundary not closed (3,5)
385	CALL logger_info("GRID GET PERIO: North boundary not closed")
386	! check pivot
387	SELECT CASE(id_pivot)
388	CASE(0)
389	! F pivot
390	grid_get_perio=3
391	CASE(1)
392	! T pivot
393	grid_get_perio=5
394	CASE DEFAULT
395	CALL logger_error("GRID GET PERIO: invalid pivot")
396	END SELECT
397	ELSE
398	! North boundary closed
399	CALL logger_info("GRID GET PERIO: North boundary closed")
400	grid_get_perio=0 ! all boundary closed
401	ENDIF
402
403	ENDIF
404
405	ENDIF
406
407	ENDIF
408	ENDIF
409
410	END FUNCTION grid_get_perio
411	!> @endcode
412	!-------------------------------------------------------------------
413	!> @brief This subroutine check domain validity.
414	!
415	!> @details
416	!> If maximum latitude greater than 88°N, program will stop.
417	!> It is not able to manage north fold boundary for now.
418	!
419	!> @author J.Paul
420	!> - Nov, 2013- Initial Version
421	!
422	!> @param[in] cd_coord : coordinate file
423	!> @param[in] id_imin : i-direction lower left point indice
424	!> @param[in] id_imax : i-direction upper right point indice
425	!> @param[in] id_jmin : j-direction lower left point indice
426	!> @param[in] id_jmax : j-direction upper right point indice
427	!>
428	!> @todo
429	!> - use domain instead of start count
430	!-------------------------------------------------------------------
431	!> @code
432	SUBROUTINE grid_check_dom(td_coord, id_imin, id_imax, id_jmin, id_jmax)
433	IMPLICIT NONE
434	! Argument
435	TYPE(TFILE), INTENT(IN) :: td_coord
436	INTEGER(i4), INTENT(IN) :: id_imin
437	INTEGER(i4), INTENT(IN) :: id_imax
438	INTEGER(i4), INTENT(IN) :: id_jmin
439	INTEGER(i4), INTENT(IN) :: id_jmax
440
441	! local variable
442	TYPE(TVAR) :: tl_var
443
444	TYPE(TFILE) :: tl_coord
445
446	TYPE(TMPP) :: tl_mppcoord
447
448	TYPE(TDOM) :: tl_dom
449
450	! loop indices
451	!----------------------------------------------------------------
452
453	IF( id_jmin >= id_jmax )THEN
454
455	CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
456	& "can not create configuration with north pole.")
457
458	ELSE
459
460	IF( td_coord%i_id == 0 )THEN
461	CALL logger_error("GRID CHECK DOM: can not check domain. "//&
462	& " file "//TRIM(td_coord%c_name)//" not opened." )
463	ELSE
464
465	IF( id_imin == id_imax .AND. td_coord%i_ew < 0 )THEN
466	CALL logger_fatal("GRID CHECK DOM: invalid domain."//&
467	& " can not create east-west cyclic fine grid"//&
468	& " inside closed coarse grid")
469	ENDIF
470
471	!1- read domain
472	tl_coord=td_coord
473	CALL iom_open(tl_coord)
474
475	!1-1 compute domain
476	tl_dom=dom_init( tl_coord, &
477	& id_imin, id_imax,&
478	& id_jmin, id_jmax )
479
480	!1-2 close file
481	CALL iom_close(tl_coord)
482
483	!1-3 read variables on domain (ugly way to do it, have to work on it)
484	!1-3-1 init mpp structure
485	tl_mppcoord=mpp_init(tl_coord)
486
487	CALL file_clean(tl_coord)
488
489	!1-3-2 get processor to be used
490	CALL mpp_get_use( tl_mppcoord, tl_dom )
491
492	!1-3-3 open mpp files
493	CALL iom_mpp_open(tl_mppcoord)
494
495	!1-3-4 read variable value on domain
496	tl_var=iom_mpp_read_var(tl_mppcoord,'latitude',td_dom=tl_dom)
497
498	!1-3-5 close mpp files
499	CALL iom_mpp_close(tl_mppcoord)
500
501	!1-3-6 clean structure
502	CALL mpp_clean(tl_mppcoord)
503
504	IF( MAXVAL(tl_var%d_value(:,:,:,:), &
505	& tl_var%d_value(:,:,:,:)/= tl_var%d_fill) >= 88.0 )THEN
506
507	CALL logger_debug("GRID CHECK DOM: max latitude "//&
508	& TRIM(fct_str(MAXVAL(tl_var%d_value(:,:,:,:)))) )
509	CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
510	& "can not create configuration too close from north pole.")
511
512	ENDIF
513
514	! clean
515	CALL var_clean(tl_var)
516
517	ENDIF
518
519
520	ENDIF
521
522	END SUBROUTINE grid_check_dom
523	!> @endcode
524	!-------------------------------------------------------------------
525	!> @brief This function get closest coarse grid indices of fine grid domain.
526	!
527	!> @details
528	!>
529	!
530	!> @author J.Paul
531	!> - Nov, 2013- Initial Version
532	!
533	!> @param[in] td_coord0 : coarse grid coordinate structure
534	!> @param[in] td_coord1 : fine grid coordinate structure
535	!> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
536	!> @todo
537	!> - use domain instead of start count
538	!-------------------------------------------------------------------
539	!> @code
540	FUNCTION grid_get_coarse_index_ff( td_coord0, td_coord1, &
541	& id_rho )
542	IMPLICIT NONE
543	! Argument
544	TYPE(TFILE), INTENT(IN) :: td_coord0
545	TYPE(TFILE), INTENT(IN) :: td_coord1
546	INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
547
548	! function
549	INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_ff
550
551	! local variable
552	TYPE(TVAR) :: tl_lon0
553	TYPE(TVAR) :: tl_lat0
554	TYPE(TVAR) :: tl_lon1
555	TYPE(TVAR) :: tl_lat1
556
557	INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho
558
559	INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
560	INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
561	INTEGER(i4), DIMENSION(2) :: il_xghost0
562	INTEGER(i4), DIMENSION(2) :: il_xghost1
563
564	INTEGER(i4) :: il_imin0
565	INTEGER(i4) :: il_imax0
566	INTEGER(i4) :: il_jmin0
567	INTEGER(i4) :: il_jmax0
568
569	INTEGER(i4) :: il_imin1
570	INTEGER(i4) :: il_imax1
571	INTEGER(i4) :: il_jmin1
572	INTEGER(i4) :: il_jmax1
573
574	! loop indices
575	!----------------------------------------------------------------
576
577	! init
578	grid_get_coarse_index_ff(:,:,:)=0
579
580	ALLOCATE(il_rho(ig_ndim))
581	il_rho(:)=1
582	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
583
584	IF( td_coord0%i_id == 0 .OR. td_coord1%i_id == 0 )THEN
585	CALL logger_error("GRID GET COARSE INDEX: can not get corase "//&
586	& "grid indices. file "//TRIM(td_coord0%c_name)//" and/or "//&
587	& TRIM(td_coord1%c_name)//" not opened." )
588	ELSE
589	!1- Coarse grid
590	! read coarse longitue and latitude
591	tl_lon0=iom_read_var(td_coord0,'longitude')
592	tl_lat0=iom_read_var(td_coord0,'latitude')
593
594	! get ghost cell factor on coarse grid
595	il_xghost0(:)=grid_get_ghost( tl_lon0, tl_lat0 )
596
597	il_imin0=1+il_xghost0(1)*ig_ghost
598	il_jmin0=1+il_xghost0(2)*ig_ghost
599
600	il_imax0=tl_lon0%t_dim(1)%i_len-il_xghost0(1)*ig_ghost
601	il_jmax0=tl_lon0%t_dim(2)%i_len-il_xghost0(2)*ig_ghost
602
603	CALL var_clean(tl_lon0)
604	CALL var_clean(tl_lat0)
605
606	! read coarse longitue and latitude without ghost cell
607	il_start(:)=(/il_imin0,il_jmin0,1,1/)
608	il_count(:)=(/il_imax0-il_imin0+1, &
609	& il_jmax0-il_jmin0+1, &
610	& tl_lon0%t_dim(3)%i_len, &
611	& tl_lon0%t_dim(4)%i_len /)
612
613	tl_lon0=iom_read_var(td_coord0,'longitude',il_start(:), il_count(:))
614	tl_lat0=iom_read_var(td_coord0,'latitude' ,il_start(:), il_count(:))
615
616	!2- Fine grid
617	! read fine longitue and latitude
618	tl_lon1=iom_read_var(td_coord1,'longitude')
619	tl_lat1=iom_read_var(td_coord1,'latitude')
620
621	! get ghost cell factor on fine grid
622	il_xghost1(:)=grid_get_ghost( tl_lon1, tl_lat1 )
623
624	il_imin1=1+il_xghost1(1)*ig_ghost
625	il_jmin1=1+il_xghost1(2)*ig_ghost
626
627	il_imax1=tl_lon1%t_dim(1)%i_len-il_xghost1(1)*ig_ghost
628	il_jmax1=tl_lon1%t_dim(2)%i_len-il_xghost1(2)*ig_ghost
629
630	CALL var_clean(tl_lon1)
631	CALL var_clean(tl_lat1)
632
633	! read fine longitue and latitude without ghost cell
634	il_start(:)=(/il_imin1,il_jmin1,1,1/)
635	il_count(:)=(/il_imax1-il_imin1+1, &
636	& il_jmax1-il_jmin1+1, &
637	& tl_lon1%t_dim(3)%i_len, &
638	& tl_lon1%t_dim(4)%i_len /)
639
640	tl_lon1=iom_read_var(td_coord1,'longitude',il_start(:), il_count(:))
641
642	tl_lat1=iom_read_var(td_coord1,'latitude' ,il_start(:), il_count(:))
643
644	!3- compute
645
646	grid_get_coarse_index_ff(:,:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
647	& tl_lon1,tl_lat1,&
648	& il_rho(:) )
649
650	il_imin0=grid_get_coarse_index_ff(1,1,1)-il_xghost0(1)*ig_ghost
651	il_imax0=grid_get_coarse_index_ff(1,2,1)+il_xghost0(1)*ig_ghost
652	il_jmin0=grid_get_coarse_index_ff(2,1,1)-il_xghost0(2)*ig_ghost
653	il_jmax0=grid_get_coarse_index_ff(2,2,1)+il_xghost0(2)*ig_ghost
654
655	grid_get_coarse_index_ff(1,1,1)=il_imin0
656	grid_get_coarse_index_ff(1,2,1)=il_imax0
657	grid_get_coarse_index_ff(2,1,1)=il_jmin0
658	grid_get_coarse_index_ff(2,2,1)=il_jmax0
659
660	CALL var_clean(tl_lon0)
661	CALL var_clean(tl_lat0)
662	CALL var_clean(tl_lon1)
663	CALL var_clean(tl_lat1)
664
665	ENDIF
666
667	END FUNCTION grid_get_coarse_index_ff
668	!> @endcode
669	!-------------------------------------------------------------------
670	!> @brief This function get closest coarse grid indices of fine grid domain.
671	!
672	!> @details
673	!>
674	!
675	!> @author J.Paul
676	!> - Nov, 2013- Initial Version
677	!
678	!> @param[in] td_longitude0 : coarse grid longitude
679	!> @param[in] td_latitude0 : coarse grid latitude
680	!> @param[in] td_coord1 : fine grid coordinate structure
681	!> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
682	!-------------------------------------------------------------------
683	!> @code
684	FUNCTION grid_get_coarse_index_cf( td_lon0, td_lat0, td_coord1, &
685	& id_rho )
686	IMPLICIT NONE
687	! Argument
688	TYPE(TVAR ), INTENT(IN) :: td_lon0
689	TYPE(TVAR ), INTENT(IN) :: td_lat0
690	TYPE(TFILE), INTENT(IN) :: td_coord1
691	INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
692
693	! function
694	INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_cf
695
696	! local variable
697	TYPE(TVAR) :: tl_lon1
698	TYPE(TVAR) :: tl_lat1
699
700	INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
701
702	INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
703	INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
704	INTEGER(i4), DIMENSION(2) :: il_xghost
705
706	INTEGER(i4) :: il_imin1
707	INTEGER(i4) :: il_imax1
708	INTEGER(i4) :: il_jmin1
709	INTEGER(i4) :: il_jmax1
710
711	! loop indices
712	!----------------------------------------------------------------
713
714	! init
715	grid_get_coarse_index_cf(:,:,:)=0
716
717	ALLOCATE(il_rho(ig_ndim) )
718	il_rho(:)=1
719	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
720
721	IF( td_coord1%i_id == 0 )THEN
722	CALL logger_error("GRID GET COARSE INDEX: file "//&
723	& TRIM(td_coord1%c_name)//" not opened." )
724
725	ELSE IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
726	& .NOT. ASSOCIATED(td_lat0%d_value) )THEN
727
728	CALL logger_error("GRID GET COARSE INDEX: some coarse grid"//&
729	& " coordinate value are not associated.")
730
731	ELSE
732
733	!1- Fine grid
734	! read fine longitue and latitude
735	tl_lon1=iom_read_var(td_coord1,'longitude')
736	tl_lat1=iom_read_var(td_coord1,'latitude')
737
738	! get ghost cell factor on fine grid
739	il_xghost(:)=grid_get_ghost( tl_lon1, tl_lat1 )
740
741	il_imin1=1+il_xghost(1)*ig_ghost
742	il_jmin1=1+il_xghost(2)*ig_ghost
743
744	il_imax1=tl_lon1%t_dim(1)%i_len-il_xghost(1)*ig_ghost
745	il_jmax1=tl_lon1%t_dim(2)%i_len-il_xghost(2)*ig_ghost
746
747	CALL var_clean(tl_lon1)
748	CALL var_clean(tl_lat1)
749
750	! read fine longitue and latitude without ghost cell
751	il_start(:)=(/il_imin1,il_jmin1,1,1/)
752	il_count(:)=(/il_imax1-il_imin1+1, &
753	& il_jmax1-il_jmin1+1, &
754	& tl_lon1%t_dim(3)%i_len, &
755	& tl_lon1%t_dim(4)%i_len /)
756
757	tl_lon1=iom_read_var(td_coord1,'longitude',il_start(:), il_count(:))
758	tl_lat1=iom_read_var(td_coord1,'latitude' ,il_start(:), il_count(:))
759
760	!3- compute
761	grid_get_coarse_index_cf(:,:,:)=grid_get_coarse_index(td_lon0,td_lat0,&
762	& tl_lon1,tl_lat1,&
763	& il_rho(:) )
764
765	CALL var_clean(tl_lon1)
766	CALL var_clean(tl_lat1)
767
768	ENDIF
769
770	END FUNCTION grid_get_coarse_index_cf
771	!> @endcode
772	!-------------------------------------------------------------------
773	!> @brief This function get closest coarse grid indices of fine grid domain.
774	!
775	!> @details
776	!>
777	!> @warning use ghost cell so can not be used on extracted domain without
778	!> ghost cell
779	!
780	!> @author J.Paul
781	!> - Nov, 2013- Initial Version
782	!
783	!> @param[in] td_coord0 : coarse grid coordinate structure
784	!> @param[in] td_lon1 : fine grid longitude
785	!> @param[in] td_lat1 : fine grid latitude
786	!> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
787	!-------------------------------------------------------------------
788	!> @code
789	FUNCTION grid_get_coarse_index_fc( td_coord0, td_lon1, td_lat1, &
790	& id_rho )
791	IMPLICIT NONE
792	! Argument
793	TYPE(TFILE), INTENT(IN) :: td_coord0
794	TYPE(TVAR ), INTENT(IN) :: td_lon1
795	TYPE(TVAR ), INTENT(IN) :: td_lat1
796	INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
797
798	! function
799	INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_fc
800
801	! local variable
802	TYPE(TVAR) :: tl_lon0
803	TYPE(TVAR) :: tl_lat0
804
805	INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
806
807	INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
808	INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
809	INTEGER(i4), DIMENSION(2) :: il_xghost
810
811	INTEGER(i4) :: il_imin0
812	INTEGER(i4) :: il_imax0
813	INTEGER(i4) :: il_jmin0
814	INTEGER(i4) :: il_jmax0
815
816
817	! loop indices
818	!----------------------------------------------------------------
819
820	! init
821	grid_get_coarse_index_fc(:,:,:)=0
822
823	ALLOCATE(il_rho(ig_ndim))
824	il_rho(:)=1
825	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
826
827	IF( td_coord0%i_id == 0 )THEN
828	CALL logger_error("GRID GET COARSE INDEX: file "//&
829	& TRIM(td_coord0%c_name)//" not opened." )
830
831	ELSE IF( .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
832	& .NOT. ASSOCIATED(td_lat1%d_value) )THEN
833
834	CALL logger_error("GRID GET COARSE INDEX: some fine grid"//&
835	& " coordinate value are not associated.")
836
837	ELSE
838	! read coarse longitue and latitude
839	tl_lon0=iom_read_var(td_coord0,'longitude')
840	tl_lat0=iom_read_var(td_coord0,'latitude')
841
842	! get ghost cell factor on coarse grid
843	il_xghost(:)=grid_get_ghost( tl_lon0, tl_lat0 )
844
845	il_imin0=1+il_xghost(1)*ig_ghost
846	il_jmin0=1+il_xghost(2)*ig_ghost
847
848	il_imax0=tl_lon0%t_dim(1)%i_len-il_xghost(1)*ig_ghost
849	il_jmax0=tl_lon0%t_dim(2)%i_len-il_xghost(2)*ig_ghost
850
851	CALL var_clean(tl_lon0)
852	CALL var_clean(tl_lat0)
853
854	! read coarse longitue and latitude without ghost cell
855	il_start(:)=(/il_imin0,il_jmin0,1,1/)
856	il_count(:)=(/il_imax0-il_imin0+1, &
857	& il_jmax0-il_jmin0+1, &
858	& tl_lon0%t_dim(3)%i_len, &
859	& tl_lon0%t_dim(4)%i_len /)
860
861	tl_lon0=iom_read_var(td_coord0,'longitude',il_start(:), il_count(:))
862	tl_lat0=iom_read_var(td_coord0,'latitude' ,il_start(:), il_count(:))
863
864	grid_get_coarse_index_fc(:,:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
865	& td_lon1,td_lat1,&
866	& il_rho(:) )
867
868	! remove ghost cell
869	il_imin0=grid_get_coarse_index_fc(1,1,1)+il_xghost(1)*ig_ghost
870	il_imax0=grid_get_coarse_index_fc(1,2,1)+il_xghost(1)*ig_ghost
871	il_jmin0=grid_get_coarse_index_fc(2,1,1)+il_xghost(2)*ig_ghost
872	il_jmax0=grid_get_coarse_index_fc(2,2,1)+il_xghost(2)*ig_ghost
873
874	grid_get_coarse_index_fc(1,1,1)=il_imin0
875	grid_get_coarse_index_fc(1,2,1)=il_imax0
876	grid_get_coarse_index_fc(2,1,1)=il_jmin0
877	grid_get_coarse_index_fc(2,2,1)=il_jmax0
878
879	CALL var_clean(tl_lon0)
880	CALL var_clean(tl_lat0)
881
882	ENDIF
883
884	END FUNCTION grid_get_coarse_index_fc
885	!> @endcode
886	!-------------------------------------------------------------------
887	!> @brief This function get closest coarse grid indices of fine grid domain.
888	!
889	!> @details
890	!>
891	!> @warning use ghost cell so can not be used on extracted domain without
892	!> ghost cell
893	!
894	!> @author J.Paul
895	!> - Nov, 2013- Initial Version
896	!
897	!> @param[in] td_lon0 : coarse grid longitude
898	!> @param[in] td_lat0 : coarse grid latitude
899	!> @param[in] td_lon1 : fine grid longitude
900	!> @param[in] td_lat1 : fine grid latitude
901	!> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
902	!>
903	!-------------------------------------------------------------------
904	!> @code
905	FUNCTION grid_get_coarse_index_cc( td_lon0, td_lat0, td_lon1, td_lat1, &
906	& id_rho )
907	IMPLICIT NONE
908	! Argument
909	TYPE(TVAR) , INTENT(IN) :: td_lon0
910	TYPE(TVAR) , INTENT(IN) :: td_lat0
911	TYPE(TVAR) , INTENT(IN) :: td_lon1
912	TYPE(TVAR) , INTENT(IN) :: td_lat1
913	INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
914
915	! function
916	INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_cc
917
918	! local variable
919	REAL(dp) :: dl_lon1_ll
920	REAL(dp) :: dl_lon1_ul
921	REAL(dp) :: dl_lon1_lr
922	REAL(dp) :: dl_lon1_ur
923
924	REAL(dp) :: dl_lat1_ll
925	REAL(dp) :: dl_lat1_ul
926	REAL(dp) :: dl_lat1_lr
927	REAL(dp) :: dl_lat1_ur
928
929	REAL(dp) :: dl_dlon
930	REAL(dp) :: dl_dlat
931
932	INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
933
934	INTEGER(i4), DIMENSION(2) :: il_ill
935	INTEGER(i4), DIMENSION(2) :: il_ilr
936	INTEGER(i4), DIMENSION(2) :: il_iul
937	INTEGER(i4), DIMENSION(2) :: il_iur
938
939	INTEGER(i4) :: il_ew0
940	INTEGER(i4) :: il_imin0
941	INTEGER(i4) :: il_imax0
942	INTEGER(i4) :: il_jmin0
943	INTEGER(i4) :: il_jmax0
944
945	INTEGER(i4) :: il_ew1
946	INTEGER(i4) :: il_imin1
947	INTEGER(i4) :: il_imax1
948	INTEGER(i4) :: il_jmin1
949	INTEGER(i4) :: il_jmax1
950
951	INTEGER(i4) :: il_imin
952	INTEGER(i4) :: il_imax
953	INTEGER(i4) :: il_jmin
954	INTEGER(i4) :: il_jmax
955
956	INTEGER(i4), DIMENSION(2,2) :: il_offset
957
958	! loop indices
959	INTEGER(i4) :: ji
960	INTEGER(i4) :: jj
961	!----------------------------------------------------------------
962
963	! init
964	grid_get_coarse_index_cc(:,:,:)=0
965
966	ALLOCATE( il_rho(ig_ndim) )
967	il_rho(:)=1
968	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
969
970	IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
971	& .NOT. ASSOCIATED(td_lat0%d_value) .OR. &
972	& .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
973	& .NOT. ASSOCIATED(td_lat1%d_value) )THEN
974	CALL logger_error("GRID GET COARSE INDEX: some fine or coarse grid"//&
975	& " coordinate value not associated.")
976	ELSE
977
978	IF( grid_is_global(td_lon1, td_lat1) )THEN
979
980	IF( grid_is_global(td_lon0, td_lat0) )THEN
981	CALL logger_trace("GRID GET COARSE INDEX: fine grid is global ")
982	grid_get_coarse_index_cc(:,:,1) = 1
983	grid_get_coarse_index_cc(:,:,2) = 0
984	ELSE
985	CALL logger_error("GRID GET COARSE INDEX: fine grid is "//&
986	& "global, coarse grid not.")
987	ENDIF
988
989	ELSE
990
991	! "global" coarse grid indice
992	il_imin0=1
993	il_jmin0=1
994
995	il_imax0=td_lon0%t_dim(1)%i_len
996	il_jmax0=td_lon0%t_dim(2)%i_len
997
998	! get east west overlap for coarse grid
999	il_ew0=dom_get_ew_overlap(td_lon0)
1000	IF( il_ew0 >= 0 )THEN
1001	! last point before overlap
1002	il_imax0=il_imax0-il_ew0
1003	ENDIF
1004
1005	! "global" fine grid indice
1006	il_imin1=1
1007	il_jmin1=1
1008
1009	il_imax1=td_lon1%t_dim(1)%i_len
1010	il_jmax1=td_lon1%t_dim(2)%i_len
1011
1012	! get east west overlap for coarse grid
1013	il_ew1=dom_get_ew_overlap(td_lon1)
1014	IF( il_ew1 >= 0 )THEN
1015	! last point before overlap
1016	il_imax1=il_imax1-il_ew1
1017	ENDIF
1018
1019	! get indices for each corner
1020	!1- search lower left corner indices
1021	dl_lon1_ll=td_lon1%d_value( il_imin1, il_jmin1, 1, 1 )
1022	dl_lat1_ll=td_lat1%d_value( il_imin1, il_jmin1, 1, 1 )
1023
1024	dl_dlon=ABS(td_lon1%d_value(il_imin1+1,il_jmin1 ,1,1)-dl_lon1_ll)
1025	dl_dlat=ABS(td_lat1%d_value(il_imin1 ,il_jmin1+1,1,1)-dl_lat1_ll)
1026
1027	! CALL logger_debug("GRID GET COARSE INDEX: lon1 ll "//&
1028	! & TRIM(fct_str(dl_lon1_ll)) )
1029	! CALL logger_debug("GRID GET COARSE INDEX: lat1 ll "//&
1030	! & TRIM(fct_str(dl_lat1_ll)) )
1031	!
1032	! CALL logger_debug("GRID GET COARSE INDEX: lon0 min "//&
1033	! & TRIM(fct_str(minval(td_lon0%d_value(2:,2:,:,:)))) )
1034	! CALL logger_debug("GRID GET COARSE INDEX: lon0 max "//&
1035	! & TRIM(fct_str(maxval(td_lon0%d_value(2:,2:,:,:)))) )
1036	!
1037	! CALL logger_debug("GRID GET COARSE INDEX: lat0 min "//&
1038	! & TRIM(fct_str(minval(td_lat0%d_value(2:,2:,:,:)))) )
1039	! CALL logger_debug("GRID GET COARSE INDEX: lat0 max "//&
1040	! & TRIM(fct_str(maxval(td_lat0%d_value(2:,2:,:,:)))) )
1041
1042	! look for closest point on coarse grid
1043	il_ill(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
1044	& il_jmin0:il_jmax0, &
1045	& 1,1), &
1046	& td_lat0%d_value(il_imin0:il_imax0, &
1047	& il_jmin0:il_jmax0, &
1048	& 1,1), &
1049	& dl_lon1_ll, dl_lat1_ll )
1050
1051	! coarse grid point should be south west of fine grid domain
1052	ji = il_ill(1)
1053	jj = il_ill(2)
1054
1055	IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ll) > dl_dlon*1.e-3 )THEN
1056	IF(td_lon0%d_value(ji,jj,1,1) > dl_lon1_ll ) il_ill(1)=il_ill(1)-1
1057	ENDIF
1058	IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ll) > dl_dlat*1.e-3 )THEN
1059	IF(td_lat0%d_value(ji,jj,1,1) > dl_lat1_ll ) il_ill(2)=il_ill(2)-1
1060	ENDIF
1061
1062	!2- search upper left corner indices
1063	dl_lon1_ul=td_lon1%d_value( il_imin1, il_jmax1, 1, 1 )
1064	dl_lat1_ul=td_lat1%d_value( il_imin1, il_jmax1, 1, 1 )
1065
1066	dl_dlon=ABS(td_lon1%d_value(il_imin1+1,il_jmax1 ,1,1)-dl_lon1_ll)
1067	dl_dlat=ABS(td_lat1%d_value(il_imin1 ,il_jmax1-1,1,1)-dl_lat1_ll)
1068
1069	! look for closest point on coarse grid
1070	il_iul(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
1071	& il_jmin0:il_jmax0, &
1072	& 1,1), &
1073	& td_lat0%d_value(il_imin0:il_imax0, &
1074	& il_jmin0:il_jmax0, &
1075	& 1,1), &
1076	& dl_lon1_ul, dl_lat1_ul )
1077
1078	! coarse grid point should be north west of fine grid domain
1079	ji = il_iul(1)
1080	jj = il_iul(2)
1081
1082	IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ul) > dl_dlon*1.e-3 )THEN
1083	IF(td_lon0%d_value(ji,jj,1,1) > dl_lon1_ul ) il_iul(1)=il_iul(1)-1
1084	ENDIF
1085	IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ul) > dl_dlat*1.e-3 )THEN
1086	IF(td_lat0%d_value(ji,jj,1,1) < dl_lat1_ul ) il_iul(2)=il_iul(2)+1
1087	ENDIF
1088
1089	!3- search lower right corner indices
1090	dl_lon1_lr=td_lon1%d_value( il_imax1, il_jmin1, 1, 1 )
1091	dl_lat1_lr=td_lat1%d_value( il_imax1, il_jmin1, 1, 1 )
1092
1093	dl_dlon=ABS(td_lon1%d_value(il_imax1-1,il_jmin1 ,1,1)-dl_lon1_ll)
1094	dl_dlat=ABS(td_lat1%d_value(il_imax1 ,il_jmin1+1,1,1)-dl_lat1_ll)
1095
1096	! look for closest point on coarse grid
1097	il_ilr(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
1098	& il_jmin0:il_jmax0, &
1099	& 1,1), &
1100	& td_lat0%d_value(il_imin0:il_imax0, &
1101	& il_jmin0:il_jmax0, &
1102	& 1,1), &
1103	& dl_lon1_lr, dl_lat1_lr )
1104
1105	! coarse grid point should be south east of fine grid domain
1106	ji = il_ilr(1)
1107	jj = il_ilr(2)
1108	IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_lr) > dl_dlon*1.e-3 )THEN
1109	IF( td_lon0%d_value(ji,jj,1,1) < dl_lon1_lr ) il_ilr(1)=il_ilr(1)+1
1110	ENDIF
1111	IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_lr) > dl_dlat*1.e-3 )THEN
1112	IF( td_lat0%d_value(ji,jj,1,1) > dl_lat1_lr ) il_ilr(2)=il_ilr(2)-1
1113	ENDIF
1114
1115	!4- search upper right corner indices
1116	dl_lon1_ur=td_lon1%d_value( il_imax1, il_jmax1, 1, 1 )
1117	dl_lat1_ur=td_lat1%d_value( il_imax1, il_jmax1, 1, 1 )
1118
1119	dl_dlon=ABS(td_lon1%d_value(il_imax1-1,il_jmax1 ,1,1)-dl_lon1_ll)
1120	dl_dlat=ABS(td_lat1%d_value(il_imax1 ,il_jmax1-1,1,1)-dl_lat1_ll)
1121
1122	! look for closest point on coarse grid
1123	il_iur(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
1124	& il_jmin0:il_jmax0, &
1125	& 1,1), &
1126	& td_lat0%d_value(il_imin0:il_imax0, &
1127	& il_jmin0:il_jmax0, &
1128	& 1,1), &
1129	& dl_lon1_ur, dl_lat1_ur )
1130
1131	! coarse grid point should be north east fine grid domain
1132	ji = il_iur(1)
1133	jj = il_iur(2)
1134	IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ur) > dl_dlon*1.e-3 )THEN
1135	IF( td_lon0%d_value(ji,jj,1,1) < dl_lon1_ur ) il_iur(1)=il_iur(1)+1
1136	ENDIF
1137	IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ur) > dl_dlat*1.e-3 )THEN
1138	IF( td_lat0%d_value(ji,jj,1,1) < dl_lat1_ur ) il_iur(2)=il_iur(2)+1
1139	ENDIF
1140
1141	! coarse grid indices
1142	il_imin = il_imin0-1+MIN(il_ill(1), il_iul(1))
1143	il_imax = il_imin0-1+MAX(il_ilr(1), il_iur(1))
1144
1145	IF( il_imax <= il_ew0 )THEN
1146	il_imax = td_lon0%t_dim(1)%i_len - il_ew0 + il_imax
1147	ENDIF
1148
1149	il_jmin = il_jmin0-1+MIN(il_ill(2), il_ilr(2))
1150	il_jmax = il_jmin0-1+MAX(il_iul(2), il_iur(2))
1151
1152	il_offset(:,:)= grid_get_fine_offset( td_lon0%d_value( :,:,1,1 ), &
1153	& td_lat0%d_value( :,:,1,1 ), &
1154	& il_imin, il_jmin, &
1155	& il_imax, il_jmax, &
1156	& td_lon1%d_value( :,:,1,1 ), &
1157	& td_lat1%d_value( :,:,1,1 ), &
1158	& il_rho(:) )
1159
1160	grid_get_coarse_index_cc(1,1,2) = il_offset(1,1)
1161	grid_get_coarse_index_cc(1,2,2) = il_offset(1,2)
1162
1163	grid_get_coarse_index_cc(2,1,2) = il_offset(2,1)
1164	grid_get_coarse_index_cc(2,2,2) = il_offset(2,2)
1165
1166	! special case if east west overlap
1167	IF( il_ew1 >= 0 )THEN
1168	CALL logger_debug("GRID GET COARSE INDEX: East-West overlap "//&
1169	& "found for fine grid " )
1170
1171	il_imin = 1
1172	il_imax = 1
1173
1174	grid_get_coarse_index_cc(1,1,2) = 0
1175	grid_get_coarse_index_cc(1,2,2) = 0
1176	ENDIF
1177
1178	ENDIF
1179
1180	IF( il_imin == il_imax ) il_imax=td_lon0%t_dim(1)%i_len
1181	IF( il_jmin == il_jmax ) il_jmax=td_lon0%t_dim(2)%i_len
1182
1183	grid_get_coarse_index_cc(1,1,1) = il_imin
1184	grid_get_coarse_index_cc(1,2,1) = il_imax
1185
1186	grid_get_coarse_index_cc(2,1,1) = il_jmin
1187	grid_get_coarse_index_cc(2,2,1) = il_jmax
1188
1189	ENDIF
1190
1191	END FUNCTION grid_get_coarse_index_cc
1192	!> @endcode
1193	!-------------------------------------------------------------------
1194	!> @brief This function check if grid is global or not
1195	!
1196	!> @details
1197	!
1198	!> @author J.Paul
1199	!> - Nov, 2013- Initial Version
1200	!
1201	!> @param[in] td_lon : longitude structure
1202	!> @param[in] td_lat : latitude structure
1203	!-------------------------------------------------------------------
1204	!> @code
1205	FUNCTION grid_is_global(td_lon, td_lat)
1206	IMPLICIT NONE
1207	! Argument
1208	TYPE(TVAR), INTENT(IN) :: td_lon
1209	TYPE(TVAR), INTENT(IN) :: td_lat
1210
1211	! function
1212	LOGICAL :: grid_is_global
1213
1214	! local variable
1215	INTEGER(i4) :: il_ew
1216	INTEGER(i4) :: il_south
1217	INTEGER(i4) :: il_north
1218
1219	REAL(dp) :: dl_lat_min
1220	REAL(dp) :: dl_lat_max
1221
1222	! loop indices
1223	!----------------------------------------------------------------
1224
1225	! init
1226	grid_is_global=.FALSE.
1227
1228	IF( ANY( td_lon%t_dim(:)%i_len /= td_lat%t_dim(:)%i_len ) )THEN
1229	CALL logger_fatal("GRID IS GLOBAL: dimension of longitude and "//&
1230	& " latitude differ")
1231	ENDIF
1232
1233	IF( .NOT. ASSOCIATED(td_lon%d_value) .OR. &
1234	& .NOT. ASSOCIATED(td_lat%d_value) )THEN
1235	CALL logger_error("GRID IS GLOBAL: na value associated to "//&
1236	& " longitude or latitude strucutre")
1237	ELSE
1238
1239	il_south=1
1240	il_north=td_lon%t_dim(2)%i_len
1241
1242	dl_lat_min=MINVAL(td_lat%d_value(:,il_south,1,1))
1243	dl_lat_max=MAXVAL(td_lat%d_value(:,il_north,1,1))
1244
1245	IF( dl_lat_min < -77.0 .AND. dl_lat_max > 89.0 )THEN
1246
1247	il_ew=td_lon%i_ew
1248	IF( il_ew >= 0 )THEN
1249
1250	grid_is_global=.TRUE.
1251
1252	ENDIF
1253
1254	ENDIF
1255	ENDIF
1256
1257	END FUNCTION grid_is_global
1258	!> @endcode
1259
1260	!-------------------------------------------------------------------
1261	!> @brief This function return coarse grid indices of the closest point
1262	!> from fine grid point (lon1,lat1)
1263	!>
1264	!
1265	!> @details
1266	!
1267	!> @note overlap band should have been already removed from coarse grid table
1268	!> of longitude and latitude, before running this function
1269	!>
1270	!> @author J.Paul
1271	!> - Nov, 2013- Initial Version
1272	!
1273	!> @param[in] dd_lon0 : coarse grid table of longitude
1274	!> @param[in] dd_lat0 : coarse grid table of latitude
1275	!> @param[in] dd_lon1 : fine grid longitude
1276	!> @param[in] dd_lat1 : fine grid latitude
1277	!> @return coarse grid indices of closest point of fine grid point
1278	!>
1279	!> @todo
1280	!-------------------------------------------------------------------
1281	!> @code
1282	FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1 )
1283	IMPLICIT NONE
1284	! Argument
1285	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
1286	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
1287	REAL(dp), INTENT(IN) :: dd_lon1
1288	REAL(dp), INTENT(IN) :: dd_lat1
1289
1290	! function
1291	INTEGER(i4), DIMENSION(2) :: grid_get_closest
1292
1293	! local variable
1294	INTEGER(i4) :: il_iinf
1295	INTEGER(i4) :: il_imid
1296	INTEGER(i4) :: il_isup
1297	INTEGER(i4) :: il_jinf
1298	INTEGER(i4) :: il_jmid
1299	INTEGER(i4) :: il_jsup
1300	INTEGER(i4), DIMENSION(2) :: il_shape
1301	INTEGER(i4), DIMENSION(1) :: il_ind
1302
1303	LOGICAL :: ll_north
1304	LOGICAL :: ll_continue
1305
1306	REAL(dp) :: dl_lon1
1307	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dist
1308	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
1309
1310	! loop indices
1311	!----------------------------------------------------------------
1312
1313	IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
1314	CALL logger_fatal("GRID GET CLOSEST: dimension of longitude and "//&
1315	& " latitude differ")
1316	ENDIF
1317
1318	il_shape(:)=SHAPE(dd_lon0(:,:))
1319
1320	ALLOCATE( dl_lon0(il_shape(1),il_shape(2)) )
1321
1322	dl_lon0(:,:) = dd_lon0(:,:)
1323	WHERE(dd_lon0(:,:) < 0 ) dl_lon0(:,:) = dd_lon0(:,:) + 360.
1324
1325	dl_lon1 = dd_lon1
1326	IF( dd_lon1 < 0 ) dl_lon1 = dd_lon1 + 360.
1327
1328	!1- first, use dichotomy to reduce domain
1329	il_iinf = 1 ; il_jinf = 1
1330	il_isup = il_shape(1) ; il_jsup = il_shape(2)
1331
1332	il_shape(1)= il_isup - il_iinf + 1
1333	il_shape(2)= il_jsup - il_jinf + 1
1334
1335	ll_north=.FALSE.
1336	ll_continue=.TRUE.
1337
1338	!1-1 look for meridian 0°/360°
1339	il_jmid = il_jinf + INT(il_shape(2)/2)
1340	il_ind(:) = MAXLOC( dl_lon0(:,il_jmid), dl_lon0(:,il_jmid) <= 360._dp )
1341
1342	il_imid=il_ind(1)
1343
1344	IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
1345	& dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN
1346
1347	il_iinf = il_imid ; il_isup = il_imid
1348	il_jinf = il_jmid ; il_jsup = il_jmid
1349
1350	ll_continue=.FALSE.
1351
1352	ELSE
1353	IF( dl_lon1 < dl_lon0(il_isup,il_jmid) .AND. &
1354	& il_imid /= il_isup )THEN
1355
1356	! point east
1357	il_iinf = il_imid
1358
1359	ELSE IF( dl_lon1 > dl_lon0(il_iinf,il_jmid) .AND. &
1360	& il_imid /= il_iinf )THEN
1361
1362	! point west
1363	il_isup = il_imid
1364
1365	ENDIF
1366
1367	il_shape(1)= il_isup - il_iinf + 1
1368	il_shape(2)= il_jsup - il_jinf + 1
1369
1370	il_imid = il_iinf + INT(il_shape(1)/2)
1371	il_jmid = il_jinf + INT(il_shape(2)/2)
1372
1373	! exit if too close from north fold (safer)
1374	IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE.
1375
1376	! exit when close enough of point
1377	IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
1378	ENDIF
1379
1380	!1-2
1381	DO WHILE( ll_continue .AND. .NOT. ll_north )
1382
1383	IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
1384	& dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN
1385
1386	il_iinf = il_imid ; il_isup = il_imid
1387	il_jinf = il_jmid ; il_jsup = il_jmid
1388
1389	ll_continue=.FALSE.
1390
1391	ELSE
1392	IF( dl_lon1 > dl_lon0(il_imid,il_jmid) )THEN
1393
1394	! point east
1395	il_iinf = il_imid
1396
1397	ELSE IF(dl_lon1 < dl_lon0(il_imid,il_jmid) )THEN
1398
1399	! point west
1400	il_isup = il_imid
1401
1402	ENDIF
1403
1404
1405	IF( dd_lat1 > dd_lat0(il_imid,il_jmid) )THEN
1406
1407	! point north
1408	il_jinf = il_jmid
1409
1410	ELSE IF(dd_lat1 < dd_lat0(il_imid,il_jmid) )THEN
1411
1412	! point south
1413	il_jsup = il_jmid
1414
1415	ENDIF
1416
1417	il_shape(1)= il_isup - il_iinf + 1
1418	il_shape(2)= il_jsup - il_jinf + 1
1419
1420	il_imid = il_iinf + INT(il_shape(1)/2)
1421	il_jmid = il_jinf + INT(il_shape(2)/2)
1422
1423	! exit if too close from north fold (safer)
1424	IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE.
1425
1426	! exit when close enough of point
1427	IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
1428	ENDIF
1429
1430	ENDDO
1431
1432	!2- then find closest point by computing distances
1433	il_shape(1)= il_isup - il_iinf + 1
1434	il_shape(2)= il_jsup - il_jinf + 1
1435
1436	ALLOCATE( dl_dist(il_shape(1), il_shape(2)) )
1437
1438	dl_dist(:,:)=grid_distance(dl_lon0(il_iinf:il_isup,il_jinf:il_jsup), &
1439	& dd_lat0(il_iinf:il_isup,il_jinf:il_jsup), &
1440	& dl_lon1, dd_lat1 )
1441
1442	grid_get_closest(:)=MINLOC(dl_dist(:,:),dl_dist(:,:)/=NF90_FILL_DOUBLE)
1443
1444	grid_get_closest(1)=grid_get_closest(1)+il_iinf-1
1445	grid_get_closest(2)=grid_get_closest(2)+il_jinf-1
1446
1447	DEALLOCATE( dl_dist )
1448	DEALLOCATE( dl_lon0 )
1449
1450	END FUNCTION grid_get_closest
1451	!> @endcode
1452	!-------------------------------------------------------------------
1453	!> @brief This function compute the distance between a point A and
1454	!> points of a grid
1455	!
1456	!> @details
1457	!
1458	!> @author J.Paul
1459	!> - Nov, 2013- Initial Version
1460	!
1461	!> @param[in] dd_lon : grid longitude table
1462	!> @param[in] dd_lat : grid latitude table
1463	!> @param[in] dd_lonA : longitude of point A
1464	!> @param[in] dd_latA : latitude of point A
1465	!-------------------------------------------------------------------
1466	!> @code
1467	FUNCTION grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA)
1468	IMPLICIT NONE
1469	! Argument
1470	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon
1471	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat
1472	REAL(dp), INTENT(IN) :: dd_lonA
1473	REAL(dp), INTENT(IN) :: dd_latA
1474
1475	! function
1476	REAL(dp), DIMENSION(SIZE(dd_lon(:,:),DIM=1),&
1477	& SIZE(dd_lon(:,:),DIM=2)) :: grid_distance
1478
1479	! local variable
1480	INTEGER(i4), DIMENSION(2) :: il_shape
1481
1482	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon
1483	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lat
1484	REAL(dp) :: dl_lonA
1485	REAL(dp) :: dl_latA
1486
1487	REAL(dp) :: dl_tmp
1488
1489	! loop indices
1490	INTEGER(i4) :: ji
1491	INTEGER(i4) :: jj
1492	!----------------------------------------------------------------
1493
1494	IF( ANY( SHAPE(dd_lon(:,:)) /= SHAPE(dd_lat(:,:)) ) )THEN
1495	CALL logger_fatal("GRID DISTANCE: dimension of longitude and "//&
1496	& " latitude differ")
1497	ENDIF
1498	il_shape(:)=SHAPE(dd_lon(:,:))
1499
1500	ALLOCATE(dl_lon(il_shape(1),il_shape(2)))
1501	ALLOCATE(dl_lat(il_shape(1),il_shape(2)))
1502
1503	dl_lon(:,:) = dd_lon(:,:)
1504	dl_lonA = dd_lonA
1505
1506	WHERE(dd_lon(:,:) < 0 ) dl_lon(:,:) = dd_lon(:,:) + 360.
1507	IF( dd_lonA < 0 ) dl_lonA = dd_lonA + 360.
1508
1509	dl_lonA = dd_lonA * dg_deg2rad
1510	dl_latA = dd_latA * dg_deg2rad
1511
1512	dl_lon(:,:) = dl_lon(:,:) * dg_deg2rad
1513	dl_lat(:,:) = dd_lat(:,:) * dg_deg2rad
1514
1515	grid_distance(:,:)=NF90_FILL_DOUBLE
1516
1517	DO jj=1,il_shape(2)
1518	DO ji=1,il_shape(1)
1519	IF( dl_lon(ji,jj) == dl_lonA .AND. &
1520	& dl_lat(ji,jj) == dl_lATA )THEN
1521	grid_distance(ji,jj)=0.0
1522	ELSE
1523	dl_tmp= SIN(dl_latA)*SIN(dl_lat(ji,jj)) + &
1524	& COS(dl_latA)COS(dl_lat(ji,jj))COS(dl_lon(ji,jj)-dl_lonA)
1525	! check to avoid mistake with ACOS
1526	IF( dl_tmp < -1.0 ) dl_tmp = -1.0
1527	IF( dl_tmp > 1.0 ) dl_tmp = 1.0
1528	grid_distance(ji,jj)=ACOS(dl_tmp)*dg_rearth
1529	ENDIF
1530	ENDDO
1531	ENDDO
1532
1533	DEALLOCATE(dl_lon)
1534	DEALLOCATE(dl_lat)
1535
1536	END FUNCTION grid_distance
1537	!> @endcode
1538	!-------------------------------------------------------------------
1539	!> @brief This function get fine grid offset.
1540	!
1541	!> @details
1542	!> offset value could be 0,1,..,rho-1
1543	!
1544	!> @author J.Paul
1545	!> - Nov, 2013- Initial Version
1546	!
1547	!> @param[in] dd_lon0 : coarse grid longitude table
1548	!> @param[in] dd_lat0 : coarse grid latitude table
1549	!> @param[in] dd_lon1 : fine grid longitude table
1550	!> @param[in] dd_lat1 : fine grid latitude table
1551	!> @param[in] id_imin0 : coarse grid lower left corner i-indice of fine grid domain
1552	!> @param[in] id_jmin0 : coarse grid lower left corner j-indice of fine grid domain
1553	!> @param[in] id_imax0 : coarse grid upper right corner i-indice of fine grid domain
1554	!> @param[in] id_jmax0 : coarse grid upper right corner j-indice of fine grid domain
1555	!> @param[in] id_rhoi : i-direction refinement factor
1556	!> @param[in] id_rhoj : j-direction refinement factor
1557	!> @return offset table (/ (/i_offset_left,i_offset_right!/),(/j_offset_lower,j_offset_upper/) /)
1558	!-------------------------------------------------------------------
1559	!> @code
1560	FUNCTION grid_get_fine_offset( dd_lon0, dd_lat0, &
1561	& id_imin0, id_jmin0, id_imax0, id_jmax0, &
1562	& dd_lon1, dd_lat1, id_rho )
1563	IMPLICIT NONE
1564	! Argument
1565	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon0
1566	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat0
1567	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon1
1568	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat1
1569
1570	INTEGER(i4), INTENT(IN) :: id_imin0
1571	INTEGER(i4), INTENT(IN) :: id_jmin0
1572	INTEGER(i4), INTENT(IN) :: id_imax0
1573	INTEGER(i4), INTENT(IN) :: id_jmax0
1574
1575	INTEGER(i4), DIMENSION(:) , INTENT(IN) :: id_rho
1576
1577	! function
1578	INTEGER(i4), DIMENSION(2,2) :: grid_get_fine_offset
1579
1580	! local variable
1581	INTEGER(i4), DIMENSION(2) :: il_shape0
1582	INTEGER(i4), DIMENSION(2) :: il_shape1
1583	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
1584	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon1
1585
1586	REAL(dp) :: dl_dlon
1587	REAL(dp) :: dl_dlat
1588
1589	! loop indices
1590	INTEGER(i4) :: ji
1591	INTEGER(i4) :: jj
1592
1593	INTEGER(i4) :: ii
1594	INTEGER(i4) :: ij
1595	!----------------------------------------------------------------
1596	IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
1597	CALL logger_fatal("GRID GET FINE OFFSET: dimension of coarse "//&
1598	& "longitude and latitude differ")
1599	ENDIF
1600
1601	IF( ANY( SHAPE(dd_lon1(:,:)) /= SHAPE(dd_lat1(:,:)) ) )THEN
1602	CALL logger_fatal("GRID GET FINE OFFSET: dimension of fine "//&
1603	& "longitude and latitude differ")
1604	ENDIF
1605
1606	il_shape0(:)=SHAPE(dd_lon0(:,:))
1607	ALLOCATE( dl_lon0(il_shape0(1),il_shape0(2)) )
1608
1609	dl_lon0(:,:)=dd_lon0(:,:)
1610	WHERE( dd_lon0(:,:) < 0 ) dl_lon0(:,:)=dd_lon0(:,:)+360.
1611
1612	il_shape1(:)=SHAPE(dd_lon1(:,:))
1613	ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) )
1614
1615	dl_lon1(:,:)=dd_lon1(:,:)
1616	WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360.
1617
1618	grid_get_fine_offset(:,:)=-1
1619
1620	! look for i-direction left offset
1621	IF( dl_lon1(1,1) < dl_lon0(id_imin0+1,id_jmin0) )THEN
1622	DO ji=1,id_rho(jp_I)+2
1623	dl_dlon=ABS(dl_lon1(ji+1,1)-dl_lon1(ji,1))*1.e-3
1624	IF( dl_lon1(ji,1) > dl_lon0(id_imin0+1,id_jmin0) + dl_dlon )THEN
1625	grid_get_fine_offset(1,1)=(id_rho(jp_I)+1)-ji+MOD(id_rho(jp_I),2)
1626	EXIT
1627	ENDIF
1628	ENDDO
1629	ELSE
1630	CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
1631	& " not match fine grid lower left corner.")
1632	ENDIF
1633
1634	! look for i-direction right offset
1635	IF( dl_lon1(il_shape1(1),1) > dl_lon0(id_imax0-1,id_jmin0) )THEN
1636	DO ji=1,id_rho(jp_I)+2
1637	ii=il_shape1(1)-ji+1
1638	dl_dlon=ABS(dl_lon1(ii,1)-dl_lon1(ii-1,1))*1.e-3
1639	IF( dl_lon1(ii,1) < dl_lon0(id_imax0-1,id_jmin0) - dl_dlon )THEN
1640	grid_get_fine_offset(1,2)=(id_rho(jp_I)+1)-ji+MOD(id_rho(jp_I),2)
1641	EXIT
1642	ENDIF
1643	ENDDO
1644	ELSE
1645	CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
1646	& " not match fine grid lower right corner.")
1647	ENDIF
1648
1649	! look for j-direction lower offset
1650	IF( dd_lat1(1,1) < dd_lat0(id_imin0,id_jmin0+1) )THEN
1651	DO jj=1,id_rho(jp_J)+2
1652	dl_dlat=ABS(dd_lat1(1,jj+1)-dd_lat1(1,jj))*1.e-3
1653	IF( dd_lat1(1,jj) > dd_lat0(id_imin0,id_jmin0+1) + dl_dlat )THEN
1654	grid_get_fine_offset(2,1)=(id_rho(jp_J)+1)-jj+MOD(id_rho(jp_J),2)
1655	EXIT
1656	ENDIF
1657	ENDDO
1658	ELSE
1659	CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
1660	& " not match fine grid upper left corner.")
1661	ENDIF
1662
1663	! look for j-direction upper offset
1664	IF( dd_lat1(1,il_shape1(2)) > dd_lat0(id_imin0,id_jmax0-1) )THEN
1665	DO jj=1,id_rho(jp_J)+2
1666	ij=il_shape1(2)-jj+1
1667	dl_dlat=ABS(dd_lat1(1,ij)-dd_lat1(1,ij-1))*1.e-3
1668	IF( dd_lat1(1,ij) < dd_lat0(id_imin0,id_jmax0-1) - dl_dlat )THEN
1669	grid_get_fine_offset(2,2)=(id_rho(jp_J)+1)-jj+MOD(id_rho(jp_J),2)
1670	EXIT
1671	ENDIF
1672	ENDDO
1673	ELSE
1674	CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
1675	& " not match fine grid upper right corner.")
1676	ENDIF
1677
1678	DEALLOCATE( dl_lon0 )
1679	DEALLOCATE( dl_lon1 )
1680
1681	END FUNCTION grid_get_fine_offset
1682	!> @endcode
1683	!-------------------------------------------------------------------
1684	!> @brief This function check if ghost cell are used or not, and return ghost
1685	!> cell factor (0,1) in i- and j-direction.
1686	!
1687	!> @details
1688	!
1689	!> @author J.Paul
1690	!> - Nov, 2013- Initial Version
1691	!
1692	!> @param[in] td_lon : grid longitude sturcture
1693	!> @param[in] td_lat : grid latitude structure
1694	!-------------------------------------------------------------------
1695	!> @code
1696	FUNCTION grid__get_ghost_ll( td_lon, td_lat )
1697	IMPLICIT NONE
1698	! Argument
1699	TYPE(TVAR), INTENT(IN) :: td_lon
1700	TYPE(TVAR), INTENT(IN) :: td_lat
1701
1702	! function
1703	INTEGER(i4), DIMENSION(2) :: grid__get_ghost_ll
1704
1705	! local variable
1706	INTEGER(i4) :: il_ew
1707	! loop indices
1708	!----------------------------------------------------------------
1709	! init
1710	grid__get_ghost_ll(:)=0
1711
1712	IF( grid_is_global(td_lon, td_lat) )THEN
1713	grid__get_ghost_ll(:)=0
1714	ELSE
1715	grid__get_ghost_ll(2)=1
1716
1717	il_ew=td_lon%i_ew
1718	IF( il_ew < 0 )THEN
1719	grid__get_ghost_ll(1)=1
1720	ELSE
1721	grid__get_ghost_ll(1)=0
1722	ENDIF
1723	ENDIF
1724
1725	END FUNCTION grid__get_ghost_ll
1726	!> @endcode
1727	!-------------------------------------------------------------------
1728	!> @brief This function check if ghost cell are used or not, and return ghost
1729	!> cell factor (0,1) in i- and j-direction.
1730	!
1731	!> @details
1732	!
1733	!> @author J.Paul
1734	!> - Nov, 2013- Initial Version
1735	!
1736	!> @param[in] td_file : file sturcture
1737	!-------------------------------------------------------------------
1738	!> @code
1739	FUNCTION grid__get_ghost_f( td_file )
1740	IMPLICIT NONE
1741	! Argument
1742	TYPE(TFILE), INTENT(IN) :: td_file
1743
1744	! function
1745	INTEGER(i4), DIMENSION(2) :: grid__get_ghost_f
1746
1747	! local variable
1748	TYPE(TVAR) :: tl_lon
1749	TYPE(TVAR) :: tl_lat
1750
1751	INTEGER(i4) :: il_lonid
1752	INTEGER(i4) :: il_latid
1753	! loop indices
1754	INTEGER(i4) :: ji
1755	!----------------------------------------------------------------
1756	! init
1757	grid__get_ghost_f(:)=0
1758
1759	IF( td_file%i_id == 0 )THEN
1760	CALL logger_error("GRID GET GHOST: file "//&
1761	& TRIM(td_file%c_name)//" not opened." )
1762
1763	ELSE
1764
1765	IF( ASSOCIATED(td_file%t_var) )THEN
1766	! read coarse longitue and latitude
1767	il_lonid=var_get_id(td_file%t_var(:),'longitude')
1768	il_latid=var_get_id(td_file%t_var(:),'latitude')
1769
1770	print *,'file ',trim(td_file%c_name),td_file%i_ew
1771	DO ji=1,td_file%i_nvar
1772	print *,ji,trim(td_file%t_var(ji)%c_name),': ',td_file%t_var(ji)%i_ew
1773	ENDDO
1774	print *,'lonid ',il_lonid
1775	print *,'latid ',il_latid
1776	IF( il_lonid /=0 .AND. il_latid /= 0 )THEN
1777	tl_lon=iom_read_var(td_file,il_lonid)
1778	print *,'lon ',tl_lon%i_ew
1779	tl_lat=iom_read_var(td_file,il_latid)
1780	print *,'lat ',tl_lat%i_ew
1781	! get ghost cell factor on coarse grid
1782	grid__get_ghost_f(:)=grid_get_ghost( tl_lon, tl_lat )
1783	ELSE
1784	CALL logger_error("GRID GET GHOST: can not find "//&
1785	& "longitude or latitude "//&
1786	& "in file "//TRIM(td_file%c_name))
1787	ENDIF
1788	ELSE
1789	CALL logger_error("GRID GET GHOST: no variable "//&
1790	& "associated to file "//TRIM(td_file%c_name))
1791	ENDIF
1792
1793	ENDIF
1794
1795	END FUNCTION grid__get_ghost_f
1796	!> @endcode
1797	!-------------------------------------------------------------------
1798	!> @brief This subroutine check fine and coarse grid coincidence
1799	!
1800	!> @details
1801	!
1802	!> @author J.Paul
1803	!> - Nov, 2013- Initial Version
1804	!
1805	!> @param[in] td_coord0 : coarse grid coordinate file structure
1806	!> @param[in] td_coord1 : fine grid coordinate file structure
1807	!> @param[in] id_imin0 : coarse grid lower left corner i-indice of fine grid domain
1808	!> @param[in] id_imax0 : coarse grid upper right corner i-indice of fine grid domain
1809	!> @param[in] id_jmin0 : coarse grid lower left corner j-indice of fine grid domain
1810	!> @param[in] id_jmax0 : coarse grid upper right corner j-indice of fine grid domain
1811	!> @param[in] id_rho : table of refinement factor
1812	!-------------------------------------------------------------------
1813	!> @code
1814	SUBROUTINE grid_check_coincidence( td_coord0, td_coord1, &
1815	& id_imin0, id_imax0, &
1816	& id_jmin0, id_jmax0, &
1817	& id_rho )
1818	IMPLICIT NONE
1819
1820	! Argument
1821	TYPE(TFILE), INTENT(IN) :: td_coord0
1822	TYPE(TFILE), INTENT(IN) :: td_coord1
1823	INTEGER(i4), INTENT(IN) :: id_imin0
1824	INTEGER(i4), INTENT(IN) :: id_imax0
1825	INTEGER(i4), INTENT(IN) :: id_jmin0
1826	INTEGER(i4), INTENT(IN) :: id_jmax0
1827	INTEGER(i4), DIMENSION(:), INTENT(IN) :: id_rho
1828
1829	! local variable
1830	INTEGER(i4) :: il_imid1
1831	INTEGER(i4) :: il_jmid1
1832
1833	INTEGER(i4) :: il_ew0
1834	INTEGER(i4) :: il_ew1
1835
1836	INTEGER(i4) :: il_imin1
1837	INTEGER(i4) :: il_imax1
1838	INTEGER(i4) :: il_jmin1
1839	INTEGER(i4) :: il_jmax1
1840
1841	INTEGER(i4), DIMENSION(2) :: il_indC
1842	INTEGER(i4), DIMENSION(2) :: il_indF
1843	INTEGER(i4), DIMENSION(2) :: il_iind
1844	INTEGER(i4), DIMENSION(2) :: il_jind
1845
1846	REAL(dp) :: dl_lon0
1847	REAL(dp) :: dl_lat0
1848	REAL(dp) :: dl_lon1
1849	REAL(dp) :: dl_lat1
1850
1851	REAL(dp) :: dl_lon1p
1852	REAL(dp) :: dl_lat1p
1853
1854	REAL(dp) :: dl_dlon
1855	REAL(dp) :: dl_dlat
1856
1857	LOGICAL :: ll_coincidence
1858
1859	TYPE(TVAR) :: tl_lon0
1860	TYPE(TVAR) :: tl_lat0
1861	TYPE(TVAR) :: tl_lon1
1862	TYPE(TVAR) :: tl_lat1
1863
1864	TYPE(TFILE) :: tl_coord0
1865
1866	TYPE(TMPP) :: tl_mppcoord0
1867
1868	TYPE(TDOM) :: tl_dom0
1869
1870	! loop indices
1871	INTEGER(i4) :: ji
1872	INTEGER(i4) :: jj
1873	!----------------------------------------------------------------
1874
1875	ll_coincidence=.TRUE.
1876
1877	! read coarse longitue and latitude on domain
1878	tl_coord0=td_coord0
1879	CALL iom_open(tl_coord0)
1880
1881	!2-1 compute domain
1882	tl_dom0=dom_init( tl_coord0, &
1883	& id_imin0, id_imax0,&
1884	& id_jmin0, id_jmax0 )
1885
1886	!2-2 close file
1887	CALL iom_close(tl_coord0)
1888
1889	!2-3 read variables on domain (ugly way to do it, have to work on it)
1890	!2-3-1 init mpp structure
1891	tl_mppcoord0=mpp_init(tl_coord0)
1892
1893	CALL file_clean(tl_coord0)
1894
1895	!2-3-2 get processor to be used
1896	CALL mpp_get_use( tl_mppcoord0, tl_dom0 )
1897
1898	!2-3-3 open mpp files
1899	CALL iom_mpp_open(tl_mppcoord0)
1900
1901	!2-3-4 read variable value on domain
1902	tl_lon0=iom_mpp_read_var(tl_mppcoord0,'longitude',td_dom=tl_dom0)
1903	tl_lat0=iom_mpp_read_var(tl_mppcoord0,'latitude' ,td_dom=tl_dom0)
1904
1905	!2-3-5 close mpp files
1906	CALL iom_mpp_close(tl_mppcoord0)
1907
1908	!2-3-6 clean structure
1909	CALL mpp_clean(tl_mppcoord0)
1910
1911	! read fine longitue and latitude
1912	tl_lon1=iom_read_var(td_coord1,'longitude')
1913	tl_lat1=iom_read_var(td_coord1,'latitude')
1914
1915	CALL logger_debug("GRID CHECK COINCIDENCE:"//&
1916	& " fine grid "//TRIM(td_coord1%c_name) )
1917	CALL logger_debug("GRID CHECK COINCIDENCE:"//&
1918	& " coarse grid "//TRIM(td_coord0%c_name) )
1919
1920	!1- check domain
1921	!1-1 check global grid
1922	IF( .NOT. grid_is_global(tl_lon0, tl_lat0) )THEN
1923	IF( grid_is_global(tl_lon1, tl_lat1) )THEN
1924
1925	ll_coincidence=.FALSE.
1926	CALL logger_fatal("GRID CHECK COINCIDENCE:"//&
1927	& " fine grid is global,"//&
1928	& " coarse grid is not ")
1929
1930	ELSE
1931	!1-2 ew overlap
1932	il_ew1=tl_lon1%i_ew
1933	IF( il_ew1 >= 0 )THEN
1934
1935	il_ew0=tl_lon0%i_ew
1936	IF( il_ew0 < 0 )THEN
1937	CALL logger_fatal("GRID CHECK COINCIDENCE: "//&
1938	& "fine grid has east west overlap,"//&
1939	& " coarse grid not ")
1940	ENDIF
1941
1942	il_jmin1=1+ig_ghost
1943	il_jmax1=tl_lon1%t_dim(2)%i_len-ig_ghost
1944
1945	ll_coincidence=grid__check_lat(&
1946	& tl_lat0%d_value(1,:,1,1),&
1947	& tl_lat1%d_value(1,il_jmin1:il_jmax1,1,1),&
1948	& id_rho(jp_J) )
1949
1950	ELSE
1951	!1-3 other case
1952	il_imin1=1+ig_ghost
1953	il_jmin1=1+ig_ghost
1954
1955	il_imax1=tl_lon1%t_dim(1)%i_len-ig_ghost
1956	il_jmax1=tl_lon1%t_dim(2)%i_len-ig_ghost
1957
1958	ll_coincidence=grid__check_corner(&
1959	& tl_lon0%d_value(:,:,1,1),&
1960	& tl_lat0%d_value(:,:,1,1),&
1961	& tl_lon1%d_value(il_imin1:il_imax1, &
1962	& il_jmin1:il_jmax1, &
1963	& 1,1),&
1964	& tl_lat1%d_value(il_imin1:il_imax1, &
1965	& il_jmin1:il_jmax1, &
1966	& 1,1) )
1967
1968	ENDIF
1969
1970	ENDIF
1971
1972	IF( .NOT. ll_coincidence )THEN
1973	CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
1974	& "between fine grid and coarse grid. invalid domain" )
1975	ENDIF
1976
1977	ENDIF
1978
1979	!2- check refinement factor
1980	! select point in middle of fine grid
1981	il_imid1=INT(tl_lon1%t_dim(1)%i_len*0.5)
1982	il_jmid1=INT(tl_lon1%t_dim(2)%i_len*0.5)
1983
1984	dl_lon1=tl_lon1%d_value(il_imid1, il_jmid1,1,1)
1985	dl_lat1=tl_lat1%d_value(il_imid1, il_jmid1,1,1)
1986
1987	! select closest point on coarse grid
1988	il_indC(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),&
1989	& tl_lat0%d_value(:,:,1,1),&
1990	& dl_lon1, dl_lat1 )
1991
1992	IF( ANY(il_indC(:)==0) )THEN
1993	CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//&
1994	& "coarse grid indices. invalid domain" )
1995	ENDIF
1996
1997	dl_lon0=tl_lon0%d_value(il_indC(1),il_indC(2),1,1)
1998	dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2),1,1)
1999
2000	! look for closest fine grid point from selected coarse grid point
2001	il_iind(:)=MAXLOC( tl_lon1%d_value(:,:,1,1), &
2002	& tl_lon1%d_value(:,:,1,1) <= dl_lon0)
2003
2004	il_jind(:)=MAXLOC( tl_lat1%d_value(:,:,1,1), &
2005	& tl_lat1%d_value(:,:,1,1) <= dl_lat0 )
2006
2007	il_indF(1)=il_iind(1)
2008	il_indF(2)=il_jind(2)
2009
2010	IF( ANY(il_indF(:)==0) )THEN
2011	CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//&
2012	& "fine grid indices. invalid domain" )
2013	ENDIF
2014
2015	dl_lon1=tl_lon1%d_value(il_indF(1),il_indF(2),1,1)
2016	dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2),1,1)
2017
2018	!2-1 check i-direction refinement factor
2019	DO ji=1,MIN(3,il_imid1)
2020
2021	IF( il_indF(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN
2022	CALL logger_debug("GRID CHECK COINCIDENCE: tl_lon1%t_dim(1)%i_len "//TRIM(fct_str(tl_lon1%t_dim(1)%i_len)))
2023	CALL logger_debug("GRID CHECK COINCIDENCE: il_indF(1)+jiid_rhoi+1 "//TRIM(fct_str(il_indF(1)+jiid_rho(jp_I)+1)))
2024	CALL logger_debug("GRID CHECK COINCIDENCE: il_indF(1) "//TRIM(fct_str(il_indF(1))))
2025	CALL logger_debug("GRID CHECK COINCIDENCE: id_rhoi "//TRIM(fct_str(id_rho(jp_I))))
2026	CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
2027	& " to check i-direction refinement factor ")
2028	EXIT
2029	ELSE
2030	dl_lon1=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I),il_indF(2),1,1)
2031	dl_lon0=tl_lon0%d_value(il_indC(1)+ji,il_indC(2),1,1)
2032
2033	dl_lon1p=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I)+1,il_indF(2),1,1)
2034
2035	dl_dlon=ABS(dl_lon1p-dl_lon1)*1.e-3
2036
2037	SELECT CASE(MOD(id_rho(jp_I),2))
2038
2039	CASE(0)
2040
2041	IF( dl_lon1 >= dl_lon0 .OR. dl_lon0 >= dl_lon1p )THEN
2042	ll_coincidence=.FALSE.
2043	CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
2044	& "i-direction refinement factor ("//&
2045	& TRIM(fct_str(id_rho(jp_I)))//&
2046	& ") between fine grid and coarse grid ")
2047	ENDIF
2048
2049	CASE DEFAULT
2050
2051	IF( ABS(dl_lon1 - dl_lon0) > dl_dlon )THEN
2052	ll_coincidence=.FALSE.
2053	CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
2054	& "i-direction refinement factor ("//&
2055	& TRIM(fct_str(id_rho(jp_I)))//&
2056	& ") between fine grid and coarse grid ")
2057	ENDIF
2058
2059	END SELECT
2060	ENDIF
2061
2062	ENDDO
2063
2064	!2-2 check j-direction refinement factor
2065	DO jj=1,MIN(3,il_jmid1)
2066
2067	IF( il_indF(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN
2068	CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
2069	& " to check j-direction refinement factor ")
2070	EXIT
2071	ELSE
2072	dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J),1,1)
2073	dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2)+jj,1,1)
2074
2075	dl_lat1p=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J)+1,1,1)
2076
2077	dl_dlat=ABS(dl_lat1p-dl_lat1)*1.e-3
2078
2079	SELECT CASE(MOD(id_rho(jp_J),2))
2080
2081	CASE(0)
2082
2083	IF( dl_lat1 >= dl_lat0 .OR. dl_lat0 >= dl_lat1p )THEN
2084	ll_coincidence=.FALSE.
2085	CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
2086	& "j-direction refinement factor ("//&
2087	& TRIM(fct_str(id_rho(jp_J)))//&
2088	& ") between fine grid and coarse grid ")
2089	ENDIF
2090
2091	CASE DEFAULT
2092
2093	IF( ABS(dl_lat1-dl_lat0) > dl_dlat )THEN
2094	ll_coincidence=.FALSE.
2095	CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
2096	& "j-direction refinement factor ("//&
2097	& TRIM(fct_str(id_rho(jp_J)))//&
2098	& ") between fine grid and coarse grid ")
2099	ENDIF
2100
2101	END SELECT
2102	ENDIF
2103
2104	ENDDO
2105
2106	IF( .NOT. ll_coincidence )THEN
2107	CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
2108	& "between fine and coarse grid: "//&
2109	& "invalid refinement factor" )
2110	ENDIF
2111
2112	END SUBROUTINE grid_check_coincidence
2113	!> @endcode
2114	!-------------------------------------------------------------------
2115	!> @brief This function check that fine grid is
2116	!> inside coarse grid
2117	!
2118	!> @details
2119	!>
2120	!> @note deltalon and delatlat are used only to avoid issue due to
2121	!> cubic interpolation approximation on the firsts grid points
2122	!
2123	!> @author J.Paul
2124	!> - Nov, 2013- Initial Version
2125	!
2126	!> @param[in] dd_lon0 : table of coarse grid longitude
2127	!> @param[in] dd_lat0 : table of coarse grid latitude
2128	!> @param[in] dd_lon1 : table of fine grid longitude
2129	!> @param[in] dd_lat1 : table of fine grid latitude
2130	!> @return logical, fine grid is inside coarse grid
2131	!-------------------------------------------------------------------
2132	!> @code
2133	FUNCTION grid__check_corner(dd_lon0, dd_lat0, &
2134	& dd_lon1, dd_lat1 )
2135	IMPLICIT NONE
2136	! Argument
2137	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
2138	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
2139	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon1
2140	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat1
2141
2142	! function
2143	LOGICAL :: grid__check_corner
2144
2145	! local variable
2146	INTEGER(i4), DIMENSION(2) :: il_shape0
2147	INTEGER(i4), DIMENSION(2) :: il_shape1
2148
2149	INTEGER(i4) :: il_imin0
2150	INTEGER(i4) :: il_jmin0
2151	INTEGER(i4) :: il_imax0
2152	INTEGER(i4) :: il_jmax0
2153
2154	INTEGER(i4) :: il_imin1
2155	INTEGER(i4) :: il_jmin1
2156	INTEGER(i4) :: il_imax1
2157	INTEGER(i4) :: il_jmax1
2158
2159	REAL(dp) :: dl_lon0
2160	REAL(dp) :: dl_lat0
2161
2162	REAL(dp) :: dl_lon1
2163	REAL(dp) :: dl_lat1
2164
2165	REAL(dp) :: dl_dlon
2166	REAL(dp) :: dl_dlat
2167	! loop indices
2168	!----------------------------------------------------------------
2169
2170	! init
2171	grid__check_corner=.TRUE.
2172
2173	il_shape0=SHAPE(dd_lon0(:,:))
2174	il_shape1=SHAPE(dd_lon1(:,:))
2175
2176	!1- check if fine grid inside coarse grid domain
2177	il_imin0=1 ; il_imax0=il_shape0(1)
2178	il_jmin0=1 ; il_jmax0=il_shape0(2)
2179
2180	il_imin1=1 ; il_imax1=il_shape1(1)
2181	il_jmin1=1 ; il_jmax1=il_shape1(2)
2182
2183	! check lower left corner
2184	dl_lon0 = dd_lon0(il_imin0, il_jmin0 )
2185	dl_lat0 = dd_lat0(il_imin0, il_jmin0 )
2186
2187	dl_lon1 = dd_lon1(il_imin1, il_jmin1)
2188	dl_lat1 = dd_lat1(il_imin1, il_jmin1)
2189
2190	dl_dlon=ABS(dd_lon1(il_imin1+1,il_jmin1 )-dl_lon1)*1.e-3
2191	dl_dlat=ABS(dd_lat1(il_imin1 ,il_jmin1+1)-dl_lat1)*1.e-3
2192
2193	IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 < dl_lon0 ) .OR. &
2194	& (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 < dl_lat0 ) )THEN
2195
2196	CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower left "//&
2197	& "corner not north east of coarse grid (imin,jmin) ")
2198	CALL logger_debug(" fine grid lower left ( "//&
2199	& TRIM(fct_str(dl_lon1))//","//&
2200	& TRIM(fct_str(dl_lat1))//")" )
2201	CALL logger_debug(" coarse grid lower left ( "//&
2202	& TRIM(fct_str(dl_lon0))//","//&
2203	& TRIM(fct_str(dl_lat0))//")" )
2204	grid__check_corner=.FALSE.
2205
2206	ENDIF
2207
2208	! check upper left corner
2209	dl_lon0 = dd_lon0(il_imin0, il_jmax0 )
2210	dl_lat0 = dd_lat0(il_imin0, il_jmax0 )
2211
2212	dl_lon1 = dd_lon1(il_imin1, il_jmax1)
2213	dl_lat1 = dd_lat1(il_imin1, il_jmax1)
2214
2215	dl_dlon=ABS(dd_lon1(il_imin1+1,il_jmax1 )-dl_lon1)*1.e-3
2216	dl_dlat=ABS(dd_lat1(il_imin1 ,il_jmax1-1)-dl_lat1)*1.e-3
2217
2218	IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 < dl_lon0) .OR. &
2219	& (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 > dl_lat0) )THEN
2220
2221	CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper left "//&
2222	& "corner not south east of coarse grid (imin,jmax) ")
2223	CALL logger_debug(" fine grid upper left ("//&
2224	& TRIM(fct_str(dl_lon1))//","//&
2225	& TRIM(fct_str(dl_lat1))//")")
2226	CALL logger_debug(" coasre grid upper left ("//&
2227	& TRIM(fct_str(dl_lon0))//","//&
2228	& TRIM(fct_str(dl_lat0))//")")
2229	grid__check_corner=.FALSE.
2230
2231	ENDIF
2232
2233	! check lower right corner
2234	dl_lon0 = dd_lon0(il_imax0, il_jmin0 )
2235	dl_lat0 = dd_lat0(il_imax0, il_jmin0 )
2236
2237	dl_lon1 = dd_lon1(il_imax1, il_jmin1)
2238	dl_lat1 = dd_lat1(il_imax1, il_jmin1)
2239
2240	dl_dlon=ABS(dd_lon1(il_imax1-1,il_jmin1 )-dl_lon1)*1.e-3
2241	dl_dlat=ABS(dd_lat1(il_imax1 ,il_jmin1+1)-dl_lat1)*1.e-3
2242
2243	IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 > dl_lon0) .OR. &
2244	& (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 < dl_lat0) )THEN
2245
2246	CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower right "//&
2247	& "corner not north west west of coarse grid (imax,jmin) ")
2248	CALL logger_debug(" fine grid lower right ( "//&
2249	& TRIM(fct_str(dl_lon1))//","//&
2250	& TRIM(fct_str(dl_lat1))//")" )
2251	CALL logger_debug(" coarse grid lower right ( "//&
2252	& TRIM(fct_str(dl_lon0))//","//&
2253	& TRIM(fct_str(dl_lat0))//")" )
2254	grid__check_corner=.FALSE.
2255
2256	ENDIF
2257
2258	! check upper right corner
2259	dl_lon0 = dd_lon0(il_imax0, il_jmax0 )
2260	dl_lat0 = dd_lat0(il_imax0, il_jmax0 )
2261
2262	dl_lon1 = dd_lon1(il_imax1, il_jmax1)
2263	dl_lat1 = dd_lat1(il_imax1, il_jmax1)
2264
2265	dl_dlon=ABS(dd_lon1(il_imax1-1,il_jmax1 )-dl_lon1)*1.e-3
2266	dl_dlat=ABS(dd_lat1(il_imax1 ,il_jmax1-1)-dl_lat1)*1.e-3
2267
2268	IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 > dl_lon0) .OR. &
2269	& (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 > dl_lat0) )THEN
2270
2271	CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper right "//&
2272	& "corner not south west of coarse grid (imax,jmax) ")
2273	CALL logger_debug(" fine grid upper right ( "//&
2274	& TRIM(fct_str(dl_lon1))//","//&
2275	& TRIM(fct_str(dl_lat1))//")" )
2276	CALL logger_debug(" fine imax1 jmax1 ( "//&
2277	& TRIM(fct_str(il_imax1))//","//&
2278	& TRIM(fct_str(il_jmax1))//")" )
2279	CALL logger_debug(" coarse grid upper right ( "//&
2280	& TRIM(fct_str(dl_lon0))//","//&
2281	& TRIM(fct_str(dl_lat0))//")" )
2282	CALL logger_debug(" fine imax0 jmax0 ( "//&
2283	& TRIM(fct_str(il_imax0))//","//&
2284	& TRIM(fct_str(il_jmax0))//")" )
2285	grid__check_corner=.FALSE.
2286
2287	ENDIF
2288
2289	END FUNCTION grid__check_corner
2290	!> @endcode
2291	!-------------------------------------------------------------------
2292	!> @brief This function check that fine grid latitude are
2293	!> inside coarse grid latitude
2294	!
2295	!> @details
2296	!
2297	!> @author J.Paul
2298	!> - Nov, 2013- Initial Version
2299	!
2300	!> @param[in] dd_lat0 : table of coarse grid latitude
2301	!> @param[in] dd_lat1 : table of fine grid latitude
2302	!-------------------------------------------------------------------
2303	!> @code
2304	FUNCTION grid__check_lat(dd_lat0, dd_lat1, id_rhoj)
2305	IMPLICIT NONE
2306	! Argument
2307	REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat0
2308	REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat1
2309	INTEGER(i4) , INTENT(IN) :: id_rhoj
2310
2311	! function
2312	LOGICAL :: grid__check_lat
2313
2314	! local variable
2315	INTEGER(i4), DIMENSION(1) :: il_shape0
2316	INTEGER(i4), DIMENSION(1) :: il_shape1
2317
2318	INTEGER(i4) :: il_jmin0
2319	INTEGER(i4) :: il_jmax0
2320
2321	INTEGER(i4) :: il_jmin1
2322	INTEGER(i4) :: il_jmax1
2323
2324	REAL(dp) :: dl_dlat
2325	! loop indices
2326	!----------------------------------------------------------------
2327
2328	! init
2329	grid__check_lat=.TRUE.
2330
2331	il_shape0(:)=SHAPE(dd_lat0(:))
2332	il_shape1(:)=SHAPE(dd_lat1(:))
2333
2334	!1- check if fine grid inside coarse grid domain
2335	il_jmin0=1+1 ; il_jmax0=il_shape0(1)-1
2336
2337	il_jmin1=1+id_rhoj ; il_jmax1=il_shape1(1)-id_rhoj
2338
2339	dl_dlat=ABS(dd_lat1(il_jmin1+1)-dd_lat1(il_jmin1))*1.e-3
2340
2341	! check lower left fine grid
2342	IF( ABS(dd_lat1(il_jmin1)-dd_lat0(il_jmin0)) > dl_dlat .AND. &
2343	& dd_lat1(il_jmin1) < dd_lat0(il_jmin0) )THEN
2344
2345	CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower point"//&
2346	& " not north of coarse grid (jmin) ")
2347	CALL logger_debug(" fine grid lower point ( "//&
2348	& TRIM(fct_str(dd_lat1(il_jmin1)))//")" )
2349	CALL logger_debug(" coarse grid lower point ( "//&
2350	& TRIM(fct_str(dd_lat0(il_jmin0)))//")" )
2351	grid__check_lat=.FALSE.
2352
2353	ENDIF
2354
2355	dl_dlat=ABS(dd_lat1(il_jmax1-1)-dd_lat1(il_jmax1))*1.e-3
2356
2357	! check upper left fine grid
2358	IF( ABS(dd_lat1(il_jmax1)-dd_lat0(il_jmax0)) > dl_dlat .AND. &
2359	& dd_lat1(il_jmax1) > dd_lat0(il_jmax0) )THEN
2360
2361	CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper point"//&
2362	& " not south of coarse grid (jmax) ")
2363	CALL logger_debug(" fine grid upper point ("//&
2364	& TRIM(fct_str(dd_lat1(il_jmax1)))//")")
2365	CALL logger_debug(" coasre grid upper point ("//&
2366	& TRIM(fct_str(dd_lat0(il_jmax0)))//")")
2367	grid__check_lat=.FALSE.
2368
2369	ENDIF
2370
2371	END FUNCTION grid__check_lat
2372	!> @endcode
2373	!-------------------------------------------------------------------
2374	!> @brief
2375	!> This subroutine add ghost cell at boundaries.
2376	!>
2377	!> @author J.Paul
2378	!> - Nov, 2013-Initial version
2379	!
2380	!> @param[inout] td_var : table of variable structure
2381	!> @param[in] id_ighost : i-direction ghost cell factor
2382	!> @param[in] id_jghost : j-direction ghost cell factor
2383	!-------------------------------------------------------------------
2384	!> @code
2385	SUBROUTINE grid_add_ghost(td_var, id_ighost, id_jghost)
2386	IMPLICIT NONE
2387	! Argument
2388	TYPE(TVAR) , INTENT(INOUT) :: td_var
2389	INTEGER(i4), INTENT(IN ) :: id_ighost
2390	INTEGER(i4), INTENT(IN ) :: id_jghost
2391
2392	! local variable
2393	INTEGER(i4) :: il_imin
2394	INTEGER(i4) :: il_jmin
2395	INTEGER(i4) :: il_imax
2396	INTEGER(i4) :: il_jmax
2397
2398	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
2399
2400	TYPE(TVAR) :: tl_var
2401
2402	! loop indices
2403	!----------------------------------------------------------------
2404
2405	IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
2406
2407	CALL logger_warn( "ADD GHOST: dimension change in variable "//&
2408	& TRIM(td_var%c_name) )
2409
2410	! copy variable
2411	tl_var=td_var
2412
2413	CALL var_del_value(td_var)
2414
2415	! compute indice to fill center
2416	il_imin=1+id_ighost*ig_ghost
2417	il_jmin=1+id_jghost*ig_ghost
2418
2419	il_imax=il_imin+tl_var%t_dim(1)%i_len-1
2420	il_jmax=il_jmin+tl_var%t_dim(2)%i_len-1
2421
2422	! compute new dimension
2423	td_var%t_dim(1)%i_len = tl_var%t_dim(1)%i_len + 2id_ighostig_ghost
2424	td_var%t_dim(2)%i_len = tl_var%t_dim(2)%i_len + 2id_jghostig_ghost
2425
2426	ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
2427	& td_var%t_dim(2)%i_len, &
2428	& td_var%t_dim(3)%i_len, &
2429	& td_var%t_dim(4)%i_len) )
2430
2431	dl_value(:,:,:,:)=tl_var%d_fill
2432
2433	dl_value(il_imin:il_imax, &
2434	& il_jmin:il_jmax, &
2435	& :,:) = tl_var%d_value(:,:,:,:)
2436
2437	! add variable value
2438	CALL var_add_value(td_var,dl_value(:,:,:,:))
2439
2440	! save variable type
2441	td_var%i_type=tl_var%i_type
2442
2443	DEALLOCATE( dl_value )
2444
2445	CALL var_clean(tl_var)
2446
2447	ENDIF
2448
2449	END SUBROUTINE grid_add_ghost
2450	!> @endcode
2451	!-------------------------------------------------------------------
2452	!> @brief
2453	!> This subroutine delete ghost cell at boundaries.
2454	!>
2455	!> @author J.Paul
2456	!> - Nov, 2013-Initial version
2457	!
2458	!> @param[inout] td_var : table of variable structure
2459	!> @param[in] id_ighost : i-direction ghost cell factor
2460	!> @param[in] id_jghost : j-direction ghost cell factor
2461	!-------------------------------------------------------------------
2462	!> @code
2463	SUBROUTINE grid_del_ghost(td_var, id_ighost, id_jghost)
2464	IMPLICIT NONE
2465	! Argument
2466	TYPE(TVAR) , INTENT(INOUT) :: td_var
2467	INTEGER(i4), INTENT(IN ) :: id_ighost
2468	INTEGER(i4), INTENT(IN ) :: id_jghost
2469
2470	! local variable
2471	INTEGER(i4) :: il_imin
2472	INTEGER(i4) :: il_jmin
2473	INTEGER(i4) :: il_imax
2474	INTEGER(i4) :: il_jmax
2475
2476	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
2477
2478	TYPE(TVAR) :: tl_var
2479
2480	! loop indices
2481	!----------------------------------------------------------------
2482
2483	IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
2484
2485	CALL logger_warn( "DEL GHOST: dimension change in variable "//&
2486	& TRIM(td_var%c_name) )
2487
2488	! copy variable
2489	tl_var=td_var
2490
2491	CALL var_del_value(td_var)
2492
2493	! compute indice to get center
2494	il_imin=1+id_ighost*ig_ghost
2495	il_jmin=1+id_jghost*ig_ghost
2496
2497	il_imax=tl_var%t_dim(1)%i_len-id_ighost*ig_ghost
2498	il_jmax=tl_var%t_dim(2)%i_len-id_jghost*ig_ghost
2499
2500	! compute new dimension
2501	td_var%t_dim(1)%i_len = tl_var%t_dim(1)%i_len - 2id_ighostig_ghost
2502	td_var%t_dim(2)%i_len = tl_var%t_dim(2)%i_len - 2id_jghostig_ghost
2503
2504	ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
2505	& td_var%t_dim(2)%i_len, &
2506	& td_var%t_dim(3)%i_len, &
2507	& td_var%t_dim(4)%i_len) )
2508
2509	dl_value(:,:,:,:)=tl_var%d_fill
2510
2511	dl_value(:,:,:,:) = tl_var%d_value(il_imin:il_imax, &
2512	& il_jmin:il_jmax, &
2513	& :,:)
2514
2515	! add variable value
2516	CALL var_add_value(td_var,dl_value(:,:,:,:))
2517
2518	! save variable type
2519	td_var%i_type=tl_var%i_type
2520
2521	DEALLOCATE( dl_value )
2522
2523	CALL var_clean(tl_var)
2524
2525	ENDIF
2526
2527	END SUBROUTINE grid_del_ghost
2528	!> @endcode
2529	!-------------------------------------------------------------------
2530	!> @brief This subroutine fill small closed sea with fill value.
2531	!
2532	!> @details
2533	!> the minimum size (nbumber of point) of closed sea to be kept could be
2534	!> sepcify with id_minsize.
2535	!> By default only the biggest sea is preserve.
2536	!
2537	!> @author J.Paul
2538	!> - Nov, 2013- Initial Version
2539	!
2540	!> @param[inout] td_var : variable structure
2541	!> @param[in] id_mask : domain mask (from grid_split_domain)
2542	!> @param[in] id_minsize : minimum size of sea to be kept
2543	!-------------------------------------------------------------------
2544	!> @code
2545	SUBROUTINE grid_fill_small_dom(td_var, id_mask, id_minsize)
2546	IMPLICIT NONE
2547	! Argument
2548	TYPE(TVAR) , INTENT(INOUT) :: td_var
2549	INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_mask
2550	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_minsize
2551
2552	! local variable
2553	INTEGER(i4) :: il_ndom
2554	INTEGER(i4) :: il_minsize
2555	INTEGER(i4), DIMENSION(2) :: il_shape
2556	INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp
2557
2558	! loop indices
2559	INTEGER(i4) :: ji
2560	INTEGER(i4) :: jk
2561	INTEGER(i4) :: jl
2562	!----------------------------------------------------------------
2563
2564	il_shape(:)=SHAPE(id_mask(:,:))
2565	IF( ANY(il_shape(:) /= td_var%t_dim(1:2)%i_len) )THEN
2566	CALL logger_error("GRID FILL SMALL DOM: variable and mask "//&
2567	& "dimension differ")
2568	ELSE
2569
2570	il_ndom=MINVAL(id_mask(:,:))
2571
2572	ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
2573	il_tmp(:,:)=0
2574	DO ji=-1,il_ndom,-1
2575	WHERE( id_mask(:,:)==ji )
2576	il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji
2577	END WHERE
2578	ENDDO
2579
2580	il_minsize=MAXVAL(il_tmp(:,:))
2581	IF( PRESENT(id_minsize) ) il_minsize=id_minsize
2582
2583	DO jl=1,td_var%t_dim(4)%i_len
2584	DO jk=1,td_var%t_dim(3)%i_len
2585	WHERE( il_tmp(:,:) < il_minsize )
2586	td_var%d_value(:,:,jk,jl)=td_var%d_fill
2587	END WHERE
2588	ENDDO
2589	ENDDO
2590
2591	DEALLOCATE( il_tmp )
2592
2593	ENDIF
2594
2595	END SUBROUTINE grid_fill_small_dom
2596	!> @endcode
2597	!-------------------------------------------------------------------
2598	!> @brief This subroutine compute closed sea domain.
2599	!
2600	!> @details
2601	!> to each domain is associated a negative value id (from -1 to ...)
2602	!
2603	!> @author J.Paul
2604	!> - Nov, 2013- Initial Version
2605	!
2606	!> @param[in] td_var : variable strucutre
2607	!> @param[in] id_level : level
2608	!> @return domain mask
2609	!-------------------------------------------------------------------
2610	!> @code
2611	FUNCTION grid_split_domain(td_var, id_level)
2612	IMPLICIT NONE
2613	! Argument
2614	TYPE(TVAR) , INTENT(IN) :: td_var
2615	INTEGER(i4), INTENT(IN), OPTIONAL :: id_level
2616
2617	! function
2618	INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &
2619	& td_var%t_dim(2)%i_len ) :: grid_split_domain
2620
2621	! local variable
2622	INTEGER(i4) :: il_domid
2623	INTEGER(i4) :: il_level
2624	INTEGER(i4), DIMENSION(2) :: il_shape
2625	INTEGER(i4), DIMENSION(2) :: il_ind
2626	INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask
2627	INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp
2628
2629	LOGICAL :: ll_full
2630
2631	! loop indices
2632	INTEGER(i4) :: ji
2633	INTEGER(i4) :: jim
2634	INTEGER(i4) :: jip
2635	INTEGER(i4) :: jj
2636	INTEGER(i4) :: jjm
2637	INTEGER(i4) :: jjp
2638	!----------------------------------------------------------------
2639	il_level=1
2640	IF( PRESENT(id_level) ) il_level=id_level
2641
2642	! init
2643	il_domid=-1
2644
2645	il_shape(:)=td_var%t_dim(1:2)%i_len
2646	ALLOCATE( il_mask(il_shape(1),il_shape(2)) )
2647	il_mask(:,:)=0
2648	WHERE( td_var%d_value(:,:,il_level,1)/=td_var%d_fill ) il_mask(:,:)=1
2649
2650	il_ind(:)=MAXLOC( il_mask(:,:) )
2651	DO WHILE( il_mask(il_ind(1), il_ind(2)) == 1 )
2652
2653	il_mask(il_ind(1),il_ind(2))=il_domid
2654	ll_full=.FALSE.
2655
2656	ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
2657
2658	DO WHILE( .NOT. ll_full )
2659	il_tmp(:,:)=0
2660
2661	ll_full=.TRUE.
2662	DO jj=1,il_shape(2)
2663	DO ji=1,il_shape(1)
2664	IF( il_mask(ji,jj)==il_domid )THEN
2665	jim=MAX(1,ji-1) ; jip=MIN(il_shape(1),ji+1)
2666	jjm=MAX(1,jj-1) ; jjp=MIN(il_shape(2),jj+1)
2667
2668	WHERE( il_mask(jim:jip,jjm:jjp)==1 )
2669	il_mask(jim:jip,jjm:jjp)=il_domid
2670	il_tmp(jim:jip,jjm:jjp)=1
2671	END WHERE
2672	ENDIF
2673	ENDDO
2674	ENDDO
2675	IF( SUM(il_tmp(:,:))/=0 ) ll_full=.FALSE.
2676
2677	ENDDO
2678	DEALLOCATE( il_tmp )
2679
2680	il_ind(:)=MAXLOC( il_mask(:,:) )
2681	il_domid=il_domid-1
2682
2683	ENDDO
2684
2685	! save result
2686	grid_split_domain(:,:)=il_mask(:,:)
2687
2688	DEALLOCATE( il_mask )
2689
2690	CALL logger_info("GRID SPLIT DOMAIN: "//TRIM( fct_str(ABS(il_domid+1)) )//&
2691	& " domain found" )
2692
2693	END FUNCTION grid_split_domain
2694	!> @endcode
2695	! !-------------------------------------------------------------------
2696	! !> @brief This function
2697	! !
2698	! !> @details
2699	! !
2700	! !> @author J.Paul
2701	! !> - Nov, 2013- Initial Version
2702	! !
2703	! !> @param[in]
2704	! !-------------------------------------------------------------------
2705	! !> @code
2706	! FUNCTION grid_()
2707	! IMPLICIT NONE
2708	! ! Argument
2709	! ! function
2710	! ! local variable
2711	! ! loop indices
2712	! !----------------------------------------------------------------
2713	!
2714	! END FUNCTION grid_
2715	! !> @endcode
2716	! !-------------------------------------------------------------------
2717	! !> @brief This subroutine
2718	! !
2719	! !> @details
2720	! !
2721	! !> @author J.Paul
2722	! !> - Nov, 2013- Initial Version
2723	! !
2724	! !> @param[in]
2725	! !-------------------------------------------------------------------
2726	! !> @code
2727	! SUBROUTINE grid_()
2728	! IMPLICIT NONE
2729	! ! Argument
2730	! ! local variable
2731	! ! loop indices
2732	! !----------------------------------------------------------------
2733	!
2734	! END SUBROUTINE grid_
2735	! !> @endcode
2736	END MODULE grid
2737

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

Download in other formats: