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/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/TOOLS/SIREN/src/grid.f90 @ 5240

Last change on this file since 5240 was 5240, checked in by davestorkey, 9 years ago

Update UKMO nn_etau_revision branch with trunk changes to rev 5107.

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