source: utils/tools/SIREN/src/merge_bathy.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months ago

update nemo trunk

File size: 42.7 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! DESCRIPTION:
6!> @file
7!> @brief
8!> this program merges bathymetry file at boundaries.
9!>
10!> @details
11!> @section sec1 method
12!> coarse grid bathymetry is interpolated on fine grid
13!> (nearest interpolation method is used).<br/> 
14!> then fine bathymetry and refined coarse bathymetry are merged at boundaries.<br/>
15!>    @f[BathyFine= Weight * BathyCoarse + (1-Weight)*BathyFine@f]
16!> the weight function used is :<br/>
17!>       @f[Weight = 0.5 + 0.5*COS( \frac{\pi*dist}{width} )@f]<br/>
18!> with
19!> - dist : number of point to border
20!> - width : boundary size
21!>
22!> @section sec2 how to
23!> USAGE: merge_bathy merge_bathy.nam [-v] [-h]<br/>
24!>    - positional arguments:<br/>
25!>       - merge_bathy.nam<br/>
26!>          namelist of merge_bathy
27!>          @note
28!>             a template of the namelist could be created running (in templates directory):
29!>             @code{.sh}
30!>                python create_templates.py merge_bathy
31!>             @endcode
32!>
33!>    - optional arguments:<br/>
34!>       - -h, --help<br/>
35!>          show this help message (and exit)<br/>
36!>       - -v, --version<br/>
37!>          show Siren's version   (and exit)
38!>
39!> @section sec_merge merge_bathy.nam
40!>    merge_bathy.nam contains 7 namelists:
41!>       - **namlog** to set logger parameters
42!>       - **namcfg** to set configuration file parameters
43!>       - **namsrc** to set source/coarse grid parameters
44!>       - **namtgt** to set target/fine grid parameters
45!>       - **namnst** to set sub domain and nesting paramters
46!>       - **nambdy** to set boundary parameters
47!>       - **namout** to set output parameters
48!>
49!>    here after, each sub-namelist parameters is detailed.
50!>    @note
51!>       default values are specified between brackets
52!>
53!> @subsection sublog namlog
54!>    the logger sub-namelist parameters are :
55!>
56!>    - **cn_logfile** [@a merge_bathy.log]<br/>
57!>       logger filename
58!>
59!>    - **cn_verbosity** [@a warning]<br/>
60!>       verbosity level, choose between :
61!>          - trace
62!>          - debug
63!>          - info
64!>          - warning
65!>          - error
66!>          - fatal
67!>          - none
68!>
69!>    - **in_maxerror** [@a 5]<br/>
70!>       maximum number of error allowed
71!>
72!> @subsection subcfg namcfg
73!>    the configuration sub-namelist parameters are :
74!>
75!>    - **cn_varcfg** [@a ./cfg/variable.cfg]<br/>
76!>       path to the variable configuration file.<br/>
77!>       the variable configuration file defines standard name,
78!>       default interpolation method, axis,...
79!>       to be used for some known variables.<br/>
80!>
81!>    - **cn_dimcfg** [@a ./cfg/dimension.cfg]<br/>
82!>       path to the dimension configuration file.<br/>
83!>       the dimension configuration file defines dimensions allowed.<br/>
84!>
85!>    - **cn_dumcfg** [@a ./cfg/dummy.cfg]<br/>
86!>       path to the useless (dummy) configuration file.<br/>
87!>       the dummy configuration file defines useless
88!>       dimension or variable. these dimension(s) or variable(s) will not be
89!>       processed.<br/>
90!>
91!> @subsection subsrc namsrc
92!>    the source/coarse grid sub-namelist parameters are :
93!>
94!>    - **cn_bathy0** [@a ]<br/>
95!>       path to the bathymetry file
96!>       @warning
97!>          variable name must be __Bathymetry__ here.
98!>
99!>    - **in_perio0** [@a ]<br/>
100!>       NEMO periodicity index<br/>
101!>       the NEMO periodicity could be choose between 0 to 6:
102!>       <dl>
103!>          <dt>in_perio=0</dt>
104!>          <dd>standard regional model</dd>
105!>          <dt>in_perio=1</dt>
106!>          <dd>east-west cyclic model</dd>
107!>          <dt>in_perio=2</dt>
108!>          <dd>model with symmetric boundary condition across the equator</dd>
109!>          <dt>in_perio=3</dt>
110!>          <dd>regional model with North fold boundary and T-point pivot</dd>
111!>          <dt>in_perio=4</dt>
112!>          <dd>global model with a T-point pivot.<br/>
113!>          example: ORCA2, ORCA025, ORCA12</dd>
114!>          <dt>in_perio=5</dt>
115!>          <dd>regional model with North fold boundary and F-point pivot</dd>
116!>          <dt>in_perio=6</dt>
117!>          <dd>global model with a F-point pivot<br/>
118!>          example: ORCA05</dd>
119!>          </dd>
120!>       </dl>
121!>       @sa For more information see @ref md_src_docsrc_6_perio
122!>       and Model Boundary Condition paragraph in the
123!>       [NEMO documentation](https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/manual/pdf/NEMO_manual.pdf)
124!>
125!> @subsection subtgt namtgt
126!>    the target/fine grid sub-namelist parameters are :
127!>
128!>    - **cn_bathy1** [@a ]<br/>
129!>       path to bathymetry file
130!>       @warning
131!>          variable name must be __Bathymetry__ here.
132!>
133!>    - **in_perio1** [@a ]<br/>
134!>       NEMO periodicity index (see above)
135!>    @note if the fine/target coordinates file (cn_coord1) was created by SIREN, you do
136!>    not need to fill this parameter. SIREN will read it on the global attributes of
137!>    the coordinates file.
138!>
139!> @subsection subnst namnst
140!>    the nesting sub-namelist parameters are (default value are specified between brackets):
141!>    - **in_rhoi**  [@a 1]<br/>
142!>       refinement factor in i-direction
143!>
144!>    - **in_rhoj**  [@a 1]<br/>
145!>       refinement factor in j-direction
146!>
147!>    @note
148!>       coarse grid indices will be deduced from fine grid
149!>       coordinate file.
150!>
151!> @subsection subbdy nambdy
152!>    the boundary sub-namelist parameters are :
153!>
154!>    - **ln_north** [@a .TRUE.]<br/>
155!>       logical to use north boundary or not
156!>    - **ln_south** [@a .TRUE.]<br/>
157!>       logical to use south boundary or not
158!>    - **ln_east**  [@a .TRUE.]<br/>
159!>       logical to use east boundary or not
160!>    - **ln_west**  [@a .TRUE.]<br/>
161!>       logical to use west  boundary or not
162!>    <br/> <br/>
163!>    - **cn_north** [@a ]<br/>
164!>       north boundary indices on fine grid<br/>
165!>    - **cn_south** [@a ]<br/>
166!>       south boundary indices on fine grid<br/>
167!>    - **cn_east**  [@a ]<br/>
168!>       east  boundary indices on fine grid<br/>
169!>    - **cn_west**  [@a ]<br/>
170!>       west  boundary indices on fine grid<br/>
171!>
172!>       *cn_north* is a string character defining boundary
173!>       segmentation.<br/>
174!>       segments are separated by '|'.<br/>
175!>       each segments of the boundary is composed of:
176!>          - indice of velocity (orthogonal to boundary .ie.
177!>             for north boundary, J-indice).
178!>          - indice of segment start (I-indice for north boundary)
179!>          - indice of segment end   (I-indice for north boundary)<br/>
180!>             indices must be separated by ':' .<br/>
181!>          - optionally, boundary size could be added between '(' and ')'
182!>          in the first segment defined.
183!>             @note
184!>                boundary size is the same for all segments of one boundary.
185!>
186!>       Examples:
187!>          - cn_north='index1,first1:last1(width)'
188!>          - cn_north='index1(width),first1:last1|index2,first2:last2'
189!>
190!>       @image html  boundary_50.png
191!>       <center>@image latex boundary_50.png
192!>       </center>
193!>
194!>    - **in_ncrs**  [@a 2]<br/>
195!>       number of point(s) with coarse value save at boundaries
196!>
197!>    - **ln_oneseg** [@a .TRUE.]<br/>
198!>       logical to use only one segment for each boundary or not
199!>
200!> @subsection subout namout
201!>    the output sub-namelist parameter is :
202!>
203!>    - **cn_fileout** [@a bathy_merged.nc]<br/>
204!>       output bathymetry filename
205!>
206!> <hr>
207!> @author J.Paul
208!>
209!> @date November, 2013 - Initial Version
210!> @date Sepember, 2014
211!> - add header for user
212!> @date July, 2015
213!> - extrapolate all land points
214!> - add attributes with boundary string character (as in namelist)
215!> @date September, 2015
216!> - manage useless (dummy) variable, attributes, and dimension
217!> @date October, 2016
218!> - allow to choose the number of boundary point with coarse grid value.
219!> - dimension to be used select from configuration file
220!> @date January, 2019
221!> - add url path to global attributes of output file(s)
222!> @date February, 2019
223!> - rename sub namelist namcrs to namsrc
224!> - rename sub namelist namfin to namtgt
225!> @date May, 2019
226!> - create and clean file structure to avoid memory leaks
227!> @date Ocober, 2019
228!> - add help and version optional arguments
229!>
230!>
231!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
232!----------------------------------------------------------------------
233PROGRAM merge_bathy
234
235   USE netcdf                          ! nf90 library
236   USE global                          ! global variable
237   USE phycst                          ! physical constant
238   USE kind                            ! F90 kind parameter
239   USE logger                          ! log file manager
240   USE fct                             ! basic useful function
241   USE date                            ! date manager
242   USE att                             ! attribute manager
243   USE dim                             ! dimension manager
244   USE var                             ! variable manager
245   USE file                            ! file manager
246   USE boundary                        ! boundary manager
247   USE iom                             ! I/O manager
248   USE grid                            ! grid manager
249   USE extrap                          ! extrapolation manager
250   USE interp                          ! interpolation manager
251   USE filter                          ! filter manager
252   USE mpp                             ! MPP manager
253   USE dom                             ! domain manager
254   USE iom_mpp                         ! MPP I/O manager
255   USE iom_dom                         ! DOM I/O manager
256
257   IMPLICIT NONE
258
259   ! parameters
260   CHARACTER(LEN=lc), PARAMETER  :: cp_myname = "merge_bathy"
261
262   ! local variable
263   CHARACTER(LEN=lc)                                  :: cl_arg
264   CHARACTER(LEN=lc)                                  :: cl_namelist
265   CHARACTER(LEN=lc)                                  :: cl_date
266   CHARACTER(LEN=lc)                                  :: cl_tmp
267   CHARACTER(LEN=lc)                                  :: cl_url
268   CHARACTER(LEN=lc)                                  :: cl_errormsg
269
270   INTEGER(i4)                                        :: il_narg
271   INTEGER(i4)                                        :: il_status
272   INTEGER(i4)                                        :: il_fileid
273   INTEGER(i4)                                        :: il_attind
274   INTEGER(i4)                                        :: il_imin0
275   INTEGER(i4)                                        :: il_imax0
276   INTEGER(i4)                                        :: il_jmin0
277   INTEGER(i4)                                        :: il_jmax0
278   INTEGER(i4)                                        :: il_shift
279   INTEGER(i4)      , DIMENSION(ip_maxdim)            :: il_rho
280   INTEGER(i4)      , DIMENSION(2,2)                  :: il_ind
281
282   LOGICAL                                            :: ll_exist
283
284   REAL(dp)                                           :: dl_fill
285   REAL(dp)         , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_refined
286   REAL(dp)         , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_weight
287
288   TYPE(TMPP)                                         :: tl_bathy0
289   TYPE(TMPP)                                         :: tl_bathy1
290
291   TYPE(TFILE)                                        :: tl_file
292   TYPE(TFILE)                                        :: tl_fileout
293   
294   TYPE(TATT)                                         :: tl_att
295   
296   TYPE(TVAR)                                         :: tl_var
297   TYPE(TVAR)                                         :: tl_lon
298   TYPE(TVAR)                                         :: tl_lat
299   
300   TYPE(TDIM)       , DIMENSION(ip_maxdim)            :: tl_dim
301   
302   TYPE(TBDY)       , DIMENSION(ip_ncard)             :: tl_bdy
303
304   ! loop indices
305   INTEGER(i4) :: ji
306   INTEGER(i4) :: jk
307   INTEGER(i4) :: jl
308
309   ! namelist variable
310   ! namlog
311   CHARACTER(LEN=lc)                       :: cn_logfile = 'merge_bathy.log' 
312   CHARACTER(LEN=lc)                       :: cn_verbosity = 'warning' 
313   INTEGER(i4)                             :: in_maxerror = 5
314
315   ! namcfg
316   CHARACTER(LEN=lc)                       :: cn_varcfg = './cfg/variable.cfg' 
317   CHARACTER(LEN=lc)                       :: cn_dimcfg = './cfg/dimension.cfg'
318   CHARACTER(LEN=lc)                       :: cn_dumcfg = './cfg/dummy.cfg'
319
320   ! namsrc
321   CHARACTER(LEN=lc)                       :: cn_bathy0 = '' 
322   INTEGER(i4)                             :: in_perio0 = -1
323
324   ! namtgt
325   CHARACTER(LEN=lc)                       :: cn_bathy1 = '' 
326   INTEGER(i4)                             :: in_perio1 = -1
327
328!   ! namvar
329!   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
330
331   ! namnst
332   INTEGER(i4)                             :: in_rhoi  = 1
333   INTEGER(i4)                             :: in_rhoj  = 1
334
335   ! nambdy
336   LOGICAL                                 :: ln_north = .TRUE.
337   LOGICAL                                 :: ln_south = .TRUE.
338   LOGICAL                                 :: ln_east  = .TRUE.
339   LOGICAL                                 :: ln_west  = .TRUE.
340   CHARACTER(LEN=lc)                       :: cn_north = ''
341   CHARACTER(LEN=lc)                       :: cn_south = ''
342   CHARACTER(LEN=lc)                       :: cn_east  = ''
343   CHARACTER(LEN=lc)                       :: cn_west  = ''
344   INTEGER(i4)                             :: in_ncrs  = 2
345   LOGICAL                                 :: ln_oneseg= .TRUE.
346
347   ! namout
348   CHARACTER(LEN=lc)                       :: cn_fileout = 'bathy_merged.nc' 
349   !-------------------------------------------------------------------
350
351   NAMELIST /namlog/ &   !< logger namelist
352   &  cn_logfile,    &   !< log file
353   &  cn_verbosity,  &   !< log verbosity
354   &  in_maxerror        !< logger maximum error
355
356   NAMELIST /namcfg/ &   !< configuration namelist
357   &  cn_varcfg, &       !< variable configuration file
358   &  cn_dimcfg, &       !< dimension configuration file
359   &  cn_dumcfg          !< dummy configuration file
360
361   NAMELIST /namsrc/ &   !< source/coarse grid namelist
362   &  cn_bathy0,  &      !< bathymetry file
363   &  in_perio0          !< periodicity index
364   
365   NAMELIST /namtgt/ &   !< target/fine grid namelist
366   &  cn_bathy1,     &   !< bathymetry file
367   &  in_perio1          !< periodicity index
368 
369!   NAMELIST /namvar/ &  !< variable namelist
370!   &  cn_varinfo        !< list of variable and interpolation
371!                        !< method to be used.
372!                        !< (ex: 'votemper|linear','vosaline|cubic' )
373   
374   NAMELIST /namnst/ &  !< nesting namelist
375   &  in_rhoi,    &     !< refinement factor in i-direction
376   &  in_rhoj           !< refinement factor in j-direction
377
378   NAMELIST /nambdy/ &  !< boundary namelist
379   &  ln_north,   &     !< use north boundary
380   &  ln_south,   &     !< use south boundary
381   &  ln_east ,   &     !< use east  boundary
382   &  ln_west ,   &     !< use west  boundary
383   &  cn_north,   &     !< north boundary indices on fine grid
384   &  cn_south,   &     !< south boundary indices on fine grid
385   &  cn_east ,   &     !< east  boundary indices on fine grid
386   &  cn_west ,   &     !< west  boundary indices on fine grid
387   &  in_ncrs,    &     !< number of point with coarse value save at boundaries
388   &  ln_oneseg         !< use only one segment for each boundary or not
389
390   NAMELIST /namout/ &  !< output namelist
391   &  cn_fileout        !< fine grid merged bathymetry file   
392   !-------------------------------------------------------------------
393
394   !
395   ! Initialisation
396   ! --------------
397   !
398   il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
399
400   ! Traitement des arguments fournis
401   ! --------------------------------
402   IF( il_narg /= 1 )THEN
403      WRITE(cl_errormsg,*) ' ERROR : one argument is needed '
404      CALL fct_help(cp_myname,cl_errormsg) 
405      CALL EXIT(1)
406   ELSE
407
408      CALL GET_COMMAND_ARGUMENT(1,cl_arg) !f03 intrinsec
409      SELECT CASE (cl_arg)
410         CASE ('-v', '--version')
411
412            CALL fct_version(cp_myname)
413            CALL EXIT(0)
414
415         CASE ('-h', '--help')
416
417            CALL fct_help(cp_myname)
418            CALL EXIT(0)
419
420         CASE DEFAULT
421
422            cl_namelist=cl_arg
423
424            ! read namelist
425            INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
426            IF( ll_exist )THEN
427
428               il_fileid=fct_getunit()
429
430               OPEN( il_fileid, FILE=TRIM(cl_namelist),  &
431               &                FORM='FORMATTED',        &
432               &                ACCESS='SEQUENTIAL',     &
433               &                STATUS='OLD',            &
434               &                ACTION='READ',           &
435               &                IOSTAT=il_status)
436               CALL fct_err(il_status)
437               IF( il_status /= 0 )THEN
438                  WRITE(cl_errormsg,*) " ERROR : error opening "//TRIM(cl_namelist)
439                  CALL fct_help(cp_myname,cl_errormsg) 
440                  CALL EXIT(1)
441               ENDIF
442
443               READ( il_fileid, NML = namlog )
444               ! define log file
445               CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
446               CALL logger_header()
447
448               READ( il_fileid, NML = namcfg )
449               ! get variable extra information
450               CALL var_def_extra(TRIM(cn_varcfg))
451
452               ! get dimension allowed
453               CALL dim_def_extra(TRIM(cn_dimcfg))
454
455               ! get dummy variable
456               CALL var_get_dummy(TRIM(cn_dumcfg))
457               ! get dummy dimension
458               CALL dim_get_dummy(TRIM(cn_dumcfg))
459               ! get dummy attribute
460               CALL att_get_dummy(TRIM(cn_dumcfg))
461
462               READ( il_fileid, NML = namsrc )
463               READ( il_fileid, NML = namtgt )
464!               READ( il_fileid, NML = namvar )
465!               ! add user change in extra information
466!               CALL var_chg_extra(cn_varinfo)
467
468               READ( il_fileid, NML = namnst )
469               READ( il_fileid, NML = nambdy )
470
471               READ( il_fileid, NML = namout )
472
473               CLOSE( il_fileid, IOSTAT=il_status )
474               CALL fct_err(il_status)
475               IF( il_status /= 0 )THEN
476                  CALL logger_error("MERGE BATHY: ERROR closing "//TRIM(cl_namelist))
477               ENDIF
478
479            ELSE
480
481               WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cl_namelist)
482               CALL fct_help(cp_myname,cl_errormsg) 
483               CALL EXIT(1)
484
485            ENDIF
486
487      END SELECT
488   ENDIF
489
490   ! open files
491   IF( TRIM(cn_bathy0) /= '' )THEN
492      tl_file=file_init(TRIM(cn_bathy0))
493      tl_bathy0=mpp_init( tl_file, id_perio=in_perio0)
494      ! clean
495      CALL file_clean(tl_file)
496      CALL grid_get_info(tl_bathy0)
497   ELSE
498      CALL logger_fatal("MERGE BATHY: can not find coarse grid bathymetry "//&
499      &  "file. check namelist")
500   ENDIF
501
502   IF( TRIM(cn_bathy1) /= '' )THEN
503      tl_file=file_init(TRIM(cn_bathy1))
504      tl_bathy1=mpp_init( tl_file, id_perio=in_perio1)
505      ! clean
506      CALL file_clean(tl_file)
507      CALL grid_get_info(tl_bathy1)
508   ELSE
509      CALL logger_fatal("MERGE BATHY: can not find fine grid bathymetry "//&
510      &  "file. check namelist")
511   ENDIF
512
513   ! check
514   ! check output file do not already exist
515   INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
516   IF( ll_exist )THEN
517      CALL logger_fatal("CREATE BATHY: output file "//TRIM(cn_fileout)//&
518      &  " already exist.")
519   ENDIF
520
521   ! check namelist
522   ! check refinament factor
523   il_rho(:)=1
524   IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
525      CALL logger_error("MERGE BATHY: invalid refinement factor."//&
526      &  " check namelist "//TRIM(cl_namelist))
527   ELSE
528      il_rho(jp_I)=in_rhoi
529      il_rho(jp_J)=in_rhoj
530   ENDIF
531
532   ! check domain indices
533   ! compute coarse grid indices around fine grid
534   il_ind(:,:)=grid_get_coarse_index(tl_bathy0, tl_bathy1, &
535   &                                 id_rho=il_rho(:) )
536
537   il_imin0=il_ind(1,1) ; il_imax0=il_ind(1,2)
538   il_jmin0=il_ind(2,1) ; il_jmax0=il_ind(2,2)
539
540   ! check domain validity
541   CALL grid_check_dom(tl_bathy0, il_imin0, il_imax0, il_jmin0, il_jmax0)
542
543   ! check coincidence between coarse and fine grid
544   CALL grid_check_coincidence( tl_bathy0, tl_bathy1, &
545   &                            il_imin0, il_imax0, &
546   &                            il_jmin0, il_jmax0, &
547   &                            il_rho(:) )
548
549   ! open mpp files
550   CALL iom_mpp_open(tl_bathy1)
551
552   ! read or compute boundary
553   tl_var=iom_mpp_read_var(tl_bathy1,'Bathymetry')
554
555   ! close mpp files
556   CALL iom_mpp_close(tl_bathy1)
557
558   tl_bdy(:)=boundary_init(tl_var, ln_north, ln_south, ln_east, ln_west, &
559   &                               cn_north, cn_south, cn_east, cn_west, &
560   &                               ln_oneseg ) 
561
562   ! get boundary on coarse grid
563   ! define refined bathymetry array (for coarse grid)
564   dl_fill=tl_var%d_fill
565   ALLOCATE( dl_refined(tl_var%t_dim(1)%i_len, &
566   &                    tl_var%t_dim(2)%i_len, &
567   &                    tl_var%t_dim(3)%i_len, &
568   &                    tl_var%t_dim(4)%i_len) )
569
570   dl_refined(:,:,:,:)=dl_fill
571
572   ! define weight array
573   ALLOCATE( dl_weight(tl_var%t_dim(1)%i_len, &
574   &                   tl_var%t_dim(2)%i_len, &
575   &                   1,1) )
576
577   dl_weight(:,:,:,:)=dl_fill 
578
579   ! compute coarse grid refined bathymetry on boundary.
580   DO jk=1,ip_ncard
581      CALL merge_bathy_get_boundary(tl_bathy0, tl_bathy1, tl_bdy(jk), &
582      &                             il_rho(:), in_ncrs,                      &
583      &                             dl_refined(:,:,:,:), dl_weight(:,:,:,:), &
584      &                             dl_fill)
585
586   ENDDO
587
588   ! merge bathy on boundary
589   DO jl=1,tl_var%t_dim(4)%i_len
590      DO jk=1,tl_var%t_dim(3)%i_len
591         WHERE(    dl_refined(:,:,jk,jl) /= dl_fill .AND. &
592         &      tl_var%d_value(:,:,jk,jl)/= tl_var%d_fill )
593            tl_var%d_value(:,:,jk,jl)=  &
594            &           dl_weight(:,:,1,1) *    dl_refined(:,:,jk,jl) + &
595            &        (1-dl_weight(:,:,1,1))*tl_var%d_value(:,:,jk,jl)
596         ELSE WHERE( dl_refined(:,:,jk,jl)== dl_fill .AND. &
597         &           dl_weight(:,:,1,1)  /= dl_fill )
598         ! to keep coarse grid mask on boundary
599            tl_var%d_value(:,:,jk,jl)=dl_fill
600         END WHERE
601      ENDDO
602   ENDDO
603
604   DEALLOCATE(dl_refined)
605
606   ! create file
607   tl_fileout=file_init(TRIM(cn_fileout),id_perio=in_perio1)
608
609   ! add dimension
610   tl_dim(:)=dim_copy(tl_var%t_dim(:))
611
612   DO ji=1,ip_maxdim
613      IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
614   ENDDO
615
616   ! add variables
617   IF( ALL( tl_dim(1:2)%l_use ) )THEN
618      ! open mpp files
619      CALL iom_mpp_open(tl_bathy1)
620
621      ! add longitude
622      tl_lon=iom_mpp_read_var(tl_bathy1,'longitude')
623      CALL file_add_var(tl_fileout, tl_lon)
624      CALL var_clean(tl_lon)
625
626      ! add latitude
627      tl_lat=iom_mpp_read_var(tl_bathy1,'latitude')
628      CALL file_add_var(tl_fileout, tl_lat)
629      CALL var_clean(tl_lat)
630
631      ! close mpp files
632      CALL iom_mpp_close(tl_bathy1)     
633   ENDIF
634
635   CALL file_add_var(tl_fileout, tl_var)
636   CALL var_clean(tl_var)
637
638   ! only 2 first dimension to be used
639   tl_dim(:)=dim_copy(tl_fileout%t_dim(:))
640   tl_dim(3:4)%l_use=.FALSE.
641   tl_var=var_init('weight',dl_weight(:,:,:,:),td_dim=tl_dim(:),dd_fill=dl_fill)
642   CALL file_add_var(tl_fileout, tl_var)
643   CALL var_clean(tl_var)
644
645   ! add some attribute
646   tl_att=att_init("Created_by","SIREN merge_bathy")
647   CALL file_add_att(tl_fileout, tl_att)
648
649   !add source url
650   cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
651   tl_att=att_init("SIREN_url",cl_url)
652   CALL file_add_att(tl_fileout, tl_att)
653
654   ! add date of creation
655   cl_date=date_print(date_now())
656   tl_att=att_init("Creation_date",cl_date)
657   CALL file_add_att(tl_fileout, tl_att)
658
659   tl_att=att_init("coarse_grid_source_file",TRIM(fct_basename(tl_bathy0%c_name)))
660   CALL file_add_att(tl_fileout, tl_att)
661
662   tl_att=att_init("fine_grid_source_file",TRIM(fct_basename(tl_bathy1%c_name)))
663   CALL file_add_att(tl_fileout, tl_att)
664
665   ! add attribute periodicity
666   il_attind=0
667   IF( ASSOCIATED(tl_fileout%t_att) )THEN
668      il_attind=att_get_index(tl_fileout%t_att(:),'periodicity')
669   ENDIF
670   IF( tl_bathy1%i_perio >= 0 .AND. il_attind == 0 )THEN
671      tl_att=att_init('periodicity',tl_bathy1%i_perio)
672      CALL file_add_att(tl_fileout,tl_att)
673   ENDIF
674
675   il_attind=0
676   IF( ASSOCIATED(tl_fileout%t_att) )THEN
677      il_attind=att_get_index(tl_fileout%t_att(:),'ew_overlap')
678   ENDIF
679   IF( tl_bathy1%i_ew >= 0 .AND. il_attind == 0 )THEN
680      tl_att=att_init('ew_overlap',tl_bathy1%i_ew)
681      CALL file_add_att(tl_fileout,tl_att)
682   ENDIF
683
684
685   IF( tl_bdy(jp_north)%l_use )THEN
686      ! add shift on north boundary
687      ! boundary compute on T point but express on U or V point
688      il_shift=1
689
690      cl_tmp=TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_index-il_shift))//','//&
691         &   TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_first))//':'//&
692         &   TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_last))//&
693         &   '('//TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_width))//')'
694      DO ji=2,tl_bdy(jp_north)%i_nseg
695         cl_tmp=TRIM(cl_tmp)//'|'//&
696            &   TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_index-il_shift))//','//&
697            &   TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_first))//':'//&
698            &   TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_last))
699      ENDDO
700      tl_att=att_init("bdy_north",TRIM(cl_tmp))
701      CALL file_add_att(tl_fileout, tl_att)
702   ENDIF
703
704   IF( tl_bdy(jp_south)%l_use )THEN
705     
706      cl_tmp=TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_index))//','//&
707         &   TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_first))//':'//&
708         &   TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_last))//&
709         &   '('//TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_width))//')'
710      DO ji=2,tl_bdy(jp_south)%i_nseg
711         cl_tmp=TRIM(cl_tmp)//'|'//&
712            &   TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_index))//','//&
713            &   TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_first))//':'//&
714            &   TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_last))
715      ENDDO
716
717      tl_att=att_init("bdy_south",TRIM(cl_tmp))
718      CALL file_add_att(tl_fileout, tl_att)
719   ENDIF
720
721   IF( tl_bdy(jp_east)%l_use )THEN
722      ! add shift on east boundary
723      ! boundary compute on T point but express on U or V point
724      il_shift=1
725
726      cl_tmp=TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_index-il_shift))//','//&
727         &   TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_first))//':'//&
728         &   TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_last))//&
729         &   '('//TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_width))//')'
730      DO ji=2,tl_bdy(jp_east)%i_nseg
731         cl_tmp=TRIM(cl_tmp)//'|'//&
732            &   TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_index-il_shift))//','//&
733            &   TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_first))//':'//&
734            &   TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_last))
735      ENDDO
736
737      tl_att=att_init("bdy_east",TRIM(cl_tmp))
738      CALL file_add_att(tl_fileout, tl_att)
739   ENDIF
740
741   IF( tl_bdy(jp_west)%l_use )THEN
742
743      cl_tmp=TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_index))//','//&
744         &   TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_first))//':'//&
745         &   TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_last))//&
746         &   '('//TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_width))//')'
747      DO ji=2,tl_bdy(jp_west)%i_nseg
748         cl_tmp=TRIM(cl_tmp)//'|'//&
749            &   TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_index))//','//&
750            &   TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_first))//':'//&
751            &   TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_last))
752      ENDDO
753
754      tl_att=att_init("bdy_west",TRIM(cl_tmp))
755      CALL file_add_att(tl_fileout, tl_att)
756   ENDIF
757
758   ! create file
759   CALL iom_create(tl_fileout)
760
761   ! write file
762   CALL iom_write_file(tl_fileout)
763
764   ! close file
765   CALL iom_close(tl_fileout)
766
767   ! clean
768   CALL att_clean(tl_att)
769   CALL file_clean(tl_fileout)
770   CALL mpp_clean(tl_bathy1)
771   CALL mpp_clean(tl_bathy0)
772   DEALLOCATE(dl_weight)
773   CALL boundary_clean(tl_bdy(:))
774   CALL var_clean_extra()
775
776   ! close log file
777   CALL logger_footer()
778   CALL logger_close()
779
780CONTAINS
781   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
782   SUBROUTINE merge_bathy_get_boundary(td_bathy0, td_bathy1, td_bdy, &
783         &                             id_rho, id_ncrs,              &
784         &                             dd_refined, dd_weight, dd_fill)
785   !-------------------------------------------------------------------
786   !> @brief
787   !> This subroutine compute refined bathymetry on boundary from coarse grid.
788   !>
789   !> @author J.Paul
790   !> @date November, 2013 - Initial Version
791   !>
792   !> @param[in] td_bathy0       coarse grid bathymetry file structure
793   !> @param[in] td_bathy1       fine grid bathymetry file structure
794   !> @param[in] td_bdy          boundary structure
795   !> @param[in] id_rho          array of refinement factor
796   !> @param[in] id_ncrs         number of point with coarse value save at boundaries
797   !> @param[inout] dd_refined   array of refined bathymetry
798   !> @param[inout] dd_weight    array of weight
799   !> @param[in] dd_fill         fillValue
800   !>
801   !-------------------------------------------------------------------
802
803      IMPLICIT NONE
804
805      ! Argument
806      TYPE(TMPP)                     , INTENT(IN   ) :: td_bathy0
807      TYPE(TMPP)                     , INTENT(IN   ) :: td_bathy1
808      TYPE(TBDY)                     , INTENT(IN   ) :: td_bdy
809      INTEGER(i4), DIMENSION(:)      , INTENT(IN   ) :: id_rho
810      INTEGER(i4)                    , INTENT(IN   ) :: id_ncrs
811      REAL(dp)   , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_refined
812      REAL(dp)   , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_weight
813      REAL(dp)                       , INTENT(IN   ) :: dd_fill 
814
815      ! local variable
816      INTEGER(i4) :: il_imin1
817      INTEGER(i4) :: il_imax1
818      INTEGER(i4) :: il_jmin1
819      INTEGER(i4) :: il_jmax1
820     
821      INTEGER(i4) :: il_imin0
822      INTEGER(i4) :: il_imax0
823      INTEGER(i4) :: il_jmin0
824      INTEGER(i4) :: il_jmax0
825
826      INTEGER(i4) :: il_width
827
828      INTEGER(i4), DIMENSION(2,2)         :: il_offset
829      INTEGER(i4), DIMENSION(2,2)         :: il_ind
830
831      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_tmp1d
832      REAL(dp)   , DIMENSION(:,:)    , ALLOCATABLE :: dl_tmp2d
833      REAL(dp)   , DIMENSION(:,:)    , ALLOCATABLE :: dl_wseg
834
835      TYPE(TVAR) :: tl_var0
836      TYPE(TVAR) :: tl_lon1
837      TYPE(TVAR) :: tl_lat1
838
839      TYPE(TMPP)  :: tl_bathy1
840      TYPE(TMPP)  :: tl_bathy0
841
842      TYPE(TDOM)  :: tl_dom1
843      TYPE(TDOM)  :: tl_dom0
844
845      ! loop indices
846      INTEGER(i4) :: ji
847      INTEGER(i4) :: jl
848      !----------------------------------------------------------------
849
850      IF( td_bdy%l_use )THEN
851         DO jl=1,td_bdy%i_nseg
852            ! get boundary definition
853            SELECT CASE(TRIM(td_bdy%c_card))
854            CASE('north')
855
856               il_imin1=td_bdy%t_seg(jl)%i_first
857               il_imax1=td_bdy%t_seg(jl)%i_last 
858               il_jmin1=td_bdy%t_seg(jl)%i_index-(td_bdy%t_seg(jl)%i_width-1)
859               il_jmax1=td_bdy%t_seg(jl)%i_index
860
861               ! do not used grid point to compute
862               ! boundaries indices (cf create_boundary)
863               ! as Bathymetry always on T point
864
865            CASE('south')
866
867               il_imin1=td_bdy%t_seg(jl)%i_first
868               il_imax1=td_bdy%t_seg(jl)%i_last 
869               il_jmin1=td_bdy%t_seg(jl)%i_index
870               il_jmax1=td_bdy%t_seg(jl)%i_index+(td_bdy%t_seg(jl)%i_width-1)
871
872            CASE('east')
873
874               il_imin1=td_bdy%t_seg(jl)%i_index-(td_bdy%t_seg(jl)%i_width-1)
875               il_imax1=td_bdy%t_seg(jl)%i_index
876               il_jmin1=td_bdy%t_seg(jl)%i_first
877               il_jmax1=td_bdy%t_seg(jl)%i_last 
878
879               ! do not used grid point to compute
880               ! boundaries indices (cf create_boundary)
881               ! as Bathymetry always on T point
882
883            CASE('west')
884
885               il_imin1=td_bdy%t_seg(jl)%i_index
886               il_imax1=td_bdy%t_seg(jl)%i_index+(td_bdy%t_seg(jl)%i_width-1)
887               il_jmin1=td_bdy%t_seg(jl)%i_first
888               il_jmax1=td_bdy%t_seg(jl)%i_last 
889
890            END SELECT
891
892            ! -read fine grid domain
893            tl_bathy1=mpp_copy(td_bathy1)
894
895            ! compute domain
896            tl_dom1=dom_init( tl_bathy1,         &
897            &                 il_imin1, il_imax1,&
898            &                 il_jmin1, il_jmax1,&
899            &                 TRIM(td_bdy%c_card))
900
901            ! add extra band to fine grid domain (if possible)
902            ! to avoid dimension of one and so be able to compute offset
903            CALL dom_add_extra(tl_dom1, id_rho(jp_I), id_rho(jp_J))
904
905            ! open mpp files over domain
906            CALL iom_dom_open(tl_bathy1, tl_dom1)
907
908            ! read variable value on domain
909            tl_lon1=iom_dom_read_var(tl_bathy1,'longitude',tl_dom1)
910            tl_lat1=iom_dom_read_var(tl_bathy1,'latitude' ,tl_dom1)
911
912            ! close mpp files
913            CALL iom_dom_close(tl_bathy1)
914
915            ! clean structure
916            CALL mpp_clean(tl_bathy1)
917
918            ! get coarse grid indices
919            il_ind(:,:)=grid_get_coarse_index(td_bathy0, tl_lon1, tl_lat1, &
920            &                                 id_rho=id_rho(:))
921
922            il_imin0=il_ind(1,1)
923            il_imax0=il_ind(1,2)
924
925            il_jmin0=il_ind(2,1)
926            il_jmax0=il_ind(2,2)
927
928            ! read coarse grid bathymetry on domain
929            tl_bathy0=mpp_copy(td_bathy0)
930
931            ! compute domain
932            tl_dom0=dom_init( tl_bathy0,         &
933            &                 il_imin0, il_imax0,&
934            &                 il_jmin0, il_jmax0 )
935
936            il_offset(:,:)= grid_get_fine_offset(tl_bathy0,         &
937            &                                    il_imin0, il_jmin0,&
938            &                                    il_imax0, il_jmax0,&
939            &                                    tl_lon1%d_value(:,:,1,1), &
940            &                                    tl_lat1%d_value(:,:,1,1), &
941            &                                    id_rho=id_rho(:))
942
943            ! clean
944            CALL var_clean(tl_lon1)
945            CALL var_clean(tl_lat1)
946
947            ! add extra band (if possible) to compute interpolation
948            CALL dom_add_extra(tl_dom0)
949
950            ! open mpp files over domain
951            CALL iom_dom_open(tl_bathy0, tl_dom0)
952
953            ! read variable value on domain
954            tl_var0=iom_dom_read_var(tl_bathy0,'Bathymetry',tl_dom0)
955
956            ! force to use nearest interpolation
957            tl_var0%c_interp(1)='nearest'
958
959            ! close mpp files
960            CALL iom_dom_close(tl_bathy0)
961
962            ! clean structure
963            CALL mpp_clean(tl_bathy0)
964
965            ! interpolate variable
966            CALL merge_bathy_interp( tl_var0,         &
967            &                        id_rho(:),       &
968            &                        il_offset(:,:) )
969
970            ! remove extraband added to domain
971            CALL dom_del_extra( tl_var0, tl_dom0, id_rho(:) )
972
973            ! remove extraband added to domain
974            CALL dom_clean_extra( tl_dom0 )
975
976            ! remove extraband added to fine grid domain
977            CALL dom_del_extra( tl_var0, tl_dom1 )
978
979            ! remove extraband added to fine grid domain
980            CALL dom_clean_extra( tl_dom1 )
981
982            ! fill refined array
983            dd_refined( il_imin1:il_imax1, &
984            &           il_jmin1:il_jmax1, &
985            &           :,: )= tl_var0%d_value(:,:,:,:)
986
987            ! clean
988            CALL var_clean(tl_var0)
989
990            ! compute weight function
991            ALLOCATE( dl_tmp1d(td_bdy%t_seg(jl)%i_width) )
992
993            SELECT CASE(TRIM(td_bdy%c_card))
994            CASE('north')
995
996
997               ! save n coarse point
998               il_width=td_bdy%t_seg(jl)%i_width-id_ncrs
999               ! compute "distance"
1000               dl_tmp1d(:)=(/(ji,ji=il_width,1,-1),(0,ji=1,id_ncrs)/)
1001
1002               ! compute weight on segment
1003               dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / &
1004               &                           (il_width) )
1005
1006               ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, &
1007               &                 tl_dom1%t_dim(2)%i_len) )
1008               dl_wseg(:,:)=dd_fill
1009               dl_wseg(:,:)=SPREAD( dl_tmp1d(:), &
1010               &                    DIM=1,       &
1011               &                    NCOPIES=tl_dom1%t_dim(1)%i_len )
1012
1013            CASE('south')
1014
1015               ! save n coarse point
1016               il_width=td_bdy%t_seg(jl)%i_width-id_ncrs
1017               ! compute "distance"
1018               dl_tmp1d(:)=(/(0,ji=1,id_ncrs),(ji,ji=1,il_width)/)
1019
1020               ! compute weight on segment
1021               dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / &
1022               &                           (il_width) )
1023
1024               ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, &
1025               &                 tl_dom1%t_dim(2)%i_len) )
1026               dl_wseg(:,:)=dd_fill
1027               dl_wseg(:,:)=SPREAD( dl_tmp1d(:), &
1028               &                    DIM=1,       &
1029               &                    NCOPIES=tl_dom1%t_dim(1)%i_len )
1030
1031            CASE('east')
1032
1033               ! save n coarse point
1034               il_width=td_bdy%t_seg(jl)%i_width-id_ncrs
1035               ! compute "distance"
1036               dl_tmp1d(:)=(/(ji,ji=il_width,1,-1),(0,ji=1,id_ncrs)/)
1037
1038               ! compute weight on segment
1039               dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / &
1040               &                           (il_width) )
1041
1042               ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, &
1043               &                 tl_dom1%t_dim(2)%i_len) )
1044               dl_wseg(:,:)=dd_fill
1045               dl_wseg(:,:)=SPREAD( dl_tmp1d(:), &
1046               &                    DIM=2,       &
1047               &                    NCOPIES=tl_dom1%t_dim(2)%i_len )
1048
1049            CASE('west')
1050
1051               ! save n coarse point
1052               il_width=td_bdy%t_seg(jl)%i_width-id_ncrs
1053               ! compute "distance"
1054               dl_tmp1d(:)=(/(0,ji=1,id_ncrs),(ji,ji=1,il_width)/)
1055
1056               ! compute weight on segment
1057               dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / &
1058               &                           (il_width) )
1059
1060               ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, &
1061               &                 tl_dom1%t_dim(2)%i_len) )
1062               dl_wseg(:,:)=dd_fill
1063               dl_wseg(:,:)=SPREAD( dl_tmp1d(:), &
1064               &                    DIM=2,       &
1065               &                    NCOPIES=tl_dom1%t_dim(2)%i_len )
1066
1067            END SELECT
1068
1069            DEALLOCATE( dl_tmp1d )
1070
1071            ! fill weight array
1072            ALLOCATE( dl_tmp2d( tl_dom1%t_dim(1)%i_len, &
1073            &                   tl_dom1%t_dim(2)%i_len) )
1074
1075            dl_tmp2d(:,:)=dd_weight( il_imin1:il_imax1, &
1076            &                        il_jmin1:il_jmax1, &
1077            &                        1,1 )
1078
1079            WHERE( dl_tmp2d(:,:) == dd_fill )
1080               dl_tmp2d(:,:)= dl_wseg(:,:)
1081            ELSE WHERE
1082               dl_tmp2d(:,:)= MAX( dl_tmp2d(:,:), dl_wseg(:,:) )
1083            END WHERE
1084
1085            dd_weight( il_imin1:il_imax1, &
1086            &          il_jmin1:il_jmax1, &
1087            &          1,1 ) = dl_tmp2d(:,:)
1088
1089            DEALLOCATE( dl_wseg )
1090            DEALLOCATE( dl_tmp2d )
1091
1092         ENDDO
1093      ENDIF
1094   END SUBROUTINE merge_bathy_get_boundary
1095   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1096   SUBROUTINE merge_bathy_interp(td_var, id_rho, id_offset, id_iext, id_jext)
1097   !-------------------------------------------------------------------
1098   !> @brief
1099   !> This subroutine interpolate variable.
1100   !>
1101   !> @author J.Paul
1102   !> @date November, 2013 - Initial Version
1103   !>
1104   !> @param[inout] td_var variable structure
1105   !> @param[in] id_rho    array of refinment factor
1106   !> @param[in] id_offset array of offset between fine and coarse grid
1107   !> @param[in] id_iext   i-direction size of extra bands (default=im_minext)
1108   !> @param[in] id_jext   j-direction size of extra bands (default=im_minext)
1109   !-------------------------------------------------------------------
1110
1111      IMPLICIT NONE
1112
1113      ! Argument
1114      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
1115      INTEGER(i4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
1116      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_offset
1117      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_iext
1118      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_jext
1119
1120      ! local variable
1121      TYPE(TVAR)  :: tl_mask
1122
1123      INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
1124
1125      INTEGER(i4) :: il_iext
1126      INTEGER(i4) :: il_jext
1127
1128      ! loop indices
1129      !----------------------------------------------------------------
1130
1131      !WARNING: two extrabands are required for cubic interpolation
1132      il_iext=3
1133      IF( PRESENT(id_iext) ) il_iext=id_iext
1134
1135      il_jext=3
1136      IF( PRESENT(id_jext) ) il_jext=id_jext
1137
1138      IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1139         CALL logger_warn("MERGE BATHY INTERP: at least extrapolation "//&
1140         &  "on two points are required with cubic interpolation ")
1141         il_iext=2
1142      ENDIF
1143
1144      IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1145         CALL logger_warn("MERGE BATHY INTERP: at least extrapolation "//&
1146         &  "on two points are required with cubic interpolation ")
1147         il_jext=2
1148      ENDIF
1149
1150      ! work on mask
1151      ! create mask
1152      ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
1153      &                td_var%t_dim(2)%i_len, &
1154      &                td_var%t_dim(3)%i_len, &
1155      &                td_var%t_dim(4)%i_len) )
1156
1157      bl_mask(:,:,:,:)=1
1158      WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0     
1159
1160      SELECT CASE(TRIM(td_var%c_point))
1161      CASE DEFAULT ! 'T'
1162         tl_mask=var_init('tmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
1163         &                id_ew=td_var%i_ew )
1164      CASE('U','V','F')
1165         CALL logger_fatal("MERGE BATHY INTERP: can not computed "//&
1166         &                 "interpolation on "//TRIM(td_var%c_point)//&
1167         &                 " grid point (variable "//TRIM(td_var%c_name)//&
1168         &                 "). check namelist.")
1169      END SELECT
1170
1171      DEALLOCATE(bl_mask)
1172
1173      ! interpolate mask
1174      CALL interp_fill_value( tl_mask, id_rho(:),  &
1175      &                       id_offset=id_offset(:,:) )
1176
1177      ! work on variable
1178      ! add extraband
1179      CALL extrap_add_extrabands(td_var, il_iext, il_iext)
1180
1181      ! extrapolate variable
1182      CALL extrap_fill_value( td_var )
1183
1184      ! interpolate Bathymetry
1185      CALL interp_fill_value( td_var, id_rho(:), &
1186      &                       id_offset=id_offset(:,:) )
1187
1188      ! remove extraband
1189      CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J))
1190
1191      ! keep original mask
1192      WHERE( tl_mask%d_value(:,:,:,:) == 0 )
1193         td_var%d_value(:,:,:,:)=td_var%d_fill
1194      END WHERE
1195
1196   END SUBROUTINE merge_bathy_interp
1197   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1198END PROGRAM merge_bathy
Note: See TracBrowser for help on using the repository browser.