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/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90 @ 4213

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

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

Download in other formats: