source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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