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/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid.f90 @ 7153

Last change on this file since 7153 was 7153, checked in by jpaul, 7 years ago

see ticket #1781

File size: 214.6 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! MODULE: grid
6!
7! DESCRIPTION:
8!> @brief This module is grid manager.
9!>
10!> @details
11!>    to get NEMO pivot point index:<br/>
12!> @code
13!>    il_pivot=grid_get_pivot(td_file)
14!> @endcode
15!>       - il_pivot is NEMO pivot point index F(0), T(1)
16!>       - td_file is mpp structure
17!>
18!>    to get NEMO periodicity index:<br/>
19!> @code
20!>    il_perio=grid_get_perio(td_file)
21!> @endcode
22!>       - il_perio is NEMO periodicity index (0,1,2,3,4,5,6)
23!>       - td_file is mpp structure
24!>
25!>    to check domain validity:<br/>
26!> @code
27!>    CALL grid_check_dom(td_coord, id_imin, id_imax, id_jmin, id_jmax)
28!> @endcode
29!>       - td_coord is coordinates mpp structure
30!>       - id_imin is i-direction lower left  point indice
31!>       - id_imax is i-direction upper right point indice
32!>       - id_jmin is j-direction lower left  point indice
33!>       - id_jmax is j-direction upper right point indice
34!>
35!>    to get closest coarse grid indices of fine grid domain:<br/>
36!> @code
37!>    il_index(:,:)=grid_get_coarse_index(td_coord0, td_coord1,
38!>                                      [id_rho,] [cd_point])
39!> @endcode
40!>    or
41!> @code
42!>    il_index(:,:)=grid_get_coarse_index(td_lon0, td_lat0, td_coord1,
43!>                                      [id_rho,] [cd_point])
44!> @endcode
45!>    or
46!> @code
47!>    il_index(:,:)=grid_get_coarse_index(td_coord0, td_lon1, td_lat1,
48!>                                      [id_rho,] [cd_point])
49!> @endcode
50!>    or
51!> @code
52!>    il_index(:,:)=grid_get_coarse_index(td_lon0, td_lat0, td_lon1, td_lat1,
53!>                                      [id_rho,] [cd_point])
54!> @endcode
55!>       - il_index(:,:) is  coarse grid indices (/ (/ imin0, imax0 /),
56!> (/ jmin0, jmax0 /) /)
57!>       - td_coord0 is coarse grid coordinate mpp structure
58!>       - td_coord1 is fine grid coordinate mpp structure
59!>       - td_lon0 is coarse grid longitude variable structure
60!>       - td_lat0 is coarse grid latitude  variable structure
61!>       - td_lon1 is fine   grid longitude variable structure
62!>       - td_lat1 is fine   grid latitude  variable structure
63!>       - id_rho is array of refinment factor (default 1)
64!>       - cd_point is Arakawa grid point (default 'T')
65!>
66!>    to know if grid is global:<br/>
67!> @code
68!>    ll_global=grid_is_global(td_lon, td_lat)
69!> @endcode
70!>       - td_lon is longitude variable structure
71!>       - td_lat is latitude variable structure
72!>   
73!>    to know if grid contains north fold:<br/>
74!> @code
75!>    ll_north=grid_is_north_fold(td_lat)
76!> @endcode
77!>       - td_lat is latitude variable structure   
78!>
79!>    to get coarse grid indices of the closest point from one fine grid
80!> point:<br/>
81!> @code
82!>    il_index(:)=grid_get_closest(dd_lon0(:,:), dd_lat0(:,:), dd_lon1, dd_lat1
83!>                                 [,dd_fill] [,cd_pos])
84!> @endcode
85!>       - il_index(:) is  coarse grid indices (/ i0, j0 /)
86!>       - dd_lon0 is coarse grid array of longitude value (real(8))
87!>       - dd_lat0 is coarse grid array of latitude  value (real(8))
88!>       - dd_lon1 is fine grid longitude value (real(8))
89!>       - dd_lat1 is fine grid latitude  value (real(8))
90!>       - dd_fill
91!>       - cd_pos
92!>
93!>    to compute distance between a point A and grid points:<br/>
94!> @code
95!>    il_dist(:,:)=grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA)
96!> @endcode
97!>       - il_dist(:,:) is array of distance between point A and grid points
98!>       - dd_lon is array of longitude value (real(8))
99!>       - dd_lat is array of longitude value (real(8))
100!>       - dd_lonA is longitude of point A (real(8))
101!>       - dd_latA is latitude  of point A (real(8))
102!>
103!>    to get offset between fine grid and coarse grid:<br/>
104!> @code
105!>    il_offset(:,:)=grid_get_fine_offset(td_coord0,
106!>                                        id_imin0, id_jmin0, id_imax0, id_jmax0,
107!>                                        td_coord1
108!>                                        [,id_rho] [,cd_point])
109!> @endcode
110!>    or
111!> @code
112!>    il_offset(:,:)=grid_get_fine_offset(dd_lon0, dd_lat0,
113!>                                        id_imin0, id_jmin0,id_imax0, id_jmax0,
114!>                                        td_coord1
115!>                                        [,id_rho] [,cd_point])
116!> @endcode
117!>    or
118!> @code
119!>    il_offset(:,:)=grid_get_fine_offset(td_coord0,
120!>                                        id_imin0, id_jmin0, id_imax0, id_jmax0,
121!>                                        dd_lon1, dd_lat1
122!>                                        [,id_rho] [,cd_point])
123!> @endcode
124!>    or
125!> @code
126!>    il_offset(:,:)=grid_get_fine_offset(dd_lon0, dd_lat0,
127!>                                        id_imin0, id_jmin0, id_imax0, id_jmax0,
128!>                                        dd_lon1, dd_lat1
129!>                                        [,id_rho] [,cd_point])
130!> @endcode
131!>       - il_offset(:,:) is offset array
132!>    (/ (/ i_offset_left, i_offset_right /), (/ j_offset_lower, j_offset_upper /) /)
133!>       - td_coord0 is coarse grid coordinate mpp structure
134!>       - dd_lon0  is coarse grid longitude array (real(8))
135!>       - dd_lat0  is coarse grid latitude  array (real(8))
136!>       - id_imin0 is coarse grid lower left  corner i-indice of fine grid
137!> domain
138!>       - id_jmin0 is coarse grid lower left  corner j-indice of fine grid
139!> domain
140!>       - id_imax0 is coarse grid upper right corner i-indice of fine grid
141!> domain
142!>       - id_jmax0 is coarse grid upper right corner j-indice of fine grid
143!> domain
144!>       - td_coord1 is fine grid coordinate mpp structure
145!>       - dd_lon1  is fine   grid longitude array (real(8))
146!>       - dd_lat1  is fine   grid latitude  array (real(8))
147!>       - id_rho is array of refinment factor (default 1)
148!>       - cd_point is Arakawa grid point (default 'T')
149!>
150!>    to check fine and coarse grid coincidence:<br/>
151!> @code
152!>    CALL grid_check_coincidence(td_coord0, td_coord1,
153!>                                id_imin0, id_imax0, id_jmin0, id_jmax0
154!>                                ,id_rho)
155!> @endcode
156!>       - td_coord0 is coarse grid coordinate mpp structure
157!>       - td_coord1 is fine   grid coordinate mpp structure
158!>       - id_imin0  is coarse grid lower left  corner i-indice of fine grid
159!> domain
160!>       - id_imax0  is coarse grid upper right corner i-indice of fine grid
161!> domain
162!>       - id_jmin0  is coarse grid lower left  corner j-indice of fine grid
163!> domain
164!>       - id_jmax0  is coarse grid upper right corner j-indice of fine grid
165!> domain
166!>       - id_rho    is array of refinement factor
167!>
168!>    to add ghost cell at boundaries:<br/>
169!> @code
170!>    CALL grid_add_ghost(td_var, id_ghost)
171!> @endcode
172!>       - td_var is array of variable structure
173!>       - id_ghost is 2D array of ghost cell factor
174!>
175!>    to delete ghost cell at boundaries:<br/>
176!> @code
177!>    CALL grid_del_ghost(td_var, id_ghost)
178!> @endcode
179!>       - td_var is array of variable structure
180!>       - id_ghost is 2D array of ghost cell factor
181!>
182!>    to get ghost cell factor (use or not):<br/>
183!> @code
184!>    il_factor(:)= grid_get_ghost( td_var )
185!> @endcode
186!>    or
187!> @code
188!>    il_factor(:)= grid_get_ghost( td_mpp )
189!> @endcode
190!>       - il_factor(:) is  array of ghost cell factor (0 or 1)
191!>       - td_var  is variable structure
192!>       - td_mpp is mpp sturcture
193!>
194!>    to compute closed sea domain:<br/>
195!> @code
196!>    il_mask(:,:)=grid_split_domain(td_var, [id_level])
197!> @endcode
198!>       - il_mask(:,:) is domain mask
199!>       - td_var is variable strucutre
200!>       - id_level is level to be used [optional]
201!>
202!>    to fill small closed sea with _FillValue:<br/>
203!> @code
204!>    CALL grid_fill_small_dom(td_var, id_mask, [id_minsize])
205!> @endcode
206!>       - td_var  is variable structure
207!>       - id_mask is domain mask (from grid_split_domain)
208!>       - id_minsize is minimum size of sea to be kept [optional]
209!>
210!> @author
211!> J.Paul
212! REVISION HISTORY:
213!> @date November, 2013 - Initial Version
214!> @date September, 2014
215!> - add header
216!> @date October, 2014
217!> - use mpp file structure instead of file
218!> @date February, 2015
219!> - add function grid_fill_small_msk to fill small domain inside bigger one
220!> @date February, 2016
221!> - improve way to check coincidence (bug fix)
222!> - manage grid cases for T,U,V or F point, with even or odd refinment (bug fix)
223!> @date April, 2016
224!> - add function to get closest grid point using coarse grid coordinates strucutre
225!
226!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
227!----------------------------------------------------------------------
228MODULE grid
229   USE netcdf
230   USE kind                            ! F90 kind parameter
231   USE fct                             ! basic usefull function
232   USE global                          ! global parameter
233   USE phycst                          ! physical constant
234   USE logger                          ! log file manager
235   USE file                            ! file manager
236   USE att                             ! attribute manager
237   USE var                             ! variable manager
238   USE dim                             ! dimension manager
239   USE iom                             ! I/O manager
240   USE mpp                             ! MPP manager
241   USE dom                             ! domain manager
242   USE iom_mpp                         ! MPP I/O manager
243   USE iom_dom                         ! DOM I/O manager
244   IMPLICIT NONE
245   ! NOTE_avoid_public_variables_if_possible
246
247   ! type and variable
248
249   ! function and subroutine
250   PUBLIC :: grid_get_info             !< get information about mpp global domain (pivot, perio, ew)
251   PUBLIC :: grid_get_pivot            !< get NEMO pivot point index
252   PUBLIC :: grid_get_perio            !< get NEMO periodicity index
253   PUBLIC :: grid_get_ew_overlap       !< get East West overlap
254   PUBLIC :: grid_check_dom            !< check domain validity
255   PUBLIC :: grid_get_coarse_index     !< get closest coarse grid indices of fine grid domain.
256   PUBLIC :: grid_is_global            !< check if grid is global or not
257   PUBLIC :: grid_is_north_fold
258   PUBLIC :: grid_get_closest          !< return closest coarse grid point from another point
259   PUBLIC :: grid_distance             !< compute grid distance to a point
260   PUBLIC :: grid_get_fine_offset      !< get fine grid offset
261   PUBLIC :: grid_check_coincidence    !< check fine and coarse grid coincidence
262   PUBLIC :: grid_add_ghost            !< add ghost cell at boundaries.
263   PUBLIC :: grid_del_ghost            !< delete ghost cell at boundaries.
264   PUBLIC :: grid_get_ghost            !< return ghost cell factor
265   PUBLIC :: grid_split_domain         !< compute closed sea domain
266   PUBLIC :: grid_fill_small_dom       !< fill small closed sea with fill value
267   PUBLIC :: grid_fill_small_msk       !< fill small domain inside bigger one
268
269                                     ! get closest coarse grid indices of fine grid domain
270   PRIVATE :: grid__get_coarse_index_ff ! - using coarse and fine grid coordinates files
271   PRIVATE :: grid__get_coarse_index_cf ! - using coarse grid array of lon,lat and fine grid coordinates files
272   PRIVATE :: grid__get_coarse_index_fc ! - using coarse grid coordinates files, and fine grid array of lon,lat
273   PRIVATE :: grid__get_coarse_index_cc ! - using coarse and fine grid array of lon,lat
274
275                                     ! return closest coarse grid point from another point
276   PRIVATE :: grid__get_closest_str    ! - using coarse grid coordinates strucutre
277   PRIVATE :: grid__get_closest_arr    ! - using coarse grid array of lon,lat
278
279                                     ! get offset between fine and coarse grid
280   PRIVATE :: grid__get_fine_offset_ff ! - using coarse and fine grid coordinates files
281   PRIVATE :: grid__get_fine_offset_cf ! - using coarse grid array of lon,lat and fine grid coordinates files
282   PRIVATE :: grid__get_fine_offset_fc ! - using coarse grid coordinates files, and fine grid array of lon,lat
283   PRIVATE :: grid__get_fine_offset_cc ! - using coarse and fine grid array of lon,lat
284
285                                     ! get information about global domain (pivot, perio, ew)
286   PRIVATE :: grid__get_info_mpp      ! - using mpp files structure
287   PRIVATE :: grid__get_info_file     ! - using files structure
288
289                                     ! get NEMO pivot point index
290   PRIVATE :: grid__get_pivot_mpp      ! - using mpp files structure
291   PRIVATE :: grid__get_pivot_file     ! - using files structure
292   PRIVATE :: grid__get_pivot_var      ! - using variable structure
293   PRIVATE :: grid__get_pivot_varT   ! compute NEMO pivot point index for variable on grid T
294   PRIVATE :: grid__get_pivot_varU   ! compute NEMO pivot point index for variable on grid U
295   PRIVATE :: grid__get_pivot_varV   ! compute NEMO pivot point index for variable on grid V
296   PRIVATE :: grid__get_pivot_varF   ! compute NEMO pivot point index for variable on grid F
297
298                                     ! get NEMO periodicity index
299   PRIVATE :: grid__get_perio_mpp      ! - using mpp files structure
300   PRIVATE :: grid__get_perio_file     ! - using files structure
301   PRIVATE :: grid__get_perio_var      ! - using variable structure
302
303                                     ! get East West overlap
304   PRIVATE :: grid__get_ew_overlap_mpp  ! - using mpp files structure
305   PRIVATE :: grid__get_ew_overlap_file ! - using files structure
306   PRIVATE :: grid__get_ew_overlap_var  ! - using longitude variable structure
307
308                                    ! return ghost cell factor
309   PRIVATE :: grid__get_ghost_mpp      ! - using mpp files structure
310   PRIVATE :: grid__get_ghost_var      ! - using array of lon,lat
311   PRIVATE :: grid__check_corner    ! check that fine grid is inside coarse grid
312   PRIVATE :: grid__check_lat       ! check that fine grid latitude are inside coarse grid latitude
313   
314   INTERFACE  grid_get_info
315      MODULE PROCEDURE grid__get_info_mpp
316      MODULE PROCEDURE grid__get_info_file
317   END INTERFACE grid_get_info
318
319   INTERFACE  grid_get_pivot
320      MODULE PROCEDURE grid__get_pivot_mpp
321      MODULE PROCEDURE grid__get_pivot_file
322      MODULE PROCEDURE grid__get_pivot_var
323   END INTERFACE grid_get_pivot
324
325   INTERFACE  grid_get_perio
326      MODULE PROCEDURE grid__get_perio_mpp
327      MODULE PROCEDURE grid__get_perio_file
328      MODULE PROCEDURE grid__get_perio_var
329   END INTERFACE grid_get_perio
330
331   INTERFACE  grid_get_ew_overlap
332      MODULE PROCEDURE grid__get_ew_overlap_mpp
333      MODULE PROCEDURE grid__get_ew_overlap_file
334      MODULE PROCEDURE grid__get_ew_overlap_var
335   END INTERFACE grid_get_ew_overlap
336
337   INTERFACE  grid_get_ghost
338      MODULE PROCEDURE grid__get_ghost_var
339      MODULE PROCEDURE grid__get_ghost_mpp
340   END INTERFACE  grid_get_ghost
341
342   INTERFACE  grid_get_closest
343      MODULE PROCEDURE grid__get_closest_str
344      MODULE PROCEDURE grid__get_closest_arr
345   END INTERFACE  grid_get_closest
346
347   INTERFACE  grid_get_coarse_index
348      MODULE PROCEDURE grid__get_coarse_index_ff
349      MODULE PROCEDURE grid__get_coarse_index_cf
350      MODULE PROCEDURE grid__get_coarse_index_fc
351      MODULE PROCEDURE grid__get_coarse_index_cc
352   END INTERFACE grid_get_coarse_index
353
354   INTERFACE  grid_get_fine_offset
355      MODULE PROCEDURE grid__get_fine_offset_ff
356      MODULE PROCEDURE grid__get_fine_offset_fc
357      MODULE PROCEDURE grid__get_fine_offset_cf
358      MODULE PROCEDURE grid__get_fine_offset_cc
359   END INTERFACE grid_get_fine_offset   
360
361CONTAINS
362   !-------------------------------------------------------------------
363   !> @brief This subroutine get information about global domain, given file
364   !> strucutre.
365   !>
366   !> @details
367   !> open edge files then:
368   !> - compute NEMO pivot point
369   !> - compute NEMO periodicity
370   !> - compute East West overlap
371   !>
372   !> @note need all processor files to be there
373   !> @author J.Paul
374   !> @date October, 2014 - Initial Version
375   !>
376   !> @param[inout] td_file file structure
377   !-------------------------------------------------------------------
378   SUBROUTINE grid__get_info_file(td_file)
379      IMPLICIT NONE
380      ! Argument     
381      TYPE(TFILE), INTENT(INOUT) :: td_file
382
383      ! local variable
384      INTEGER(i4) :: il_ew
385      INTEGER(i4) :: il_pivot
386      INTEGER(i4) :: il_perio
387      INTEGER(i4) :: il_attid
388
389      TYPE(TATT)  :: tl_att
390
391      TYPE(TFILE) :: tl_file
392
393      ! loop indices
394      INTEGER(i4) :: ji
395      !----------------------------------------------------------------
396      ! intialise
397      il_pivot=-1
398      il_perio=-1
399      il_ew   =-1 
400
401      ! copy structure
402      tl_file=file_copy(td_file)
403
404      ! open file to be used
405      CALL iom_open(tl_file)
406
407      IF( td_file%i_perio >= 0 .AND. td_file%i_perio <= 6 )THEN
408         il_perio=td_file%i_perio
409      ELSE
410         ! look for attribute in file
411         il_attid=att_get_index(tl_file%t_att(:),'periodicity')
412         IF( il_attid /= 0 )THEN
413            il_perio=INT(tl_file%t_att(il_attid)%d_value(1),i4)
414         ENDIF
415      ENDIF
416
417      IF( td_file%i_ew >= 0 )THEN
418         il_ew=td_file%i_ew
419      ELSE
420         ! look for attribute in file
421         il_attid=att_get_index(tl_file%t_att(:),'ew_overlap')
422         IF( il_attid /= 0 )THEN
423            il_ew=INT(tl_file%t_att(il_attid)%d_value(1),i4)
424         ENDIF
425      ENDIF
426
427      SELECT CASE(il_perio)
428      CASE(3,4)
429         il_pivot=0
430      CASE(5,6)
431         il_pivot=1
432      CASE(0,1,2)
433         il_pivot=1
434      END SELECT
435     
436      IF( il_pivot < 0 .OR. il_pivot > 1 )THEN
437         ! get pivot
438         il_pivot=grid_get_pivot(tl_file)
439      ENDIF
440
441      IF( il_perio < 0 .OR. il_perio > 6 )THEN
442         ! get periodicity
443         il_perio=grid_get_perio(tl_file, il_pivot)
444      ENDIF
445
446      IF( il_ew < 0 )THEN
447         ! get periodicity
448         il_ew=grid_get_ew_overlap(tl_file)
449      ENDIF
450
451      ! close
452      CALL iom_close(tl_file)
453
454      !save in file structure
455      td_file%i_ew=il_ew
456      td_file%i_pivot=il_pivot
457      td_file%i_perio=il_perio
458
459      ! save in variable of file structure
460      tl_att=att_init("ew_overlap",il_ew)
461      DO ji=1,td_file%i_nvar
462         IF( td_file%t_var(ji)%t_dim(jp_I)%l_use )THEN
463            CALL var_move_att(td_file%t_var(ji),tl_att)
464         ENDIF
465      ENDDO
466
467      ! clean
468      CALL file_clean(tl_file)
469      CALL att_clean(tl_att)
470
471      IF( td_file%i_perio == -1 )THEN
472         CALL logger_fatal("GRID GET INFO: can not read or compute "//&
473         &  "domain periodicity from file "//TRIM(td_file%c_name)//"."//&
474         &  " you have to inform periodicity in namelist.")
475      ENDIF
476
477   END SUBROUTINE grid__get_info_file
478   !-------------------------------------------------------------------
479   !> @brief This subroutine get information about global domain, given mpp
480   !> structure.
481   !>
482   !> @details
483   !> open edge files then:
484   !> - compute NEMO pivot point
485   !> - compute NEMO periodicity
486   !> - compute East West overlap
487   !>
488   !> @note need all processor files
489   !> @author J.Paul
490   !> @date October, 2014 - Initial Version
491   !>
492   !> @param[in] td_mpp mpp structure
493   !-------------------------------------------------------------------
494   SUBROUTINE grid__get_info_mpp(td_mpp)
495      IMPLICIT NONE
496      ! Argument     
497      TYPE(TMPP) , INTENT(INOUT) :: td_mpp
498
499      ! local variable
500      INTEGER(i4) :: il_ew
501      INTEGER(i4) :: il_pivot
502      INTEGER(i4) :: il_perio
503      INTEGER(i4) :: il_attid
504
505      TYPE(TATT)  :: tl_att
506
507      TYPE(TMPP)  :: tl_mpp
508
509      ! loop indices
510      INTEGER(i4) :: ji
511      INTEGER(i4) :: jj
512      !----------------------------------------------------------------
513      ! intialise
514      il_pivot=-1
515      il_perio=-1
516      il_ew   =-1
517
518      CALL logger_info("GRID GET INFO: look for "//TRIM(td_mpp%c_name))
519      ! copy structure
520      tl_mpp=mpp_copy(td_mpp)
521      ! select edge files
522      CALL mpp_get_contour(tl_mpp)
523      ! open mpp file to be used
524      CALL iom_mpp_open(tl_mpp)
525
526      IF( td_mpp%i_perio >= 0 .AND. td_mpp%i_perio <= 6 )THEN
527         il_perio=td_mpp%i_perio
528      ELSE
529         ! look for attribute in mpp files
530         il_attid=att_get_index(tl_mpp%t_proc(1)%t_att(:),'periodicity')
531         IF( il_attid /= 0 )THEN
532            il_perio=INT(tl_mpp%t_proc(1)%t_att(il_attid)%d_value(1),i4)
533         ENDIF
534      ENDIF
535     
536      IF( td_mpp%i_ew >= 0 )THEN
537         il_ew=td_mpp%i_ew
538      ELSE
539         ! look for attribute in mpp files
540         il_attid=att_get_index(tl_mpp%t_proc(1)%t_att(:),'ew_overlap')
541         IF( il_attid /= 0 )THEN
542            il_ew=INT(tl_mpp%t_proc(1)%t_att(il_attid)%d_value(1),i4)
543         ENDIF
544      ENDIF
545
546      CALL logger_info("GRID GET INFO: perio "//TRIM(fct_str(il_perio)))
547
548      SELECT CASE(il_perio)
549         CASE(3,4)
550            il_pivot=1
551         CASE(5,6)
552            il_pivot=0
553         CASE(0,1,2)
554            il_pivot=1
555      END SELECT
556
557      IF( il_pivot < 0 .OR. il_pivot > 1 )THEN
558         ! get pivot
559         CALL logger_info("GRID GET INFO: look for pivot ")
560         il_pivot=grid_get_pivot(tl_mpp)
561      ENDIF
562
563      IF( il_perio < 0 .OR. il_perio > 6 )THEN
564         ! get periodicity
565         CALL logger_info("GRID GET INFO: look for perio ")
566         il_perio=grid_get_perio(tl_mpp, il_pivot)
567      ENDIF
568
569      IF( il_ew < 0 )THEN
570         ! get periodicity
571         CALL logger_info("GRID GET INFO: look for overlap ")
572         il_ew=grid_get_ew_overlap(tl_mpp)
573      ENDIF
574
575      ! close
576      CALL iom_mpp_close(tl_mpp)
577
578      !save in mpp structure
579      td_mpp%i_ew=il_ew
580      td_mpp%i_pivot=il_pivot
581      td_mpp%i_perio=il_perio
582
583      ! save in variable of mpp structure
584      IF( ASSOCIATED(td_mpp%t_proc) )THEN
585         tl_att=att_init("ew_overlap",il_ew)
586         DO jj=1,td_mpp%i_nproc
587            DO ji=1,td_mpp%t_proc(jj)%i_nvar
588               IF( td_mpp%t_proc(jj)%t_var(ji)%t_dim(jp_I)%l_use )THEN
589                  CALL var_move_att(td_mpp%t_proc(jj)%t_var(ji),tl_att)
590               ENDIF
591            ENDDO
592         ENDDO
593      ENDIF
594
595      ! clean
596      CALL mpp_clean(tl_mpp)
597      CALL att_clean(tl_att)
598
599      IF( td_mpp%i_perio == -1 )THEN
600         CALL logger_fatal("GRID GET INFO: can not read or compute "//&
601         &  "domain periodicity from mpp "//TRIM(td_mpp%c_name)//"."//&
602         &  " you have to inform periodicity in namelist.")
603      ENDIF
604
605   END SUBROUTINE grid__get_info_mpp
606   !-------------------------------------------------------------------
607   !> @brief
608   !> This function compute NEMO pivot point index of the input variable.
609   !> - F-point : 0
610   !> - T-point : 1
611   !>
612   !> @details
613   !> check north points of latitude grid (indices jpj to jpj-3) depending on which grid point
614   !> (T,F,U,V) variable is defined
615   !>
616   !> @note variable must be at least 2D variable, and should not be coordinate
617   !> variable (i.e lon, lat)
618   !>
619   !> @warning
620   !> - do not work with ORCA2 grid (T-point)
621   !>
622   !> @author J.Paul
623   !> @date November, 2013 - Initial version
624   !> @date September, 2014
625   !> - add dummy loop in case variable not over right point.
626   !> @date October, 2014
627   !> - work on variable structure instead of file structure
628   !
629   !> @param[in] td_lat  latitude variable structure
630   !> @param[in] td_var  variable structure
631   !> @return pivot point index
632   !-------------------------------------------------------------------
633   FUNCTION grid__get_pivot_var(td_var)
634      IMPLICIT NONE
635      ! Argument     
636      TYPE(TVAR), INTENT(IN) :: td_var
637
638      ! function
639      INTEGER(i4) :: grid__get_pivot_var
640
641      ! local variable
642      INTEGER(i4), DIMENSION(ip_maxdim)            :: il_dim
643
644      REAL(dp)   , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
645
646      ! loop indices
647      INTEGER(i4) :: jj
648      !----------------------------------------------------------------
649      ! intitalise
650      grid__get_pivot_var=-1
651
652      IF( .NOT. ASSOCIATED(td_var%d_value) .OR. &
653      &   .NOT. ALL(td_var%t_dim(1:2)%l_use) )THEN
654         CALL logger_error("GRID GET PIVOT: can not compute pivot point"//&
655         &  " with variable "//TRIM(td_var%c_name)//"."//&
656         &  " no value associated or missing dimension.")
657      ELSE
658         il_dim(:)=td_var%t_dim(:)%i_len
659
660         ALLOCATE(dl_value(il_dim(1),4,1,1))
661         ! extract value
662         dl_value(:,:,:,:)=td_var%d_value( 1:il_dim(1),          &
663         &                                 il_dim(2)-3:il_dim(2),&
664         &                                 1:1,                  &
665         &                                 1:1 )
666
667         SELECT CASE(TRIM(td_var%c_point))
668         CASE('T')
669            grid__get_pivot_var=grid__get_pivot_varT(dl_value)
670         CASE('U')
671            grid__get_pivot_var=grid__get_pivot_varU(dl_value)
672         CASE('V')
673            grid__get_pivot_var=grid__get_pivot_varV(dl_value)
674         CASE('F')
675            grid__get_pivot_var=grid__get_pivot_varF(dl_value)
676         END SELECT
677
678         ! dummy loop in case variable not over right point
679         ! (ex: nav_lon over U-point)
680         IF( grid__get_pivot_var == -1 )THEN
681           
682            ! no pivot point found
683            CALL logger_warn("GRID GET PIVOT: something wrong "//&
684            &  "when computing pivot point with variable "//&
685            &  TRIM(td_var%c_name))
686
687            DO jj=1,ip_npoint
688               SELECT CASE(TRIM(cp_grid_point(jj)))
689               CASE('T')
690                  CALL logger_debug("GRID GET PIVOT: check variable on point T")
691                  grid__get_pivot_var=grid__get_pivot_varT(dl_value)
692               CASE('U')
693                  CALL logger_debug("GRID GET PIVOT: check variable on point U")
694                  grid__get_pivot_var=grid__get_pivot_varU(dl_value)
695               CASE('V')
696                  CALL logger_debug("GRID GET PIVOT: check variable on point V")
697                  grid__get_pivot_var=grid__get_pivot_varV(dl_value)
698               CASE('F')
699                  CALL logger_debug("GRID GET PIVOT: check variable on point F")
700                  grid__get_pivot_var=grid__get_pivot_varF(dl_value)
701               END SELECT
702
703               IF( grid__get_pivot_var /= -1 )THEN
704                  CALL logger_info("GRID GET PIVOT: variable "//&
705                  &  TRIM(td_var%c_name)//" seems to be on grid point "//&
706                  &  TRIM(cp_grid_point(jj)) )
707                  EXIT
708               ENDIF
709
710            ENDDO
711         ENDIF
712
713         IF( grid__get_pivot_var == -1 )THEN
714            CALL logger_warn("GRID GET PIVOT: not able to found pivot point. "//&
715            &  "Force to use pivot point T.")
716            grid__get_pivot_var = 1
717         ENDIF
718
719         ! clean
720         DEALLOCATE(dl_value)
721
722      ENDIF
723
724   END FUNCTION grid__get_pivot_var
725   !-------------------------------------------------------------------
726   !> @brief
727   !> This function compute NEMO pivot point index for variable on grid T.
728   !>
729   !> @details
730   !> - F-point : 0
731   !> - T-point : 1
732   !>
733   !> @note array of value must be only the top border of the domain.
734   !>
735   !> @author J.Paul
736   !> @date October, 2014 - Initial version
737   !
738   !> @param[in] dd_value array of value
739   !> @return pivot point index
740   !-------------------------------------------------------------------
741   FUNCTION grid__get_pivot_varT(dd_value)
742      IMPLICIT NONE
743      ! Argument     
744      REAL(dp), DIMENSION(:,:,:,:), INTENT(IN) :: dd_value
745
746      ! function
747      INTEGER(i4) :: grid__get_pivot_varT
748
749      ! local variable
750      INTEGER(i4)                       :: il_midT
751      INTEGER(i4)                       :: il_midF
752
753      INTEGER(i4)                       :: it1
754      INTEGER(i4)                       :: it2
755      INTEGER(i4)                       :: jt1
756      INTEGER(i4)                       :: jt2
757
758      INTEGER(i4)                       :: if1
759      INTEGER(i4)                       :: if2
760      INTEGER(i4)                       :: jf1
761      INTEGER(i4)                       :: jf2
762
763      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
764
765      LOGICAL                           :: ll_check
766
767      ! loop indices
768      INTEGER(i4) :: ji
769      !----------------------------------------------------------------
770      ! intitalise
771      grid__get_pivot_varT=-1
772
773      il_dim(:)=SHAPE(dd_value(:,:,:,:))
774
775      ! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
776      jt1=4  ; jt2=2 
777      il_midT=il_dim(1)/2+1
778
779      ! F-point pivot !case of ORCA05 grid
780      jf1=4 ; jf2=3
781      il_midF=il_dim(1)/2
782
783      ! check T-point pivot
784      DO ji=2,il_midT
785         ll_check=.TRUE.
786         it1=ji
787         it2=il_dim(1)-(ji-2)
788         IF( dd_value(it1,jt1,1,1) /= dd_value(it2,jt2,1,1) )THEN
789            ll_check=.FALSE.
790            EXIT
791         ENDIF
792      ENDDO
793
794      IF( ll_check )THEN
795         CALL logger_info("GRID GET PIVOT: T-pivot")
796         grid__get_pivot_varT=1
797      ELSE
798
799         ! check F-point pivot
800         DO ji=1,il_midF
801            ll_check=.TRUE.
802            if1=ji
803            if2=il_dim(1)-(ji-1)
804            IF( dd_value(if1,jf1,1,1) /= dd_value(if2,jf2,1,1) )THEN
805               ll_check=.FALSE.
806               EXIT
807            ENDIF
808         ENDDO
809
810         IF( ll_check )THEN
811            CALL logger_info("GRID GET PIVOT: F-pivot")
812            grid__get_pivot_varT=0
813         ENDIF
814
815      ENDIF
816
817   END FUNCTION grid__get_pivot_varT
818   !-------------------------------------------------------------------
819   !> @brief
820   !> This function compute NEMO pivot point index for variable on grid U.
821   !>
822   !> @details
823   !> - F-point : 0
824   !> - T-point : 1
825   !>
826   !> @note array of value must be only the top border of the domain.
827   !>
828   !> @author J.Paul
829   !> @date October, 2014 - Initial version
830   !
831   !> @param[in] dd_value array of value
832   !> @return pivot point index
833   !-------------------------------------------------------------------
834   FUNCTION grid__get_pivot_varU(dd_value)
835      IMPLICIT NONE
836      ! Argument     
837      REAL(dp), DIMENSION(:,:,:,:), INTENT(IN) :: dd_value
838
839      ! function
840      INTEGER(i4) :: grid__get_pivot_varU
841
842      ! local variable
843      INTEGER(i4)                       :: il_midT
844      INTEGER(i4)                       :: il_midF
845
846      INTEGER(i4)                       :: it1
847      INTEGER(i4)                       :: it2
848      INTEGER(i4)                       :: jt1
849      INTEGER(i4)                       :: jt2
850
851      INTEGER(i4)                       :: if1
852      INTEGER(i4)                       :: if2
853      INTEGER(i4)                       :: jf1
854      INTEGER(i4)                       :: jf2
855
856      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
857
858      LOGICAL                           :: ll_check
859
860      ! loop indices
861      INTEGER(i4) :: ji
862      !----------------------------------------------------------------
863      ! intitalise
864      grid__get_pivot_varU=-1
865
866      il_dim(:)=SHAPE(dd_value(:,:,:,:))
867
868      ! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
869      jt1=4 ; jt2=2 
870      il_midT=il_dim(1)/2+1
871
872      ! F-point pivot !case of ORCA05 grid
873      jf1=4 ; jf2=3
874      il_midF=il_dim(1)/2
875
876      ! check T-point pivot
877      DO ji=1,il_midT
878         ll_check=.TRUE.
879         it1=ji
880         it2=il_dim(1)-(ji-2)
881         IF( dd_value(it1,jt1,1,1) /= dd_value(it2-1,jt2,1,1) )THEN
882            ll_check=.FALSE.
883            EXIT
884         ENDIF
885      ENDDO
886
887      IF( ll_check )THEN
888         CALL logger_info("GRID GET PIVOT: T-pivot")
889         grid__get_pivot_varU=1
890      ELSE
891
892         ! check F-point pivot
893         DO ji=1,il_midF
894            ll_check=.TRUE.
895            if1=ji
896            if2=il_dim(1)-(ji-1)
897            IF( dd_value(if1,jf1,1,1) /= dd_value(if2-1,jf2,1,1) )THEN
898               ll_check=.FALSE.
899               EXIT
900            ENDIF
901         ENDDO
902
903         IF( ll_check )THEN
904            CALL logger_info("GRID GET PIVOT: F-pivot")
905            grid__get_pivot_varU=0
906         ENDIF
907
908      ENDIF
909
910   END FUNCTION grid__get_pivot_varU
911   !-------------------------------------------------------------------
912   !> @brief
913   !> This function compute NEMO pivot point index for variable on grid V.
914   !>
915   !> @details
916   !> - F-point : 0
917   !> - T-point : 1
918   !>
919   !> @note array of value must be only the top border of the domain.
920   !>
921   !> @author J.Paul
922   !> @date October, 2014 - Initial version
923   !
924   !> @param[in] dd_value array of value
925   !> @return pivot point index
926   !-------------------------------------------------------------------
927   FUNCTION grid__get_pivot_varV(dd_value)
928      IMPLICIT NONE
929      ! Argument     
930      REAL(dp), DIMENSION(:,:,:,:), INTENT(IN) :: dd_value
931
932      ! function
933      INTEGER(i4) :: grid__get_pivot_varV
934
935      ! local variable
936      INTEGER(i4)                       :: il_midT
937      INTEGER(i4)                       :: il_midF
938
939      INTEGER(i4)                       :: it1
940      INTEGER(i4)                       :: it2
941      INTEGER(i4)                       :: jt1
942      INTEGER(i4)                       :: jt2
943
944      INTEGER(i4)                       :: if1
945      INTEGER(i4)                       :: if2
946      INTEGER(i4)                       :: jf1
947      INTEGER(i4)                       :: jf2
948
949      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
950
951      LOGICAL                           :: ll_check
952
953      ! loop indices
954      INTEGER(i4) :: ji
955      !----------------------------------------------------------------
956      ! intitalise
957      grid__get_pivot_varV=-1
958
959      il_dim(:)=SHAPE(dd_value(:,:,:,:))
960
961      ! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
962      jt1=4 ; jt2=2 
963      il_midT=il_dim(1)/2+1
964
965      ! F-point pivot !case of ORCA05 grid
966      jf1=4 ; jf2=3
967      il_midF=il_dim(1)/2
968
969      ! check T-point pivot
970      DO ji=2,il_midT
971         ll_check=.TRUE.
972         it1=ji
973         it2=il_dim(1)-(ji-2)
974         IF( dd_value(it1,jt1,1,1) /= dd_value(it2,jt2-1,1,1) )THEN
975            ll_check=.FALSE.
976            EXIT
977         ENDIF
978      ENDDO
979
980      IF( ll_check )THEN
981         CALL logger_info("GRID GET PIVOT: T-pivot")
982         grid__get_pivot_varV=1
983      ELSE
984
985         ! check F-point pivot
986         DO ji=1,il_midF
987            ll_check=.TRUE.
988            if1=ji
989            if2=il_dim(1)-(ji-1)
990            IF( dd_value(if1,jf1,1,1) /= dd_value(if2,jf2-1,1,1) )THEN
991               ll_check=.FALSE.
992               EXIT
993            ENDIF
994         ENDDO
995
996         IF( ll_check )THEN
997            CALL logger_info("GRID GET PIVOT: F-pivot")
998            grid__get_pivot_varV=0
999         ENDIF
1000
1001      ENDIF
1002
1003   END FUNCTION grid__get_pivot_varV
1004   !-------------------------------------------------------------------
1005   !> @brief
1006   !> This function compute NEMO pivot point index for variable on grid F.
1007   !>
1008   !> @details
1009   !> - F-point : 0
1010   !> - T-point : 1
1011   !>
1012   !> @note array of value must be only the top border of the domain.
1013   !>
1014   !> @author J.Paul
1015   !> @date October, 2014 - Initial version
1016   !
1017   !> @param[in] dd_value array of value
1018   !> @return pivot point index
1019   !-------------------------------------------------------------------
1020   FUNCTION grid__get_pivot_varF(dd_value)
1021      IMPLICIT NONE
1022      ! Argument     
1023      REAL(dp), DIMENSION(:,:,:,:), INTENT(IN) :: dd_value
1024
1025      ! function
1026      INTEGER(i4) :: grid__get_pivot_varF
1027
1028      ! local variable
1029      INTEGER(i4)                       :: il_midT
1030      INTEGER(i4)                       :: il_midF
1031
1032      INTEGER(i4)                       :: it1
1033      INTEGER(i4)                       :: it2
1034      INTEGER(i4)                       :: jt1
1035      INTEGER(i4)                       :: jt2
1036
1037      INTEGER(i4)                       :: if1
1038      INTEGER(i4)                       :: if2
1039      INTEGER(i4)                       :: jf1
1040      INTEGER(i4)                       :: jf2
1041
1042      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1043
1044      LOGICAL                           :: ll_check
1045
1046      ! loop indices
1047      INTEGER(i4) :: ji
1048      !----------------------------------------------------------------
1049      ! intitalise
1050      grid__get_pivot_varF=-1
1051
1052      il_dim(:)=SHAPE(dd_value(:,:,:,:))
1053
1054      ! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
1055      jt1=4 ; jt2=2 
1056      il_midT=il_dim(1)/2+1
1057
1058      ! F-point pivot !case of ORCA05 grid
1059      jf1=4 ; jf2=3
1060      il_midF=il_dim(1)/2
1061
1062      ! check T-point pivot
1063      DO ji=1,il_midT
1064         ll_check=.TRUE.
1065         it1=ji
1066         it2=il_dim(1)-(ji-2)
1067         IF( dd_value(it1,jt1,1,1) /= dd_value(it2-1,jt2-1,1,1) )THEN
1068            ll_check=.FALSE.
1069            EXIT
1070         ENDIF
1071      ENDDO
1072
1073      IF( ll_check )THEN
1074         CALL logger_info("GRID GET PIVOT: T-pivot")
1075         grid__get_pivot_varF=1
1076      ELSE
1077
1078         ! check F-point pivot
1079         DO ji=1,il_midF
1080            ll_check=.TRUE.
1081            if1=ji
1082            if2=il_dim(1)-(ji-1)
1083            IF( dd_value(if1,jf1,1,1) /= dd_value(if2-1,jf2-1,1,1) )THEN
1084               ll_check=.FALSE.
1085               EXIT
1086            ENDIF
1087         ENDDO
1088
1089         IF( ll_check )THEN
1090            CALL logger_info("GRID GET PIVOT: F-pivot")
1091            grid__get_pivot_varF=0
1092         ENDIF
1093
1094      ENDIF
1095
1096   END FUNCTION grid__get_pivot_varF
1097   !-------------------------------------------------------------------
1098   !> @brief
1099   !> This function compute NEMO pivot point index from input file variable.
1100   !> - F-point : 0
1101   !> - T-point : 1
1102   !>
1103   !> @details
1104   !> check north points of latitude grid (indices jpj to jpj-3) depending on which grid point
1105   !> (T,F,U,V) variable is defined
1106   !>
1107   !> @warning
1108   !> - do not work with ORCA2 grid (T-point)
1109   !>
1110   !> @author J.Paul
1111   !> @date Ocotber, 2014 - Initial version
1112   !
1113   !> @param[in] td_file file structure
1114   !> @return pivot point index
1115   !-------------------------------------------------------------------
1116   FUNCTION grid__get_pivot_file(td_file)
1117      IMPLICIT NONE
1118      ! Argument     
1119      TYPE(TFILE), INTENT(IN) :: td_file
1120
1121      ! function
1122      INTEGER(i4) :: grid__get_pivot_file
1123
1124      ! local variable
1125      INTEGER(i4)                       :: il_varid
1126      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1127
1128      LOGICAL                           :: ll_north
1129
1130      TYPE(TVAR)                        :: tl_var
1131      TYPE(TVAR)                        :: tl_lat
1132
1133      ! loop indices
1134      INTEGER(i4) :: ji
1135      !----------------------------------------------------------------
1136      ! intitalise
1137      grid__get_pivot_file=-1
1138
1139      ! look for north fold
1140      il_varid=var_get_index(td_file%t_var(:), 'latitude')
1141      IF( il_varid == 0 )THEN
1142         CALL logger_error("GRID GET PIVOT: no variable with name "//&
1143         &  "or standard name latitude in file structure "//&
1144         &  TRIM(td_file%c_name))
1145      ENDIF
1146      IF( ASSOCIATED(td_file%t_var(il_varid)%d_value) )THEN
1147         tl_lat=var_copy(td_file%t_var(il_varid))
1148      ELSE
1149         tl_lat=iom_read_var(td_file, 'latitude')
1150      ENDIF
1151
1152      ll_north=grid_is_north_fold(tl_lat)
1153      ! clean
1154      CALL var_clean(tl_lat)
1155
1156      IF( ll_north )THEN     
1157         ! look for suitable variable
1158         DO ji=1,td_file%i_nvar
1159            IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
1160
1161            IF( ASSOCIATED(td_file%t_var(ji)%d_value) )THEN
1162               tl_var=var_copy(td_file%t_var(ji))
1163            ELSE
1164               il_dim(:)=td_file%t_var(ji)%t_dim(:)%i_len
1165               tl_var=iom_read_var(td_file, &
1166               &                   td_file%t_var(ji)%c_name, &
1167               &                   id_start=(/1,il_dim(2)-3,1,1/), &
1168               &                   id_count=(/il_dim(1),4,1,1/) )
1169            ENDIF
1170         ENDDO
1171
1172         IF( ASSOCIATED(tl_var%d_value) )THEN
1173
1174            grid__get_pivot_file=grid_get_pivot(tl_var)
1175
1176         ENDIF
1177
1178         ! clean
1179         CALL var_clean(tl_var)
1180      ELSE
1181         CALL logger_warn("GRID GET PIVOT: no north fold. force to use T-PIVOT")
1182         grid__get_pivot_file=1
1183      ENDIF
1184
1185   END FUNCTION grid__get_pivot_file
1186   !-------------------------------------------------------------------
1187   !> @brief
1188   !> This function compute NEMO pivot point index from input mpp variable.
1189   !> - F-point : 0
1190   !> - T-point : 1
1191   !>
1192   !> @details
1193   !> check north points of latitude grid (indices jpj to jpj-3) depending on which grid point
1194   !> (T,F,U,V) variable is defined
1195   !>
1196   !> @warning
1197   !> - do not work with ORCA2 grid (T-point)
1198   !>
1199   !> @author J.Paul
1200   !> @date October, 2014 - Initial version
1201   !
1202   !> @param[in] td_mpp   mpp file structure
1203   !> @return pivot point index
1204   !-------------------------------------------------------------------
1205   FUNCTION grid__get_pivot_mpp(td_mpp)
1206      IMPLICIT NONE
1207      ! Argument     
1208      TYPE(TMPP), INTENT(IN) :: td_mpp
1209
1210      ! function
1211      INTEGER(i4) :: grid__get_pivot_mpp
1212
1213      ! local variable
1214      INTEGER(i4)                       :: il_varid
1215      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1216
1217      LOGICAL                           :: ll_north
1218
1219      TYPE(TVAR)                        :: tl_var
1220      TYPE(TVAR)                        :: tl_lat
1221 
1222      ! loop indices
1223      INTEGER(i4) :: ji
1224      !----------------------------------------------------------------
1225      ! intitalise
1226      grid__get_pivot_mpp=-1
1227
1228      ! look for north fold
1229      il_varid=var_get_index(td_mpp%t_proc(1)%t_var(:), 'latitude')
1230      IF( il_varid == 0 )THEN
1231         CALL logger_error("GRID GET PIVOT: no variable with name "//&
1232         &  "or standard name latitude in mpp structure "//&
1233         &  TRIM(td_mpp%c_name)//". Assume there is north fold and "//&
1234         &  "try to get pivot point")
1235
1236         ll_north=.TRUE.
1237      ELSE
1238         IF( ASSOCIATED(td_mpp%t_proc(1)%t_var(il_varid)%d_value) )THEN
1239            !
1240            tl_lat=mpp_recombine_var(td_mpp, 'latitude')
1241         ELSE
1242            tl_lat=iom_mpp_read_var(td_mpp, 'latitude')
1243         ENDIF     
1244
1245         ll_north=grid_is_north_fold(tl_lat)
1246      ENDIF
1247
1248      IF( ll_north )THEN
1249
1250         IF( ASSOCIATED(tl_lat%d_value) )THEN
1251            grid__get_pivot_mpp=grid_get_pivot(tl_lat)
1252         ELSE
1253            ! look for suitable variable
1254            DO ji=1,td_mpp%t_proc(1)%i_nvar
1255               IF(.NOT. ALL(td_mpp%t_proc(1)%t_var(ji)%t_dim(1:2)%l_use)) CYCLE
1256
1257               IF( ASSOCIATED(td_mpp%t_proc(1)%t_var(ji)%d_value) )THEN
1258                  CALL logger_debug("GRID GET PIVOT: mpp_recombine_var"//&
1259                  &         TRIM(td_mpp%t_proc(1)%t_var(ji)%c_name))
1260                  tl_var=mpp_recombine_var(td_mpp, &
1261                  &              TRIM(td_mpp%t_proc(1)%t_var(ji)%c_name))
1262               ELSE
1263                  CALL logger_debug("GRID GET PIVOT: iom_mpp_read_var "//&
1264                  &        TRIM(td_mpp%t_proc(1)%t_var(ji)%c_name))
1265                  il_dim(:)=td_mpp%t_dim(:)%i_len
1266
1267                  ! read variable
1268                  tl_var=iom_mpp_read_var(td_mpp, &
1269                  &                       td_mpp%t_proc(1)%t_var(ji)%c_name, &
1270                  &                       id_start=(/1,il_dim(2)-3,1,1/), &
1271                  &                       id_count=(/il_dim(1),4,1,1/) )
1272               ENDIF
1273               EXIT
1274            ENDDO
1275
1276            IF( ASSOCIATED(tl_var%d_value) )THEN
1277
1278               grid__get_pivot_mpp=grid_get_pivot(tl_var)
1279
1280           ELSE
1281               CALL logger_warn("GRID GET PIVOT: force to use T-PIVOT")
1282               grid__get_pivot_mpp=1
1283            ENDIF
1284
1285            ! clean
1286            CALL var_clean(tl_var)
1287         ENDIF
1288      ELSE
1289         CALL logger_warn("GRID GET PIVOT: no north fold. force to use T-PIVOT")
1290         grid__get_pivot_mpp=1
1291      ENDIF
1292
1293      CALL var_clean(tl_lat)
1294   END FUNCTION grid__get_pivot_mpp
1295   !-------------------------------------------------------------------
1296   !> @brief
1297   !> This subroutine search NEMO periodicity index given variable structure and
1298   !> pivot point index.
1299   !> @details
1300   !> The variable must be on T point.
1301   !>
1302   !> 0: closed boundaries
1303   !> 1: cyclic east-west boundary
1304   !> 2: symmetric boundary condition across the equator
1305   !> 3: North fold boundary (with a T-point pivot)
1306   !> 4: North fold boundary (with a T-point pivot) and cyclic east-west boundary
1307   !> 5: North fold boundary (with a F-point pivot)
1308   !> 6: North fold boundary (with a F-point pivot) and cyclic east-west boundary
1309   !>
1310   !> @warning pivot point should have been computed before run this script. see grid_get_pivot.
1311   !>
1312   !> @author J.Paul
1313   !> @date November, 2013 - Initial version
1314   !> @date October, 2014
1315   !> - work on variable structure instead of file structure
1316   !
1317   !> @param[in] td_var   variable structure
1318   !> @param[in] id_pivot pivot point index
1319   !-------------------------------------------------------------------
1320   FUNCTION grid__get_perio_var(td_var, id_pivot)
1321      IMPLICIT NONE
1322
1323      ! Argument     
1324      TYPE(TVAR) , INTENT(IN) :: td_var
1325      INTEGER(i4), INTENT(IN) :: id_pivot
1326
1327      ! function
1328      INTEGER(i4) :: grid__get_perio_var
1329
1330      ! local variable
1331      INTEGER(i4)                       :: il_perio
1332
1333      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1334
1335      ! loop indices
1336      !----------------------------------------------------------------
1337      ! intitalise
1338      grid__get_perio_var=-1
1339
1340      IF( id_pivot < 0 .OR. id_pivot > 1 )THEN
1341         CALL logger_error("GRID GET PERIO: invalid pivot point index. "//&
1342         &  "you should use grid_get_pivot to compute it")
1343      ENDIF
1344
1345      IF( .NOT. ASSOCIATED(td_var%d_value) .OR. &
1346      &   .NOT. ALL(td_var%t_dim(1:2)%l_use) )THEN
1347         CALL logger_error("GRID GET PERIO: can not compute periodicity"//&
1348         &  " with variable "//TRIM(td_var%c_name)//"."//&
1349         &  " no value associated or missing dimension.")
1350      ELSE
1351
1352         il_dim(:)=td_var%t_dim(:)%i_len
1353
1354         CALL logger_debug("GRID GET PERIO: use varibale "//TRIM(td_var%c_name))
1355         CALL logger_debug("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_fill)))
1356         CALL logger_debug("GRID GET PERIO: first value "//TRIM(fct_str(td_var%d_value(1,1,1,1))))
1357
1358         IF(ALL(td_var%d_value(    1    ,    :    ,1,1)/=td_var%d_fill).AND.&
1359         &  ALL(td_var%d_value(il_dim(1),    :    ,1,1)/=td_var%d_fill).AND.&
1360         &  ALL(td_var%d_value(    :    ,    1    ,1,1)/=td_var%d_fill).AND.&
1361         &  ALL(td_var%d_value(    :    ,il_dim(2),1,1)/=td_var%d_fill))THEN
1362         ! no boundary closed
1363            CALL logger_error("GRID GET PERIO: can't determined periodicity. "//&
1364            &             "there is no boundary closed for variable "//&
1365            &              TRIM(td_var%c_name) )
1366            ! check pivot
1367            SELECT CASE(id_pivot)
1368               CASE(0)
1369                  ! F pivot
1370                  CALL logger_warn("GRID GET PERIO: assume domain is global")
1371                  grid__get_perio_var=6
1372               CASE(1)
1373                  ! T pivot
1374                  CALL logger_warn("GRID GET PERIO: assume domain is global")
1375                  grid__get_perio_var=4
1376            END SELECT
1377         ELSE
1378            il_perio=-1
1379            ! check periodicity
1380            IF(ANY(td_var%d_value(   1     ,:,1,1)/=td_var%d_fill).OR.&
1381            &  ANY(td_var%d_value(il_dim(1),:,1,1)/=td_var%d_fill))THEN
1382            ! East-West cyclic (1,4,6)
1383
1384               IF( ANY(td_var%d_value(:, 1, 1, 1) /= td_var%d_fill) )THEN
1385               ! South boundary not closed
1386
1387                  CALL logger_debug("GRID GET PERIO: East_West cyclic")
1388                  CALL logger_debug("GRID GET PERIO: South boundary not closed")
1389                  CALL logger_error("GRID GET PERIO: should have been an "//&
1390                  &              "impossible case")
1391
1392               ELSE
1393               ! South boundary closed (1,4,6)
1394                  CALL logger_info("GRID GET PERIO: South boundary closed")
1395
1396                  IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill) )THEN
1397                  ! North boundary not closed (4,6)
1398                     CALL logger_info("GRID GET PERIO: North boundary not closed")
1399                     ! check pivot
1400                     SELECT CASE(id_pivot)
1401                        CASE(0)
1402                           ! F pivot
1403                           il_perio=6
1404                        CASE(1)
1405                           ! T pivot
1406                           il_perio=4
1407                        CASE DEFAULT
1408                           CALL logger_error("GRID GET PERIO: invalid pivot ")
1409                     END SELECT
1410                  ELSE
1411                  ! North boundary closed
1412                     CALL logger_info("GRID GET PERIO: North boundary closed")
1413                     il_perio=1 ! North and South boundaries closed
1414                  ENDIF
1415
1416               ENDIF
1417
1418            ELSE
1419            ! East-West boundaries closed (0,2,3,5)
1420               CALL logger_info("GRID GET PERIO: East West boundaries closed")
1421
1422               IF( ANY(td_var%d_value(:, 1, 1, 1) /= td_var%d_fill) )THEN
1423               ! South boundary not closed (2)
1424                  CALL logger_info("GRID GET PERIO: South boundary not closed")
1425
1426                  IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill))THEN
1427                  ! North boundary not closed
1428                     CALL logger_debug("GRID GET PERIO: East West boundaries "//&
1429                     &              "closed")
1430                     CALL logger_debug("GRID GET PERIO: South boundary not closed")
1431                     CALL logger_debug("GRID GET PERIO: North boundary not closed")
1432                     CALL logger_error("GRID GET PERIO: should have been "//&
1433                     &              "an impossible case")
1434                  ELSE
1435                  ! North boundary closed
1436                     il_perio=2   ! East-West and North boundaries closed
1437                  ENDIF
1438
1439               ELSE
1440               ! South boundary closed (0,3,5)
1441                  CALL logger_info("GRID GET PERIO: South boundary closed")
1442
1443                  IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill))THEN
1444                  ! North boundary not closed (3,5)
1445                     CALL logger_info("GRID GET PERIO: North boundary not closed")
1446                     ! check pivot
1447                     SELECT CASE(id_pivot)
1448                        CASE(0)
1449                           ! F pivot
1450                           il_perio=5
1451                        CASE(1)
1452                           ! T pivot
1453                           il_perio=3
1454                        CASE DEFAULT
1455                           CALL logger_error("GRID GET PERIO: invalid pivot")
1456                     END SELECT
1457                  ELSE
1458                  ! North boundary closed   
1459                     CALL logger_info("GRID GET PERIO: North boundary closed")
1460                     il_perio=0   ! all boundary closed
1461                  ENDIF
1462
1463               ENDIF
1464
1465            ENDIF
1466
1467            grid__get_perio_var=il_perio
1468
1469         ENDIF
1470
1471      ENDIF
1472
1473   END FUNCTION grid__get_perio_var
1474   !-------------------------------------------------------------------
1475   !> @brief
1476   !> This subroutine search NEMO periodicity index given file structure, and
1477   !> optionaly pivot point index.
1478   !> @details
1479   !> The variable used must be on T point.
1480   !>
1481   !> 0: closed boundaries
1482   !> 1: cyclic east-west boundary
1483   !> 2: symmetric boundary condition across the equator
1484   !> 3: North fold boundary (with a F-point pivot)
1485   !> 4: North fold boundary (with a F-point pivot) and cyclic east-west boundary
1486   !> 5: North fold boundary (with a T-point pivot)
1487   !> 6: North fold boundary (with a T-point pivot) and cyclic east-west boundary
1488   !>
1489   !> @warning pivot point should have been computed before run this script. see grid_get_pivot.
1490   !>
1491   !> @author J.Paul
1492   !> @date October, 2014 - Initial version
1493   !>
1494   !> @param[in] td_file   file structure
1495   !> @param[in] id_pivot  pivot point index
1496   !-------------------------------------------------------------------
1497   FUNCTION grid__get_perio_file(td_file, id_pivot)
1498      IMPLICIT NONE
1499
1500      ! Argument     
1501      TYPE(TFILE), INTENT(IN) :: td_file
1502      INTEGER(i4), INTENT(IN), OPTIONAL :: id_pivot
1503
1504      ! function
1505      INTEGER(i4) :: grid__get_perio_file
1506
1507      ! local variable
1508      INTEGER(i4)                       :: il_varid
1509      INTEGER(i4)                       :: il_pivot
1510      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1511
1512      TYPE(TVAR)                        :: tl_var
1513
1514      ! loop indices
1515      INTEGER(i4) :: ji
1516      !----------------------------------------------------------------
1517      !initialise
1518      grid__get_perio_file=-1
1519
1520      IF(PRESENT(id_pivot) )THEN
1521         il_pivot=id_pivot
1522      ELSE
1523         il_pivot=grid_get_pivot(td_file)
1524      ENDIF
1525
1526      IF( il_pivot < 0 .OR. il_pivot > 1 )THEN
1527         CALL logger_error("GRID GET PERIO: invalid pivot point index. "//&
1528         &  "you should use grid_get_pivot to compute it")
1529      ENDIF
1530
1531      ! look for suitable variable
1532      il_varid=0
1533      DO ji=1,td_file%i_nvar
1534         IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
1535         SELECT CASE(TRIM(fct_lower(td_file%t_var(ji)%c_stdname)) )
1536            CASE('longitude','latitude')
1537            CASE DEFAULT
1538               il_varid=ji
1539               EXIT
1540         END SELECT
1541      ENDDO
1542
1543      IF( il_varid==0 )THEN
1544     
1545         CALL logger_error("GRID GET PERIO: no suitable variable to compute "//&
1546         &              " periodicity in file "//TRIM(td_file%c_name))
1547
1548      ELSE
1549
1550         il_dim(:)= td_file%t_var(il_varid)%t_dim(:)%i_len
1551
1552         ! read variable
1553         tl_var=iom_read_var(td_file, &
1554         &                   td_file%t_var(il_varid)%c_name, &
1555         &                   id_start=(/1,1,1,1/), &
1556         &                   id_count=(/il_dim(1),il_dim(2),1,1/) )
1557
1558
1559         grid__get_perio_file=grid_get_perio(tl_var,il_pivot)
1560 
1561         ! clean
1562         CALL var_clean(tl_var)
1563
1564      ENDIF
1565
1566   END FUNCTION grid__get_perio_file
1567   !-------------------------------------------------------------------
1568   !> @brief
1569   !> This subroutine search NEMO periodicity given mpp structure and optionaly
1570   !> pivot point index.
1571   !> @details
1572   !> The variable used must be on T point.
1573   !>
1574   !> 0: closed boundaries
1575   !> 1: cyclic east-west boundary
1576   !> 2: symmetric boundary condition across the equator
1577   !> 3: North fold boundary (with a T-point pivot)
1578   !> 4: North fold boundary (with a T-point pivot) and cyclic east-west boundary
1579   !> 5: North fold boundary (with a F-point pivot)
1580   !> 6: North fold boundary (with a F-point pivot) and cyclic east-west boundary
1581   !>
1582   !> @warning pivot point should have been computed before run this script. see grid_get_pivot.
1583   !>
1584   !> @author J.Paul
1585   !> @date October, 2014 - Initial version
1586   !
1587   !> @param[in] td_mpp   mpp file structure
1588   !> @param[in] id_pivot pivot point index
1589   !-------------------------------------------------------------------
1590   FUNCTION grid__get_perio_mpp(td_mpp, id_pivot)
1591      IMPLICIT NONE
1592
1593      ! Argument     
1594      TYPE(TMPP) , INTENT(IN) :: td_mpp
1595      INTEGER(i4), INTENT(IN), OPTIONAL :: id_pivot
1596
1597      ! function
1598      INTEGER(i4) :: grid__get_perio_mpp
1599
1600      ! local variable
1601      INTEGER(i4)                       :: il_varid
1602      INTEGER(i4)                       :: il_pivot
1603      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
1604
1605      TYPE(TVAR)                        :: tl_var
1606
1607      ! loop indices
1608      INTEGER(i4) :: ji
1609      !----------------------------------------------------------------
1610      ! initialise
1611      grid__get_perio_mpp=-1
1612
1613      IF(PRESENT(id_pivot) )THEN
1614         il_pivot=id_pivot
1615      ELSE
1616         il_pivot=grid_get_pivot(td_mpp)
1617      ENDIF
1618
1619      IF( il_pivot < 0 .OR. il_pivot > 1 )THEN
1620         CALL logger_error("GRID GET PERIO: invalid pivot point index. "//&
1621         &  "you should use grid_get_pivot to compute it")
1622      ENDIF
1623
1624      ! look for suitable variable
1625      il_varid=0
1626      DO ji=1,td_mpp%t_proc(1)%i_nvar
1627         IF( .NOT. ALL(td_mpp%t_proc(1)%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
1628         SELECT CASE(TRIM(fct_lower(td_mpp%t_proc(1)%t_var(ji)%c_stdname)) )
1629            CASE('longitude','latitude')
1630            CASE DEFAULT
1631               il_varid=ji
1632               EXIT
1633         END SELECT
1634      ENDDO
1635
1636      IF( il_varid==0 )THEN
1637     
1638         CALL logger_error("GRID GET PERIO: no suitable variable to compute "//&
1639         &              " periodicity in file "//TRIM(td_mpp%c_name))
1640      ELSE
1641 
1642         DO ji=1,ip_maxdim
1643            IF( td_mpp%t_proc(1)%t_var(il_varid)%t_dim(ji)%l_use )THEN
1644               il_dim(ji)=td_mpp%t_dim(ji)%i_len
1645            ELSE
1646               il_dim(ji)=1
1647            ENDIF
1648         ENDDO
1649
1650         ! read variable
1651         tl_var=iom_mpp_read_var(td_mpp, &
1652         &                       td_mpp%t_proc(1)%t_var(il_varid)%c_name, &
1653         &                       id_start=(/1,1,1,1/), &
1654         &                       id_count=(/il_dim(1),il_dim(2),1,1/) )
1655
1656         grid__get_perio_mpp=grid_get_perio(tl_var, il_pivot)
1657
1658         ! clean
1659         CALL var_clean(tl_var)
1660      ENDIF
1661
1662   END FUNCTION grid__get_perio_mpp
1663   !-------------------------------------------------------------------
1664   !> @brief This function get East-West overlap.
1665   !
1666   !> @details
1667   !> If no East-West wrap return -1,
1668   !> else return the size of the ovarlap band.
1669   !> East-West overlap is computed comparing longitude value of the 
1670   !> South part of the domain, to avoid  north fold boundary.
1671   !>
1672   !
1673   !> @author J.Paul
1674   !> @date November, 2013 - Initial Version
1675   !> @date October, 2014
1676   !> - work on mpp file structure instead of file structure
1677   !> @date October, 2016
1678   !> - check longitude as longname
1679   !>
1680   !> @param[in] td_lon longitude variable structure
1681   !> @return East West overlap
1682   !-------------------------------------------------------------------
1683   FUNCTION grid__get_ew_overlap_var(td_var)
1684      IMPLICIT NONE
1685      ! Argument     
1686      TYPE(TVAR), INTENT(INOUT) :: td_var
1687      ! function
1688      INTEGER(i4) :: grid__get_ew_overlap_var
1689
1690      ! local variable
1691      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_value
1692      REAL(dp), DIMENSION(:)  , ALLOCATABLE :: dl_vare
1693      REAL(dp), DIMENSION(:)  , ALLOCATABLE :: dl_varw
1694
1695      REAL(dp)    :: dl_delta
1696      REAL(dp)    :: dl_varmax
1697      REAL(dp)    :: dl_varmin
1698
1699      INTEGER(i4) :: il_east
1700      INTEGER(i4) :: il_west
1701      INTEGER(i4) :: il_jmin
1702      INTEGER(i4) :: il_jmax
1703
1704      INTEGER(i4), PARAMETER :: il_max_overlap = 5
1705
1706      ! loop indices
1707      INTEGER(i4) :: ji
1708      !----------------------------------------------------------------
1709      ! initialise
1710      grid__get_ew_overlap_var=-1
1711
1712      IF( ASSOCIATED(td_var%d_value) )THEN
1713         IF( td_var%t_dim(1)%i_len > 1 )THEN
1714            il_west=1
1715            il_east=td_var%t_dim(1)%i_len
1716
1717            ALLOCATE( dl_value(td_var%t_dim(1)%i_len, &
1718            &                  td_var%t_dim(2)%i_len) )
1719
1720            dl_value(:,:)=td_var%d_value(:,:,1,1)
1721
1722            ! we do not use jmax as dimension length due to north fold boundary
1723            il_jmin=1+ip_ghost
1724            il_jmax=(td_var%t_dim(2)%i_len-ip_ghost)/2
1725
1726            ALLOCATE( dl_vare(il_jmax-il_jmin+1) )
1727            ALLOCATE( dl_varw(il_jmax-il_jmin+1) )
1728
1729            dl_vare(:)=dl_value(il_east,il_jmin:il_jmax)
1730            dl_varw(:)=dl_value(il_west,il_jmin:il_jmax)
1731
1732            IF( .NOT.( ALL(dl_vare(:)==td_var%d_fill) .AND. &
1733            &          ALL(dl_varw(:)==td_var%d_fill) ) )THEN
1734
1735               IF( TRIM(td_var%c_stdname) == 'longitude' .OR. &
1736                 & SCAN(TRIM(td_var%c_longname), 'longitude') == 0 )THEN
1737                  WHERE( dl_value(:,:) > 180._dp .AND. &
1738                  &      dl_value(:,:) /= td_var%d_fill ) 
1739                     dl_value(:,:)=360.-dl_value(:,:)
1740                  END WHERE
1741
1742                  dl_varmax=MAXVAL(dl_value(:,il_jmin:il_jmax))
1743                  dl_varmin=MINVAL(dl_value(:,il_jmin:il_jmax))
1744
1745                  dl_delta=(dl_varmax-dl_varmin)/td_var%t_dim(1)%i_len
1746
1747                  IF( ALL(ABS(dl_vare(:)) - ABS(dl_varw(:)) == dl_delta) )THEN
1748                     grid__get_ew_overlap_var=0
1749                  ENDIF
1750               ENDIF
1751
1752               IF( grid__get_ew_overlap_var == -1 )THEN
1753                  DO ji=0,il_max_overlap
1754
1755                     IF( il_east-ji == il_west )THEN
1756                        ! case of small domain
1757                        EXIT
1758                     ELSE
1759                        dl_vare(:)=dl_value(il_east-ji,il_jmin:il_jmax)
1760
1761                        IF( ALL( dl_varw(:) == dl_vare(:) ) )THEN
1762                           grid__get_ew_overlap_var=ji+1
1763                           EXIT
1764                        ENDIF
1765                     ENDIF
1766
1767                  ENDDO
1768               ENDIF
1769            ENDIF
1770
1771         ENDIF
1772      ELSE
1773         CALL logger_error("GRID GET EW OVERLAP: input variable standard name"//&
1774         &  TRIM(td_var%c_stdname)//" can not be used to compute East West "//&
1775         &  "overalp. no value associated. ")
1776      ENDIF
1777
1778   END FUNCTION grid__get_ew_overlap_var
1779   !-------------------------------------------------------------------
1780   !> @brief This function get East-West overlap.
1781   !
1782   !> @details
1783   !> If no East-West wrap return -1,
1784   !> else return the size of the ovarlap band.
1785   !> East-West overlap is computed comparing longitude value of the 
1786   !> South part of the domain, to avoid  north fold boundary.
1787   !>
1788   !> @author J.Paul
1789   !> @date October, 2014 - Initial Version
1790   !> @date October, 2016
1791   !> - check varid for longitude_T
1792   !>
1793   !> @param[in] td_file file structure
1794   !> @return East West overlap
1795   !-------------------------------------------------------------------
1796   FUNCTION grid__get_ew_overlap_file(td_file)
1797      IMPLICIT NONE
1798      ! Argument     
1799      TYPE(TFILE), INTENT(INOUT) :: td_file
1800      ! function
1801      INTEGER(i4) :: grid__get_ew_overlap_file
1802
1803      ! local variable
1804      INTEGER(i4) :: il_varid
1805
1806      TYPE(TVAR)  :: tl_var
1807
1808      ! loop indices
1809      INTEGER(i4) :: ji
1810      !----------------------------------------------------------------
1811
1812      il_varid=var_get_id(td_file%t_var(:), 'longitude', 'longitude_T')
1813      IF( il_varid /= 0 )THEN
1814         ! read longitude on boundary
1815         tl_var=iom_read_var(td_file, il_varid)
1816      ELSE
1817         DO ji=1,td_file%i_nvar
1818            IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
1819
1820            tl_var=iom_read_var(td_file, td_file%t_var(ji)%c_name)
1821            EXIT
1822         ENDDO
1823      ENDIF
1824
1825      grid__get_ew_overlap_file=grid_get_ew_overlap(tl_var)
1826
1827      ! clean
1828      CALL var_clean(tl_var)
1829
1830   END FUNCTION grid__get_ew_overlap_file
1831   !-------------------------------------------------------------------
1832   !> @brief This function get East-West overlap.
1833   !
1834   !> @details
1835   !> If no East-West wrap return -1,
1836   !> else return the size of the ovarlap band.
1837   !> East-West overlap is computed comparing longitude value of the 
1838   !> South part of the domain, to avoid  north fold boundary.
1839   !>
1840   !
1841   !> @author J.Paul
1842   !> @date November, 2013 - Initial Version
1843   !> @date October, 2014
1844   !> - work on mpp file structure instead of file structure
1845   !> @date October, 2016
1846   !> - check varid for longitude_T
1847   !>
1848   !> @param[in] td_mpp mpp structure
1849   !> @return East West overlap
1850   !-------------------------------------------------------------------
1851   FUNCTION grid__get_ew_overlap_mpp(td_mpp)
1852      IMPLICIT NONE
1853      ! Argument     
1854      TYPE(TMPP), INTENT(INOUT) :: td_mpp
1855      ! function
1856      INTEGER(i4) :: grid__get_ew_overlap_mpp
1857
1858      ! local variable
1859      INTEGER(i4) :: il_ew
1860      INTEGER(i4) :: il_varid
1861
1862      TYPE(TVAR)  :: tl_var
1863      ! loop indices
1864      INTEGER(i4) :: ji
1865      !----------------------------------------------------------------
1866
1867      ! initialise
1868      grid__get_ew_overlap_mpp=td_mpp%i_ew
1869
1870      ! read longitude on boundary
1871      il_varid=var_get_id(td_mpp%t_proc(1)%t_var(:),'longitude', 'longitude_T')
1872      IF( il_varid /= 0 )THEN
1873         tl_var=iom_mpp_read_var(td_mpp, il_varid)
1874      ELSE
1875         DO ji=1,td_mpp%t_proc(1)%i_nvar
1876            IF( .NOT. ALL(td_mpp%t_proc(1)%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
1877
1878            tl_var=iom_mpp_read_var(td_mpp, td_mpp%t_proc(1)%t_var(ji)%c_name)
1879            EXIT
1880         ENDDO         
1881      ENDIF
1882
1883      il_ew=grid_get_ew_overlap(tl_var)
1884      IF( il_ew >= 0 )THEN
1885         grid__get_ew_overlap_mpp=il_ew
1886      ENDIF
1887
1888      ! clean
1889      CALL var_clean(tl_var)
1890
1891   END FUNCTION grid__get_ew_overlap_mpp
1892   !-------------------------------------------------------------------
1893   !> @brief This subroutine check if there is north fold.
1894   !>
1895   !> @details
1896   !> check if maximum latitude greater than 88°N
1897   !>
1898   !> @author J.Paul
1899   !> @date November, 2013 - Initial Version
1900   !>
1901   !> @param[in] td_lat latitude variable structure
1902   !-------------------------------------------------------------------
1903   LOGICAL FUNCTION grid_is_north_fold(td_lat)
1904      IMPLICIT NONE
1905      ! Argument     
1906      TYPE(TVAR), INTENT(IN) :: td_lat
1907
1908      ! local variable
1909      ! loop indices
1910      !----------------------------------------------------------------
1911   
1912      ! init
1913      grid_is_north_fold=.FALSE.
1914
1915      IF( .NOT. ASSOCIATED(td_lat%d_value) )THEN
1916         CALL logger_error("GRID IS NORTH FOLD: "//&
1917         &                 " no value associated to latitude")
1918      ELSE     
1919         IF( MAXVAL(td_lat%d_value(:,:,:,:), &
1920         &          td_lat%d_value(:,:,:,:)/= td_lat%d_fill) >= 88.0 )THEN
1921
1922            grid_is_north_fold=.TRUE.
1923           
1924         ENDIF
1925      ENDIF
1926
1927   END FUNCTION grid_is_north_fold
1928   !-------------------------------------------------------------------
1929   !> @brief This subroutine check domain validity.
1930   !
1931   !> @details
1932   !> If maximum latitude greater than 88°N, program will stop.
1933   !> @note Not able to manage north fold for now.
1934   !
1935   !> @author J.Paul
1936   !> @date November, 2013 - Initial Version
1937   !> @date October, 2014
1938   !> - work on mpp file structure instead of file structure
1939   !
1940   !> @param[in] cd_coord  coordinate file
1941   !> @param[in] id_imin   i-direction lower left  point indice 
1942   !> @param[in] id_imax   i-direction upper right point indice
1943   !> @param[in] id_jmin   j-direction lower left  point indice
1944   !> @param[in] id_jmax   j-direction upper right point indice
1945   !-------------------------------------------------------------------
1946   SUBROUTINE grid_check_dom(td_coord, id_imin, id_imax, id_jmin, id_jmax)
1947      IMPLICIT NONE
1948      ! Argument     
1949      TYPE(TMPP) , INTENT(IN) :: td_coord
1950      INTEGER(i4), INTENT(IN) :: id_imin
1951      INTEGER(i4), INTENT(IN) :: id_imax
1952      INTEGER(i4), INTENT(IN) :: id_jmin
1953      INTEGER(i4), INTENT(IN) :: id_jmax
1954
1955      ! local variable
1956      TYPE(TVAR) :: tl_var
1957
1958      TYPE(TMPP) :: tl_coord
1959
1960      TYPE(TDOM) :: tl_dom
1961      ! loop indices
1962      !----------------------------------------------------------------
1963
1964      IF( id_jmin > id_jmax .OR. id_jmax == 0 )THEN
1965
1966         CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
1967         &  "can not create configuration with north pole.")
1968
1969      ELSE
1970
1971            IF( id_imin == id_imax .AND. td_coord%i_ew < 0 )THEN
1972               CALL logger_fatal("GRID CHECK DOM: invalid domain."//&
1973               &  " can not create east-west cyclic fine grid"//&
1974               &  " inside closed coarse grid")
1975            ENDIF
1976
1977            ! copy structure
1978            tl_coord=mpp_copy(td_coord)
1979
1980            ! compute domain
1981            tl_dom=dom_init( tl_coord,        &
1982            &                id_imin, id_imax,&
1983            &                id_jmin, id_jmax )
1984           
1985            ! open mpp files to be used
1986            CALL iom_dom_open(tl_coord, tl_dom)
1987
1988            ! read variable value on domain
1989            tl_var=iom_dom_read_var(tl_coord,'latitude',tl_dom)
1990
1991            ! close mpp files
1992            CALL iom_dom_close(tl_coord)
1993
1994            ! clean structure
1995            CALL mpp_clean(tl_coord)
1996
1997            IF( MAXVAL(tl_var%d_value(:,:,:,:), &
1998            &          tl_var%d_value(:,:,:,:)/= tl_var%d_fill) >= 88.0 )THEN
1999               
2000               CALL logger_debug("GRID CHECK DOM: max latitude "//&
2001               &  TRIM(fct_str(MAXVAL(tl_var%d_value(:,:,:,:)))) )
2002               CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
2003               &  "can not create configuration too close from north pole.")
2004
2005            ENDIF
2006
2007            ! clean
2008            CALL dom_clean(tl_dom)
2009            CALL var_clean(tl_var)
2010
2011      ENDIF
2012
2013   END SUBROUTINE grid_check_dom
2014   !-------------------------------------------------------------------
2015   !> @brief This function get closest coarse grid indices of fine grid domain.
2016   !
2017   !> @details
2018   !> it use coarse and fine grid coordinates files.
2019   !> optionally, you could specify the array of refinment factor (default 1.)
2020   !> optionally, you could specify on which Arakawa grid point you want to
2021   !> work (default 'T')
2022   !>
2023   !> @author J.Paul
2024   !> @date November, 2013 - Initial Version
2025   !> @date September, 2014
2026   !> - use grid point to read coordinates variable.
2027   !> @date October, 2014
2028   !> - work on mpp file structure instead of file structure
2029   !> @date February, 2015
2030   !> - use longitude or latitude as standard name, if can not find
2031   !> longitude_T, latitude_T...
2032   !>
2033   !> @param[in] td_coord0 coarse grid coordinate mpp structure
2034   !> @param[in] td_coord1 fine grid coordinate mpp structure
2035   !> @param[in] id_rho    array of refinment factor (default 1.)
2036   !> @param[in] cd_point  Arakawa grid point (default 'T').
2037   !> @return coarse grid indices(/(/imin0, imax0/), (/jmin0, jmax0/)/)
2038   !>                                     
2039   !-------------------------------------------------------------------
2040   FUNCTION grid__get_coarse_index_ff( td_coord0, td_coord1, &
2041   &                                   id_rho, cd_point )
2042      IMPLICIT NONE
2043      ! Argument
2044      TYPE(TMPP)                    , INTENT(IN) :: td_coord0
2045      TYPE(TMPP)                    , INTENT(IN) :: td_coord1
2046      INTEGER(i4)     , DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
2047      CHARACTER(LEN=*)              , INTENT(IN), OPTIONAL :: cd_point
2048
2049      ! function
2050      INTEGER(i4), DIMENSION(2,2) :: grid__get_coarse_index_ff
2051
2052      ! local variable
2053      CHARACTER(LEN= 1)                        :: cl_point
2054      CHARACTER(LEN=lc)                        :: cl_name
2055
2056      INTEGER(i4)                              :: il_imin0
2057      INTEGER(i4)                              :: il_imax0
2058      INTEGER(i4)                              :: il_jmin0
2059      INTEGER(i4)                              :: il_jmax0
2060      INTEGER(i4)                              :: il_ind
2061
2062      INTEGER(i4), DIMENSION(2,2)              :: il_xghost0
2063      INTEGER(i4), DIMENSION(2,2)              :: il_xghost1
2064
2065      INTEGER(i4), DIMENSION(:)  , ALLOCATABLE :: il_rho
2066
2067      TYPE(TVAR)                               :: tl_lon0
2068      TYPE(TVAR)                               :: tl_lat0
2069      TYPE(TVAR)                               :: tl_lon1
2070      TYPE(TVAR)                               :: tl_lat1
2071
2072      TYPE(TMPP)                               :: tl_coord0
2073      TYPE(TMPP)                               :: tl_coord1
2074
2075      ! loop indices
2076      !----------------------------------------------------------------
2077
2078      ! init
2079      grid__get_coarse_index_ff(:,:)=0
2080
2081      ALLOCATE(il_rho(ip_maxdim))
2082      il_rho(:)=1
2083      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
2084
2085      cl_point='T'
2086      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
2087
2088      ! copy structure
2089      tl_coord0=mpp_copy(td_coord0)
2090      tl_coord1=mpp_copy(td_coord1)
2091
2092      IF( .NOT. ASSOCIATED(tl_coord0%t_proc) .OR. &
2093      &   .NOT. ASSOCIATED(tl_coord1%t_proc) )THEN
2094         CALL logger_error("GRID GET COARSE INDEX: can not get coarse "//&
2095         &  "grid indices. decompsition of mpp file "//TRIM(tl_coord0%c_name)//&
2096         &  " and/or "//TRIM(tl_coord1%c_name)//" not defined." )
2097      ELSE
2098         ! Coarse grid
2099         ! get ghost cell factor on coarse grid
2100         il_xghost0(:,:)=grid_get_ghost( tl_coord0 )
2101
2102         ! open mpp files
2103         CALL iom_mpp_open(tl_coord0)
2104
2105         ! read coarse longitue and latitude
2106         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
2107         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
2108         IF( il_ind == 0 )THEN
2109            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2110            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
2111            &  try to use longitude.")
2112            WRITE(cl_name,*) 'longitude'
2113         ENDIF
2114         tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
2115
2116         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
2117         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
2118         IF( il_ind == 0 )THEN
2119            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2120            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
2121            &  try to use latitude.")
2122            WRITE(cl_name,*) 'latitude'
2123         ENDIF
2124         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
2125         
2126         CALL grid_del_ghost(tl_lon0, il_xghost0(:,:))
2127         CALL grid_del_ghost(tl_lat0, il_xghost0(:,:))
2128
2129         ! close mpp files
2130         CALL iom_mpp_close(tl_coord0)
2131
2132         ! Fine grid
2133
2134         ! get ghost cell factor on fine grid
2135         il_xghost1(:,:)=grid_get_ghost( tl_coord1 )
2136
2137         ! open mpp files
2138         CALL iom_mpp_open(tl_coord1)
2139
2140         ! read fine longitue and latitude
2141         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
2142         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
2143         IF( il_ind == 0 )THEN
2144            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2145            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
2146            &  try to use longitude.")
2147            WRITE(cl_name,*) 'longitude'
2148         ENDIF
2149         tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 
2150
2151         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
2152         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
2153         IF( il_ind == 0 )THEN
2154            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2155            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
2156            &  try to use latitude.")
2157            WRITE(cl_name,*) 'latitude'
2158         ENDIF
2159         tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
2160 
2161         CALL grid_del_ghost(tl_lon1, il_xghost1(:,:))
2162         CALL grid_del_ghost(tl_lat1, il_xghost1(:,:))
2163
2164         ! close mpp files
2165         CALL iom_mpp_close(tl_coord1)
2166
2167         ! compute
2168         grid__get_coarse_index_ff(:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
2169         &                                                    tl_lon1,tl_lat1,&
2170         &                                                    il_rho(:) )
2171
2172         ! add ghost cell to indices
2173         il_imin0=grid__get_coarse_index_ff(1,1)+il_xghost0(jp_I,1)*ip_ghost
2174         il_imax0=grid__get_coarse_index_ff(1,2)+il_xghost0(jp_I,1)*ip_ghost
2175
2176         il_jmin0=grid__get_coarse_index_ff(2,1)+il_xghost0(jp_J,1)*ip_ghost
2177         il_jmax0=grid__get_coarse_index_ff(2,2)+il_xghost0(jp_J,1)*ip_ghost
2178
2179         grid__get_coarse_index_ff(jp_I,1)=il_imin0
2180         grid__get_coarse_index_ff(jp_I,2)=il_imax0
2181         grid__get_coarse_index_ff(jp_J,1)=il_jmin0
2182         grid__get_coarse_index_ff(jp_J,2)=il_jmax0
2183
2184         CALL var_clean(tl_lon0)
2185         CALL var_clean(tl_lat0)         
2186         CALL var_clean(tl_lon1)
2187         CALL var_clean(tl_lat1)         
2188
2189      ENDIF
2190
2191      ! clean
2192      CALL mpp_clean(tl_coord0)
2193      CALL mpp_clean(tl_coord1)
2194      DEALLOCATE(il_rho)
2195
2196   END FUNCTION grid__get_coarse_index_ff
2197   !-------------------------------------------------------------------
2198   !> @brief This function get closest coarse grid indices of fine grid domain.
2199   !
2200   !> @details
2201   !> it use coarse array of longitude and latitude and fine grid coordinates file.
2202   !> optionaly, you could specify the array of refinment factor (default 1.)
2203   !> optionally, you could specify on which Arakawa grid point you want to
2204   !> work (default 'T')
2205   !>
2206   !> @author J.Paul
2207   !> @date November, 2013 - Initial Version
2208   !> @date September, 2014
2209   !> - use grid point to read coordinates variable.
2210   !> @date October, 2014
2211   !> - work on mpp file structure instead of file structure
2212   !> @date February, 2015
2213   !> - use longitude or latitude as standard name, if can not find
2214   !> longitude_T, latitude_T...
2215   !>
2216   !> @param[in] td_longitude0   coarse grid longitude
2217   !> @param[in] td_latitude0    coarse grid latitude
2218   !> @param[in] td_coord1       fine grid coordinate mpp structure
2219   !> @param[in] id_rho          array of refinment factor
2220   !> @param[in] cd_point        Arakawa grid point (default 'T')
2221   !> @return coarse grid indices (/(/imin0, imax0/), (/jmin0, jmax0/)/)
2222   !-------------------------------------------------------------------
2223   FUNCTION grid__get_coarse_index_cf( td_lon0, td_lat0, td_coord1, &
2224   &                                   id_rho, cd_point )
2225      IMPLICIT NONE
2226      ! Argument
2227      TYPE(TVAR )                   , INTENT(IN) :: td_lon0
2228      TYPE(TVAR )                   , INTENT(IN) :: td_lat0
2229      TYPE(TMPP )                   , INTENT(IN) :: td_coord1
2230      INTEGER(i4)     , DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
2231      CHARACTER(LEN=*)              , INTENT(IN), OPTIONAL :: cd_point
2232
2233      ! function
2234      INTEGER(i4), DIMENSION(2,2) :: grid__get_coarse_index_cf
2235
2236      ! local variable
2237      CHARACTER(LEN= 1)                        :: cl_point
2238      CHARACTER(LEN=lc)                        :: cl_name
2239
2240      INTEGER(i4)                              :: il_ind
2241
2242      INTEGER(i4), DIMENSION(:)  , ALLOCATABLE :: il_rho
2243
2244      INTEGER(i4), DIMENSION(2,2)              :: il_xghost
2245
2246      TYPE(TVAR)                               :: tl_lon1
2247      TYPE(TVAR)                               :: tl_lat1
2248
2249      TYPE(TMPP)                               :: tl_coord1
2250
2251      ! loop indices
2252      !----------------------------------------------------------------
2253
2254      ! init
2255      grid__get_coarse_index_cf(:,:)=0
2256
2257      ALLOCATE(il_rho(ip_maxdim) )
2258      il_rho(:)=1
2259      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
2260
2261      ! copy structure
2262      tl_coord1=mpp_copy(td_coord1)
2263     
2264      cl_point='T'
2265      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
2266
2267      IF( .NOT. ASSOCIATED(tl_coord1%t_proc) )THEN
2268         CALL logger_error("GRID GET COARSE INDEX: decompsition of mpp "//&
2269         &  "file "//TRIM(tl_coord1%c_name)//" not defined." )
2270
2271      ELSE IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
2272      &        .NOT. ASSOCIATED(td_lat0%d_value) )THEN
2273
2274         CALL logger_error("GRID GET COARSE INDEX: some coarse grid"//&
2275         &                 " coordinate value are not associated.")         
2276
2277      ELSE
2278         
2279         IF( TRIM(td_lon0%c_point)/='' )THEN
2280            cl_point=TRIM(td_lon0%c_point)
2281         ELSEIF( TRIM(td_lat0%c_point)/='' )THEN
2282            cl_point=TRIM(td_lat0%c_point)
2283         ENDIF
2284
2285         ! Fine grid
2286         ! get ghost cell factor on fine grid
2287         il_xghost(:,:)=grid_get_ghost( tl_coord1 )
2288
2289         ! open mpp files
2290         CALL iom_mpp_open(tl_coord1)
2291
2292         ! read fine longitue and latitude
2293         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
2294         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
2295         IF( il_ind == 0 )THEN
2296            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2297            &  TRIM(cl_name)//"in file "//TRIM(tl_coord1%c_name)//". &
2298            &  try to use longitude.")
2299            WRITE(cl_name,*) 'longitude'
2300         ENDIF
2301         tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
2302
2303         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
2304         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
2305         IF( il_ind == 0 )THEN
2306            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2307            &  TRIM(cl_name)//"in file "//TRIM(tl_coord1%c_name)//". &
2308            &  try to use longitude.")
2309            WRITE(cl_name,*) 'latitude'
2310         ENDIF
2311         tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
2312         
2313         CALL grid_del_ghost(tl_lon1, il_xghost(:,:))
2314         CALL grid_del_ghost(tl_lat1, il_xghost(:,:))
2315
2316         ! close mpp files
2317         CALL iom_mpp_close(tl_coord1)
2318
2319         ! compute
2320         grid__get_coarse_index_cf(:,:)=grid_get_coarse_index(td_lon0,td_lat0,&
2321         &                                                    tl_lon1,tl_lat1,&
2322         &                                                    il_rho(:), cl_point )
2323
2324         CALL var_clean(tl_lon1)
2325         CALL var_clean(tl_lat1)         
2326
2327      ENDIF
2328
2329      DEALLOCATE(il_rho)
2330      CALL mpp_clean(tl_coord1)
2331
2332   END FUNCTION grid__get_coarse_index_cf
2333   !-------------------------------------------------------------------
2334   !> @brief This function get closest coarse grid indices of fine grid domain.
2335   !
2336   !> @details
2337   !> it use coarse grid coordinates file and fine grid array of longitude and latitude.
2338   !> optionaly, you could specify the array of refinment factor (default 1.)
2339   !> optionally, you could specify on which Arakawa grid point you want to
2340   !> work (default 'T')
2341   !>
2342   !> @author J.Paul
2343   !> @date November, 2013 - Initial Version
2344   !> @date September, 2014
2345   !> - use grid point to read coordinates variable.
2346   !> @date October, 2014
2347   !> - work on mpp file structure instead of file structure
2348   !> @date February, 2015
2349   !> - use longitude or latitude as standard name, if can not find
2350   !> longitude_T, latitude_T...
2351   !>
2352   !> @param[in] td_coord0 coarse grid coordinate mpp structure
2353   !> @param[in] td_lon1   fine grid longitude
2354   !> @param[in] td_lat1   fine grid latitude
2355   !> @param[in] id_rho    array of refinment factor (default 1.)
2356   !> @param[in] cd_point  Arakawa grid point (default 'T')
2357   !> @return coarse grid indices (/(/imin0, imax0/), (/jmin0, jmax0/)/)
2358   !-------------------------------------------------------------------
2359   FUNCTION grid__get_coarse_index_fc( td_coord0, td_lon1, td_lat1, &
2360   &                                  id_rho, cd_point )
2361      IMPLICIT NONE
2362      ! Argument
2363      TYPE(TMPP )                   , INTENT(IN) :: td_coord0
2364      TYPE(TVAR )                   , INTENT(IN) :: td_lon1
2365      TYPE(TVAR )                   , INTENT(IN) :: td_lat1
2366      INTEGER(i4)     , DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
2367      CHARACTER(LEN=*)              , INTENT(IN), OPTIONAL :: cd_point
2368
2369      ! function
2370      INTEGER(i4), DIMENSION(2,2) :: grid__get_coarse_index_fc
2371
2372      ! local variable
2373      CHARACTER(LEN= 1)                        :: cl_point
2374      CHARACTER(LEN=lc)                        :: cl_name     
2375
2376      INTEGER(i4)                              :: il_imin0
2377      INTEGER(i4)                              :: il_imax0
2378      INTEGER(i4)                              :: il_jmin0
2379      INTEGER(i4)                              :: il_jmax0
2380      INTEGER(i4)                              :: il_ind
2381
2382      INTEGER(i4), DIMENSION(:), ALLOCATABLE   :: il_rho
2383
2384      INTEGER(i4), DIMENSION(2,2)              :: il_xghost
2385
2386      TYPE(TVAR)                               :: tl_lon0
2387      TYPE(TVAR)                               :: tl_lat0
2388
2389      TYPE(TMPP)                               :: tl_coord0
2390
2391      ! loop indices
2392      !----------------------------------------------------------------
2393
2394      ! init
2395      grid__get_coarse_index_fc(:,:)=0
2396
2397      ALLOCATE(il_rho(ip_maxdim))
2398      il_rho(:)=1
2399      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
2400
2401      cl_point='T'
2402      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
2403
2404      ! copy structure
2405      tl_coord0=mpp_copy(td_coord0)
2406
2407      IF( .NOT. ASSOCIATED(tl_coord0%t_proc) )THEN
2408         CALL logger_error("GRID GET COARSE INDEX: decompsition of mpp "//&
2409         &  "file "//TRIM(tl_coord0%c_name)//" not defined." )
2410
2411      ELSE IF( .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
2412      &        .NOT. ASSOCIATED(td_lat1%d_value) )THEN
2413
2414         CALL logger_error("GRID GET COARSE INDEX: some fine grid"//&
2415         &                 " coordinate value are not associated.")
2416
2417      ELSE
2418
2419         IF( TRIM(td_lon1%c_point)/='' )THEN
2420            cl_point=TRIM(td_lon1%c_point)
2421         ELSEIF( TRIM(td_lat1%c_point)/='' )THEN
2422            cl_point=TRIM(td_lat1%c_point)
2423         ENDIF
2424
2425         ! get ghost cell factor on coarse grid
2426         il_xghost(:,:)=grid_get_ghost( tl_coord0 )
2427
2428         ! open mpp files
2429         CALL iom_mpp_open(tl_coord0)
2430
2431         ! read coarse longitue and latitude
2432         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
2433         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
2434         IF( il_ind == 0 )THEN
2435            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2436            &  TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". &
2437            &  try to use longitude.")
2438            WRITE(cl_name,*) 'longitude'
2439         ENDIF
2440         tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
2441         
2442         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
2443         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
2444         IF( il_ind == 0 )THEN
2445            CALL logger_warn("GRID GET COARSE INDEX: no variable "//&
2446            &  TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". &
2447            &  try to use latitude.")
2448            WRITE(cl_name,*) 'latitude'
2449         ENDIF
2450         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
2451
2452         CALL grid_del_ghost(tl_lon0, il_xghost(:,:))
2453         CALL grid_del_ghost(tl_lat0, il_xghost(:,:))
2454
2455         ! close mpp files
2456         CALL iom_mpp_close(tl_coord0)
2457
2458         grid__get_coarse_index_fc(:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
2459         &                                                    td_lon1,td_lat1,&
2460         &                                                    il_rho(:), cl_point )
2461
2462         ! remove ghost cell
2463         il_imin0=grid__get_coarse_index_fc(1,1)+il_xghost(jp_I,1)*ip_ghost
2464         il_imax0=grid__get_coarse_index_fc(1,2)+il_xghost(jp_I,1)*ip_ghost
2465
2466         il_jmin0=grid__get_coarse_index_fc(2,1)+il_xghost(jp_J,1)*ip_ghost
2467         il_jmax0=grid__get_coarse_index_fc(2,2)+il_xghost(jp_J,1)*ip_ghost
2468
2469         grid__get_coarse_index_fc(1,1)=il_imin0
2470         grid__get_coarse_index_fc(1,2)=il_imax0
2471         grid__get_coarse_index_fc(2,1)=il_jmin0
2472         grid__get_coarse_index_fc(2,2)=il_jmax0
2473
2474         CALL var_clean(tl_lon0)
2475         CALL var_clean(tl_lat0)
2476
2477      ENDIF
2478
2479      CALL mpp_clean(tl_coord0)
2480      DEALLOCATE(il_rho)
2481
2482   END FUNCTION grid__get_coarse_index_fc
2483   !-------------------------------------------------------------------
2484   !> @brief This function get closest coarse grid indices of fine grid domain.
2485   !
2486   !> @details
2487   !> it use coarse and fine grid array of longitude and latitude.
2488   !> optionaly, you could specify the array of refinment factor (default 1.)
2489   !> optionally, you could specify on which Arakawa grid point you want to
2490   !> work (default 'T')
2491   !>
2492   !> @note do not use ghost cell
2493   !>
2494   !> @author J.Paul
2495   !> @date November, 2013 - Initial Version
2496   !> @date September, 2014
2497   !> - check grid point
2498   !> - take into account EW overlap
2499   !> @date February, 2016
2500   !> - use delta (lon or lat)
2501   !> - manage cases for T,U,V or F point, with even or odd refinment
2502   !>
2503   !> @param[in] td_lon0   coarse grid longitude
2504   !> @param[in] td_lat0   coarse grid latitude
2505   !> @param[in] td_lon1   fine grid longitude
2506   !> @param[in] td_lat1   fine grid latitude
2507   !> @param[in] id_rho    array of refinment factor
2508   !> @param[in] cd_point  Arakawa grid point ('T','U','V','F')
2509   !> @return coarse grid indices (/(/imin0, imax0/), (/jmin0, jmax0/)/)
2510   !>
2511   !> @todo
2512   !> -check case boundary domain on overlap band
2513   !-------------------------------------------------------------------
2514   FUNCTION grid__get_coarse_index_cc( td_lon0, td_lat0, td_lon1, td_lat1, &
2515   &                                   id_rho, cd_point )
2516      IMPLICIT NONE
2517      ! Argument
2518      TYPE(TVAR)                    , INTENT(IN) :: td_lon0
2519      TYPE(TVAR)                    , INTENT(IN) :: td_lat0
2520      TYPE(TVAR)                    , INTENT(IN) :: td_lon1
2521      TYPE(TVAR)                    , INTENT(IN) :: td_lat1
2522      INTEGER(i4)     , DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
2523      CHARACTER(LEN=*)              , INTENT(IN), OPTIONAL :: cd_point
2524
2525      ! function
2526      INTEGER(i4), DIMENSION(2,2) :: grid__get_coarse_index_cc
2527
2528      ! local variable
2529      CHARACTER(LEN= 1)                      :: cl_point0
2530      CHARACTER(LEN= 1)                      :: cl_point1
2531
2532      LOGICAL    , DIMENSION(2)              :: ll_even
2533
2534      REAL(dp)                               :: dl_lon1
2535      REAL(dp)                               :: dl_dlon
2536      REAL(dp)                               :: dl_lat1
2537      REAL(dp)                               :: dl_dlat
2538
2539      INTEGER(i4)                            :: il_ew0 
2540      INTEGER(i4)                            :: il_imin0
2541      INTEGER(i4)                            :: il_imax0
2542      INTEGER(i4)                            :: il_jmin0
2543      INTEGER(i4)                            :: il_jmax0
2544
2545      INTEGER(i4)                            :: il_ew1 
2546      INTEGER(i4)                            :: il_imin1
2547      INTEGER(i4)                            :: il_imax1
2548      INTEGER(i4)                            :: il_jmin1
2549      INTEGER(i4)                            :: il_jmax1
2550
2551      INTEGER(i4)                            :: il_imin
2552      INTEGER(i4)                            :: il_imax
2553      INTEGER(i4)                            :: il_jmin
2554      INTEGER(i4)                            :: il_jmax     
2555
2556      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
2557
2558      INTEGER(i4), DIMENSION(2)              :: il_ill
2559      INTEGER(i4), DIMENSION(2)              :: il_ilr
2560      INTEGER(i4), DIMENSION(2)              :: il_iul
2561      INTEGER(i4), DIMENSION(2)              :: il_iur
2562
2563      INTEGER(i4), DIMENSION(2,2)            :: il_xghost0
2564      INTEGER(i4), DIMENSION(2,2)            :: il_yghost0
2565      INTEGER(i4), DIMENSION(2,2)            :: il_xghost1
2566      INTEGER(i4), DIMENSION(2,2)            :: il_yghost1
2567
2568      TYPE(TVAR)                             :: tl_lon0
2569      TYPE(TVAR)                             :: tl_lat0
2570      TYPE(TVAR)                             :: tl_lon1
2571      TYPE(TVAR)                             :: tl_lat1
2572
2573      ! loop indices
2574      !----------------------------------------------------------------
2575      ! init
2576      grid__get_coarse_index_cc(:,:)=0
2577
2578      ALLOCATE( il_rho(ip_maxdim) )
2579      il_rho(:)=1
2580      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
2581
2582      ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /)
2583
2584      cl_point0='T'
2585      cl_point1='T'
2586      IF( PRESENT(cd_point) )THEN
2587         cl_point0=TRIM(fct_upper(cd_point))
2588         cl_point1=TRIM(fct_upper(cd_point))
2589      ENDIF
2590     
2591      IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
2592      &   .NOT. ASSOCIATED(td_lat0%d_value) .OR. &
2593      &   .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
2594      &   .NOT. ASSOCIATED(td_lat1%d_value) )THEN
2595         CALL logger_error("GRID GET COARSE INDEX: some fine or coarse grid"//&
2596         &                 " coordinate value not associated.")
2597      ELSE
2598
2599         IF( TRIM(td_lon0%c_point)/='' )THEN
2600            cl_point0=TRIM(td_lon0%c_point)
2601         ELSEIF( TRIM(td_lat0%c_point)/='' )THEN
2602            cl_point0=TRIM(td_lat0%c_point)
2603         ENDIF
2604         IF( TRIM(td_lon1%c_point)/='' )THEN
2605            cl_point1=TRIM(td_lon1%c_point)
2606         ELSEIF( TRIM(td_lat1%c_point)/='' )THEN
2607            cl_point1=TRIM(td_lat1%c_point)
2608         ENDIF
2609         IF( cl_point0 /= cl_point1 )THEN
2610            CALL logger_error("GRID GET COARSE INDEX: fine and coarse grid"//&
2611         &                 " coordinate not on same grid point.")
2612         ENDIF
2613
2614         IF( grid_is_global(td_lon1, td_lat1) )THEN
2615
2616            IF( grid_is_global(td_lon0, td_lat0) )THEN
2617               CALL logger_trace("GRID GET COARSE INDEX: fine grid is global ")
2618               grid__get_coarse_index_cc(:,:) = 1
2619               grid__get_coarse_index_cc(:,:) = 0
2620            ELSE
2621               CALL logger_error("GRID GET COARSE INDEX: fine grid is "//&
2622               &                 "global, coarse grid not.")
2623            ENDIF
2624
2625         ELSE
2626
2627            il_xghost0(:,:)=grid_get_ghost( td_lon0 )
2628            il_yghost0(:,:)=grid_get_ghost( td_lat0 )
2629            IF( ANY(il_xghost0(:,:) /= il_yghost0(:,:)) )THEN
2630               CALL logger_error("GRID GET COARSE INDEX: coarse grid "//&
2631               &        "coordinate do not share same ghost cell")
2632            ENDIF
2633
2634            tl_lon0=var_copy(td_lon0)
2635            tl_lat0=var_copy(td_lat0)
2636            CALL grid_del_ghost(tl_lon0, il_xghost0(:,:))
2637            CALL grid_del_ghost(tl_lat0, il_xghost0(:,:))
2638 
2639            ! "global" coarse grid indice
2640            il_imin0=1
2641            il_jmin0=1
2642
2643            il_imax0=tl_lon0%t_dim(1)%i_len
2644            il_jmax0=tl_lon0%t_dim(2)%i_len
2645
2646            ! get east west overlap for coarse grid
2647            il_ew0=tl_lon0%i_ew
2648            IF( il_ew0 >= 0 )THEN
2649               ! last point before overlap
2650               il_imax0=il_imax0-il_ew0
2651            ENDIF
2652
2653            il_xghost1(:,:)=grid_get_ghost( td_lon1 )
2654            il_yghost1(:,:)=grid_get_ghost( td_lat1 )
2655            IF( ANY(il_xghost1(:,:) /= il_yghost1(:,:)) )THEN
2656               CALL logger_error("GRID GET COARSE INDEX: fine grid "//&
2657               &        "coordinate do not share same ghost cell")
2658            ENDIF
2659
2660            tl_lon1=var_copy(td_lon1)
2661            tl_lat1=var_copy(td_lat1)
2662            CALL grid_del_ghost(tl_lon1, il_xghost1(:,:))
2663            CALL grid_del_ghost(tl_lat1, il_xghost1(:,:))
2664           
2665            ! "global" fine grid indice
2666            il_imin1=1
2667            il_jmin1=1
2668
2669            il_imax1=tl_lon1%t_dim(1)%i_len
2670            il_jmax1=tl_lon1%t_dim(2)%i_len
2671
2672            ! get east west overlap for fine grid
2673            il_ew1=tl_lon1%i_ew
2674            IF( il_ew1 >= 0 )THEN
2675               ! last point before overlap
2676               il_imax1=il_imax1-il_ew1
2677            ENDIF
2678
2679            ! get indices for each corner
2680            !1- search lower left corner indices
2681            dl_lon1=tl_lon1%d_value( il_imin1, il_jmin1, 1, 1 )
2682            dl_lat1=tl_lat1%d_value( il_imin1, il_jmin1, 1, 1 )
2683
2684            IF( dl_lon1 == tl_lon1%d_fill .OR. &
2685            &   dl_lat1 == tl_lat1%d_fill )THEN
2686               CALL logger_error("GRID GET COARSE INDEX: lower left corner "//&
2687               &                 "point is FillValue. remove ghost cell "//&
2688               &                 "before running grid_get_coarse_index.")
2689            ENDIF
2690
2691            !!!!! i-direction !!!!!
2692            IF( ll_even(jp_I) )THEN
2693               ! even
2694               SELECT CASE(TRIM(cl_point1))
2695                  CASE('F','U')
2696                     dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) -   &
2697                        &       tl_lon1%d_value(il_imin1  ,il_jmin1,1,1) ) / &
2698                        &     2.
2699                  CASE DEFAULT
2700                     dl_dlon=0
2701               END SELECT
2702            ELSE
2703               ! odd
2704               dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) -   &
2705                  &       tl_lon1%d_value(il_imin1  ,il_jmin1,1,1) ) / &
2706                  &     2.
2707            ENDIF
2708
2709            !!!!! j-direction !!!!!
2710            IF( ll_even(jp_J) )THEN
2711               ! even
2712               SELECT CASE(TRIM(cl_point1))
2713                  CASE('F','V')
2714                     dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) -   &
2715                        &       tl_lat1%d_value(il_imin1,il_jmin1  ,1,1) ) / &
2716                        &     2.
2717                  CASE DEFAULT
2718                     dl_dlat=0
2719               END SELECT
2720            ELSE
2721               ! odd
2722               dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) -   &
2723                  &       tl_lat1%d_value(il_imin1,il_jmin1  ,1,1) ) / &
2724                  &     2.
2725            ENDIF
2726
2727            dl_lon1 = dl_lon1 + dl_dlon
2728            dl_lat1 = dl_lat1 + dl_dlat
2729
2730            ! look for closest point on coarse grid
2731            il_ill(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, &
2732            &                                           il_jmin0:il_jmax0, &
2733            &                                           1,1), &
2734            &                           tl_lat0%d_value(il_imin0:il_imax0, &
2735            &                                           il_jmin0:il_jmax0, &
2736            &                                           1,1), &
2737            &                           dl_lon1, dl_lat1, 'll'   )
2738
2739
2740            !2- search upper left corner indices
2741            dl_lon1=tl_lon1%d_value( il_imin1, il_jmax1, 1, 1 )
2742            dl_lat1=tl_lat1%d_value( il_imin1, il_jmax1, 1, 1 )
2743
2744            IF( dl_lon1 == tl_lon1%d_fill .OR. &
2745            &   dl_lat1 == tl_lat1%d_fill )THEN
2746               CALL logger_error("GRID GET COARSE INDEX: upper left corner "//&
2747               &                 "point is FillValue. remove ghost cell "//&
2748               &                 "running grid_get_coarse_index.")
2749            ENDIF           
2750
2751            !!!!! i-direction !!!!!
2752            IF( ll_even(jp_I) )THEN
2753               ! even
2754               SELECT CASE(TRIM(cl_point1))
2755                  CASE('F','U')
2756                     dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) -   &
2757                        &       tl_lon1%d_value(il_imin1  ,il_jmax1,1,1) ) / &
2758                        &     2.
2759                  CASE DEFAULT
2760                     dl_dlon=0
2761               END SELECT
2762            ELSE
2763               ! odd
2764               dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) -   &
2765                  &       tl_lon1%d_value(il_imin1  ,il_jmax1,1,1) ) / &
2766                  &     2.
2767            ENDIF
2768
2769            !!!!! j-direction !!!!!
2770            IF( ll_even(jp_J) )THEN
2771               ! even
2772               SELECT CASE(TRIM(cl_point1))
2773                  CASE('F','V')
2774                     dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1  ,1,1) -   &
2775                        &       tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / &
2776                        &     2.
2777                  CASE DEFAULT
2778                     dl_dlat=0
2779               END SELECT
2780            ELSE
2781               ! odd
2782               dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1  ,1,1) -   &
2783                  &       tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / &
2784                  &     2.
2785            ENDIF
2786
2787            dl_lon1 = dl_lon1 + dl_dlon
2788            dl_lat1 = dl_lat1 - dl_dlat
2789
2790            ! look for closest point on coarse grid
2791            il_iul(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, &
2792            &                                           il_jmin0:il_jmax0, &
2793            &                                           1,1), &
2794            &                           tl_lat0%d_value(il_imin0:il_imax0, &
2795            &                                           il_jmin0:il_jmax0, &
2796            &                                           1,1), &
2797            &                           dl_lon1, dl_lat1, 'ul' )
2798
2799            !3- search lower right corner indices
2800            dl_lon1=tl_lon1%d_value( il_imax1, il_jmin1, 1, 1 )
2801            dl_lat1=tl_lat1%d_value( il_imax1, il_jmin1, 1, 1 )
2802
2803            IF( dl_lon1 == tl_lon1%d_fill .OR. &
2804            &   dl_lat1 == tl_lat1%d_fill )THEN
2805               CALL logger_error("GRID GET COARSE INDEX: lower right corner "//&
2806               &                 "point is FillValue. remove ghost cell "//&
2807               &                 "running grid_get_coarse_index.")
2808            ENDIF           
2809
2810            !!!!! i-direction !!!!!
2811            IF( ll_even(jp_I) )THEN
2812               ! even
2813               SELECT CASE(TRIM(cl_point1))
2814                  CASE('F','U')
2815                     dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmin1,1,1) -   &
2816                        &       tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / &
2817                        &     2.
2818                  CASE DEFAULT
2819                     dl_dlon=0
2820               END SELECT
2821            ELSE
2822               ! odd
2823               dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmin1,1,1) -   &
2824                  &       tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / &
2825                  &     2.
2826            ENDIF
2827
2828            !!!!! j-direction !!!!!
2829            IF( ll_even(jp_J) )THEN
2830               ! even
2831               SELECT CASE(TRIM(cl_point1))
2832                  CASE('F','V')
2833                     dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) -   &
2834                        &       tl_lat1%d_value(il_imax1,il_jmin1  ,1,1) ) / &
2835                        &     2.
2836                  CASE DEFAULT
2837                     dl_dlat=0
2838               END SELECT
2839            ELSE
2840               ! odd
2841               dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) -   &
2842                  &       tl_lat1%d_value(il_imax1,il_jmin1  ,1,1) ) / &
2843                  &     2.
2844            ENDIF
2845
2846            dl_lon1 = dl_lon1 - dl_dlon
2847            dl_lat1 = dl_lat1 + dl_dlat
2848
2849            ! look for closest point on coarse grid
2850            il_ilr(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, &
2851            &                                           il_jmin0:il_jmax0, &
2852            &                                           1,1), &
2853            &                           tl_lat0%d_value(il_imin0:il_imax0, &
2854            &                                           il_jmin0:il_jmax0, &
2855            &                                           1,1), &
2856            &                           dl_lon1, dl_lat1, 'lr' )
2857
2858            !4- search upper right corner indices
2859            dl_lon1=tl_lon1%d_value( il_imax1, il_jmax1, 1, 1 )
2860            dl_lat1=tl_lat1%d_value( il_imax1, il_jmax1, 1, 1 )
2861
2862            IF( dl_lon1 == tl_lon1%d_fill .OR. &
2863            &   dl_lat1 == tl_lat1%d_fill )THEN
2864               CALL logger_error("GRID GET COARSE INDEX: upper right corner "//&
2865               &                 "point is FillValue. remove ghost cell "//&
2866               &                 "before running grid_get_coarse_index.")
2867            ENDIF           
2868
2869            !!!!! i-direction !!!!!
2870            IF( ll_even(jp_I) )THEN
2871               ! even
2872               SELECT CASE(TRIM(cl_point1))
2873                  CASE('F','U')
2874                     dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmax1,1,1) -   &
2875                        &       tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / &
2876                        &     2.
2877                  CASE DEFAULT
2878                     dl_dlon=0
2879               END SELECT
2880            ELSE
2881               ! odd
2882               dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmax1,1,1) -   &
2883                  &       tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / &
2884                  &     2.
2885            ENDIF
2886
2887            !!!!! j-direction !!!!!
2888            IF( ll_even(jp_J) )THEN
2889               ! even
2890               SELECT CASE(TRIM(cl_point1))
2891                  CASE('F','V')
2892                     dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1  ,1,1) -   &
2893                        &       tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / &
2894                        &     2.
2895                  CASE DEFAULT
2896                     dl_dlat=0
2897               END SELECT
2898            ELSE
2899               ! odd
2900               dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1  ,1,1) -   &
2901                  &       tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / &
2902                  &     2.
2903            ENDIF
2904
2905            dl_lon1 = dl_lon1 - dl_dlon
2906            dl_lat1 = dl_lat1 - dl_dlat
2907
2908            ! look for closest point on coarse grid
2909            il_iur(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, &
2910            &                                           il_jmin0:il_jmax0, &
2911            &                                           1,1), &
2912            &                           tl_lat0%d_value(il_imin0:il_imax0, &
2913            &                                           il_jmin0:il_jmax0, &
2914            &                                           1,1), &
2915            &                           dl_lon1, dl_lat1, 'ur' )
2916
2917            ! coarse grid indices
2918            il_imin = il_imin0-1+MIN(il_ill(1), il_iul(1))
2919            il_imax = il_imin0-1+MAX(il_ilr(1), il_iur(1))
2920
2921            IF( il_imax <= il_ew0 )THEN
2922               !il_imin = 1
2923               il_imax = tl_lon0%t_dim(1)%i_len - il_ew0 + il_imax 
2924            ENDIF
2925
2926            il_jmin = il_jmin0-1+MIN(il_ill(2), il_ilr(2))
2927            il_jmax = il_jmin0-1+MAX(il_iul(2), il_iur(2))
2928
2929            ! special case if east west overlap
2930            IF( il_ew1 >= 0 )THEN
2931               CALL logger_debug("GRID GET COARSE INDEX: East-West overlap "//&
2932               &                 "found for fine grid " )
2933
2934               il_imin = 1
2935               il_imax = tl_lon0%t_dim(1)%i_len
2936
2937            ENDIF
2938         ENDIF
2939
2940         grid__get_coarse_index_cc(1,1) = il_imin
2941         grid__get_coarse_index_cc(1,2) = il_imax
2942
2943         grid__get_coarse_index_cc(2,1) = il_jmin
2944         grid__get_coarse_index_cc(2,2) = il_jmax
2945 
2946         ! clean
2947         CALL var_clean(tl_lon1)
2948         CALL var_clean(tl_lat1)
2949         CALL var_clean(tl_lon0)
2950         CALL var_clean(tl_lat0)
2951      ENDIF
2952
2953      DEALLOCATE( il_rho )
2954
2955   END FUNCTION grid__get_coarse_index_cc
2956   !-------------------------------------------------------------------
2957   !> @brief This function check if grid is global or not
2958   !
2959   !> @details
2960   !
2961   !> @author J.Paul
2962   !> @date November, 2013 - Initial Version
2963   !
2964   !> @param[in] td_lon longitude structure
2965   !> @param[in] td_lat latitude structure
2966   !-------------------------------------------------------------------
2967   FUNCTION grid_is_global(td_lon, td_lat)
2968      IMPLICIT NONE
2969      ! Argument     
2970      TYPE(TVAR), INTENT(IN) :: td_lon
2971      TYPE(TVAR), INTENT(IN) :: td_lat
2972
2973      ! function
2974      LOGICAL :: grid_is_global
2975     
2976      ! local variable
2977      INTEGER(i4)               :: il_ew
2978      INTEGER(i4)               :: il_south
2979      INTEGER(i4)               :: il_north
2980
2981      REAL(dp)                  :: dl_lat_min
2982      REAL(dp)                  :: dl_lat_max
2983
2984      ! loop indices
2985      !----------------------------------------------------------------
2986
2987      ! init
2988      grid_is_global=.FALSE.
2989
2990      IF( ANY( td_lon%t_dim(:)%i_len /= td_lat%t_dim(:)%i_len )  )THEN
2991         CALL logger_fatal("GRID IS GLOBAL: dimension of longitude and "//&
2992         &              " latitude differ")
2993      ENDIF
2994
2995      IF( .NOT. ASSOCIATED(td_lon%d_value) .OR. &
2996      &   .NOT. ASSOCIATED(td_lat%d_value) )THEN
2997         CALL logger_error("GRID IS GLOBAL: no value associated to "//&
2998         &              " longitude or latitude strucutre")     
2999      ELSE
3000
3001         il_south=1
3002         il_north=td_lon%t_dim(2)%i_len
3003
3004         dl_lat_min=MINVAL(td_lat%d_value(:,il_south,1,1))
3005         dl_lat_max=MAXVAL(td_lat%d_value(:,il_north,1,1))
3006
3007         IF( dl_lat_min < -77.0 .AND. dl_lat_max > 89.0 )THEN
3008
3009            il_ew=td_lon%i_ew
3010            IF( il_ew >= 0 )THEN
3011
3012               grid_is_global=.TRUE.
3013
3014            ENDIF
3015
3016         ENDIF
3017      ENDIF
3018
3019   END FUNCTION grid_is_global
3020   !-------------------------------------------------------------------
3021   !> @brief This function return grid indices of the closest point
3022   !> from point (lon1,lat1)
3023   !>
3024   !> @details
3025   !>
3026   !> @note overlap band should have been already removed from coarse grid array
3027   !> of longitude and latitude, before running this function
3028   !>
3029   !> if you add cd_pos argument, you could choice to return closest point at
3030   !> - lower left  (ll) of the point
3031   !> - lower right (lr) of the point
3032   !> - upper left  (ul) of the point
3033   !> - upper right (ur) of the point
3034   !> - lower       (lo) of the point
3035   !> - upper       (up) of the point
3036   !> -       left  (le) of the point
3037   !> -       right (ri) of the point
3038   !>
3039   !> @author J.Paul
3040   !> @date April, 2016 - Initial Version
3041   !> @date October, 2016
3042   !> - use max of zero and east-west overlap instead of east-west overlap
3043   !>
3044   !> @param[in] td_coord0 coarse grid coordinate mpp structure
3045   !> @param[in] dd_lon1   fine   grid longitude
3046   !> @param[in] dd_lat1   fine   grid latitude
3047   !> @param[in] cd_pos    relative position of grid point from point
3048   !> @param[in] dd_fill   fill value
3049   !> @return coarse grid indices of closest point of fine grid point
3050   !-------------------------------------------------------------------
3051   FUNCTION grid__get_closest_str( td_coord0, dd_lon1, dd_lat1, cd_pos, dd_fill ) &
3052   &  RESULT(id_res)
3053
3054      IMPLICIT NONE
3055      ! Argument
3056      TYPE(TMPP )     , INTENT(IN) :: td_coord0
3057      REAL(dp),         INTENT(IN) :: dd_lon1
3058      REAL(dp),         INTENT(IN) :: dd_lat1
3059      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: cd_pos
3060      REAL(dp),         INTENT(IN), OPTIONAL :: dd_fill
3061
3062      ! function
3063      INTEGER(i4), DIMENSION(2) :: id_res
3064
3065      ! local variable
3066      CHARACTER(LEN=lc)                        :: cl_point
3067      CHARACTER(LEN=lc)                        :: cl_name
3068
3069      INTEGER(i4)                              :: il_ind
3070      INTEGER(i4)                              :: il_ew
3071
3072      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
3073      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat0
3074
3075      TYPE(TVAR)                               :: tl_lon0
3076      TYPE(TVAR)                               :: tl_lat0
3077      TYPE(TMPP)                               :: tl_coord0
3078      !----------------------------------------------------------------
3079
3080      id_res(:)=-1
3081      cl_point='T'
3082
3083      ! copy structure
3084      tl_coord0=mpp_copy(td_coord0)
3085
3086      IF( .NOT. ASSOCIATED(tl_coord0%t_proc) )THEN
3087
3088         CALL logger_error("GRID GET CLOSEST: decompsition of mpp "//&
3089         &  "file "//TRIM(tl_coord0%c_name)//" not defined." )
3090
3091      ELSE
3092
3093         ! open mpp files
3094         CALL iom_mpp_open(tl_coord0)
3095 
3096         ! read coarse longitue and latitude
3097         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
3098         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3099         IF( il_ind == 0 )THEN
3100            CALL logger_warn("GRID GET CLOSEST: no variable "//&
3101            &  TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". &
3102            &  try to use longitude.")
3103            WRITE(cl_name,*) 'longitude'
3104         ENDIF
3105         tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3106 
3107         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
3108         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3109         IF( il_ind == 0 )THEN
3110            CALL logger_warn("GRID GET CLOSEST: no variable "//&
3111            &  TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". &
3112            &  try to use latitude.")
3113            WRITE(cl_name,*) 'latitude'
3114         ENDIF
3115         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3116
3117         ! close mpp files
3118         CALL iom_mpp_close(tl_coord0)
3119
3120         il_ew=MAX(0,tl_coord0%i_ew)
3121         ALLOCATE(dl_lon0(tl_coord0%t_dim(jp_I)%i_len-il_ew, &
3122            &             tl_coord0%t_dim(jp_J)%i_len) )             
3123         ALLOCATE(dl_lat0(tl_coord0%t_dim(jp_I)%i_len-il_ew, &
3124            &             tl_coord0%t_dim(jp_J)%i_len) )
3125
3126         dl_lon0(:,:)=tl_lon0%d_value(il_ew+1:,:,1,1)
3127         dl_lat0(:,:)=tl_lat0%d_value(il_ew+1:,:,1,1)
3128
3129         id_res(:)=grid_get_closest( dl_lon0, dl_lat0, dd_lon1, dd_lat1, cd_pos, dd_fill )
3130
3131         DEALLOCATE(dl_lon0, dl_lat0)
3132         CALL var_clean(tl_lon0)
3133         CALL var_clean(tl_lat0)
3134         CALL mpp_clean(tl_coord0)
3135
3136      ENDIF
3137
3138   END FUNCTION  grid__get_closest_str
3139   !-------------------------------------------------------------------
3140   !> @brief This function return grid indices of the closest point
3141   !> from point (lon1,lat1)
3142   !>
3143   !> @details
3144   !>
3145   !> @note overlap band should have been already removed from coarse grid array
3146   !> of longitude and latitude, before running this function
3147   !>
3148   !> if you add cd_pos argument, you could choice to return closest point at
3149   !> - lower left  (ll) of the point
3150   !> - lower right (lr) of the point
3151   !> - upper left  (ul) of the point
3152   !> - upper right (ur) of the point
3153   !> - lower       (lo) of the point
3154   !> - upper       (up) of the point
3155   !> -       left  (le) of the point
3156   !> -       right (ri) of the point
3157   !>
3158   !> @author J.Paul
3159   !> @date November, 2013 - Initial Version
3160   !> @date February, 2015
3161   !> - change dichotomy method to manage ORCA grid
3162   !> @date February, 2016
3163   !> - add optional use of relative position
3164   !
3165   !> @param[in] dd_lon0   coarse grid array of longitude
3166   !> @param[in] dd_lat0   coarse grid array of latitude
3167   !> @param[in] dd_lon1   fine   grid longitude
3168   !> @param[in] dd_lat1   fine   grid latitude
3169   !> @param[in] cd_pos    relative position of grid point from point
3170   !> @param[in] dd_fill   fill value
3171   !> @return coarse grid indices of closest point of fine grid point
3172   !-------------------------------------------------------------------
3173   FUNCTION grid__get_closest_arr( dd_lon0, dd_lat0, dd_lon1, dd_lat1, cd_pos, dd_fill )
3174      IMPLICIT NONE
3175      ! Argument
3176      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
3177      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
3178      REAL(dp),                 INTENT(IN) :: dd_lon1
3179      REAL(dp),                 INTENT(IN) :: dd_lat1
3180      CHARACTER(LEN=*),         INTENT(IN), OPTIONAL :: cd_pos
3181      REAL(dp),                 INTENT(IN), OPTIONAL :: dd_fill
3182
3183      ! function
3184      INTEGER(i4), DIMENSION(2) :: grid__get_closest_arr
3185
3186      ! local variable
3187      INTEGER(i4)                              :: il_iinf
3188      INTEGER(i4)                              :: il_imid
3189      INTEGER(i4)                              :: il_isup
3190      INTEGER(i4)                              :: il_jinf
3191      INTEGER(i4)                              :: il_jmid
3192      INTEGER(i4)                              :: il_jsup
3193      INTEGER(i4), DIMENSION(2)                :: il_shape
3194      INTEGER(i4), DIMENSION(1)                :: il_ind
3195   
3196      LOGICAL                                  :: ll_north
3197      LOGICAL                                  :: ll_continue
3198
3199      REAL(dp)                                 :: dl_lon1
3200      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dist
3201      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
3202
3203      ! loop indices
3204      !----------------------------------------------------------------
3205
3206      IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
3207         CALL logger_fatal("GRID GET CLOSEST: dimension of longitude and "//&
3208         &              " latitude differ")
3209      ENDIF
3210
3211      il_shape(:)=SHAPE(dd_lon0(:,:))
3212     
3213      ALLOCATE( dl_lon0(il_shape(1),il_shape(2)) ) 
3214     
3215      dl_lon0(:,:) = dd_lon0(:,:)
3216      WHERE(dd_lon0(:,:) < 0 ) dl_lon0(:,:) = dd_lon0(:,:) + 360.
3217
3218      dl_lon1 = dd_lon1
3219      IF( dd_lon1 < 0 ) dl_lon1 = dd_lon1 + 360.
3220
3221      ! first, use dichotomy to reduce domain
3222      il_iinf = 1              ; il_jinf = 1
3223      il_isup = il_shape(1)    ; il_jsup = il_shape(2)
3224
3225      il_shape(1)= il_isup - il_iinf + 1
3226      il_shape(2)= il_jsup - il_jinf + 1
3227
3228      ll_north=.FALSE.
3229      ll_continue=.FALSE.
3230
3231      ! avoid to use fillvalue for reduce domain on first time
3232      IF( PRESENT(dd_fill) )THEN
3233         DO WHILE( ALL(dl_lon0(il_isup,:) == dd_fill) )
3234            il_isup=il_isup-1
3235         ENDDO
3236         DO WHILE( ALL(dl_lon0(il_iinf,:) == dd_fill) )
3237            il_iinf=il_iinf+1
3238         ENDDO
3239         DO WHILE( ALL(dd_lat0(:,il_jsup) == dd_fill) )
3240            il_jsup=il_jsup-1
3241         ENDDO
3242         DO WHILE( ALL(dd_lat0(:,il_jinf) == dd_fill) )
3243            il_jinf=il_jinf+1
3244         ENDDO
3245
3246         il_shape(1)= il_isup - il_iinf + 1
3247         il_shape(2)= il_jsup - il_jinf + 1
3248
3249      ENDIF
3250
3251      ! special case for north ORCA grid
3252      IF( dd_lat1 > 19. .AND. dl_lon1 < 74.  )THEN
3253         ll_north=.TRUE.
3254      ENDIF
3255
3256      IF( .NOT. ll_north )THEN
3257         ! look for meridian 0°/360°
3258         il_jmid = il_jinf + INT(il_shape(2)/2)
3259         il_ind(:) = MAXLOC( dl_lon0(il_iinf:il_isup,il_jmid), &
3260         &                   dl_lon0(il_iinf:il_isup,il_jmid) <= 360._dp )
3261
3262         il_imid=il_ind(1)
3263
3264         IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
3265         &   dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN
3266
3267            il_iinf = il_imid ;  il_isup = il_imid
3268            il_jinf = il_jmid ;  il_jsup = il_jmid
3269
3270         ELSE
3271            IF( ALL(dl_lon0(il_isup,il_jinf:il_jsup) >  dl_lon1 ) .AND. &
3272            &   il_imid /= il_isup )THEN
3273               ! 0 < lon1 < lon0(isup)
3274               ! point east
3275               il_iinf = il_imid+1
3276               ll_continue=.TRUE.
3277         
3278            ELSE IF( ALL(dl_lon0(il_iinf,il_jinf:il_jsup) <  dl_lon1 ) .AND. &
3279            &        il_imid /= il_iinf )THEN
3280               ! lon0(iinf) < lon1 < 360
3281               ! point west
3282               il_isup = il_imid
3283               ll_continue=.TRUE.
3284
3285            ENDIF
3286
3287            il_shape(1)= il_isup - il_iinf + 1
3288            il_shape(2)= il_jsup - il_jinf + 1
3289
3290            il_imid = il_iinf + INT(il_shape(1)/2) 
3291            il_jmid = il_jinf + INT(il_shape(2)/2)
3292
3293            ! exit when close enough of point
3294            IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
3295         ENDIF
3296      ENDIF
3297
3298      !
3299      DO WHILE( ll_continue .AND. .NOT. ll_north )
3300
3301         ll_continue=.FALSE.
3302         IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
3303         &   dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN
3304
3305            il_iinf = il_imid ;  il_isup = il_imid
3306            il_jinf = il_jmid ;  il_jsup = il_jmid
3307
3308         ELSE
3309            IF( ALL(dl_lon0(il_imid,il_jinf:il_jsup) <  dl_lon1) )THEN   
3310
3311               ! point east
3312               il_iinf = il_imid
3313               ll_continue=.TRUE.
3314       
3315            ELSE IF( ALL(dl_lon0(il_imid,il_jinf:il_jsup) >  dl_lon1) )THEN   
3316
3317               ! point west
3318               il_isup = il_imid
3319               ll_continue=.TRUE.
3320
3321            ENDIF
3322
3323            IF( ALL(dd_lat0(il_iinf:il_isup,il_jmid) <  dd_lat1) )THEN   
3324               
3325               ! point north
3326               il_jinf = il_jmid
3327               ll_continue=.TRUE.
3328
3329            ELSE IF( ALL(dd_lat0(il_iinf:il_isup,il_jmid) > dd_lat1) )THEN   
3330
3331               ! point south
3332               il_jsup = il_jmid
3333               ll_continue=.TRUE.
3334           
3335            ENDIF
3336
3337            il_shape(1)= il_isup - il_iinf + 1
3338            il_shape(2)= il_jsup - il_jinf + 1
3339
3340            il_imid = il_iinf + INT(il_shape(1)/2) 
3341            il_jmid = il_jinf + INT(il_shape(2)/2)
3342
3343            ! exit when close enough of point
3344            IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
3345         ENDIF
3346         
3347      ENDDO
3348
3349      ! then find closest point by computing distances
3350      il_shape(1)= il_isup - il_iinf + 1
3351      il_shape(2)= il_jsup - il_jinf + 1
3352
3353      ALLOCATE( dl_dist(il_shape(1), il_shape(2)) )
3354
3355      dl_dist(:,:)=grid_distance(dl_lon0(il_iinf:il_isup,il_jinf:il_jsup), &
3356      &                          dd_lat0(il_iinf:il_isup,il_jinf:il_jsup), &
3357      &                          dl_lon1, dd_lat1 )
3358
3359      IF( PRESENT(cd_pos) )THEN
3360         !
3361         SELECT CASE(TRIM(cd_pos))
3362            CASE('le')
3363               WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 )
3364                  dl_dist(:,:)=NF90_FILL_DOUBLE
3365               END WHERE
3366            CASE('ri')
3367               WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 )
3368                  dl_dist(:,:)=NF90_FILL_DOUBLE
3369               END WHERE
3370            CASE('up')
3371               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 )
3372                  dl_dist(:,:)=NF90_FILL_DOUBLE
3373               END WHERE
3374            CASE('lo')
3375               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 )
3376                  dl_dist(:,:)=NF90_FILL_DOUBLE
3377               END WHERE
3378            CASE('ll')
3379               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. &
3380                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 )
3381                  dl_dist(:,:)=NF90_FILL_DOUBLE
3382               END WHERE
3383            CASE('lr')
3384               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. &
3385                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 )
3386                  dl_dist(:,:)=NF90_FILL_DOUBLE
3387               END WHERE               
3388            CASE('ul')
3389               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. &
3390                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 )
3391                  dl_dist(:,:)=NF90_FILL_DOUBLE
3392               END WHERE               
3393            CASE('ur')
3394               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. &
3395                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 )
3396                  dl_dist(:,:)=NF90_FILL_DOUBLE
3397               END WHERE
3398         END SELECT
3399      ENDIF
3400      grid__get_closest_arr(:)=MINLOC(dl_dist(:,:),dl_dist(:,:)/=NF90_FILL_DOUBLE)
3401
3402      grid__get_closest_arr(1)=grid__get_closest_arr(1)+il_iinf-1
3403      grid__get_closest_arr(2)=grid__get_closest_arr(2)+il_jinf-1
3404
3405      DEALLOCATE( dl_dist )
3406      DEALLOCATE( dl_lon0 )
3407
3408   END FUNCTION grid__get_closest_arr
3409   !-------------------------------------------------------------------
3410   !> @brief This function compute the distance between a point A and grid points. 
3411   !
3412   !> @details
3413   !
3414   !> @author J.Paul
3415   !> @date November, 2013 - Initial Version
3416   !
3417   !> @param[in] dd_lon    grid longitude array
3418   !> @param[in] dd_lat    grid latitude  array
3419   !> @param[in] dd_lonA   longitude of point A
3420   !> @param[in] dd_latA   latitude  of point A
3421   !> @param[in] dd_fill
3422   !> @return array of distance between point A and grid points.
3423   !-------------------------------------------------------------------
3424   FUNCTION grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA )
3425      IMPLICIT NONE
3426      ! Argument     
3427      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon
3428      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat
3429      REAL(dp),                 INTENT(IN) :: dd_lonA
3430      REAL(dp),                 INTENT(IN) :: dd_latA
3431
3432      ! function
3433      REAL(dp), DIMENSION(SIZE(dd_lon(:,:),DIM=1),&
3434      &                   SIZE(dd_lon(:,:),DIM=2)) :: grid_distance
3435
3436      ! local variable
3437      INTEGER(i4), DIMENSION(2) :: il_shape
3438
3439      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon
3440      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat
3441      REAL(dp)                                 :: dl_lonA
3442      REAL(dp)                                 :: dl_latA
3443
3444      REAL(dp)                                 :: dl_tmp
3445
3446      ! loop indices
3447      INTEGER(i4) :: ji
3448      INTEGER(i4) :: jj
3449      !----------------------------------------------------------------
3450
3451      IF( ANY( SHAPE(dd_lon(:,:)) /= SHAPE(dd_lat(:,:)) ) )THEN
3452         CALL logger_fatal("GRID DISTANCE: dimension of longitude and "//&
3453         &              " latitude differ")
3454      ENDIF
3455      il_shape(:)=SHAPE(dd_lon(:,:))
3456     
3457      ALLOCATE(dl_lon(il_shape(1),il_shape(2)))
3458      ALLOCATE(dl_lat(il_shape(1),il_shape(2)))
3459
3460      dl_lon(:,:) = dd_lon(:,:)
3461      dl_lonA     = dd_lonA
3462
3463      WHERE(dd_lon(:,:) < 0 ) dl_lon(:,:) = dd_lon(:,:) + 360.
3464      IF(   dd_lonA     < 0 ) dl_lonA     = dd_lonA     + 360.
3465     
3466      dl_lonA = dd_lonA * dp_deg2rad
3467      dl_latA = dd_latA * dp_deg2rad
3468
3469      dl_lon(:,:) = dl_lon(:,:) * dp_deg2rad
3470      dl_lat(:,:) = dd_lat(:,:) * dp_deg2rad
3471
3472      grid_distance(:,:)=NF90_FILL_DOUBLE
3473
3474      DO jj=1,il_shape(2)
3475         DO ji=1,il_shape(1)
3476            IF( dl_lon(ji,jj) == dl_lonA .AND. &
3477            &   dl_lat(ji,jj) == dl_latA )THEN
3478               grid_distance(ji,jj)=0.0
3479            ELSE
3480               dl_tmp= SIN(dl_latA)*SIN(dl_lat(ji,jj)) + &
3481               &       COS(dl_latA)*COS(dl_lat(ji,jj)) * &
3482               &       COS(dl_lon(ji,jj)-dl_lonA)
3483               ! check to avoid mistake with ACOS
3484               IF( dl_tmp < -1.0 ) dl_tmp = -1.0
3485               IF( dl_tmp >  1.0 ) dl_tmp =  1.0
3486               grid_distance(ji,jj)=ACOS(dl_tmp)*dp_rearth
3487            ENDIF
3488         ENDDO
3489      ENDDO
3490
3491      DEALLOCATE(dl_lon)
3492      DEALLOCATE(dl_lat)
3493
3494   END FUNCTION grid_distance
3495   !-------------------------------------------------------------------
3496   !> @brief This function get offset between fine grid and coarse grid.
3497   !
3498   !> @details
3499   !> optionally, you could specify on which Arakawa grid point you want to
3500   !> work (default 'T')
3501   !> offset value could be 0,1,..,rho-1
3502   !
3503   !> @author J.Paul
3504   !> @date September, 2014 - Initial Version
3505   !> @date October, 2014
3506   !> - work on mpp file structure instead of file structure
3507   !
3508   !> @param[in] td_coord0 coarse grid coordinate
3509   !> @param[in] id_imin0  coarse grid lower left corner i-indice of fine grid domain
3510   !> @param[in] id_jmin0  coarse grid lower left corner j-indice of fine grid domain
3511   !> @param[in] id_imax0  coarse grid upper right corner i-indice of fine grid domain
3512   !> @param[in] id_jmax0  coarse grid upper right corner j-indice of fine grid domain
3513   !> @param[in] td_coord1 fine   grid coordinate
3514   !> @param[in] id_rho    array of refinement factor
3515   !> @param[in] cd_point  Arakawa grid point
3516   !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /)
3517   !-------------------------------------------------------------------
3518   FUNCTION grid__get_fine_offset_ff( td_coord0, &
3519   &                                  id_imin0, id_jmin0, id_imax0, id_jmax0, &
3520   &                                  td_coord1, id_rho, cd_point )
3521      IMPLICIT NONE
3522      ! Argument
3523      TYPE(TMPP)                    , INTENT(IN) :: td_coord0
3524      TYPE(TMPP)                    , INTENT(IN) :: td_coord1
3525
3526      INTEGER(i4)                   , INTENT(IN) :: id_imin0
3527      INTEGER(i4)                   , INTENT(IN) :: id_jmin0
3528      INTEGER(i4)                   , INTENT(IN) :: id_imax0
3529      INTEGER(i4)                   , INTENT(IN) :: id_jmax0
3530
3531      INTEGER(i4)     , DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
3532      CHARACTER(LEN=*)              , INTENT(IN), OPTIONAL :: cd_point
3533
3534      ! function
3535      INTEGER(i4), DIMENSION(2,2) :: grid__get_fine_offset_ff
3536
3537      ! local variable
3538      INTEGER(i4)                              :: il_imin0
3539      INTEGER(i4)                              :: il_jmin0
3540      INTEGER(i4)                              :: il_imax0
3541      INTEGER(i4)                              :: il_jmax0
3542      INTEGER(i4)                              :: il_ind
3543     
3544      INTEGER(i4), DIMENSION(:), ALLOCATABLE   :: il_rho
3545     
3546      INTEGER(i4), DIMENSION(2,2)              :: il_xghost0
3547      INTEGER(i4), DIMENSION(2,2)              :: il_xghost1
3548
3549      CHARACTER(LEN= 1)                        :: cl_point
3550      CHARACTER(LEN=lc)                        :: cl_name
3551
3552      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
3553      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat0
3554      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon1
3555      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat1
3556
3557      TYPE(TVAR)                               :: tl_lon0
3558      TYPE(TVAR)                               :: tl_lat0
3559      TYPE(TVAR)                               :: tl_lon1
3560      TYPE(TVAR)                               :: tl_lat1
3561
3562      TYPE(TMPP)                               :: tl_coord0
3563      TYPE(TMPP)                               :: tl_coord1
3564     
3565      ! loop indices
3566      !----------------------------------------------------------------
3567      ! init
3568      grid__get_fine_offset_ff(:,:)=-1
3569
3570      ALLOCATE(il_rho(ip_maxdim))
3571      il_rho(:)=1
3572      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
3573
3574      cl_point='T'
3575      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
3576
3577      ! copy structure
3578      tl_coord0=mpp_copy(td_coord0)
3579      tl_coord1=mpp_copy(td_coord1)
3580
3581      IF( .NOT. ASSOCIATED(tl_coord0%t_proc) .OR. &
3582      &   .NOT. ASSOCIATED(tl_coord1%t_proc) )THEN
3583         CALL logger_error("GRID GET FINE OFFSET: can not get coarse "//&
3584         &  "grid indices. decompsition of mpp file "//TRIM(tl_coord0%c_name)//&
3585         &  " and/or "//TRIM(tl_coord1%c_name)//" not defined." )
3586      ELSE     
3587         !1- Coarse grid
3588         ! get ghost cell factor on coarse grid
3589         il_xghost0(:,:)=grid_get_ghost( tl_coord0 )
3590
3591         ! open mpp files
3592         CALL iom_mpp_open(tl_coord0)
3593
3594         ! read coarse longitue and latitude
3595         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
3596         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3597         IF( il_ind == 0 )THEN
3598            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3599            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
3600            &  try to use longitude.")
3601            WRITE(cl_name,*) 'longitude'
3602         ENDIF
3603         tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3604
3605         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
3606         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3607         IF( il_ind == 0 )THEN
3608            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3609            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
3610            &  try to use latitude.")
3611            WRITE(cl_name,*) 'latitude'
3612         ENDIF
3613         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3614
3615         ! close mpp files
3616         CALL iom_mpp_close(tl_coord0)
3617
3618         CALL grid_del_ghost(tl_lon0, il_xghost0(:,:))
3619         CALL grid_del_ghost(tl_lat0, il_xghost0(:,:))
3620
3621         ALLOCATE(dl_lon0(tl_lon0%t_dim(jp_I)%i_len, &
3622         &                tl_lon0%t_dim(jp_J)%i_len ))
3623
3624         dl_lon0(:,:)=tl_lon0%d_value(:,:,1,1)
3625
3626         ALLOCATE(dl_lat0(tl_lat0%t_dim(jp_I)%i_len, &
3627         &                tl_lat0%t_dim(jp_J)%i_len ))
3628
3629         dl_lat0(:,:)=tl_lat0%d_value(:,:,1,1)
3630
3631         ! clean
3632         CALL var_clean(tl_lon0)
3633         CALL var_clean(tl_lat0)
3634
3635         ! adjust coarse grid indices
3636         il_imin0=id_imin0-il_xghost0(jp_I,1)
3637         il_imax0=id_imax0-il_xghost0(jp_I,1)
3638
3639         il_jmin0=id_jmin0-il_xghost0(jp_J,1)
3640         il_jmax0=id_jmax0-il_xghost0(jp_J,1)
3641
3642         !2- Fine grid
3643         ! get ghost cell factor on fine grid
3644         il_xghost1(:,:)=grid_get_ghost( tl_coord1 )
3645
3646         ! open mpp files
3647         CALL iom_mpp_open(tl_coord1)
3648
3649         ! read fine longitue and latitude
3650         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
3651         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
3652         IF( il_ind == 0 )THEN
3653            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3654            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
3655            &  try to use longitude.")
3656            WRITE(cl_name,*) 'longitude'
3657         ENDIF
3658         tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
3659
3660         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
3661         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
3662         IF( il_ind == 0 )THEN
3663            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3664            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
3665            &  try to use latitude.")
3666            WRITE(cl_name,*) 'latitude'
3667         ENDIF
3668         tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
3669 
3670         ! close mpp files
3671         CALL iom_mpp_close(tl_coord1)
3672
3673         CALL grid_del_ghost(tl_lon1, il_xghost1(:,:))
3674         CALL grid_del_ghost(tl_lat1, il_xghost1(:,:))
3675
3676         ALLOCATE(dl_lon1(tl_lon1%t_dim(jp_I)%i_len, &
3677         &                tl_lon1%t_dim(jp_J)%i_len ))
3678
3679         dl_lon1(:,:)=tl_lon1%d_value(:,:,1,1)
3680
3681         ALLOCATE(dl_lat1(tl_lat1%t_dim(jp_I)%i_len, &
3682         &                tl_lat1%t_dim(jp_J)%i_len ))
3683
3684         dl_lat1(:,:)=tl_lat1%d_value(:,:,1,1)
3685
3686         ! clean
3687         CALL var_clean(tl_lon1)
3688         CALL var_clean(tl_lat1)
3689 
3690         !3- compute
3691         grid__get_fine_offset_ff(:,:)=grid_get_fine_offset( &
3692         &                                         dl_lon0(:,:), dl_lat0(:,:),&
3693         &                                         il_imin0, il_jmin0, &
3694         &                                         il_imax0, il_jmax0, &
3695         &                                         dl_lon1(:,:), dl_lat1(:,:),&
3696         &                                         id_rho(:), cl_point )
3697 
3698         DEALLOCATE(dl_lon0, dl_lat0)
3699         DEALLOCATE(dl_lon1, dl_lat1)
3700      ENDIF
3701
3702      ! clean
3703      CALL mpp_clean(tl_coord0)
3704      CALL mpp_clean(tl_coord1)
3705      DEALLOCATE(il_rho)
3706
3707   END FUNCTION grid__get_fine_offset_ff
3708   !-------------------------------------------------------------------
3709   !> @brief This function get offset between fine grid and coarse grid.
3710   !
3711   !> @details
3712   !> optionally, you could specify on which Arakawa grid point you want to
3713   !> work (default 'T')
3714   !> offset value could be 0,1,..,rho-1
3715   !
3716   !> @author J.Paul
3717   !> @date September, 2014 - Initial Version
3718   !> @date October, 2014
3719   !> - work on mpp file structure instead of file structure
3720   !
3721   !> @param[in] dd_lon0   coarse grid longitude array
3722   !> @param[in] dd_lat0   coarse grid latitude  array
3723   !> @param[in] id_imin0  coarse grid lower left corner i-indice of fine grid domain
3724   !> @param[in] id_jmin0  coarse grid lower left corner j-indice of fine grid domain
3725   !> @param[in] id_imax0  coarse grid upper right corner i-indice of fine grid domain
3726   !> @param[in] id_jmax0  coarse grid upper right corner j-indice of fine grid domain
3727   !> @param[in] td_coord1 fine   grid coordinate
3728   !> @param[in] id_rho    array of refinement factor
3729   !> @param[in] cd_point  Arakawa grid point
3730   !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /)
3731   !-------------------------------------------------------------------
3732   FUNCTION grid__get_fine_offset_cf( dd_lon0, dd_lat0, &
3733   &                                  id_imin0, id_jmin0, id_imax0, id_jmax0, &
3734   &                                  td_coord1, id_rho, cd_point )
3735      IMPLICIT NONE
3736      ! Argument
3737      REAL(dp)       , DIMENSION(:,:), INTENT(IN) :: dd_lon0
3738      REAL(dp)       , DIMENSION(:,:), INTENT(IN) :: dd_lat0
3739      TYPE(TMPP)                     , INTENT(IN) :: td_coord1
3740
3741      INTEGER(i4)                    , INTENT(IN) :: id_imin0
3742      INTEGER(i4)                    , INTENT(IN) :: id_jmin0
3743      INTEGER(i4)                    , INTENT(IN) :: id_imax0
3744      INTEGER(i4)                    , INTENT(IN) :: id_jmax0
3745
3746      INTEGER(i4)     , DIMENSION(:) , INTENT(IN), OPTIONAL :: id_rho
3747      CHARACTER(LEN=*)               , INTENT(IN), OPTIONAL :: cd_point
3748
3749      ! function
3750      INTEGER(i4), DIMENSION(2,2) :: grid__get_fine_offset_cf
3751
3752      ! local variable
3753      INTEGER(i4)                              :: il_ind
3754      INTEGER(i4), DIMENSION(2,2)              :: il_xghost1
3755      INTEGER(i4), DIMENSION(:), ALLOCATABLE   :: il_rho
3756     
3757      CHARACTER(LEN= 1)                        :: cl_point
3758      CHARACTER(LEN=lc)                        :: cl_name
3759
3760      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon1
3761      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat1
3762
3763      TYPE(TVAR)                               :: tl_lon1
3764      TYPE(TVAR)                               :: tl_lat1
3765
3766      TYPE(TMPP)                               :: tl_coord1
3767      ! loop indices
3768      !----------------------------------------------------------------
3769      ! init
3770      grid__get_fine_offset_cf(:,:)=-1
3771
3772      ALLOCATE(il_rho(ip_maxdim))
3773      il_rho(:)=1
3774      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
3775
3776      cl_point='T'
3777      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
3778
3779      ! copy structure
3780      tl_coord1=mpp_copy(td_coord1)
3781
3782      IF( .NOT. ASSOCIATED(tl_coord1%t_proc) )THEN
3783         CALL logger_error("GRID GET FINE OFFSET: decompsition of mpp "//&
3784         &  "file "//TRIM(tl_coord1%c_name)//" not defined." )
3785      ELSE     
3786
3787         ! Fine grid
3788         ! get ghost cell factor on fine grid
3789         il_xghost1(:,:)=grid_get_ghost( tl_coord1 )
3790
3791         ! open mpp files
3792         CALL iom_mpp_open(tl_coord1)
3793
3794         ! read fine longitue and latitude
3795         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
3796         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
3797         IF( il_ind == 0 )THEN
3798            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3799            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
3800            &  try to use longitude.")
3801            WRITE(cl_name,*) 'longitude'
3802         ENDIF
3803         tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
3804
3805         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
3806         il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name)
3807         IF( il_ind == 0 )THEN
3808            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3809            &  TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". &
3810            &  try to use latitude.")
3811            WRITE(cl_name,*) 'latitude'
3812         ENDIF
3813         tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name))
3814 
3815         ! close mpp files
3816         CALL iom_mpp_close(tl_coord1)
3817
3818         CALL grid_del_ghost(tl_lon1, il_xghost1(:,:))
3819         CALL grid_del_ghost(tl_lat1, il_xghost1(:,:))
3820
3821         ALLOCATE(dl_lon1(tl_lon1%t_dim(jp_I)%i_len, &
3822         &                tl_lon1%t_dim(jp_J)%i_len ))
3823
3824         dl_lon1(:,:)=tl_lon1%d_value(:,:,1,1)
3825
3826         ALLOCATE(dl_lat1(tl_lat1%t_dim(jp_I)%i_len, &
3827         &                tl_lat1%t_dim(jp_J)%i_len ))
3828
3829         dl_lat1(:,:)=tl_lat1%d_value(:,:,1,1)
3830
3831         ! clean
3832         CALL var_clean(tl_lon1)
3833         CALL var_clean(tl_lat1)
3834     
3835         ! compute
3836         grid__get_fine_offset_cf(:,:)=grid_get_fine_offset( &
3837         &                                         dd_lon0(:,:), dd_lat0(:,:),&
3838         &                                         id_imin0, id_jmin0, &
3839         &                                         id_imax0, id_jmax0, &
3840         &                                         dl_lon1(:,:), dl_lat1(:,:),&
3841         &                                         id_rho(:), cl_point )
3842         
3843         DEALLOCATE(dl_lon1, dl_lat1)
3844      ENDIF
3845
3846      ! clean
3847      CALL mpp_clean(tl_coord1)
3848      DEALLOCATE(il_rho)
3849
3850   END FUNCTION grid__get_fine_offset_cf
3851   !-------------------------------------------------------------------
3852   !> @brief This function get offset between fine grid and coarse grid.
3853   !
3854   !> @details
3855   !> optionally, you could specify on which Arakawa grid point you want to
3856   !> work (default 'T')
3857   !> offset value could be 0,1,..,rho-1
3858   !
3859   !> @author J.Paul
3860   !> @date September, 2014 - Initial Version
3861   !> @date October, 2014
3862   !> - work on mpp file structure instead of file structure
3863   !
3864   !> @param[in] td_coord0 coarse grid coordinate
3865   !> @param[in] id_imin0  coarse grid lower left corner i-indice of fine grid domain
3866   !> @param[in] id_jmin0  coarse grid lower left corner j-indice of fine grid domain
3867   !> @param[in] id_imax0  coarse grid upper right corner i-indice of fine grid domain
3868   !> @param[in] id_jmax0  coarse grid upper right corner j-indice of fine grid domain
3869   !> @param[in] dd_lon1   fine   grid longitude array
3870   !> @param[in] dd_lat1   fine   grid latitude  array
3871   !> @param[in] id_rho    array of refinement factor
3872   !> @param[in] cd_point  Arakawa grid point
3873   !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /)
3874   !-------------------------------------------------------------------
3875   FUNCTION grid__get_fine_offset_fc( td_coord0, &
3876   &                                  id_imin0, id_jmin0, id_imax0, id_jmax0, &
3877   &                                  dd_lon1, dd_lat1, &
3878   &                                  id_rho, cd_point )
3879      IMPLICIT NONE
3880      ! Argument
3881      TYPE(TMPP)                      , INTENT(IN) :: td_coord0
3882      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lon1
3883      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lat1
3884
3885      INTEGER(i4)                     , INTENT(IN) :: id_imin0
3886      INTEGER(i4)                     , INTENT(IN) :: id_jmin0
3887      INTEGER(i4)                     , INTENT(IN) :: id_imax0
3888      INTEGER(i4)                     , INTENT(IN) :: id_jmax0
3889
3890      INTEGER(i4)     , DIMENSION(:)  , INTENT(IN), OPTIONAL :: id_rho
3891      CHARACTER(LEN=*)                , INTENT(IN), OPTIONAL :: cd_point
3892
3893      ! function
3894      INTEGER(i4), DIMENSION(2,2) :: grid__get_fine_offset_fc
3895
3896      ! local variable
3897      INTEGER(i4)                              :: il_imin0
3898      INTEGER(i4)                              :: il_jmin0
3899      INTEGER(i4)                              :: il_imax0
3900      INTEGER(i4)                              :: il_jmax0
3901      INTEGER(i4)                              :: il_ind
3902     
3903      INTEGER(i4), DIMENSION(:), ALLOCATABLE   :: il_rho
3904     
3905      INTEGER(i4), DIMENSION(2,2)              :: il_xghost0
3906
3907      CHARACTER(LEN= 1)                        :: cl_point
3908      CHARACTER(LEN=lc)                        :: cl_name
3909
3910      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
3911      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat0
3912
3913      TYPE(TVAR)                               :: tl_lon0
3914      TYPE(TVAR)                               :: tl_lat0
3915
3916      TYPE(TMPP)                               :: tl_coord0
3917      ! loop indices
3918      !----------------------------------------------------------------
3919      ! init
3920      grid__get_fine_offset_fc(:,:)=-1
3921      ALLOCATE(il_rho(ip_maxdim))
3922      il_rho(:)=1
3923      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
3924
3925      cl_point='T'
3926      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
3927
3928      ! copy structure
3929      tl_coord0=mpp_copy(td_coord0)
3930
3931      IF( .NOT. ASSOCIATED(tl_coord0%t_proc) )THEN
3932         CALL logger_error("GRID GET FINE OFFSET: decompsition of mpp "//&
3933         &  "file "//TRIM(tl_coord0%c_name)//" not defined." )
3934      ELSE     
3935         !1- Coarse grid
3936         ! get ghost cell factor on coarse grid
3937         il_xghost0(:,:)=grid_get_ghost( tl_coord0 )
3938
3939         ! open mpp files
3940         CALL iom_mpp_open(tl_coord0)
3941
3942         ! read coarse longitude and latitude
3943         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point)
3944         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3945         IF( il_ind == 0 )THEN
3946            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3947            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
3948            &  try to use longitude.")
3949            WRITE(cl_name,*) 'longitude'
3950         ENDIF
3951         tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3952
3953         WRITE(cl_name,*) 'latitude_'//TRIM(cl_point)
3954         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name)
3955         IF( il_ind == 0 )THEN
3956            CALL logger_warn("GRID GET FINE OFFSET: no variable "//&
3957            &  TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". &
3958            &  try to use latitude.")
3959            WRITE(cl_name,*) 'latitude'
3960         ENDIF
3961         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name))
3962 
3963         ! close mpp files
3964         CALL iom_mpp_close(tl_coord0)
3965
3966         CALL grid_del_ghost(tl_lon0, il_xghost0(:,:))
3967         CALL grid_del_ghost(tl_lat0, il_xghost0(:,:))
3968
3969
3970         ALLOCATE(dl_lon0(tl_lon0%t_dim(jp_I)%i_len, &
3971         &                tl_lon0%t_dim(jp_J)%i_len ))
3972
3973         dl_lon0(:,:)=tl_lon0%d_value(:,:,1,1)
3974
3975         ALLOCATE(dl_lat0(tl_lat0%t_dim(jp_I)%i_len, &
3976         &                tl_lat0%t_dim(jp_J)%i_len ))
3977
3978         dl_lat0(:,:)=tl_lat0%d_value(:,:,1,1)
3979
3980         ! clean
3981         CALL var_clean(tl_lon0)
3982         CALL var_clean(tl_lat0)
3983
3984         ! adjust coarse grid indices
3985         il_imin0=id_imin0-il_xghost0(jp_I,1)
3986         il_imax0=id_imax0-il_xghost0(jp_I,1)
3987
3988         il_jmin0=id_jmin0-il_xghost0(jp_J,1)
3989         il_jmax0=id_jmax0-il_xghost0(jp_J,1)
3990
3991         !3- compute
3992         grid__get_fine_offset_fc(:,:)=grid_get_fine_offset(&
3993         &                                         dl_lon0(:,:), dl_lat0(:,:),&
3994         &                                         il_imin0, il_jmin0, &
3995         &                                         il_imax0, il_jmax0, &
3996         &                                         dd_lon1(:,:), dd_lat1(:,:),&
3997         &                                         id_rho(:), cl_point )
3998         
3999         DEALLOCATE(dl_lon0, dl_lat0)
4000      ENDIF
4001
4002      ! clean
4003      CALL mpp_clean(tl_coord0)
4004      DEALLOCATE(il_rho)
4005
4006   END FUNCTION grid__get_fine_offset_fc
4007   !-------------------------------------------------------------------
4008   !> @brief This function get offset between fine grid and coarse grid.
4009   !
4010   !> @details
4011   !> offset value could be 0,1,..,rho-1
4012   !
4013   !> @author J.Paul
4014   !> @date November, 2013 - Initial Version
4015   !> @date September, 2014
4016   !> - rename from grid_get_fine_offset
4017   !> @date May, 2015
4018   !> - improve way to find offset
4019   !> @date July, 2015
4020   !> - manage case close to greenwich meridian
4021   !> @date February, 2016
4022   !> - use grid_get_closest to assess offset
4023   !> - use delta (lon or lat)
4024   !> - manage cases for T,U,V or F point, with even or odd refinment
4025   !> - check lower left(upper right) fine grid point inside lower left(upper
4026   !> right) coarse grid cell.
4027   !>
4028   !> @todo check case close from North fold.
4029   !>
4030   !> @param[in] dd_lon0   coarse grid longitude array
4031   !> @param[in] dd_lat0   coarse grid latitude  array
4032   !> @param[in] id_imin0  coarse grid lower left corner i-indice of fine grid domain
4033   !> @param[in] id_jmin0  coarse grid lower left corner j-indice of fine grid domain
4034   !> @param[in] id_imax0  coarse grid upper right corner i-indice of fine grid domain
4035   !> @param[in] id_jmax0  coarse grid upper right corner j-indice of fine grid domain
4036   !> @param[in] dd_lon1   fine   grid longitude array
4037   !> @param[in] dd_lat1   fine   grid latitude  array
4038   !> @param[in] id_rho    array of refinement factor
4039   !> @param[in] cd_point  Arakawa grid point
4040   !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /)
4041   !-------------------------------------------------------------------
4042   FUNCTION grid__get_fine_offset_cc( dd_lon0, dd_lat0, &
4043   &                                  id_imin0, id_jmin0, id_imax0, id_jmax0, &
4044   &                                  dd_lon1, dd_lat1, id_rho, cd_point )
4045      IMPLICIT NONE
4046      ! Argument
4047      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lon0
4048      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lat0
4049      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lon1
4050      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lat1
4051
4052      INTEGER(i4)     ,                 INTENT(IN) :: id_imin0
4053      INTEGER(i4)     ,                 INTENT(IN) :: id_jmin0
4054      INTEGER(i4)     ,                 INTENT(IN) :: id_imax0
4055      INTEGER(i4)     ,                 INTENT(IN) :: id_jmax0
4056
4057      INTEGER(i4)     , DIMENSION(:)  , INTENT(IN) :: id_rho
4058      CHARACTER(LEN=*)                , INTENT(IN), OPTIONAL :: cd_point
4059
4060      ! function
4061      INTEGER(i4), DIMENSION(2,2) :: grid__get_fine_offset_cc
4062
4063      ! local variable
4064      CHARACTER(LEN= 1)                        :: cl_point
4065
4066      INTEGER(i4)                              :: i1
4067      INTEGER(i4)                              :: i2
4068      INTEGER(i4)                              :: j1
4069      INTEGER(i4)                              :: j2
4070
4071      INTEGER(i4), DIMENSION(2)                :: il_shape0
4072      INTEGER(i4), DIMENSION(2)                :: il_shape1
4073
4074      INTEGER(i4), DIMENSION(2)                :: il_ind
4075
4076      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
4077      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon1
4078
4079      REAL(dp)                                 :: dl_lonmax0
4080      REAL(dp)                                 :: dl_latmax0
4081      REAL(dp)                                 :: dl_lonmin0
4082      REAL(dp)                                 :: dl_latmin0
4083
4084      REAL(dp)                                 :: dl_lon0F
4085      REAL(dp)                                 :: dl_lat0F
4086      REAL(dp)                                 :: dl_dlon
4087      REAL(dp)                                 :: dl_dlat
4088
4089      LOGICAL    , DIMENSION(2)                :: ll_even
4090      LOGICAL                                  :: ll_greenwich
4091     
4092      ! loop indices
4093      INTEGER(i4) :: ii
4094      INTEGER(i4) :: ij
4095      !----------------------------------------------------------------
4096      IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
4097         CALL logger_fatal("GRID GET FINE OFFSET: dimension of coarse "//&
4098         &              "longitude and latitude differ")
4099      ENDIF
4100
4101      IF( ANY( SHAPE(dd_lon1(:,:)) /= SHAPE(dd_lat1(:,:)) ) )THEN
4102         CALL logger_fatal("GRID GET FINE OFFSET: dimension of fine "//&
4103         &              "longitude and latitude differ")
4104      ENDIF
4105
4106      ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /)
4107
4108      cl_point='T'
4109      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point))
4110
4111      il_shape0(:)=SHAPE(dd_lon0(:,:))
4112      ALLOCATE( dl_lon0(il_shape0(1),il_shape0(2)) )
4113
4114      il_shape1(:)=SHAPE(dd_lon1(:,:))
4115      ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) )
4116
4117      dl_lon0(:,:)=dd_lon0(:,:)
4118      WHERE( dd_lon0(:,:) < 0 ) dl_lon0(:,:)=dd_lon0(:,:)+360.
4119
4120      dl_lon1(:,:)=dd_lon1(:,:)
4121      WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360.
4122
4123      ! init
4124      grid__get_fine_offset_cc(:,:)=-1
4125      ll_greenwich=.FALSE.
4126
4127      IF( il_shape1(jp_J) == 1 )THEN
4128 
4129         grid__get_fine_offset_cc(jp_J,:)=((id_rho(jp_J)-1)/2)
4130
4131         !!! work on i-direction
4132         !!! look for i-direction left offset
4133         i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I))
4134         j1=1 ; j2=1
4135
4136         ! check if cross greenwich meridien
4137         IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))<5. .OR. &
4138           & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))>355. )THEN
4139            ! close to greenwich meridien
4140            ll_greenwich=.TRUE.
4141            ! 0:360 => -180:180
4142            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) > 180. )
4143               dl_lon0(id_imin0:id_imin0+1,id_jmin0) = &
4144                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0)-360.
4145            END WHERE
4146
4147            WHERE( dl_lon1(i1:i2,j1:j2) > 180. )
4148               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360.
4149            END WHERE
4150         ENDIF
4151
4152         ! max lognitude of the left cell
4153         dl_lonmax0=dl_lon0(id_imin0+1,id_jmin0)
4154         IF( dl_lon1(1,1) < dl_lonmax0 )THEN
4155
4156            !!!!! i-direction !!!!!
4157            IF( ll_even(jp_I) )THEN
4158               ! even
4159               SELECT CASE(TRIM(cl_point))
4160                  CASE('F','U')
4161                     dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) -   &
4162                        &       dl_lon0(id_imin0  ,id_jmin0) ) / &
4163                        &     ( 2.*id_rho(jp_I) )
4164                  CASE DEFAULT
4165                     dl_dlon=0
4166               END SELECT
4167            ELSE
4168               ! odd
4169               dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) -   &
4170                  &       dl_lon0(id_imin0  ,id_jmin0) ) / &
4171                  &     ( 2.*id_rho(jp_I) )
4172            ENDIF
4173
4174            dl_lon0F= dl_lon0(id_imin0+1,id_jmin0) + dl_dlon
4175            dl_lat0F= dd_lat0(id_imin0+1,id_jmin0)
4176
4177            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4178            &                           dl_lon0F, dl_lat0F, 'le' )
4179     
4180            ii=il_ind(1)
4181
4182            !!!!! i-direction !!!!!
4183            IF( ll_even(jp_I) )THEN
4184               ! even
4185               SELECT CASE(TRIM(cl_point))
4186                  CASE('T','V')
4187                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii
4188                  CASE DEFAULT !'F','U'
4189                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4190               END SELECT
4191            ELSE
4192               ! odd
4193               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4194            ENDIF
4195
4196         ELSE
4197            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
4198            &                 " not match fine grid left corner.")
4199         ENDIF
4200
4201         IF( ll_greenwich )THEN
4202            ! close to greenwich meridien
4203            ll_greenwich=.FALSE.
4204            ! -180:180 => 0:360
4205            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) < 0. )
4206               dl_lon0(id_imin0:id_imin0+1,id_jmin0) = &
4207                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0)+360.
4208            END WHERE
4209
4210            WHERE( dl_lon1(i1:i2,j1:j2) < 0. )
4211               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360.
4212            END WHERE
4213         ENDIF
4214
4215         !!!!!! look for i-direction right offset !!!!!!
4216         i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I)
4217         j1=1                                         ; j2=1
4218
4219         ! check if cross greenwich meridien
4220         IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))<5. .OR. &
4221           & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))>355. )THEN
4222            ! close to greenwich meridien
4223            ll_greenwich=.TRUE.
4224            ! 0:360 => -180:180
4225            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) > 180. )
4226               dl_lon0(id_imax0-1:id_imax0,id_jmin0) = &
4227                  & dl_lon0(id_imax0-1:id_imax0,id_jmin0)-360.
4228            END WHERE
4229
4230            WHERE( dl_lon1(i1:i2,j1:j2) > 180. )
4231               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360.
4232            END WHERE
4233         ENDIF
4234
4235         ! min lognitude of the right cell
4236         dl_lonmin0=dl_lon0(id_imax0-1,id_jmin0)
4237         IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 )THEN
4238
4239            !!!!! i-direction !!!!!
4240            IF( ll_even(jp_I) )THEN
4241               ! even
4242               SELECT CASE(TRIM(cl_point))
4243                  CASE('F','U')
4244                     dl_dlon= ( dl_lon0(id_imax0  ,id_jmin0) -   &
4245                        &       dl_lon0(id_imax0-1,id_jmin0) ) / &
4246                        &     ( 2.*id_rho(jp_I) )
4247                  CASE DEFAULT
4248                     dl_dlon=0
4249               END SELECT
4250            ELSE
4251               ! odd
4252               dl_dlon= ( dl_lon0(id_imax0  ,id_jmin0) -   &
4253                  &       dl_lon0(id_imax0-1,id_jmin0) ) / &
4254                  &     ( 2.*id_rho(jp_I) )
4255            ENDIF
4256
4257            dl_lon0F= dl_lon0(id_imax0-1,id_jmin0) - dl_dlon
4258            dl_lat0F= dd_lat0(id_imax0-1,id_jmin0)
4259
4260            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4261            &                           dl_lon0F, dl_lat0F, 'ri' )
4262
4263            ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1)
4264
4265            !!!!! i-direction !!!!!
4266            IF( ll_even(jp_I) )THEN
4267               ! even
4268               SELECT CASE(TRIM(cl_point))
4269                  CASE('T','V')
4270                     grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii
4271                  CASE DEFAULT !'F','U'
4272                     grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii
4273               END SELECT
4274            ELSE
4275               ! odd
4276               grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii
4277            ENDIF
4278
4279         ELSE
4280            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
4281            &                 " not match fine grid right corner.")
4282         ENDIF
4283
4284         IF( ll_greenwich )THEN
4285            ! close to greenwich meridien
4286            ll_greenwich=.FALSE.
4287            ! -180:180 => 0:360
4288            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) < 0. )
4289               dl_lon0(id_imax0-1:id_imax0,id_jmin0) = &
4290                  & dl_lon0(id_imax0-1:id_imax0,id_jmin0)+360.
4291            END WHERE
4292
4293            WHERE( dl_lon1(i1:i2,j1:j2) < 0. )
4294               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360.
4295            END WHERE
4296         ENDIF
4297
4298      ELSEIF( il_shape1(jp_I) == 1 )THEN
4299         
4300         grid__get_fine_offset_cc(jp_I,:)=((id_rho(jp_I)-1)/2)
4301         
4302         !!! work on j-direction
4303         !!! look for j-direction lower offset
4304         i1=1 ; i2=1
4305         j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J))
4306
4307
4308         ! max latitude of the lower cell
4309         dl_latmax0=dd_lat0(id_imin0,id_jmin0+1)
4310         IF( dd_lat1(1,1) < dl_latmax0 )THEN
4311
4312            IF( ll_even(jp_J) )THEN
4313               ! even
4314               SELECT CASE(TRIM(cl_point))
4315                  CASE('F','V')
4316                     dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) -   &
4317                        &       dd_lat0(id_imin0,id_jmin0  ) ) / &
4318                        &     ( 2.*id_rho(jp_J) )
4319                  CASE DEFAULT
4320                     dl_dlat=0
4321               END SELECT
4322            ELSE
4323               ! odd
4324               dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) -   &
4325                  &       dd_lat0(id_imin0,id_jmin0  ) ) / &
4326                  &     ( 2.*id_rho(jp_J) )
4327            ENDIF
4328
4329            dl_lon0F= dl_lon0(id_imin0,id_jmin0+1)
4330            dl_lat0F= dd_lat0(id_imin0,id_jmin0+1) + dl_dlat 
4331           
4332            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4333            &                           dl_lon0F, dl_lat0F, 'lo' )
4334
4335            ij=il_ind(2)
4336
4337            !!!!! i-direction !!!!!
4338            IF( ll_even(jp_I) )THEN
4339               ! even
4340               SELECT CASE(TRIM(cl_point))
4341                  CASE('T','V')
4342                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii
4343                  CASE DEFAULT !'F','U'
4344                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4345               END SELECT
4346            ELSE
4347               ! odd
4348               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4349            ENDIF
4350
4351         ELSE
4352            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
4353            &                 " not match fine grid lower corner.")
4354         ENDIF
4355
4356         !!! look for j-direction upper offset
4357         i1=1                                         ; i2=1
4358         j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J)
4359
4360         ! min latitude of the upper cell
4361         dl_latmin0=dd_lat0(id_imin0,id_jmax0-1)
4362         IF( dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN
4363
4364            IF( ll_even(jp_J) )THEN
4365               ! even
4366               SELECT CASE(TRIM(cl_point))
4367                  CASE('F','V')
4368                     dl_dlat= ( dd_lat0(id_imin0,id_jmax0  ) -   &
4369                        &       dd_lat0(id_imin0,id_jmax0-1) ) / &
4370                        &     ( 2.*id_rho(jp_J) )
4371                  CASE DEFAULT
4372                     dl_dlat=0
4373               END SELECT
4374            ELSE
4375               ! odd
4376               dl_dlat= ( dd_lat0(id_imin0,id_jmax0  ) -   &
4377                  &       dd_lat0(id_imin0,id_jmax0-1) ) / &
4378                  &     ( 2*id_rho(jp_J) )
4379            ENDIF
4380
4381            dl_lon0F= dl_lon0(id_imin0,id_jmax0-1) 
4382            dl_lat0F= dd_lat0(id_imin0,id_jmax0-1) - dl_dlat
4383           
4384            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4385            &                           dl_lon0F, dl_lat0F, 'up' )
4386
4387            ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1)
4388
4389            !!!!! j-direction !!!!!
4390            IF( ll_even(jp_J) )THEN
4391               ! even
4392               SELECT CASE(TRIM(cl_point))
4393                  CASE('T','U')
4394                     grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij
4395                  CASE DEFAULT !'F','V'
4396                     grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij
4397               END SELECT
4398            ELSE
4399               ! odd
4400               grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij
4401            ENDIF
4402
4403         ELSE
4404            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
4405            &                 " not match fine grid upper corner.")
4406         ENDIF
4407
4408      ELSE ! il_shape1(1) > 1 .AND. il_shape1(2) > 1
4409
4410         !!!!!! look for lower left offset !!!!!!
4411         i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I))
4412         j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J))
4413
4414         ! check if cross greenwich meridien
4415         IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))<5. .OR. &
4416           & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))>355. )THEN
4417            ! close to greenwich meridien
4418            ll_greenwich=.TRUE.
4419            ! 0:360 => -180:180
4420            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) > 180. )
4421               dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = &
4422                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)-360.
4423            END WHERE
4424
4425            WHERE( dl_lon1(i1:i2,j1:j2) > 180. )
4426               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360.
4427            END WHERE
4428         ENDIF
4429
4430         ! max longitude of the lower left cell
4431         dl_lonmax0=MAX(dl_lon0(id_imin0+1,id_jmin0),dl_lon0(id_imin0+1,id_jmin0+1))
4432         ! max latitude of the lower left cell
4433         dl_latmax0=MAX(dd_lat0(id_imin0,id_jmin0+1),dd_lat0(id_imin0+1,id_jmin0+1))
4434         IF( dl_lon1(1,1) < dl_lonmax0 .AND. &
4435           & dd_lat1(1,1) < dl_latmax0 )THEN
4436
4437            !!!!! i-direction !!!!!
4438            IF( ll_even(jp_I) )THEN
4439               ! even
4440               SELECT CASE(TRIM(cl_point))
4441                  CASE('F','U')
4442                     dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) -   &
4443                        &       dl_lon0(id_imin0  ,id_jmin0+1) ) / &
4444                        &     ( 2.*id_rho(jp_I) )
4445                  CASE DEFAULT
4446                     dl_dlon=0
4447               END SELECT
4448            ELSE
4449               ! odd
4450               dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) -   &
4451                  &       dl_lon0(id_imin0  ,id_jmin0+1) ) / &
4452                  &     ( 2.*id_rho(jp_I) )
4453            ENDIF
4454
4455            !!!!! j-direction !!!!!
4456            IF( ll_even(jp_J) )THEN
4457               ! even
4458               SELECT CASE(TRIM(cl_point))
4459                  CASE('F','V')
4460                     dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) -   &
4461                        &       dd_lat0(id_imin0+1,id_jmin0  ) ) / &
4462                        &     ( 2.*id_rho(jp_J) )
4463                  CASE DEFAULT
4464                     dl_dlat=0
4465               END SELECT
4466            ELSE
4467               ! odd
4468               dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) -   &
4469                  &       dd_lat0(id_imin0+1,id_jmin0  ) ) / &
4470                  &     ( 2.*id_rho(jp_J) )
4471            ENDIF
4472
4473            dl_lon0F= dl_lon0(id_imin0+1,id_jmin0+1) + dl_dlon
4474            dl_lat0F= dd_lat0(id_imin0+1,id_jmin0+1) + dl_dlat
4475
4476            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4477            &                           dl_lon0F, dl_lat0F, 'll' )
4478
4479            ii=il_ind(1)
4480            ij=il_ind(2)
4481
4482            !!!!! i-direction !!!!!
4483            IF( ll_even(jp_I) )THEN
4484               ! even
4485               SELECT CASE(TRIM(cl_point))
4486                  CASE('T','V')
4487                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii
4488                  CASE DEFAULT !'F','U'
4489                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4490               END SELECT
4491            ELSE
4492               ! odd
4493               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii
4494            ENDIF
4495
4496            !!!!! j-direction !!!!!
4497            IF( ll_even(jp_J) )THEN
4498               ! even
4499               SELECT CASE(TRIM(cl_point))
4500                  CASE('T','U')
4501                     grid__get_fine_offset_cc(jp_J,1)=id_rho(jp_J)-ij
4502                  CASE DEFAULT !'F','V'
4503                     grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij
4504               END SELECT
4505            ELSE
4506               ! odd
4507               grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij
4508            ENDIF
4509
4510         ELSE
4511            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//&
4512            &                 " not match fine grid lower left corner.")
4513         ENDIF
4514
4515         IF( ll_greenwich )THEN
4516            ! close to greenwich meridien
4517            ll_greenwich=.FALSE.
4518            ! -180:180 => 0:360
4519            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) < 0. )
4520               dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = &
4521                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)+360.
4522            END WHERE
4523
4524            WHERE( dl_lon1(i1:i2,j1:j2) < 0. )
4525               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360.
4526            END WHERE
4527         ENDIF
4528
4529         !!!!!! look for upper right offset !!!!!!
4530         i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I)
4531         j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J)
4532
4533         ! check if cross greenwich meridien
4534         IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))<5. .OR. &
4535           & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))>355. )THEN
4536            ! close to greenwich meridien
4537            ll_greenwich=.TRUE.
4538            ! 0:360 => -180:180
4539            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) > 180. )
4540               dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = &
4541                  & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)-360.
4542            END WHERE
4543
4544            WHERE( dl_lon1(i1:i2,j1:j2) > 180. )
4545               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360.
4546            END WHERE
4547         ENDIF
4548
4549         ! min latitude of the upper right cell
4550         dl_lonmin0=MIN(dl_lon0(id_imax0-1,id_jmax0-1),dl_lon0(id_imax0-1,id_jmax0))
4551         ! min latitude of the upper right cell
4552         dl_latmin0=MIN(dd_lat0(id_imax0-1,id_jmax0-1),dd_lat0(id_imax0,id_jmax0-1))
4553         IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 .AND. &
4554           & dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN
4555
4556            !!!!! i-direction !!!!!
4557            IF( ll_even(jp_I) )THEN
4558               ! even
4559               SELECT CASE(TRIM(cl_point))
4560                  CASE('F','U')
4561                     dl_dlon= ( dl_lon0(id_imax0  ,id_jmax0-1) -   &
4562                        &       dl_lon0(id_imax0-1,id_jmax0-1) ) / &
4563                        &     ( 2.*id_rho(jp_I) )
4564                  CASE DEFAULT
4565                     dl_dlon=0
4566               END SELECT               
4567            ELSE
4568               ! odd
4569               dl_dlon= ( dl_lon0(id_imax0  ,id_jmax0-1) -   &
4570                  &       dl_lon0(id_imax0-1,id_jmax0-1) ) / &
4571                  &     ( 2*id_rho(jp_I) )
4572            ENDIF
4573
4574            !!!!! j-direction !!!!!
4575            IF( ll_even(jp_J) )THEN
4576               ! even
4577               SELECT CASE(TRIM(cl_point))
4578                  CASE('F','V')
4579                     dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0  ) -   &
4580                        &       dd_lat0(id_imax0-1,id_jmax0-1) ) / &
4581                        &     ( 2.*id_rho(jp_J) )
4582                  CASE DEFAULT
4583                     dl_dlat=0
4584               END SELECT
4585            ELSE
4586               ! odd
4587               dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0  ) -   &
4588                  &       dd_lat0(id_imax0-1,id_jmax0-1) ) / &
4589                  &     ( 2*id_rho(jp_J) )
4590            ENDIF
4591
4592            dl_lon0F= dl_lon0(id_imax0-1,id_jmax0-1) - dl_dlon
4593            dl_lat0F= dd_lat0(id_imax0-1,id_jmax0-1) - dl_dlat
4594
4595            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), &
4596            &                           dl_lon0F, dl_lat0F, 'ur' )
4597
4598            ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1)
4599            ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1)
4600
4601            !!!!! i-direction !!!!!
4602            IF( ll_even(jp_I) )THEN
4603               ! even
4604               SELECT CASE(TRIM(cl_point))
4605                  CASE('T','V')
4606                     grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii
4607                  CASE DEFAULT !'F','U'
4608                     grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii
4609               END SELECT
4610            ELSE
4611               ! odd
4612               grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii
4613            ENDIF
4614
4615            !!!!! j-direction !!!!!
4616            IF( ll_even(jp_J) )THEN
4617               ! even
4618               SELECT CASE(TRIM(cl_point))
4619                  CASE('T','U')
4620                     grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij
4621                  CASE DEFAULT !'F','V'
4622                     grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij
4623               END SELECT
4624            ELSE
4625               ! odd
4626               grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij
4627            ENDIF
4628
4629         ELSE
4630            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//&
4631            &                 " not match fine grid upper right corner.")
4632         ENDIF
4633
4634         IF( ll_greenwich )THEN
4635            ! close to greenwich meridien
4636            ll_greenwich=.FALSE.
4637            ! -180:180 => 0:360
4638            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) < 0. )
4639               dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = &
4640                  & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)+360.
4641            END WHERE
4642
4643            WHERE( dl_lon1(i1:i2,j1:j2) < 0. )
4644               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360.
4645            END WHERE
4646         ENDIF
4647
4648      ENDIF
4649
4650      DEALLOCATE( dl_lon0 )
4651      DEALLOCATE( dl_lon1 )
4652
4653      IF( ANY(grid__get_fine_offset_cc(:,:)==-1) )THEN
4654         CALL logger_fatal("GRID GET FINE OFFSET: can not found "//&
4655         &                 " offset between coarse and fine grid.")
4656      ENDIF
4657
4658   END FUNCTION grid__get_fine_offset_cc
4659   !-------------------------------------------------------------------
4660   !> @brief This subroutine check fine and coarse grid coincidence.
4661   !
4662   !> @details
4663   !
4664   !> @author J.Paul
4665   !> @date November, 2013- Initial Version
4666   !> @date October, 2014
4667   !> - work on mpp file structure instead of file structure
4668   !> @date February, 2016
4669   !> - use F-point to check coincidence for even refinment
4670   !> - use F-point estimation, if can not read it.
4671   !>
4672   !> @param[in] td_coord0 coarse grid coordinate file structure
4673   !> @param[in] td_coord1 fine   grid coordinate file structure
4674   !> @param[in] id_imin0  coarse grid lower left  corner i-indice of fine grid domain
4675   !> @param[in] id_imax0  coarse grid upper right corner i-indice of fine grid domain
4676   !> @param[in] id_jmin0  coarse grid lower left  corner j-indice of fine grid domain
4677   !> @param[in] id_jmax0  coarse grid upper right corner j-indice of fine grid domain 
4678   !> @param[in] id_rho    array of refinement factor
4679   !-------------------------------------------------------------------
4680   SUBROUTINE grid_check_coincidence( td_coord0, td_coord1, &
4681   &                                  id_imin0, id_imax0, &
4682   &                                  id_jmin0, id_jmax0, &
4683   &                                  id_rho )
4684      IMPLICIT NONE
4685     
4686      ! Argument     
4687      TYPE(TMPP)               , INTENT(IN) :: td_coord0
4688      TYPE(TMPP)               , INTENT(IN) :: td_coord1
4689      INTEGER(i4)              , INTENT(IN) :: id_imin0
4690      INTEGER(i4)              , INTENT(IN) :: id_imax0
4691      INTEGER(i4)              , INTENT(IN) :: id_jmin0
4692      INTEGER(i4)              , INTENT(IN) :: id_jmax0
4693      INTEGER(i4), DIMENSION(:), INTENT(IN) :: id_rho
4694
4695      ! local variable
4696      INTEGER(i4)               :: il_imid1
4697      INTEGER(i4)               :: il_jmid1
4698     
4699      INTEGER(i4)               :: il_ew0
4700      INTEGER(i4)               :: il_ew1
4701
4702      INTEGER(i4)               :: il_ind
4703
4704      INTEGER(i4)               :: il_imin1
4705      INTEGER(i4)               :: il_imax1
4706      INTEGER(i4)               :: il_jmin1
4707      INTEGER(i4)               :: il_jmax1
4708
4709      INTEGER(i4), DIMENSION(2) :: il_ind0
4710      INTEGER(i4), DIMENSION(2) :: il_ind1
4711
4712      INTEGER(i4), DIMENSION(2) :: il_ill1
4713      INTEGER(i4), DIMENSION(2) :: il_ilr1
4714      INTEGER(i4), DIMENSION(2) :: il_iul1
4715      INTEGER(i4), DIMENSION(2) :: il_iur1
4716
4717      REAL(dp)                  :: dl_lon0F
4718      REAL(dp)                  :: dl_lat0F
4719      REAL(dp)                  :: dl_lon0
4720      REAL(dp)                  :: dl_lat0
4721      REAL(dp)                  :: dl_lon1F
4722      REAL(dp)                  :: dl_lat1F
4723      REAL(dp)                  :: dl_lon1
4724      REAL(dp)                  :: dl_lat1
4725
4726      REAL(dp)                  :: dl_delta
4727
4728      LOGICAL                   :: ll_coincidence
4729      LOGICAL                   :: ll_even
4730      LOGICAL                   :: ll_grid0F
4731      LOGICAL                   :: ll_grid1F
4732
4733      TYPE(TVAR)                :: tl_lon0
4734      TYPE(TVAR)                :: tl_lat0
4735      TYPE(TVAR)                :: tl_lon0F
4736      TYPE(TVAR)                :: tl_lat0F
4737      TYPE(TVAR)                :: tl_lon1
4738      TYPE(TVAR)                :: tl_lat1
4739      TYPE(TVAR)                :: tl_lon1F
4740      TYPE(TVAR)                :: tl_lat1F
4741
4742      TYPE(TMPP)                :: tl_coord0
4743      TYPE(TMPP)                :: tl_coord1
4744
4745      TYPE(TDOM)                :: tl_dom0
4746
4747      ! loop indices
4748      INTEGER(i4) :: ji
4749      INTEGER(i4) :: jj
4750      !----------------------------------------------------------------
4751
4752      ll_coincidence=.TRUE.
4753
4754      ll_even=.FALSE.
4755      IF( MOD(id_rho(jp_I)*id_rho(jp_J),2) == 0 )THEN
4756         ll_even=.TRUE.
4757      ENDIF
4758
4759      ! copy structure
4760      tl_coord0=mpp_copy(td_coord0)
4761
4762      ! compute domain
4763      tl_dom0=dom_init( tl_coord0,         &
4764      &                 id_imin0, id_imax0,&
4765      &                 id_jmin0, id_jmax0 )
4766
4767      ! open mpp files
4768      CALL iom_dom_open(tl_coord0, tl_dom0)
4769
4770      ! read variable value on domain
4771      il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_T')
4772      IF( il_ind /= 0 )THEN
4773         tl_lon0=iom_dom_read_var(tl_coord0,'longitude_T',tl_dom0)
4774      ELSE
4775         tl_lon0=iom_dom_read_var(tl_coord0,'longitude',tl_dom0)
4776      ENDIF
4777
4778      il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_T')
4779      IF( il_ind /= 0 )THEN
4780         tl_lat0=iom_dom_read_var(tl_coord0,'latitude_T' ,tl_dom0)
4781      ELSE
4782         tl_lat0=iom_dom_read_var(tl_coord0,'latitude' ,tl_dom0)
4783      ENDIF
4784
4785      IF( ll_even )THEN
4786
4787         ! look for variable value on domain for F point
4788         il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_F')
4789         IF( il_ind /= 0 )THEN
4790            tl_lon0F=iom_dom_read_var(tl_coord0,'longitude_F',tl_dom0)
4791         ENDIF
4792
4793         il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_F')
4794         IF( il_ind /= 0 )THEN
4795            tl_lat0F=iom_dom_read_var(tl_coord0,'latitude_F' ,tl_dom0)
4796         ENDIF
4797
4798         ll_grid0F=.FALSE.
4799         IF( ASSOCIATED(tl_lon0F%d_value) .AND. &
4800         &   ASSOCIATED(tl_lat0F%d_value) )THEN
4801            ll_grid0F=.TRUE.
4802         ENDIF
4803
4804      ENDIF
4805
4806      ! close mpp files
4807      CALL iom_dom_close(tl_coord0)
4808
4809      ! clean structure
4810      CALL mpp_clean(tl_coord0)
4811      CALL dom_clean(tl_dom0)
4812
4813      ! copy structure
4814      tl_coord1=mpp_copy(td_coord1)
4815
4816      ! open mpp files
4817      CALL iom_mpp_open(tl_coord1)
4818
4819      ! read fine longitue and latitude
4820      il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lon0%c_longname))
4821      IF( il_ind /= 0 )THEN
4822         tl_lon1=iom_mpp_read_var(tl_coord1,TRIM(tl_lon0%c_longname))
4823      ELSE
4824         tl_lon1=iom_mpp_read_var(tl_coord1,'longitude')
4825      ENDIF
4826      il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lat0%c_longname))
4827      IF( il_ind /= 0 )THEN
4828         tl_lat1=iom_mpp_read_var(tl_coord1,TRIM(tl_lat0%c_longname))
4829      ELSE
4830         tl_lat1=iom_mpp_read_var(tl_coord1,'latitude')
4831      ENDIF
4832     
4833      IF( ll_even )THEN
4834
4835         ! look for variable value on domain for F point
4836         il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'longitude_F')
4837         IF( il_ind /= 0 )THEN
4838            tl_lon1F=iom_mpp_read_var(tl_coord1,'longitude_F')
4839         ENDIF
4840
4841         il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'latitude_F')
4842         IF( il_ind /= 0 )THEN
4843            tl_lat1F=iom_mpp_read_var(tl_coord1,'latitude_F')
4844         ENDIF
4845
4846         ll_grid1F=.FALSE.
4847         IF( ASSOCIATED(tl_lon1F%d_value) .AND. &
4848         &   ASSOCIATED(tl_lat1F%d_value) )THEN
4849            ll_grid1F=.TRUE.
4850         ENDIF
4851
4852      ENDIF
4853
4854      ! close mpp files
4855      CALL iom_mpp_close(tl_coord1)
4856      ! clean structure
4857      CALL mpp_clean(tl_coord1)
4858
4859      CALL logger_debug("GRID CHECK COINCIDENCE:"//&
4860      &        " fine   grid "//TRIM(td_coord1%c_name) )
4861      CALL logger_debug("GRID CHECK COINCIDENCE:"//&
4862      &        " coarse grid "//TRIM(td_coord0%c_name) )
4863
4864      ! check domain
4865      ! check global grid
4866      IF( .NOT. grid_is_global(tl_lon0, tl_lat0) )THEN
4867         IF( grid_is_global(tl_lon1, tl_lat1) )THEN
4868
4869            ll_coincidence=.FALSE.
4870            CALL logger_fatal("GRID CHECK COINCIDENCE:"//&
4871            &        " fine   grid is global,"//&
4872            &        " coarse grid is not ")
4873
4874         ELSE
4875            il_ew1=tl_lon1%i_ew
4876            IF( il_ew1 >= 0 )THEN
4877               ! ew overlap
4878
4879               il_ew0=tl_lon0%i_ew
4880               IF( il_ew0 < 0 )THEN
4881                  CALL logger_fatal("GRID CHECK COINCIDENCE: "//&
4882                  &        "fine grid has east west overlap,"//&
4883                  &        " coarse grid not ")
4884               ENDIF
4885
4886               il_jmin1=1+ip_ghost
4887               il_jmax1=tl_lon1%t_dim(2)%i_len-ip_ghost
4888
4889               ll_coincidence=grid__check_lat(&
4890               &                     tl_lat0%d_value(1,:,1,1),&
4891               &                     tl_lat1%d_value(1,il_jmin1:il_jmax1,1,1))
4892
4893            ELSE
4894               ! other case
4895               il_imin1=1+ip_ghost
4896               il_jmin1=1+ip_ghost
4897
4898               il_imax1=tl_lon1%t_dim(1)%i_len-ip_ghost
4899               il_jmax1=tl_lon1%t_dim(2)%i_len-ip_ghost
4900
4901               ll_coincidence=grid__check_corner(&
4902               &                      tl_lon0%d_value(:,:,1,1),&
4903               &                      tl_lat0%d_value(:,:,1,1),&
4904               &                      tl_lon1%d_value(il_imin1:il_imax1, &
4905               &                                      il_jmin1:il_jmax1, &
4906               &                                      1,1),&
4907               &                      tl_lat1%d_value(il_imin1:il_imax1, &
4908               &                                      il_jmin1:il_jmax1, &
4909               &                                      1,1) )
4910
4911            ENDIF
4912 
4913         ENDIF
4914
4915         IF( .NOT. ll_coincidence )THEN
4916            CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
4917            &              "between fine grid and coarse grid: invalid domain." )
4918         ENDIF
4919
4920      ENDIF
4921 
4922      ! check refinement factor
4923      ! select point in middle of fine grid
4924      il_imid1=INT(tl_lon1%t_dim(1)%i_len*0.5)
4925      il_jmid1=INT(tl_lon1%t_dim(2)%i_len*0.5)
4926 
4927      dl_lon1=tl_lon1%d_value(il_imid1, il_jmid1,1,1)
4928      dl_lat1=tl_lat1%d_value(il_imid1, il_jmid1,1,1)
4929
4930      ! select closest point on coarse grid
4931      il_ind0(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),&
4932      &                           tl_lat0%d_value(:,:,1,1),&
4933      &                           dl_lon1, dl_lat1   )
4934
4935      IF( ANY(il_ind0(:)==0) )THEN
4936         CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//&
4937         &              "coarse grid indices: invalid domain." )
4938      ENDIF
4939
4940      IF( .NOT. ll_even )THEN
4941         ! case odd refinment in both direction
4942         ! work on T-point
4943
4944         dl_lon0=tl_lon0%d_value(il_ind0(1),il_ind0(2),1,1)
4945         dl_lat0=tl_lat0%d_value(il_ind0(1),il_ind0(2),1,1)
4946
4947         il_ind1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),&
4948         &                           tl_lat1%d_value(:,:,1,1),&
4949         &                           dl_lon0, dl_lat0 )
4950
4951         ! check i-direction refinement factor
4952         DO ji=0,MIN(3,il_imid1)
4953
4954            IF( il_ind1(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN
4955               CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
4956               &  " to check i-direction refinement factor ")
4957               EXIT
4958            ELSE
4959               dl_lon0=tl_lon0%d_value(il_ind0(1)+ji             ,il_ind0(2),1,1)
4960               dl_lon1=tl_lon1%d_value(il_ind1(1)+ji*id_rho(jp_I),il_ind1(2),1,1)
4961
4962               ! assume there could be little difference due to interpolation
4963               IF( ABS(dl_lon1 - dl_lon0) > dp_delta )THEN
4964                  ll_coincidence=.FALSE.
4965                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
4966                  &  "i-direction refinement factor ("//&
4967                  &   TRIM(fct_str(id_rho(jp_I)))//&
4968                  &  ") between fine grid and coarse grid ")
4969               ENDIF
4970            ENDIF
4971
4972         ENDDO
4973
4974         ! check j-direction refinement factor
4975         DO jj=0,MIN(3,il_jmid1)
4976
4977            IF( il_ind1(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN
4978               CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
4979               &  " to check j-direction refinement factor ")
4980               EXIT
4981            ELSE     
4982               dl_lat0=tl_lat0%d_value(il_ind0(1),il_ind0(2)+jj             ,1,1)
4983               dl_lat1=tl_lat1%d_value(il_ind1(1),il_ind1(2)+jj*id_rho(jp_J),1,1)
4984
4985               ! assume there could be little difference due to interpolation
4986               IF( ABS(dl_lat1-dl_lat0) > dp_delta )THEN
4987                  ll_coincidence=.FALSE.
4988                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
4989                  &  "j-direction refinement factor ("//&
4990                  &   TRIM(fct_str(id_rho(jp_J)))//&
4991                  &  ") between fine grid and coarse grid ")
4992               ENDIF
4993            ENDIF
4994
4995         ENDDO
4996
4997      ELSE
4998         ! case even refinment at least in one direction
4999         ! work on F-point
5000
5001         dl_delta=dp_delta
5002         ! look for lower left fine point in coarse cell.
5003         IF( ll_grid0F )THEN
5004     
5005            ! lower left corner of coarse cell
5006            dl_lon0F=tl_lon0F%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1)
5007            dl_lat0F=tl_lat0F%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1)
5008
5009         ELSE
5010
5011            ! approximate lower left corner of coarse cell (with T point)
5012            dl_lon0F=( tl_lon0%d_value(il_ind0(1)  ,il_ind0(2)  ,1,1) + &
5013            &          tl_lon0%d_value(il_ind0(1)  ,il_ind0(2)-1,1,1) + &
5014            &          tl_lon0%d_value(il_ind0(1)-1,il_ind0(2)  ,1,1) + &
5015            &          tl_lon0%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) ) * 0.25
5016
5017            dl_lat0F=( tl_lat0%d_value(il_ind0(1)  ,il_ind0(2)  ,1,1) + &
5018            &          tl_lat0%d_value(il_ind0(1)  ,il_ind0(2)-1,1,1) + &
5019            &          tl_lat0%d_value(il_ind0(1)-1,il_ind0(2)  ,1,1) + &
5020            &          tl_lat0%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) ) * 0.25
5021
5022            ! as we use approximation of F-point we relax condition
5023            dl_delta=100*dp_delta
5024
5025         ENDIF
5026
5027         IF( ll_grid1F )THEN
5028     
5029            il_ind1(:)=grid_get_closest(tl_lon1F%d_value(:,:,1,1),&
5030            &                           tl_lat1F%d_value(:,:,1,1),&
5031            &                           dl_lon0F, dl_lat0F )
5032
5033         ELSE
5034
5035            il_ill1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),&
5036            &                           tl_lat1%d_value(:,:,1,1),&
5037            &                           dl_lon0F, dl_lat0F, 'll' )
5038
5039            il_ilr1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),&
5040            &                           tl_lat1%d_value(:,:,1,1),&
5041            &                           dl_lon0F, dl_lat0F, 'lr' )
5042
5043            il_iul1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),&
5044            &                           tl_lat1%d_value(:,:,1,1),&
5045            &                           dl_lon0F, dl_lat0F, 'ul' )
5046
5047            il_iur1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),&
5048            &                           tl_lat1%d_value(:,:,1,1),&
5049            &                           dl_lon0F, dl_lat0F, 'ur' )
5050
5051            ! as we use approximation of F-point we relax condition
5052            dl_delta=100*dp_delta
5053
5054         ENDIF
5055
5056         ! check i-direction refinement factor
5057         DO ji=0,MIN(3,il_imid1)
5058
5059            IF( il_ind1(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN
5060               CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
5061               &  " to check i-direction refinement factor ")
5062               EXIT
5063            ELSE
5064               IF( ll_grid0F )THEN
5065                  dl_lon0F=tl_lon0F%d_value(il_ind0(1)+ji-1, il_ind0(2)-1,1,1)
5066               ELSE
5067                  dl_lon0F= 0.25 * &
5068                  & ( tl_lon0%d_value(il_ind0(1)+ji  , il_ind0(2)  ,1,1) + &
5069                  &   tl_lon0%d_value(il_ind0(1)+ji-1, il_ind0(2)  ,1,1) + &
5070                  &   tl_lon0%d_value(il_ind0(1)+ji  , il_ind0(2)-1,1,1) + &
5071                  &   tl_lon0%d_value(il_ind0(1)+ji-1, il_ind0(2)-1,1,1) )
5072               ENDIF
5073
5074               IF( ll_grid1F )THEN
5075                  dl_lon1F= tl_lon1F%d_value( il_ind1(1)+ji*id_rho(jp_I), &
5076                                            & il_ind1(2),1,1)
5077               ELSE
5078                  dl_lon1F= 0.25 * &
5079                  & ( tl_lon1%d_value( il_ill1(1)+ji*id_rho(jp_I), &
5080                                     & il_ill1(2),1,1) + &
5081                  &   tl_lon1%d_value( il_ilr1(1)+ji*id_rho(jp_I), &
5082                                     & il_ilr1(2),1,1) + &
5083                  &   tl_lon1%d_value( il_iul1(1)+ji*id_rho(jp_I), &
5084                                     & il_iul1(2),1,1) + &
5085                  &   tl_lon1%d_value( il_iur1(1)+ji*id_rho(jp_I), &
5086                                     & il_iur1(2),1,1) )
5087                 
5088               ENDIF
5089
5090               ! assume there could be little difference due to interpolation
5091               IF( ABS(dl_lon1F - dl_lon0F) > dl_delta )THEN
5092                  ll_coincidence=.FALSE.
5093                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
5094                  &  "i-direction refinement factor ("//&
5095                  &   TRIM(fct_str(id_rho(jp_I)))//&
5096                  &  ") between fine grid and coarse grid ")
5097               ENDIF
5098            ENDIF
5099
5100         ENDDO
5101
5102         ! check j-direction refinement factor
5103         DO jj=0,MIN(3,il_jmid1)
5104
5105            IF( il_ind1(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN
5106               CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
5107               &  " to check j-direction refinement factor ")
5108               EXIT
5109            ELSE     
5110               IF( ll_grid0F )THEN
5111                  dl_lat0F=tl_lat0F%d_value(il_ind0(1)-1, il_ind0(2)+jj-1,1,1)
5112               ELSE
5113                  dl_lat0F= 0.25 * &
5114                  & ( tl_lat0%d_value(il_ind0(1)  , il_ind0(2)+jj  ,1,1) + &
5115                  &   tl_lat0%d_value(il_ind0(1)-1, il_ind0(2)+jj  ,1,1) + &
5116                  &   tl_lat0%d_value(il_ind0(1)  , il_ind0(2)+jj-1,1,1) + &
5117                  &   tl_lat0%d_value(il_ind0(1)-1, il_ind0(2)+jj-1,1,1) )
5118               ENDIF
5119
5120               IF( ll_grid1F )THEN
5121                  dl_lat1F= tl_lat1F%d_value( il_ind1(1), &
5122                                            & il_ind1(2)+jj*id_rho(jp_J),1,1)
5123               ELSE
5124                  dl_lat1F= 0.25 * &
5125                  & ( tl_lat1%d_value( il_ill1(1), &
5126                                     & il_ill1(2)+jj*id_rho(jp_J),1,1) + &
5127                  &   tl_lat1%d_value( il_ilr1(1), &
5128                                     & il_ilr1(2)+jj*id_rho(jp_J),1,1) + &
5129                  &   tl_lat1%d_value( il_iul1(1), &
5130                                     & il_iul1(2)+jj*id_rho(jp_J),1,1) + &
5131                  &   tl_lat1%d_value( il_iur1(1), &
5132                                     & il_iur1(2)+jj*id_rho(jp_J),1,1) )
5133                 
5134               ENDIF
5135
5136               ! assume there could be little difference due to interpolation
5137               IF( ABS(dl_lat1F - dl_lat0F) > dl_delta )THEN
5138                  ll_coincidence=.FALSE.
5139                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
5140                  &  "i-direction refinement factor ("//&
5141                  &   TRIM(fct_str(id_rho(jp_I)))//&
5142                  &  ") between fine grid and coarse grid ")
5143               ENDIF
5144            ENDIF
5145
5146         ENDDO
5147      ENDIF
5148
5149      ! clean
5150      CALL var_clean(tl_lon1)
5151      CALL var_clean(tl_lat1)
5152      CALL var_clean(tl_lon0)
5153      CALL var_clean(tl_lat0)
5154
5155      IF( .NOT. ll_coincidence )THEN
5156         CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
5157         &              "between fine and coarse grid: "//&
5158         &              "invalid refinement factor" )
5159      ENDIF
5160
5161   END SUBROUTINE grid_check_coincidence
5162   !-------------------------------------------------------------------
5163   !> @brief This function check that fine grid is
5164   !> inside coarse grid
5165   !
5166   !> @details
5167   !>
5168   !> @author J.Paul
5169   !> @date November, 2013 - Initial Version
5170   !
5171   !> @param[in] dd_lon0   array of coarse grid longitude
5172   !> @param[in] dd_lat0   array of coarse grid latitude
5173   !> @param[in] dd_lon1   array of fine   grid longitude
5174   !> @param[in] dd_lat1   array of fine   grid latitude
5175   !> @return logical, fine grid is inside coarse grid
5176   !-------------------------------------------------------------------
5177   FUNCTION grid__check_corner(dd_lon0, dd_lat0, &
5178   &                           dd_lon1, dd_lat1 )
5179   IMPLICIT NONE
5180      ! Argument     
5181      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
5182      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
5183      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon1
5184      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat1
5185
5186      ! function
5187      LOGICAL :: grid__check_corner
5188
5189      ! local variable
5190      INTEGER(i4), DIMENSION(2) :: il_shape0
5191      INTEGER(i4), DIMENSION(2) :: il_shape1
5192
5193      INTEGER(i4) :: il_imin0
5194      INTEGER(i4) :: il_jmin0
5195      INTEGER(i4) :: il_imax0
5196      INTEGER(i4) :: il_jmax0
5197
5198      INTEGER(i4) :: il_imin1
5199      INTEGER(i4) :: il_jmin1
5200      INTEGER(i4) :: il_imax1
5201      INTEGER(i4) :: il_jmax1
5202
5203      REAL(dp)    :: dl_lon0
5204      REAL(dp)    :: dl_lat0
5205
5206      REAL(dp)    :: dl_lon1
5207      REAL(dp)    :: dl_lat1
5208      ! loop indices
5209      !----------------------------------------------------------------
5210
5211      ! init
5212      grid__check_corner=.TRUE.
5213
5214      il_shape0=SHAPE(dd_lon0(:,:))
5215      il_shape1=SHAPE(dd_lon1(:,:))
5216
5217      !1- check if fine grid inside coarse grid domain
5218      il_imin0=1 ; il_imax0=il_shape0(1)
5219      il_jmin0=1 ; il_jmax0=il_shape0(2)
5220
5221      il_imin1=1 ; il_imax1=il_shape1(1)
5222      il_jmin1=1 ; il_jmax1=il_shape1(2)
5223
5224      ! check lower left corner
5225      dl_lon0 = dd_lon0(il_imin0, il_jmin0)
5226      dl_lat0 = dd_lat0(il_imin0, il_jmin0)
5227
5228      dl_lon1 = dd_lon1(il_imin1, il_jmin1)
5229      dl_lat1 = dd_lat1(il_imin1, il_jmin1)
5230
5231      IF( (ABS(dl_lon1-dl_lon0)>dp_delta) .AND. (dl_lon1 < dl_lon0 ) .OR. & 
5232      &   (ABS(dl_lat1-dl_lat0)>dp_delta) .AND. (dl_lat1 < dl_lat0 ) )THEN
5233
5234         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower left "//&
5235         &     "corner  not north east of coarse grid (imin,jmin) ")
5236         CALL logger_debug(" fine   grid lower left ( "//&
5237         &              TRIM(fct_str(dl_lon1))//","//&
5238         &              TRIM(fct_str(dl_lat1))//")" )
5239         CALL logger_debug(" coarse grid lower left ( "//&
5240         &              TRIM(fct_str(dl_lon0))//","//&
5241         &              TRIM(fct_str(dl_lat0))//")" )
5242         grid__check_corner=.FALSE.
5243
5244      ENDIF
5245
5246      ! check upper left corner
5247      dl_lon0 = dd_lon0(il_imin0, il_jmax0)
5248      dl_lat0 = dd_lat0(il_imin0, il_jmax0)
5249
5250      dl_lon1 = dd_lon1(il_imin1, il_jmax1)
5251      dl_lat1 = dd_lat1(il_imin1, il_jmax1)
5252
5253
5254      IF( (ABS(dl_lon1-dl_lon0)>dp_delta) .AND. (dl_lon1 < dl_lon0) .OR. &
5255      &   (ABS(dl_lat1-dl_lat0)>dp_delta) .AND. (dl_lat1 > dl_lat0) )THEN
5256
5257         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper left "//&
5258         &     "corner not south east of coarse grid (imin,jmax) ")
5259         CALL logger_debug(" fine   grid upper left ("//&
5260         &              TRIM(fct_str(dl_lon1))//","//&
5261         &              TRIM(fct_str(dl_lat1))//")")
5262         CALL logger_debug(" coasre grid upper left ("//&
5263         &              TRIM(fct_str(dl_lon0))//","//&
5264         &              TRIM(fct_str(dl_lat0))//")")
5265         grid__check_corner=.FALSE.
5266
5267      ENDIF
5268
5269      ! check lower right corner
5270      dl_lon0 = dd_lon0(il_imax0, il_jmin0)
5271      dl_lat0 = dd_lat0(il_imax0, il_jmin0)
5272
5273      dl_lon1 = dd_lon1(il_imax1, il_jmin1)
5274      dl_lat1 = dd_lat1(il_imax1, il_jmin1)
5275
5276
5277      IF( (ABS(dl_lon1-dl_lon0)>dp_delta) .AND. (dl_lon1 > dl_lon0) .OR. &
5278      &   (ABS(dl_lat1-dl_lat0)>dp_delta) .AND. (dl_lat1 < dl_lat0) )THEN
5279
5280         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower right "//&
5281         &     "corner not north west of coarse grid (imax,jmin) ")
5282         CALL logger_debug(" fine   grid lower right ( "//&
5283         &              TRIM(fct_str(dl_lon1))//","//&
5284         &              TRIM(fct_str(dl_lat1))//")" )
5285         CALL logger_debug(" coarse grid lower right ( "//&
5286         &              TRIM(fct_str(dl_lon0))//","//&
5287         &              TRIM(fct_str(dl_lat0))//")" )   
5288         grid__check_corner=.FALSE.
5289
5290      ENDIF
5291
5292      ! check upper right corner
5293      dl_lon0 = dd_lon0(il_imax0, il_jmax0)
5294      dl_lat0 = dd_lat0(il_imax0, il_jmax0)
5295
5296      dl_lon1 = dd_lon1(il_imax1, il_jmax1)
5297      dl_lat1 = dd_lat1(il_imax1, il_jmax1)
5298
5299      IF( (ABS(dl_lon1-dl_lon0)>dp_delta) .AND. (dl_lon1 > dl_lon0) .OR. &
5300      &   (ABS(dl_lat1-dl_lat0)>dp_delta) .AND. (dl_lat1 > dl_lat0) )THEN
5301
5302         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper right "//&
5303         &     "corner not south west of coarse grid (imax,jmax) ")
5304         CALL logger_debug(" fine   grid upper right ( "//&
5305         &              TRIM(fct_str(dl_lon1))//","//&
5306         &              TRIM(fct_str(dl_lat1))//")" )
5307         CALL logger_debug(" fine   imax1 jmax1 ( "//&
5308         &              TRIM(fct_str(il_imax1))//","//&
5309         &              TRIM(fct_str(il_jmax1))//")" )
5310         CALL logger_debug(" coarse grid upper right ( "//&
5311         &              TRIM(fct_str(dl_lon0))//","//&
5312         &              TRIM(fct_str(dl_lat0))//")" )   
5313         CALL logger_debug(" fine   imax0 jmax0 ( "//&
5314         &              TRIM(fct_str(il_imax0))//","//&
5315         &              TRIM(fct_str(il_jmax0))//")" )
5316         grid__check_corner=.FALSE.
5317
5318      ENDIF     
5319
5320   END FUNCTION grid__check_corner
5321   !-------------------------------------------------------------------
5322   !> @brief This function check that fine grid latitude are
5323   !> inside coarse grid latitude
5324   !
5325   !> @details
5326   !
5327   !> @author J.Paul
5328   !> @date November, 2013 - Initial Version
5329   !
5330   !> @param[in] dd_lat0   array of coarse grid latitude
5331   !> @param[in] dd_lat1   array of fine grid latitude
5332   !-------------------------------------------------------------------
5333   FUNCTION grid__check_lat(dd_lat0, dd_lat1)
5334      IMPLICIT NONE
5335      ! Argument
5336      REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat0
5337      REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat1
5338
5339      ! function
5340      LOGICAL :: grid__check_lat
5341
5342      ! local variable
5343      INTEGER(i4), DIMENSION(1) :: il_shape0
5344      INTEGER(i4), DIMENSION(1) :: il_shape1
5345
5346      INTEGER(i4) :: il_jmin0
5347      INTEGER(i4) :: il_jmax0
5348
5349      INTEGER(i4) :: il_jmin1
5350      INTEGER(i4) :: il_jmax1
5351      ! loop indices
5352      !----------------------------------------------------------------
5353
5354      ! init
5355      grid__check_lat=.TRUE.
5356
5357      il_shape0(:)=SHAPE(dd_lat0(:))
5358      il_shape1(:)=SHAPE(dd_lat1(:))
5359
5360      !1- check if fine grid inside coarse grid domain
5361      il_jmin0=1 ; il_jmax0=il_shape0(1)
5362      il_jmin1=1 ; il_jmax1=il_shape1(1)
5363
5364      ! check lower left fine grid
5365      IF( ABS(dd_lat1(il_jmin1)-dd_lat0(il_jmin0)) > dp_delta .AND. &
5366      &   dd_lat1(il_jmin1) < dd_lat0(il_jmin0) )THEN
5367
5368         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower point"//&
5369         &     " not north of coarse grid (jmin) ")
5370         CALL logger_debug(" fine grid lower point ( "//&
5371         &              TRIM(fct_str(dd_lat1(il_jmin1)))//")" )
5372         CALL logger_debug(" coarse grid lower point ( "//&
5373         &              TRIM(fct_str(dd_lat0(il_jmin0)))//")" )
5374         grid__check_lat=.FALSE.
5375
5376      ENDIF
5377
5378      ! check upper left fine grid
5379      IF( ABS(dd_lat1(il_jmax1)-dd_lat0(il_jmax0)) > dp_delta .AND. &
5380      &   dd_lat1(il_jmax1) > dd_lat0(il_jmax0) )THEN
5381
5382         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper point"//&
5383         &     " not south of coarse grid (jmax) ")
5384         CALL logger_debug(" fine grid upper point ("//&
5385         &              TRIM(fct_str(dd_lat1(il_jmax1)))//")")
5386         CALL logger_debug(" coasre grid upper point ("//&
5387         &              TRIM(fct_str(dd_lat0(il_jmax0)))//")")
5388         grid__check_lat=.FALSE.
5389
5390      ENDIF
5391     
5392   END FUNCTION grid__check_lat
5393   !-------------------------------------------------------------------
5394   !> @brief
5395   !> This subroutine add ghost cell at boundaries.
5396   !>
5397   !> @author J.Paul
5398   !> @date November, 2013 - Initial version
5399   !
5400   !> @param[inout] td_var array of variable structure
5401   !> @param[in] id_ghost  array of ghost cell factor
5402   !-------------------------------------------------------------------
5403   SUBROUTINE grid_add_ghost(td_var, id_ghost)
5404      IMPLICIT NONE
5405      ! Argument
5406      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
5407      INTEGER(i4), DIMENSION(2,2), INTENT(IN   ) :: id_ghost
5408
5409      ! local variable
5410      INTEGER(i4) :: il_imin
5411      INTEGER(i4) :: il_jmin
5412      INTEGER(i4) :: il_imax
5413      INTEGER(i4) :: il_jmax
5414
5415      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
5416     
5417      TYPE(TVAR) :: tl_var
5418
5419      ! loop indices
5420      !----------------------------------------------------------------
5421
5422      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
5423
5424         CALL logger_warn( "ADD GHOST: dimension change in variable "//&
5425         &              TRIM(td_var%c_name) )
5426
5427         ! copy variable
5428         tl_var=var_copy(td_var)
5429
5430         CALL var_del_value(td_var)
5431
5432         ! compute indice to fill center
5433         il_imin=1+id_ghost(jp_I,1)*ip_ghost
5434         il_jmin=1+id_ghost(jp_J,1)*ip_ghost
5435
5436         il_imax=tl_var%t_dim(1)%i_len+id_ghost(jp_I,1)*ip_ghost
5437         il_jmax=tl_var%t_dim(2)%i_len+id_ghost(jp_J,1)*ip_ghost
5438
5439         ! compute new dimension
5440         td_var%t_dim(1)%i_len = tl_var%t_dim(1)%i_len + &
5441         &                             SUM(id_ghost(jp_I,:))*ip_ghost
5442         td_var%t_dim(2)%i_len = tl_var%t_dim(2)%i_len + &
5443         &                             SUM(id_ghost(jp_J,:))*ip_ghost
5444
5445         ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
5446         &                 td_var%t_dim(2)%i_len, &
5447         &                 td_var%t_dim(3)%i_len, &
5448         &                 td_var%t_dim(4)%i_len) )
5449
5450         dl_value(:,:,:,:)=tl_var%d_fill
5451
5452         dl_value(il_imin:il_imax, &
5453         &        il_jmin:il_jmax, &
5454         &              :,:)  =  tl_var%d_value(:,:,:,:)
5455
5456         ! add variable value
5457         CALL var_add_value(td_var,dl_value(:,:,:,:))
5458
5459         ! save variable type
5460         td_var%i_type=tl_var%i_type
5461         
5462         DEALLOCATE( dl_value )
5463
5464         CALL var_clean(tl_var)
5465
5466      ENDIF
5467
5468   END SUBROUTINE grid_add_ghost
5469   !-------------------------------------------------------------------
5470   !> @brief
5471   !> This subroutine delete ghost cell at boundaries.
5472   !>
5473   !> @author J.Paul
5474   !> @date November, 2013 - Initial version
5475   !
5476   !> @param[inout] td_var array of variable structure
5477   !> @param[in] id_ghost  array of ghost cell factor
5478   !-------------------------------------------------------------------
5479   SUBROUTINE grid_del_ghost(td_var, id_ghost)
5480      IMPLICIT NONE
5481      ! Argument
5482      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
5483      INTEGER(i4), DIMENSION(2,2), INTENT(IN   ) :: id_ghost
5484
5485      ! local variable
5486      INTEGER(i4) :: il_imin
5487      INTEGER(i4) :: il_jmin
5488      INTEGER(i4) :: il_imax
5489      INTEGER(i4) :: il_jmax
5490
5491      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
5492     
5493      TYPE(TVAR) :: tl_var
5494
5495      ! loop indices
5496      !----------------------------------------------------------------
5497
5498      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
5499
5500         IF( ANY(id_ghost(:,:)/=0) )THEN
5501            CALL logger_warn( "GRID DEL GHOST: dimension change in variable "//&
5502            &              TRIM(td_var%c_name) )
5503         ENDIF
5504
5505         ! copy variable
5506         tl_var=var_copy(td_var)
5507
5508         CALL var_del_value(td_var)
5509
5510         ! compute indice to get center
5511         il_imin=1+id_ghost(jp_I,1)*ip_ghost
5512         il_jmin=1+id_ghost(jp_J,1)*ip_ghost
5513
5514         il_imax=tl_var%t_dim(1)%i_len-id_ghost(jp_I,2)*ip_ghost
5515         il_jmax=tl_var%t_dim(2)%i_len-id_ghost(jp_J,2)*ip_ghost
5516
5517         ! compute new dimension
5518         td_var%t_dim(1)%i_len = il_imax - il_imin +1 
5519         td_var%t_dim(2)%i_len = il_jmax - il_jmin +1 
5520
5521         ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
5522         &                 td_var%t_dim(2)%i_len, &
5523         &                 td_var%t_dim(3)%i_len, &
5524         &                 td_var%t_dim(4)%i_len) )
5525
5526         dl_value(:,:,:,:)=tl_var%d_fill
5527
5528         dl_value(:,:,:,:)  =  tl_var%d_value(il_imin:il_imax, &
5529         &                                    il_jmin:il_jmax, &
5530         &                                    :,:)
5531
5532         ! add variable value
5533         CALL var_add_value(td_var,dl_value(:,:,:,:))
5534
5535         ! save variable type
5536         td_var%i_type=tl_var%i_type
5537         
5538         DEALLOCATE( dl_value )
5539
5540         CALL var_clean(tl_var)
5541
5542      ENDIF
5543
5544   END SUBROUTINE grid_del_ghost
5545   !-------------------------------------------------------------------
5546   !> @brief This function check if ghost cell are used or not, and return ghost
5547   !> cell factor (0,1) in horizontal plan.
5548   !
5549   !> @details
5550   !> check if domain is global, and if there is an East-West overlap.
5551   !>
5552   !> @author J.Paul
5553   !> @date September, 2014 - Initial Version
5554   !
5555   !> @param[in] td_var variable sturcture
5556   !> @return array of ghost cell factor
5557   !-------------------------------------------------------------------
5558   FUNCTION grid__get_ghost_var( td_var )
5559      IMPLICIT NONE
5560      ! Argument
5561      TYPE(TVAR), INTENT(IN) :: td_var
5562
5563      ! function
5564      INTEGER(i4), DIMENSION(2,2) :: grid__get_ghost_var
5565
5566      ! local variable
5567      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim
5568
5569      ! loop indices
5570      !----------------------------------------------------------------
5571      ! init
5572      grid__get_ghost_var(:,:)=0
5573
5574      IF( .NOT. ALL(td_var%t_dim(1:2)%l_use) )THEN
5575         CALL logger_error("GRID GET GHOST: "//TRIM(td_var%c_name)//" is not a suitable"//&
5576         &                 " variable to look for ghost cell (not 2D).")
5577      ELSE
5578         IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
5579            CALL logger_error("GRID GET GHOST: no value associated to "//TRIM(td_var%c_name)//&
5580            &                 ". can't look for ghost cell.")
5581         ELSE
5582            il_dim(:)=td_var%t_dim(:)%i_len
5583
5584            IF(ALL(td_var%d_value(    1    ,    :    ,1,1)/=td_var%d_fill).AND.&
5585            &  ALL(td_var%d_value(il_dim(1),    :    ,1,1)/=td_var%d_fill).AND.&
5586            &  ALL(td_var%d_value(    :    ,    1    ,1,1)/=td_var%d_fill).AND.&
5587            &  ALL(td_var%d_value(    :    ,il_dim(2),1,1)/=td_var%d_fill))THEN
5588            ! no boundary closed
5589               CALL logger_warn("GRID GET GHOST: can't determined ghost cell. "//&
5590               &             "there is no boundary closed for variable "//&
5591               &              TRIM(td_var%c_name))
5592
5593            ELSE
5594               ! check periodicity
5595               IF(ANY(td_var%d_value(   1     ,:,1,1)/=td_var%d_fill).OR.&
5596               &  ANY(td_var%d_value(il_dim(1),:,1,1)/=td_var%d_fill))THEN
5597               ! East-West cyclic (1,4,6)
5598                  CALL logger_info("GRID GET GHOST: East West cyclic")
5599                  grid__get_ghost_var(jp_I,:)=0
5600
5601                  IF( ANY(td_var%d_value(:, 1, 1, 1) /= td_var%d_fill) )THEN
5602                  ! South boundary not closed
5603
5604                     CALL logger_debug("GRID GET GHOST: East_West cyclic")
5605                     CALL logger_debug("GRID GET GHOST: South boundary not closed")
5606                     CALL logger_error("GRID GET GHOST: should have been an "//&
5607                     &              "impossible case")
5608
5609                  ELSE
5610                  ! South boundary closed (1,4,6)
5611                     CALL logger_info("GRID GET GHOST: South boundary closed")
5612                     grid__get_ghost_var(jp_J,1)=1
5613
5614                     IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill) )THEN
5615                     ! North boundary not closed (4,6)
5616                        CALL logger_info("GRID GET GHOST: North boundary not closed")
5617                        grid__get_ghost_var(jp_J,2)=0
5618                     ELSE
5619                     ! North boundary closed
5620                        CALL logger_info("GRID GET GHOST: North boundary closed")
5621                        grid__get_ghost_var(jp_J,2)=1
5622                     ENDIF
5623
5624                  ENDIF
5625
5626               ELSE
5627               ! East-West boundaries closed (0,2,3,5)
5628                  CALL logger_info("GRID GET GHOST: East West boundaries closed")
5629                  grid__get_ghost_var(jp_I,:)=1
5630
5631                  IF( ANY(td_var%d_value(:, 1, 1, 1) /= td_var%d_fill) )THEN
5632                  ! South boundary not closed (2)
5633                     CALL logger_info("GRID GET GHOST: South boundary not closed")
5634                     grid__get_ghost_var(jp_J,1)=0
5635
5636                     IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill))THEN
5637                     ! North boundary not closed
5638                        CALL logger_debug("GRID GET GHOST: East West boundaries "//&
5639                        &              "closed")
5640                        CALL logger_debug("GRID GET GHOST: South boundary not closed")
5641                        CALL logger_debug("GRID GET GHOST: North boundary not closed")
5642                        CALL logger_error("GRID GET GHOST: should have been "//&
5643                        &              "an impossible case")
5644                     ELSE
5645                     ! North boundary closed
5646                        grid__get_ghost_var(jp_J,2)=1
5647                     ENDIF
5648
5649                  ELSE
5650                  ! South boundary closed (0,3,5)
5651                     CALL logger_info("GRID GET GHOST: South boundary closed")
5652                     grid__get_ghost_var(jp_J,1)=1
5653
5654                     IF(ANY(td_var%d_value(:,il_dim(2),1,1)/=td_var%d_fill))THEN
5655                     ! North boundary not closed (3,5)
5656                        CALL logger_info("GRID GET GHOST: North boundary not closed")
5657                        grid__get_ghost_var(jp_J,2)=0
5658                     ELSE
5659                     ! North boundary closed   
5660                        CALL logger_info("GRID GET GHOST: North boundary closed")
5661                        grid__get_ghost_var(jp_J,2)=1
5662                     ENDIF
5663
5664                  ENDIF
5665
5666               ENDIF
5667
5668            ENDIF
5669
5670         ENDIF
5671      ENDIF
5672
5673   END FUNCTION grid__get_ghost_var
5674   !-------------------------------------------------------------------
5675   !> @brief This function check if ghost cell are used or not, and return ghost
5676   !> cell factor (0,1) in i- and j-direction.
5677   !
5678   !> @details
5679   !> get longitude an latitude array, then
5680   !> check if domain is global, and if there is an East-West overlap
5681   !>
5682   !> @author J.Paul
5683   !> @date September, 2014 - Initial Version
5684   !> @date October, 2014
5685   !> - work on mpp file structure instead of file structure
5686   !
5687   !> @param[in] td_file   file sturcture
5688   !> @return array of ghost cell factor
5689   !-------------------------------------------------------------------
5690   FUNCTION grid__get_ghost_mpp( td_mpp )
5691      IMPLICIT NONE
5692      ! Argument
5693      TYPE(TMPP), INTENT(IN) :: td_mpp
5694
5695      ! function
5696      INTEGER(i4), DIMENSION(2,2) :: grid__get_ghost_mpp
5697
5698      ! local variable
5699      !TYPE(TVAR)  :: tl_lon
5700      !TYPE(TVAR)  :: tl_lat
5701
5702      TYPE(TMPP) :: tl_mpp
5703
5704      !INTEGER(i4) :: il_lonid
5705      !INTEGER(i4) :: il_latid
5706      ! loop indices
5707      !----------------------------------------------------------------
5708      ! init
5709      grid__get_ghost_mpp(:,:)=0
5710
5711      IF( .NOT. ASSOCIATED(td_mpp%t_proc) )THEN
5712         CALL logger_error("GRID GET GHOST: decomposition of mpp file "//&
5713         &                 TRIM(td_mpp%c_name)//" not defined." )
5714
5715      ELSE
5716
5717         ! copy structure
5718         tl_mpp=mpp_copy(td_mpp)
5719
5720         CALL logger_info("GRID GET FINE GHOST perio"//TRIM(fct_str(tl_mpp%i_perio)))
5721         IF( tl_mpp%i_perio < 0 )THEN
5722            ! compute NEMO periodicity index
5723            CALL grid_get_info(tl_mpp)
5724         ENDIF
5725
5726         SELECT CASE(tl_mpp%i_perio)
5727            CASE(0)
5728               grid__get_ghost_mpp(:,:)=1
5729            CASE(1)
5730               grid__get_ghost_mpp(jp_J,:)=1
5731            CASE(2)
5732               grid__get_ghost_mpp(jp_I,:)=1
5733               grid__get_ghost_mpp(jp_J,2)=1
5734            CASE(3,5)
5735               grid__get_ghost_mpp(jp_I,:)=1
5736               grid__get_ghost_mpp(jp_J,1)=1
5737            CASE(4,6)
5738               grid__get_ghost_mpp(jp_J,1)=1
5739            CASE DEFAULT
5740         END SELECT
5741
5742         ! clean
5743         CALL mpp_clean(tl_mpp)
5744
5745      ENDIF
5746
5747   END FUNCTION grid__get_ghost_mpp
5748   !-------------------------------------------------------------------
5749   !> @brief This subroutine compute closed sea domain.
5750   !
5751   !> @details
5752   !> to each domain is associated a negative value id (from -1 to ...)<br/>
5753   !> optionaly you could specify which level use (default 1)
5754   !>
5755   !> @author J.Paul
5756   !> @date November, 2013 - Initial Version
5757   !
5758   !> @param[in] td_var    variable strucutre
5759   !> @param[in] id_level  level
5760   !> @return domain mask 
5761   !-------------------------------------------------------------------
5762   FUNCTION grid_split_domain(td_var, id_level)
5763      IMPLICIT NONE
5764      ! Argument
5765      TYPE(TVAR) , INTENT(IN) :: td_var
5766      INTEGER(i4), INTENT(IN), OPTIONAL :: id_level
5767
5768      ! function
5769      INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &
5770      &                      td_var%t_dim(2)%i_len ) :: grid_split_domain
5771
5772      ! local variable
5773      INTEGER(i4)                              :: il_domid
5774      INTEGER(i4)                              :: il_level
5775      INTEGER(i4), DIMENSION(2)                :: il_shape
5776      INTEGER(i4), DIMENSION(2)                :: il_ind
5777      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask
5778      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp
5779
5780      LOGICAL                                  :: ll_full
5781
5782      ! loop indices
5783      INTEGER(i4) :: ji
5784      INTEGER(i4) :: jim
5785      INTEGER(i4) :: jip
5786      INTEGER(i4) :: jj
5787      INTEGER(i4) :: jjm
5788      INTEGER(i4) :: jjp
5789      !----------------------------------------------------------------
5790      il_level=1
5791      IF( PRESENT(id_level) ) il_level=id_level
5792
5793      ! init
5794      il_domid=-1
5795
5796      il_shape(:)=td_var%t_dim(1:2)%i_len
5797      ALLOCATE( il_mask(il_shape(1),il_shape(2)) )
5798      il_mask(:,:)=0
5799      WHERE( td_var%d_value(:,:,il_level,1)/=td_var%d_fill ) il_mask(:,:)=1
5800
5801      il_ind(:)=MAXLOC( il_mask(:,:) )
5802      DO WHILE( il_mask(il_ind(1), il_ind(2)) == 1 )
5803
5804         il_mask(il_ind(1),il_ind(2))=il_domid
5805         ll_full=.FALSE.
5806
5807         ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
5808
5809         DO WHILE( .NOT. ll_full )
5810            il_tmp(:,:)=0
5811
5812            ll_full=.TRUE.
5813            DO jj=1,il_shape(2)
5814               DO ji=1,il_shape(1)
5815                  IF( il_mask(ji,jj)==il_domid )THEN
5816                     jim=MAX(1,ji-1)   ;  jip=MIN(il_shape(1),ji+1)
5817                     jjm=MAX(1,jj-1)   ;  jjp=MIN(il_shape(2),jj+1)
5818                     
5819                     WHERE( il_mask(jim:jip,jjm:jjp)==1 )
5820                        il_mask(jim:jip,jjm:jjp)=il_domid
5821                        il_tmp(jim:jip,jjm:jjp)=1
5822                     END WHERE
5823
5824                  ENDIF
5825               ENDDO
5826            ENDDO
5827            IF( SUM(il_tmp(:,:))/=0 ) ll_full=.FALSE.
5828
5829         ENDDO
5830         DEALLOCATE( il_tmp )
5831
5832         il_ind(:)=MAXLOC( il_mask(:,:) )
5833         il_domid=il_domid-1
5834
5835      ENDDO
5836
5837      ! save result
5838      grid_split_domain(:,:)=il_mask(:,:)
5839
5840      DEALLOCATE( il_mask )
5841
5842      CALL logger_info("GRID SPLIT DOMAIN: "//TRIM( fct_str(ABS(il_domid+1)) )//&
5843      &             " domain found" ) 
5844
5845   END FUNCTION grid_split_domain
5846   !-------------------------------------------------------------------
5847   !> @brief This subroutine fill small closed sea with fill value.
5848   !>
5849   !> @details
5850   !> the minimum size (number of point) of closed sea to be kept could be
5851   !> sepcify with id_minsize.
5852   !> By default only the biggest sea is preserve.
5853   !>
5854   !> @author J.Paul
5855   !> @date November, 2013 - Initial Version
5856   !>
5857   !> @param[inout] td_var    variable structure
5858   !> @param[in] id_mask      domain mask (from grid_split_domain)
5859   !> @param[in] id_minsize   minimum size of sea to be kept
5860   !-------------------------------------------------------------------
5861   SUBROUTINE grid_fill_small_dom(td_var, id_mask, id_minsize)
5862      IMPLICIT NONE
5863      ! Argument     
5864      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
5865      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_mask
5866      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_minsize
5867
5868      ! local variable
5869      INTEGER(i4)                              :: il_ndom
5870      INTEGER(i4)                              :: il_minsize
5871      INTEGER(i4), DIMENSION(2)                :: il_shape
5872      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp
5873
5874      ! loop indices
5875      INTEGER(i4) :: ji
5876      INTEGER(i4) :: jk
5877      INTEGER(i4) :: jl
5878      !----------------------------------------------------------------
5879
5880      il_shape(:)=SHAPE(id_mask(:,:))
5881      IF( ANY(il_shape(:) /= td_var%t_dim(1:2)%i_len) )THEN
5882         CALL logger_error("GRID FILL SMALL DOM: variable and mask "//&
5883         &              "dimension differ")
5884      ELSE
5885
5886         il_ndom=MINVAL(id_mask(:,:))
5887
5888         ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
5889         il_tmp(:,:)=0
5890         DO ji=-1,il_ndom,-1
5891            WHERE( id_mask(:,:)==ji ) 
5892               il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji
5893            END WHERE
5894         ENDDO
5895
5896         il_minsize=MAXVAL(il_tmp(:,:))
5897         IF( PRESENT(id_minsize) ) il_minsize=id_minsize
5898
5899         DO jl=1,td_var%t_dim(4)%i_len
5900            DO jk=1,td_var%t_dim(3)%i_len
5901               WHERE( il_tmp(:,:) < il_minsize ) 
5902                  td_var%d_value(:,:,jk,jl)=td_var%d_fill
5903               END WHERE
5904            ENDDO
5905         ENDDO
5906
5907         DEALLOCATE( il_tmp )
5908
5909      ENDIF
5910
5911   END SUBROUTINE grid_fill_small_dom
5912   !-------------------------------------------------------------------
5913   !> @brief This subroutine fill small domain inside bigger one.
5914   !>
5915   !> @details
5916   !> the minimum size (number of point) of domain sea to be kept could be
5917   !> is sepcified with id_minsize.
5918   !> smaller domain are included in the one they are embedded.
5919   !>
5920   !> @author J.Paul
5921   !> @date Ferbruay, 2015 - Initial Version
5922   !>
5923   !> @param[inout] id_mask      domain mask (from grid_split_domain)
5924   !> @param[in] id_minsize   minimum size of sea to be kept
5925   !-------------------------------------------------------------------
5926   SUBROUTINE grid_fill_small_msk(id_mask, id_minsize)
5927      IMPLICIT NONE
5928      ! Argument     
5929      INTEGER(i4), DIMENSION(:,:), INTENT(INOUT) :: id_mask
5930      INTEGER(i4),                 INTENT(IN   ) :: id_minsize
5931
5932      ! local variable
5933      INTEGER(i4)                              :: il_ndom
5934      INTEGER(i4)                              :: il_minsize
5935      INTEGER(i4)                              :: il_msk
5936     
5937      INTEGER(i4)                              :: jim
5938      INTEGER(i4)                              :: jjm
5939      INTEGER(i4)                              :: jip
5940      INTEGER(i4)                              :: jjp
5941
5942      INTEGER(i4), DIMENSION(2)                :: il_shape
5943      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp
5944
5945      ! loop indices
5946      INTEGER(i4) :: ii
5947      INTEGER(i4) :: ij
5948
5949      INTEGER(i4) :: ji
5950      INTEGER(i4) :: jj
5951      !----------------------------------------------------------------
5952
5953      il_shape(:)=SHAPE(id_mask(:,:))
5954      il_ndom=MINVAL(id_mask(:,:))
5955
5956      ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
5957      il_tmp(:,:)=0
5958      DO ji=-1,il_ndom,-1
5959         WHERE( id_mask(:,:)==ji ) 
5960            il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji
5961         END WHERE
5962      ENDDO
5963
5964      DO WHILE( id_minsize > MINVAL(il_tmp(:,:)) )
5965
5966         DO jj=1,il_shape(2)
5967            DO ji=1,il_shape(1)
5968
5969               IF( il_tmp(ji,jj) < il_minsize )THEN
5970                  jim=MAX(1,ji-1)   ;  jip=MIN(il_shape(1),ji+1)
5971                  jjm=MAX(1,jj-1)   ;  jjp=MIN(il_shape(2),jj+1)
5972                 
5973                  il_msk=0
5974                  DO ij=jjm,jjp
5975                     DO ii=jim,jip
5976                        IF( id_mask(ii,ij) /= id_mask(ji,jj) )THEN
5977                           IF( il_msk == 0 )THEN
5978                              il_msk=id_mask(ii,ij)
5979                           ELSEIF( il_msk /= id_mask(ii,ij) )THEN
5980                              CALL logger_error("GRID FILL SMALL MSK: "//&
5981                              &  "small domain not embedded in bigger one"//&
5982                              &  ". point should be between two different"//&
5983                              &  " domain.")
5984                           ENDIF
5985                        ENDIF
5986                     ENDDO
5987                  ENDDO
5988                  IF( il_msk /= 0 ) id_mask(ji,jj)=il_msk
5989
5990               ENDIF
5991
5992            ENDDO
5993         ENDDO
5994
5995
5996         il_tmp(:,:)=0
5997         DO ji=-1,il_ndom,-1
5998            WHERE( id_mask(:,:)==ji ) 
5999               il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji
6000            END WHERE
6001         ENDDO           
6002
6003      ENDDO
6004
6005      DEALLOCATE( il_tmp )
6006
6007
6008   END SUBROUTINE grid_fill_small_msk
6009END MODULE grid
6010
Note: See TracBrowser for help on using the repository browser.