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/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/TOOLS/SIREN/src/grid.f90 @ 5972

Last change on this file since 5972 was 5972, checked in by timgraham, 8 years ago

Upgraded to head of trunk (r5936)

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