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.

create_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/create_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: 32.6 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	!
6	! PROGRAM: create_bathy
7	!
8	! DESCRIPTION:
9	!> @brief
10	!> This program create bathymetry file.
11	!>
12	!> @details
13	!> Bathymetry could be extracted from fine grid Bathymetry file, or interpolated
14	!> from coarse grid Bathymetry file.
15	!>
16	!> @author
17	!> J.Paul
18	! REVISION HISTORY:
19	!> @date Nov, 2013 - Initial Version
20	!
21	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
22	!>
23	!> @todo
24	!> - add attributes indices and refinement in output file
25	!----------------------------------------------------------------------
26	!> @code
27	PROGRAM create_bathy
28
29	! USE netcdf ! nf90 library
30	USE global ! global variable
31	USE kind ! F90 kind parameter
32	USE logger ! log file manager
33	USE fct ! basic useful function
34	USE date ! date manager
35	USE att ! attribute manager
36	USE dim ! dimension manager
37	USE var ! variable manager
38	USE file ! file manager
39	USE multi ! multi file manager
40	USE iom ! I/O manager
41	USE dom ! domain manager
42	USE grid ! grid manager
43	USE extrap ! extrapolation manager
44	USE interp ! interpolation manager
45	USE filter ! filter manager
46	USE mpp ! MPP manager
47	USE iom_mpp ! MPP I/O manager
48
49	IMPLICIT NONE
50
51	! local variable
52	CHARACTER(LEN=lc) :: cl_namelist
53	CHARACTER(LEN=lc) :: cl_date
54	CHARACTER(LEN=lc) :: cl_data
55
56	INTEGER(i4) :: il_narg
57	INTEGER(i4) :: il_status
58	INTEGER(i4) :: il_fileid
59	INTEGER(i4) :: il_attid
60	INTEGER(i4) :: il_imin
61	INTEGER(i4) :: il_imax
62	INTEGER(i4) :: il_jmin
63	INTEGER(i4) :: il_jmax
64	INTEGER(i4) , DIMENSION(ip_maxdim) :: il_rho
65	INTEGER(i4) , DIMENSION(2,2) :: il_offset
66	INTEGER(i4) , DIMENSION(2,2,2) :: il_ind
67	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_mask
68
69	LOGICAL :: ll_exist
70
71	TYPE(TFILE) :: tl_coord0
72	TYPE(TFILE) :: tl_coord1
73	TYPE(TFILE) :: tl_file
74	TYPE(TFILE) :: tl_fileout
75
76	TYPE(TATT) :: tl_att
77
78	TYPE(TVAR) :: tl_lon
79	TYPE(TVAR) :: tl_lat
80	TYPE(TVAR) :: tl_depth
81	TYPE(TVAR) :: tl_time
82
83	TYPE(TVAR) :: tl_tmp
84	TYPE(TVAR) , DIMENSION(:), ALLOCATABLE :: tl_var
85
86	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
87
88	TYPE(TMULTI) :: tl_multi
89
90	! loop indices
91	INTEGER(i4) :: ji
92	INTEGER(i4) :: jj
93	INTEGER(i4) :: jk
94
95	! namelist variable
96	CHARACTER(LEN=lc) :: cn_logfile = 'create_bathy.log'
97	CHARACTER(LEN=lc) :: cn_verbosity = 'warning'
98
99	CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg'
100
101	CHARACTER(LEN=lc) :: cn_coord0 = ''
102	INTEGER(i4) :: in_perio0 = -1
103
104	CHARACTER(LEN=lc) :: cn_coord1 = ''
105	INTEGER(i4) :: in_perio1 = -1
106	LOGICAL :: ln_fillclosed = .TRUE.
107
108	CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varinfo = ''
109	CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varfile = ''
110
111	INTEGER(i4) :: in_imin0 = 0
112	INTEGER(i4) :: in_imax0 = 0
113	INTEGER(i4) :: in_jmin0 = 0
114	INTEGER(i4) :: in_jmax0 = 0
115	INTEGER(i4) :: in_rhoi = 1
116	INTEGER(i4) :: in_rhoj = 1
117
118	CHARACTER(LEN=lc) :: cn_fileout = 'bathy_fine.nc'
119	!-------------------------------------------------------------------
120
121	NAMELIST /namlog/ & !< logger namelist
122	& cn_logfile, & !< log file
123	& cn_verbosity !< log verbosity
124
125	NAMELIST /namcfg/ & !< configuration namelist
126	& cn_varcfg !< variable configuration file
127
128	NAMELIST /namcrs/ & !< coarse grid namelist
129	& cn_coord0, & !< coordinate file
130	& in_perio0 !< periodicity index
131
132	NAMELIST /namfin/ & !< fine grid namelist
133	& cn_coord1, & !< coordinate file
134	& in_perio1, & !< periodicity index
135	& ln_fillclosed !< fill closed sea
136
137	NAMELIST /namvar/ & !< variable namelist
138	& cn_varinfo, & !< list of variable and interpolation method to be used. (ex: 'votemper:linear','vosaline:cubic' )
139	& cn_varfile !< list of variable file
140
141	NAMELIST /namnst/ & !< nesting namelist
142	& in_imin0, & !< i-direction lower left point indice
143	& in_imax0, & !< i-direction upper right point indice
144	& in_jmin0, & !< j-direction lower left point indice
145	& in_jmax0, & !< j-direction upper right point indice
146	& in_rhoi, & !< refinement factor in i-direction
147	& in_rhoj !< refinement factor in j-direction
148
149	NAMELIST /namout/ & !< output namlist
150	& cn_fileout !< fine grid bathymetry file
151	!-------------------------------------------------------------------
152
153	!1- namelist
154	!1-1 get namelist
155	il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
156	IF( il_narg/=1 )THEN
157	PRINT *,"ERROR in create_bathy: need a namelist"
158	STOP
159	ELSE
160	CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
161	ENDIF
162
163	!1-2 read namelist
164	INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
165	IF( ll_exist )THEN
166
167	il_fileid=fct_getunit()
168
169	OPEN( il_fileid, FILE=TRIM(cl_namelist), &
170	& FORM='FORMATTED', &
171	& ACCESS='SEQUENTIAL', &
172	& STATUS='OLD', &
173	& ACTION='READ', &
174	& IOSTAT=il_status)
175	CALL fct_err(il_status)
176	IF( il_status /= 0 )THEN
177	PRINT *,"ERROR in create_bathy: error opening "//TRIM(cl_namelist)
178	STOP
179	ENDIF
180
181	READ( il_fileid, NML = namlog )
182	!1-2-1 define log file
183	CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity))
184	CALL logger_header()
185
186	READ( il_fileid, NML = namcfg )
187	!1-2-2 get variable extra information
188	CALL var_def_extra(TRIM(cn_varcfg))
189
190	READ( il_fileid, NML = namcrs )
191	READ( il_fileid, NML = namfin )
192	READ( il_fileid, NML = namvar )
193	!1-2-3 add user change in extra information
194	CALL var_chg_extra(cn_varinfo)
195	!1-2-4 match variable with file
196	tl_multi=multi_init(cn_varfile)
197
198	READ( il_fileid, NML = namnst )
199	READ( il_fileid, NML = namout )
200
201	CLOSE( il_fileid, IOSTAT=il_status )
202	CALL fct_err(il_status)
203	IF( il_status /= 0 )THEN
204	CALL logger_error("CREATE BATHY: closing "//TRIM(cl_namelist))
205	ENDIF
206
207	ELSE
208
209	PRINT *,"ERROR in create_bathy: can't find "//TRIM(cl_namelist)
210
211	ENDIF
212
213	!2- open files
214	IF( cn_coord0 /= '' )THEN
215	tl_coord0=file_init(TRIM(cn_coord0),id_perio=in_perio0)
216	CALL iom_open(tl_coord0)
217	ELSE
218	CALL logger_fatal("CREATE BATHY: no coarse grid coordinate found. "//&
219	& "check namelist")
220	ENDIF
221
222	IF( TRIM(cn_coord1) /= '' )THEN
223	tl_coord1=file_init(TRIM(cn_coord1),id_perio=in_perio1)
224	CALL iom_open(tl_coord1)
225	ELSE
226	CALL logger_fatal("CREATE BATHY: no fine grid coordinate found. "//&
227	& "check namelist")
228	ENDIF
229
230	!3- check
231	!3-1 check output file do not already exist
232	INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
233	IF( ll_exist )THEN
234	CALL logger_fatal("CREATE BATHY: output file "//TRIM(cn_fileout)//&
235	& " already exist.")
236	ENDIF
237
238	!3-2 check namelist
239	!3-2-1 check refinement factor
240	il_rho(:)=1
241	IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
242	CALL logger_error("CREATE BATHY: invalid refinement factor."//&
243	& " check namelist "//TRIM(cl_namelist))
244	ELSE
245	il_rho(jp_I)=in_rhoi
246	il_rho(jp_J)=in_rhoj
247	ENDIF
248
249	!3-2-2 check domain indices
250	IF( in_imin0 < 1 .OR. in_imax0 < 1 .OR. &
251	& in_jmin0 < 1 .OR. in_jmax0 < 1)THEN
252	! compute coarse grid indices around fine grid
253	IF( cn_coord0 /= '' )THEN
254	il_ind(:,:,:)=grid_get_coarse_index( tl_coord0, tl_coord1, &
255	& id_rho=il_rho(:) )
256	ENDIF
257
258	il_imin=il_ind(1,1,1) ; il_imax=il_ind(1,2,1)
259	il_jmin=il_ind(2,1,1) ; il_jmax=il_ind(2,2,1)
260
261	il_offset(:,:)=il_ind(:,:,2)
262	ELSE
263	il_imin=in_imin0 ; il_imax=in_imax0
264	il_jmin=in_jmin0 ; il_jmax=in_jmax0
265
266	il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5)
267	il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5)
268	ENDIF
269
270	!3-2-3 check domain validity
271	IF( cn_coord0 /= '' )THEN
272	CALL grid_check_dom(tl_coord0, il_imin, il_imax, il_jmin, il_jmax)
273	ENDIF
274
275	!3-2-4 check coincidence between coarse and fine grid
276	IF( cn_coord0 /= '' )THEN
277	CALL grid_check_coincidence( tl_coord0, tl_coord1, &
278	& il_imin, il_imax, &
279	& il_jmin, il_jmax, &
280	& il_rho(:) )
281	ENDIF
282
283	IF( .NOT. ASSOCIATED(tl_multi%t_file) )THEN
284	CALL logger_error("CREATE BATHY: no file to work on. "//&
285	& "check cn_varfile in namelist.")
286	ELSE
287	ALLOCATE( tl_var( tl_multi%i_nvar ) )
288	jk=0
289	DO ji=1,tl_multi%i_nfile
290	WRITE(cl_data,'(a,i2.2)') 'data_',jk+1
291	IF( .NOT. ASSOCIATED(tl_multi%t_file(ji)%t_var) )THEN
292	CALL logger_error("CREATE BATHY: no variable to work on for "//&
293	& "file"//TRIM(tl_multi%t_file(ji)%c_name)//&
294	& ". check cn_varfile in namelist.")
295	ELSEIF( TRIM(tl_multi%t_file(ji)%c_name) == TRIM(cl_data) )THEN
296	DO jj=1,tl_multi%t_file(ji)%i_nvar
297	jk=jk+1
298	tl_tmp=tl_multi%t_file(ji)%t_var(jj)
299	!- use input matrix to initialise variable
300	tl_var(jk)=create_bathy_matrix(tl_tmp, tl_coord1)
301	ENDDO
302	ELSE
303	! open file
304	tl_file=file_init(TRIM(tl_multi%t_file(ji)%c_name))
305	CALL iom_open(tl_file)
306
307	! get or check depth value
308	IF( tl_file%i_depthid /= 0 )THEN
309	IF( ASSOCIATED(tl_depth%d_value) )THEN
310	IF( ANY( tl_depth%d_value(:,:,:,:) /= &
311	& tl_tmp%d_value(:,:,:,:) ) )THEN
312	CALL logger_fatal("CREATE BATHY: depth value from "//&
313	& TRIM(tl_multi%t_file(ji)%c_name)//" not conform "//&
314	& " to those from former file(s).")
315	ENDIF
316	ELSE
317	tl_depth=iom_read_var(tl_file,tl_file%i_depthid)
318	ENDIF
319	ENDIF
320
321	! get or check time value
322	IF( tl_file%i_timeid /= 0 )THEN
323	IF( ASSOCIATED(tl_time%d_value) )THEN
324	IF( ANY( tl_time%d_value(:,:,:,:) /= &
325	& tl_tmp%d_value(:,:,:,:) ) )THEN
326	CALL logger_fatal("CREATE BATHY: time value from "//&
327	& TRIM(tl_multi%t_file(ji)%c_name)//" not conform "//&
328	& " to those from former file(s).")
329	ENDIF
330	ELSE
331	tl_time=iom_read_var(tl_file,tl_file%i_timeid)
332	ENDIF
333	ENDIF
334
335	IF( ANY( tl_file%t_dim(1:2)%i_len /= &
336	& tl_coord0%t_dim(1:2)%i_len) )THEN
337	DO jj=1,tl_multi%t_file(ji)%i_nvar
338	jk=jk+1
339	tl_tmp=tl_multi%t_file(ji)%t_var(jj)
340	!- extract bathymetry from fine grid bathymetry
341	tl_var(jk)=create_bathy_extract( tl_tmp, tl_file, &
342	& tl_coord1 )
343	ENDDO
344	ELSE
345	DO jj=1,tl_multi%t_file(ji)%i_nvar
346	jk=jk+1
347	tl_tmp=tl_multi%t_file(ji)%t_var(jj)
348	!- get bathymetry from coarse grid bathymetry
349	tl_var(jk)=create_bathy_get_var( tl_tmp, tl_file, &
350	& il_imin, il_jmin, &
351	& il_imax, il_jmax, &
352	& il_offset(:,:), &
353	& il_rho(:) )
354	ENDDO
355	ENDIF
356
357	! close file
358	CALL iom_close(tl_file)
359	! clean structure
360	CALL file_clean(tl_file)
361	ENDIF
362	ENDDO
363	ENDIF
364
365	DO jk=1,tl_multi%i_nvar
366	!6- forced min and max value
367	CALL var_limit_value(tl_var(jk))
368
369	!7- fill closed sea
370	IF( TRIM(tl_var(jk)%c_stdname) == 'bathymetry' .AND. &
371	ln_fillclosed )THEN
372	ALLOCATE( il_mask(tl_var(jk)%t_dim(1)%i_len, &
373	& tl_var(jk)%t_dim(2)%i_len) )
374
375	!7-1 split domain in N sea subdomain
376	il_mask(:,:)=grid_split_domain(tl_var(jk))
377
378	!7-2 fill smallest domain
379	CALL grid_fill_small_dom( tl_var(jk), il_mask(:,:) )
380
381	DEALLOCATE( il_mask )
382	ENDIF
383
384	!8- filter
385	CALL filter_fill_value(tl_var(jk))
386
387	!9- check bathymetry
388	IF( TRIM(tl_var(jk)%c_stdname) == 'bathymetry' .AND. &
389	& MINVAL(tl_var(jk)%d_value(:,:,:,:)) <= 0._dp )THEN
390	CALL logger_error("CREATE BATHY: Bathymetry has value <= 0")
391	ENDIF
392	ENDDO
393
394
395	!10- create file
396	tl_fileout=file_init(TRIM(cn_fileout),id_perio=in_perio1)
397
398	!10-1 add dimension
399	tl_dim(:)=var_max_dim(tl_var(:))
400
401	DO ji=1,ip_maxdim
402	IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
403	ENDDO
404
405	!10-2 add variables
406	IF( ALL( tl_dim(1:2)%l_use ) )THEN
407
408	! add longitude
409	tl_lon=iom_read_var(tl_coord1,'longitude')
410	CALL file_add_var(tl_fileout, tl_lon)
411	CALL var_clean(tl_lon)
412
413	! add latitude
414	tl_lat=iom_read_var(tl_coord1,'latitude')
415	CALL file_add_var(tl_fileout, tl_lat)
416	CALL var_clean(tl_lat)
417
418	ENDIF
419
420	IF( tl_dim(3)%l_use )THEN
421	! add depth
422	CALL file_add_var(tl_fileout, tl_depth)
423	CALL var_clean(tl_depth)
424	ENDIF
425
426	IF( tl_dim(4)%l_use )THEN
427	! add time
428	CALL file_add_var(tl_fileout, tl_time)
429	CALL var_clean(tl_time)
430	ENDIF
431
432	! add other variables
433	DO jk=1,tl_multi%i_nvar
434	CALL file_add_var(tl_fileout, tl_var(jk))
435	CALL var_clean(tl_var(jk))
436	ENDDO
437
438	!10-3 add some attribute
439	tl_att=att_init("Created_by","SIREN create_bathy")
440	CALL file_add_att(tl_fileout, tl_att)
441
442	cl_date=date_print(date_now())
443	tl_att=att_init("Creation_date",cl_date)
444	CALL file_add_att(tl_fileout, tl_att)
445
446	! add attribute periodicity
447	il_attid=0
448	IF( ASSOCIATED(tl_fileout%t_att) )THEN
449	il_attid=att_get_id(tl_fileout%t_att(:),'periodicity')
450	ENDIF
451	IF( tl_coord1%i_perio >= 0 .AND. il_attid == 0 )THEN
452	tl_att=att_init('periodicity',tl_coord1%i_perio)
453	CALL file_add_att(tl_fileout,tl_att)
454	ENDIF
455
456	il_attid=0
457	IF( ASSOCIATED(tl_fileout%t_att) )THEN
458	il_attid=att_get_id(tl_fileout%t_att(:),'ew_overlap')
459	ENDIF
460	IF( tl_coord1%i_ew >= 0 .AND. il_attid == 0 )THEN
461	tl_att=att_init('ew_overlap',tl_coord1%i_ew)
462	CALL file_add_att(tl_fileout,tl_att)
463	ENDIF
464
465	!10-4 create file
466	CALL iom_create(tl_fileout)
467
468	!10-5 write file
469	CALL iom_write_file(tl_fileout)
470
471	!10-6 close file
472	CALL iom_close(tl_fileout)
473	IF( cn_coord0 /= '' ) CALL iom_close(tl_coord0)
474
475	!11- clean
476	DEALLOCATE(tl_var)
477
478	CALL file_clean(tl_fileout)
479	CALL file_clean(tl_coord1)
480	CALL file_clean(tl_coord0)
481
482	! close log file
483	CALL logger_footer()
484	CALL logger_close()
485
486	!> @endcode
487	CONTAINS
488	!-------------------------------------------------------------------
489	!> @brief
490	!> This subroutine
491	!>
492	!> @details
493	!>
494	!> @author J.Paul
495	!> - Nov, 2013- Initial Version
496	!>
497	!> @param[in]
498	!-------------------------------------------------------------------
499	!> @code
500	FUNCTION create_bathy_matrix(td_var, td_coord)
501	IMPLICIT NONE
502	! Argument
503	TYPE(TVAR) , INTENT(IN) :: td_var
504	TYPE(TFILE), INTENT(IN) :: td_coord
505
506	! function
507	TYPE(TVAR) :: create_bathy_matrix
508
509	! local variable
510	INTEGER(i4) :: il_ighost
511	INTEGER(i4) :: il_jghost
512	INTEGER(i4) , DIMENSION(2) :: il_xghost
513	INTEGER(i4) , DIMENSION(3) :: il_dim
514	INTEGER(i4) , DIMENSION(3) :: il_size
515	INTEGER(i4) , DIMENSION(3) :: il_rest
516
517	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ishape
518	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_jshape
519	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_kshape
520
521	REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
522
523	TYPE(TVAR) :: tl_lon
524	TYPE(TVAR) :: tl_lat
525	TYPE(TVAR) :: tl_var
526	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
527
528	! loop indices
529	INTEGER(i4) :: ji
530	INTEGER(i4) :: jj
531	INTEGER(i4) :: jk
532	!----------------------------------------------------------------
533
534	!1- read output grid
535	tl_lon=iom_read_var(td_coord,'longitude')
536	tl_lat=iom_read_var(td_coord,'latitude')
537
538	!2- look for ghost cell
539	il_xghost(:)=grid_get_ghost( tl_lon, tl_lat )
540
541	il_ighost=il_xghost(1)*ig_ghost
542	il_jghost=il_xghost(2)*ig_ghost
543
544	!3- write value on grid
545	!3-1 get matrix dimension
546	il_dim(:)=td_var%t_dim(1:3)%i_len
547	!3-2 output dimension
548	tl_dim(:)=tl_lon%t_dim(:)
549
550	! remove ghost cell
551	tl_dim(1)%i_len=tl_dim(1)%i_len - 2il_xghost(1)ig_ghost
552	tl_dim(2)%i_len=tl_dim(2)%i_len - 2il_xghost(2)ig_ghost
553
554	!3-3 split output domain in N subdomain depending of matrix dimension
555	il_size(:) = tl_dim(1:3)%i_len / il_dim(:)
556	il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:))
557
558	ALLOCATE( il_ishape(il_dim(1)+1) )
559	il_ishape(:)=0
560	DO ji=2,il_dim(1)+1
561	il_ishape(ji)=il_ishape(ji-1)+il_size(1)
562	ENDDO
563	! add rest to last cell
564	il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1)
565
566
567	ALLOCATE( il_jshape(il_dim(2)+1) )
568	il_jshape(:)=0
569	DO jj=2,il_dim(2)+1
570	il_jshape(jj)=il_jshape(jj-1)+il_size(2)
571	ENDDO
572	! add rest to last cell
573	il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2)
574
575	ALLOCATE( il_kshape(il_dim(3)+1) )
576	il_kshape(:)=0
577	DO jk=2,il_dim(3)+1
578	il_kshape(jk)=il_kshape(jk-1)+il_size(3)
579	ENDDO
580	! add rest to last cell
581	il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3)
582
583	!3-3 write ouput table of value
584	ALLOCATE(dl_value( tl_dim(1)%i_len, &
585	& tl_dim(2)%i_len, &
586	& tl_dim(3)%i_len, &
587	& tl_dim(4)%i_len) )
588
589	dl_value(:,:,:,:)=0
590
591	DO jk=2,il_dim(3)+1
592	DO jj=2,il_dim(2)+1
593	DO ji=2,il_dim(1)+1
594
595	dl_value( 1+il_ishape(ji-1):il_ishape(ji), &
596	& 1+il_jshape(jj-1):il_jshape(jj), &
597	& 1+il_kshape(jk-1):il_kshape(jk), &
598	& 1 ) = td_var%d_value(ji-1,jj-1,jk-1,1)
599
600	ENDDO
601	ENDDO
602	ENDDO
603
604	!3-4 initialise variable with value
605	tl_var=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:))
606
607	DEALLOCATE(dl_value)
608
609	!4- add ghost cell
610	CALL grid_add_ghost(tl_var,il_ighost,il_jghost)
611
612	!5- save result
613	create_bathy_matrix=tl_var
614
615	END FUNCTION create_bathy_matrix
616	!> @endcode
617	!-------------------------------------------------------------------
618	!> @brief
619	!> This subroutine
620	!>
621	!> @details
622	!>
623	!> @author J.Paul
624	!> - Nov, 2013- Initial Version
625	!>
626	!> @param[in]
627	!-------------------------------------------------------------------
628	!> @code
629	FUNCTION create_bathy_extract(td_var, td_file, &
630	& td_coord)
631	IMPLICIT NONE
632	! Argument
633	TYPE(TVAR) , INTENT(IN) :: td_var
634	TYPE(TFILE), INTENT(IN) :: td_file
635	TYPE(TFILE), INTENT(IN) :: td_coord
636
637	! function
638	TYPE(TVAR) :: create_bathy_extract
639
640	! local variable
641	INTEGER(i4), DIMENSION(2,2,2) :: il_ind
642
643	INTEGER(i4) :: il_pivot
644	INTEGER(i4) :: il_perio
645
646	INTEGER(i4) :: il_imin
647	INTEGER(i4) :: il_jmin
648	INTEGER(i4) :: il_imax
649	INTEGER(i4) :: il_jmax
650
651	TYPE(TFILE) :: tl_file
652
653	TYPE(TMPP) :: tl_mpp
654
655	TYPE(TATT) :: tl_att
656
657	TYPE(TVAR) :: tl_var
658
659	TYPE(TDOM) :: tl_dom
660	! loop indices
661	!----------------------------------------------------------------
662
663	IF( td_file%i_id == 0 )THEN
664	CALL logger_error("CREATE BATHY EXTRACT: file "//&
665	& TRIM(td_file%c_name)//" not opened ")
666	ELSE
667
668	!init
669	tl_file=td_file
670
671	!1- open file
672	CALL iom_open(tl_file)
673
674	! get periodicity
675	il_pivot=grid_get_pivot(tl_file)
676	il_perio=grid_get_perio(tl_file,il_pivot)
677
678	tl_file%i_perio=il_perio
679
680	!2- compute file grid indices around coord grid
681	il_ind(:,:,:)=grid_get_coarse_index(tl_file, td_coord )
682
683	il_imin=il_ind(1,1,1) ; il_imax=il_ind(1,2,1)
684	il_jmin=il_ind(2,1,1) ; il_jmax=il_ind(2,2,1)
685
686	IF( ANY( il_ind(:,:,2) /= 0 ) )THEN
687	CALL logger_error("CREATE BATHY EXTRACT: something wrong "//&
688	& " find offset when extracting data")
689	ENDIF
690	!3- check grid coincidence
691	CALL grid_check_coincidence( tl_file, td_coord, &
692	& il_imin, il_imax, &
693	& il_jmin, il_jmax, &
694	& (/1, 1, 1/) )
695
696	!4- compute domain
697	tl_dom=dom_init(tl_file, &
698	& il_imin, il_imax, &
699	& il_jmin, il_jmax)
700
701	! close file
702	CALL iom_close(tl_file)
703
704	!5- read bathymetry on domain (ugly way to do it, have to work on it)
705	!5-1 init mpp structure
706	tl_mpp=mpp_init(tl_file)
707
708	CALL file_clean(tl_file)
709
710	!5-2 get processor to be used
711	CALL mpp_get_use( tl_mpp, tl_dom )
712
713	!5-3 open mpp files
714	CALL iom_mpp_open(tl_mpp)
715
716	!5-4 read variable on domain
717	tl_var=iom_mpp_read_var(tl_mpp,TRIM(td_var%c_name),td_dom=tl_dom)
718
719	!5-5 close mpp file
720	CALL iom_mpp_close(tl_mpp)
721
722	!6- add ghost cell
723	CALL grid_add_ghost(tl_var,tl_dom%i_ighost,tl_dom%i_jghost)
724
725	!7- check result
726	IF( ANY( tl_var%t_dim(:)%l_use .AND. &
727	& tl_var%t_dim(:)%i_len /= td_coord%t_dim(:)%i_len) )THEN
728	CALL logger_debug("CREATE BATHY EXTRACT: "//&
729	& "dimensoin of variable "//TRIM(td_var%c_name)//" "//&
730	& TRIM(fct_str(tl_var%t_dim(1)%i_len))//","//&
731	& TRIM(fct_str(tl_var%t_dim(2)%i_len))//","//&
732	& TRIM(fct_str(tl_var%t_dim(3)%i_len))//","//&
733	& TRIM(fct_str(tl_var%t_dim(4)%i_len)) )
734	CALL logger_debug("CREATE BATHY EXTRACT: "//&
735	& "dimensoin of coordinate file "//&
736	& TRIM(fct_str(td_coord%t_dim(1)%i_len))//","//&
737	& TRIM(fct_str(td_coord%t_dim(2)%i_len))//","//&
738	& TRIM(fct_str(td_coord%t_dim(3)%i_len))//","//&
739	& TRIM(fct_str(td_coord%t_dim(4)%i_len)) )
740	CALL logger_fatal("CREATE BATHY EXTRACT: "//&
741	& "dimensoin of extracted "//&
742	& "variable and coordinate file dimension differ")
743	ENDIF
744
745	!8- add attribute to variable
746	tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
747	CALL var_move_att(tl_var, tl_att)
748
749	tl_att=att_init('src_i-indices',(/tl_dom%i_imin, tl_dom%i_imax/))
750	CALL var_move_att(tl_var, tl_att)
751
752	tl_att=att_init('src_j-indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
753	CALL var_move_att(tl_var, tl_att)
754
755	!9- save result
756	create_bathy_extract=tl_var
757
758	! clean structure
759	CALL var_clean(tl_var)
760	CALL mpp_clean(tl_mpp)
761	ENDIF
762
763	END FUNCTION create_bathy_extract
764	!> @endcode
765	!-------------------------------------------------------------------
766	!> @brief
767	!> This subroutine
768	!>
769	!> @details
770	!>
771	!> @author J.Paul
772	!> - Nov, 2013- Initial Version
773	!>
774	!> @param[in] td_var : variable structure
775	!> @param[in] td_file : file structure
776	!> @param[in] id_imin : i-direction lower left corner indice
777	!> @param[in] id_imax : i-direction upper right corner indice
778	!> @param[in] id_jmin : j-direction lower left corner indice
779	!> @param[in] id_jmax : j-direction upper right corner indice
780	!> @param[in] id_rho : table of refinement factor
781	!-------------------------------------------------------------------
782	!> @code
783	FUNCTION create_bathy_get_var(td_var, td_file, &
784	& id_imin, id_jmin, &
785	& id_imax, id_jmax, &
786	& id_offset, &
787	& id_rho )
788	IMPLICIT NONE
789	! Argument
790	TYPE(TVAR) , INTENT(IN) :: td_var
791	TYPE(TFILE), INTENT(IN) :: td_file
792	INTEGER(i4), INTENT(IN) :: id_imin
793	INTEGER(i4), INTENT(IN) :: id_imax
794	INTEGER(i4), INTENT(IN) :: id_jmin
795	INTEGER(i4), INTENT(IN) :: id_jmax
796	INTEGER(i4), DIMENSION(:,:), INTENT(IN) :: id_offset
797	INTEGER(i4), DIMENSION(:) , INTENT(IN) :: id_rho
798
799	! function
800	TYPE(TVAR) :: create_bathy_get_var
801
802	! local variable
803	INTEGER(i4) :: il_pivot
804	INTEGER(i4) :: il_perio
805
806	TYPE(TFILE) :: tl_file
807
808	TYPE(TMPP) :: tl_mpp
809
810	TYPE(TATT) :: tl_att
811
812	TYPE(TVAR) :: tl_var
813
814	TYPE(TDOM) :: tl_dom
815
816	! loop indices
817	!----------------------------------------------------------------
818	IF( ANY(SHAPE(id_offset(:,:)) /= 2) )THEN
819	CALL logger_error("CREATE BATHY GET VAR: invalid dimensio of "//&
820	& "offset table")
821	ENDIF
822
823	!init
824	tl_file=td_file
825
826	!1- open file
827	CALL iom_open(tl_file)
828
829	! get periodicity
830	il_pivot=grid_get_pivot(tl_file)
831	il_perio=grid_get_perio(tl_file,il_pivot)
832
833	tl_file%i_perio=il_perio
834
835	!2- compute domain
836	tl_dom=dom_init(tl_file, &
837	& id_imin, id_imax, &
838	& id_jmin, id_jmax)
839
840	CALL dom_print(tl_dom)
841	print *,'id_offset ',id_offset(:,:)
842	!3- close file
843	CALL iom_close(tl_file)
844
845	!4- add extra band (if possible) to compute interpolation
846	CALL dom_add_extra(tl_dom)
847
848	!5- read bathymetry on domain (ugly way to do it, have to work on it)
849	!5-1 init mpp sturcutre
850	tl_mpp=mpp_init(tl_file)
851
852	CALL file_clean(tl_file)
853
854	!5-2 get processor to be used
855	CALL mpp_get_use( tl_mpp, tl_dom )
856
857	!5-3 open mpp files
858	CALL iom_mpp_open(tl_mpp)
859
860	!5-4 read variable value on domain
861	tl_var=iom_mpp_read_var(tl_mpp,TRIM(td_var%c_name),td_dom=tl_dom)
862
863	!5-5 close mpp files
864	CALL iom_mpp_close(tl_mpp)
865
866	!6- interpolate variable
867	CALL create_bathy_interp(tl_var, id_rho(:), id_offset(:,:))
868
869	!7- remove extraband added to domain
870	CALL dom_del_extra( tl_var, tl_dom, id_rho(:) )
871
872	!8- add ghost cell
873	CALL grid_add_ghost(tl_var,tl_dom%i_ighost,tl_dom%i_jghost)
874
875	!9- add attribute to variable
876	tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
877	CALL var_move_att(tl_var, tl_att)
878
879	tl_att=att_init('src_i-indices',(/tl_dom%i_imin, tl_dom%i_imax/))
880	CALL var_move_att(tl_var, tl_att)
881
882	tl_att=att_init('src_j-indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
883	CALL var_move_att(tl_var, tl_att)
884
885	!10- save result
886	create_bathy_get_var=tl_var
887
888	!11- clean structure
889	CALL mpp_clean(tl_mpp)
890
891	END FUNCTION create_bathy_get_var
892	!> @endcode
893	!-------------------------------------------------------------------
894	!> @brief
895	!> This subroutine
896	!>
897	!> @details
898	!>
899	!> @author J.Paul
900	!> - Nov, 2013- Initial Version
901	!>
902	!> @param[in]
903	!> @todo
904	!-------------------------------------------------------------------
905	!> @code
906	SUBROUTINE create_bathy_interp( td_var, &
907	& id_rho, &
908	& id_offset, &
909	& id_iext, id_jext)
910
911	IMPLICIT NONE
912
913	! Argument
914	TYPE(TVAR) , INTENT(INOUT) :: td_var
915	INTEGER(i4), DIMENSION(:) , INTENT(IN ) :: id_rho
916	INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset
917	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext
918	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext
919
920	! local variable
921	TYPE(TVAR) :: tl_var
922	TYPE(TVAR) :: tl_mask
923
924	INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
925
926	INTEGER(i4) :: il_iext
927	INTEGER(i4) :: il_jext
928
929	! loop indices
930	!----------------------------------------------------------------
931
932	! copy variable
933	tl_var=td_var
934
935	!WARNING: two extrabands are required for cubic interpolation
936	il_iext=3
937	IF( PRESENT(id_iext) ) il_iext=id_iext
938
939	il_jext=3
940	IF( PRESENT(id_jext) ) il_jext=id_jext
941
942	IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
943	CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
944	& "on two points are required with cubic interpolation ")
945	il_iext=2
946	ENDIF
947
948	IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
949	CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
950	& "on two points are required with cubic interpolation ")
951	il_jext=2
952	ENDIF
953
954	!1- work on mask
955	!1-1 create mask
956	ALLOCATE(bl_mask(tl_var%t_dim(1)%i_len, &
957	& tl_var%t_dim(2)%i_len, &
958	& tl_var%t_dim(3)%i_len, &
959	& tl_var%t_dim(4)%i_len) )
960
961	bl_mask(:,:,:,:)=1
962	WHERE(tl_var%d_value(:,:,:,:)==tl_var%d_fill) bl_mask(:,:,:,:)=0
963
964	SELECT CASE(TRIM(tl_var%c_point))
965	CASE DEFAULT ! 'T'
966	tl_mask=var_init('tmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))
967	CASE('U')
968	tl_mask=var_init('umask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))
969	CASE('V')
970	tl_mask=var_init('vmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))
971	CASE('F')
972	tl_mask=var_init('fmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))
973	END SELECT
974
975	DEALLOCATE(bl_mask)
976
977	!1-2 interpolate mask
978	CALL interp_fill_value( tl_mask, id_rho(:), &
979	& id_offset=id_offset(:,:) )
980
981	!2- work on variable
982	!2-0 add extraband
983	CALL extrap_add_extrabands(tl_var, il_iext, il_jext)
984
985	!2-1 extrapolate variable
986	CALL extrap_fill_value( tl_var, id_offset=id_offset(:,:), &
987	& id_rho=id_rho(:), &
988	& id_iext=il_iext, id_jext=il_jext )
989
990	!2-2 interpolate Bathymetry
991	CALL interp_fill_value( tl_var, id_rho(:), &
992	& id_offset=id_offset(:,:) )
993
994	!2-3 remove extraband
995	CALL extrap_del_extrabands(tl_var, il_iextid_rho(jp_I), il_jextid_rho(jp_J))
996
997	!2-2-5 keep original mask
998	WHERE( tl_mask%d_value(:,:,:,:) == 0 )
999	tl_var%d_value(:,:,:,:)=tl_var%d_fill
1000	END WHERE
1001
1002	!3- save result
1003	td_var=tl_var
1004
1005	! clean variable structure
1006	CALL var_clean(tl_mask)
1007	CALL var_clean(tl_var)
1008
1009	END SUBROUTINE create_bathy_interp
1010	!> @endcode
1011	END PROGRAM create_bathy

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

Download in other formats: