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/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90 @ 7261

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: