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 branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90 @ 5989

Last change on this file since 5989 was 5989, checked in by deazer, 8 years ago

Merging TMB and 25h diagnostics to head of trunk
added brief documentation

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