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.
merge_bathy.f90 in trunk/NEMOGCM/TOOLS/SIREN/src – NEMO

source: trunk/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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