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/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_bathy.f90 @ 7025

Last change on this file since 7025 was 7025, checked in by jpaul, 8 years ago
see ticket #1781
File size: 40.4 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! PROGRAM: create_bathy
6	!
7	! DESCRIPTION:
8	!> @file
9	!> @brief
10	!> This program creates fine grid bathymetry file.
11	!>
12	!> @details
13	!> @section sec1 method
14	!> Bathymetry could be extracted from fine grid Bathymetry file, interpolated
15	!> from coarse grid Bathymetry file, or manually written.
16	!>
17	!> @section sec2 how to
18	!> to create fine grid bathymetry file:<br/>
19	!> @code{.sh}
20	!> ./SIREN/bin/create_bathy create_bathy.nam
21	!> @endcode
22	!> <br/>
23	!> \image html bathy_40.png
24	!> <center>\image latex bathy_30.png
25	!> </center>
26	!>
27	!> @note
28	!> you could find a template of the namelist in templates directory.
29	!>
30	!> create_bathy.nam contains 7 namelists:<br/>
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	!> - output namelist (namout)
38	!>
39	!> * _logger namelist (namlog)_:<br/>
40	!> - cn_logfile : log filename
41	!> - cn_verbosity : verbosity ('trace','debug','info',
42	!> 'warning','error','fatal','none')
43	!> - in_maxerror : maximum number of error allowed
44	!>
45	!> * _config namelist (namcfg)_:<br/>
46	!> - cn_varcfg : variable configuration file
47	!> (see ./SIREN/cfg/variable.cfg)
48	!> - cn_dumcfg : useless (dummy) configuration file, for useless
49	!> dimension or variable (see ./SIREN/cfg/dummy.cfg).
50	!>
51	!> * _coarse grid namelist (namcrs)_:<br/>
52	!> - cn_coord0 : coordinate file
53	!> - in_perio0 : NEMO periodicity index (see Model Boundary Condition in
54	!> [NEMO documentation](http://www.nemo-ocean.eu/About-NEMO/Reference-manuals))
55	!>
56	!> * _fine grid namelist (namfin)_:<br/>
57	!> - cn_coord1 : coordinate file
58	!> - in_perio1 : periodicity index
59	!> - ln_fillclosed : fill closed sea or not (default is .TRUE.)
60	!>
61	!> * _variable namelist (namvar)_:<br/>
62	!> - cn_varfile : list of variable, and corresponding file.<br/>
63	!> cn_varfile is the path and filename of the file where find
64	!> variable.
65	!> @note
66	!> cn_varfile could be a matrix of value, if you want to filled
67	!> manually variable value.<br/>
68	!> the variable array of value is split into equal subdomain.<br/>
69	!> Each subdomain is filled with the corresponding value
70	!> of the matrix.<br/>
71	!> separators used to defined matrix are:
72	!> - ',' for line
73	!> - '/' for row
74	!> Example:<br/>
75	!> 3,2,3/1,4,5 => @f$ \left( \begin{array}{ccc}
76	!> 3 & 2 & 3 \\
77	!> 1 & 4 & 5 \end{array} \right) @f$
78	!>
79	!> Examples:
80	!> - 'Bathymetry:gridT.nc'
81	!> - 'Bathymetry:5000,5000,5000/5000,3000,5000/5000,5000,5000'
82	!>
83	!> - cn_varinfo : list of variable and extra information about request(s)
84	!> to be used.<br/>
85	!> each elements of cn_varinfo is a string character
86	!> (separated by ',').<br/>
87	!> it is composed of the variable name follow by ':',
88	!> then request(s) to be used on this variable.<br/>
89	!> request could be:
90	!> - int = interpolation method
91	!> - ext = extrapolation method
92	!> - flt = filter method
93	!> - min = minimum value
94	!> - max = maximum value
95	!> - unt = new units
96	!> - unf = unit scale factor (linked to new units)
97	!>
98	!> requests must be separated by ';'.<br/>
99	!> order of requests does not matter.<br/>
100	!>
101	!> informations about available method could be find in @ref interp,
102	!> @ref extrap and @ref filter modules.<br/>
103	!> Example: 'Bathymetry: flt=2*hamming(2,3); min=0'
104	!> @note
105	!> If you do not specify a method which is required,
106	!> default one is apply.
107	!> @warning
108	!> variable name must be __Bathymetry__ here.
109	!>
110	!> * _nesting namelist (namnst)_:<br/>
111	!> - in_rhoi : refinement factor in i-direction
112	!> - in_rhoj : refinement factor in j-direction
113	!> @note
114	!> coarse grid indices will be deduced from fine grid
115	!> coordinate file.
116	!>
117	!> * _output namelist (namout)_:<br/>
118	!> - cn_fileout : output bathymetry file
119	!>
120	!> @author J.Paul
121	! REVISION HISTORY:
122	!> @date November, 2013 - Initial Version
123	!> @date Sepember, 2014
124	!> - add header for user
125	!> - Bug fix, compute offset depending of grid point
126	!> @date June, 2015
127	!> - extrapolate all land points.
128	!> - allow to change unit.
129	!> @date September, 2015
130	!> - manage useless (dummy) variable, attributes, and dimension
131	!> @date January,2016
132	!> - add create_bathy_check_depth as in create_boundary
133	!> - add create_bathy_check_time as in create_boundary
134	!> @date February, 2016
135	!> - do not closed sea for east-west cyclic domain
136	!
137	!> @todo
138	!> - check tl_multi is not empty
139	!>
140	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
141	!----------------------------------------------------------------------
142	PROGRAM create_bathy
143
144	USE global ! global variable
145	USE kind ! F90 kind parameter
146	USE logger ! log file manager
147	USE fct ! basic useful function
148	USE date ! date manager
149	USE att ! attribute manager
150	USE dim ! dimension manager
151	USE var ! variable manager
152	USE file ! file manager
153	USE multi ! multi file manager
154	USE iom ! I/O manager
155	USE grid ! grid manager
156	USE extrap ! extrapolation manager
157	USE interp ! interpolation manager
158	USE filter ! filter manager
159	USE mpp ! MPP manager
160	USE dom ! domain manager
161	USE iom_mpp ! MPP I/O manager
162	USE iom_dom ! DOM I/O manager
163
164	IMPLICIT NONE
165
166	! local variable
167	CHARACTER(LEN=lc) :: cl_namelist
168	CHARACTER(LEN=lc) :: cl_date
169	CHARACTER(LEN=lc) :: cl_data
170
171	INTEGER(i4) :: il_narg
172	INTEGER(i4) :: il_status
173	INTEGER(i4) :: il_fileid
174	INTEGER(i4) :: il_attid
175	INTEGER(i4) :: il_imin0
176	INTEGER(i4) :: il_imax0
177	INTEGER(i4) :: il_jmin0
178	INTEGER(i4) :: il_jmax0
179	INTEGER(i4) , DIMENSION(ip_maxdim) :: il_rho
180	INTEGER(i4) , DIMENSION(2,2) :: il_offset
181	INTEGER(i4) , DIMENSION(2,2) :: il_ind
182	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_mask
183
184	LOGICAL :: ll_exist
185	LOGICAL :: ll_fillclosed
186
187	TYPE(TMPP) :: tl_coord0
188	TYPE(TMPP) :: tl_coord1
189	TYPE(TMPP) :: tl_mpp
190	TYPE(TFILE) :: tl_fileout
191
192	TYPE(TATT) :: tl_att
193
194	TYPE(TVAR) :: tl_lon
195	TYPE(TVAR) :: tl_lat
196	TYPE(TVAR) :: tl_depth
197	TYPE(TVAR) :: tl_time
198
199	TYPE(TVAR) :: tl_tmp
200	TYPE(TVAR) , DIMENSION(:), ALLOCATABLE :: tl_var
201
202	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
203
204	TYPE(TMULTI) :: tl_multi
205
206	REAL(dp) :: dl_minbat
207
208	! loop indices
209	INTEGER(i4) :: ji
210	INTEGER(i4) :: jj
211	INTEGER(i4) :: jk
212
213	! namelist variable
214	! namlog
215	CHARACTER(LEN=lc) :: cn_logfile = 'create_bathy.log'
216	CHARACTER(LEN=lc) :: cn_verbosity = 'warning'
217	INTEGER(i4) :: in_maxerror = 5
218
219	! namcfg
220	CHARACTER(LEN=lc) :: cn_varcfg = './cfg/variable.cfg'
221	CHARACTER(LEN=lc) :: cn_dumcfg = './cfg/dummy.cfg'
222
223	! namcrs
224	CHARACTER(LEN=lc) :: cn_coord0 = ''
225	INTEGER(i4) :: in_perio0 = -1
226
227	! namfin
228	CHARACTER(LEN=lc) :: cn_coord1 = ''
229	INTEGER(i4) :: in_perio1 = -1
230	LOGICAL :: ln_fillclosed = .TRUE.
231
232	! namvar
233	CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = ''
234	CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
235
236	! namnst
237	INTEGER(i4) :: in_rhoi = 1
238	INTEGER(i4) :: in_rhoj = 1
239
240	! namout
241	CHARACTER(LEN=lc) :: cn_fileout = 'bathy_fine.nc'
242	!-------------------------------------------------------------------
243
244	NAMELIST /namlog/ & !< logger namelist
245	& cn_logfile, & !< log file
246	& cn_verbosity, & !< log verbosity
247	& in_maxerror !< logger maximum error
248
249	NAMELIST /namcfg/ & !< configuration namelist
250	& cn_varcfg, & !< variable configuration file
251	& cn_dumcfg !< dummy configuration file
252
253	NAMELIST /namcrs/ & !< coarse grid namelist
254	& cn_coord0, & !< coordinate file
255	& in_perio0 !< periodicity index
256
257	NAMELIST /namfin/ & !< fine grid namelist
258	& cn_coord1, & !< coordinate file
259	& in_perio1, & !< periodicity index
260	& ln_fillclosed !< fill closed sea
261
262	NAMELIST /namvar/ & !< variable namelist
263	& cn_varfile, & !< list of variable file
264	& cn_varinfo !< list of variable and interpolation method to be used. (ex: 'votemper:linear','vosaline:cubic' )
265
266	NAMELIST /namnst/ & !< nesting namelist
267	& in_rhoi, & !< refinement factor in i-direction
268	& in_rhoj !< refinement factor in j-direction
269
270	NAMELIST /namout/ & !< output namlist
271	& cn_fileout !< fine grid bathymetry file
272	!-------------------------------------------------------------------
273
274	! namelist
275	! get namelist
276	il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
277	IF( il_narg/=1 )THEN
278	PRINT *,"ERROR in create_bathy: need a namelist"
279	STOP
280	ELSE
281	CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
282	ENDIF
283
284	! read namelist
285	INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
286	IF( ll_exist )THEN
287
288	il_fileid=fct_getunit()
289
290	OPEN( il_fileid, FILE=TRIM(cl_namelist), &
291	& FORM='FORMATTED', &
292	& ACCESS='SEQUENTIAL', &
293	& STATUS='OLD', &
294	& ACTION='READ', &
295	& IOSTAT=il_status)
296	CALL fct_err(il_status)
297	IF( il_status /= 0 )THEN
298	PRINT *,"ERROR in create_bathy: error opening "//TRIM(cl_namelist)
299	STOP
300	ENDIF
301
302	READ( il_fileid, NML = namlog )
303	! define log file
304	CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
305	CALL logger_header()
306
307	READ( il_fileid, NML = namcfg )
308	! get variable extra information
309	CALL var_def_extra(TRIM(cn_varcfg))
310
311	! get dummy variable
312	CALL var_get_dummy(TRIM(cn_dumcfg))
313	! get dummy dimension
314	CALL dim_get_dummy(TRIM(cn_dumcfg))
315	! get dummy attribute
316	CALL att_get_dummy(TRIM(cn_dumcfg))
317
318	READ( il_fileid, NML = namcrs )
319	READ( il_fileid, NML = namfin )
320	READ( il_fileid, NML = namvar )
321	! add user change in extra information
322	CALL var_chg_extra( cn_varinfo )
323	! match variable with file
324	tl_multi=multi_init(cn_varfile)
325
326	READ( il_fileid, NML = namnst )
327	READ( il_fileid, NML = namout )
328
329	CLOSE( il_fileid, IOSTAT=il_status )
330	CALL fct_err(il_status)
331	IF( il_status /= 0 )THEN
332	CALL logger_error("CREATE BATHY: closing "//TRIM(cl_namelist))
333	ENDIF
334
335	ELSE
336
337	PRINT *,"ERROR in create_bathy: can't find "//TRIM(cl_namelist)
338	STOP
339
340	ENDIF
341
342	CALL multi_print(tl_multi)
343
344	! open files
345	IF( cn_coord0 /= '' )THEN
346	tl_coord0=mpp_init( file_init(TRIM(cn_coord0)), id_perio=in_perio0)
347	CALL grid_get_info(tl_coord0)
348	ELSE
349	CALL logger_fatal("CREATE BATHY: no coarse grid coordinate found. "//&
350	& "check namelist")
351	ENDIF
352
353	IF( TRIM(cn_coord1) /= '' )THEN
354	tl_coord1=mpp_init( file_init(TRIM(cn_coord1)),id_perio=in_perio1)
355	CALL grid_get_info(tl_coord1)
356	ELSE
357	CALL logger_fatal("CREATE BATHY: no fine grid coordinate found. "//&
358	& "check namelist")
359	ENDIF
360
361	! do not closed sea for east-west cyclic domain
362	ll_fillclosed=ln_fillclosed
363	IF( tl_coord1%i_perio == 1 ) ll_fillclosed=.FALSE.
364
365	! check
366	! check output file do not already exist
367	INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
368	IF( ll_exist )THEN
369	CALL logger_fatal("CREATE BATHY: output file "//TRIM(cn_fileout)//&
370	& " already exist.")
371	ENDIF
372
373	! check namelist
374	! check refinement factor
375	il_rho(:)=1
376	IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
377	CALL logger_error("CREATE BATHY: invalid refinement factor."//&
378	& " check namelist "//TRIM(cl_namelist))
379	ELSE
380	il_rho(jp_I)=in_rhoi
381	il_rho(jp_J)=in_rhoj
382	ENDIF
383
384	! check domain indices
385	! compute coarse grid indices around fine grid
386	il_ind(:,:)=grid_get_coarse_index( tl_coord0, tl_coord1, &
387	& id_rho=il_rho(:) )
388
389	il_imin0=il_ind(jp_I,1) ; il_imax0=il_ind(jp_I,2)
390	il_jmin0=il_ind(jp_J,1) ; il_jmax0=il_ind(jp_J,2)
391
392	! check domain validity
393	CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0)
394
395	! check coincidence between coarse and fine grid
396	CALL grid_check_coincidence( tl_coord0, tl_coord1, &
397	& il_imin0, il_imax0, &
398	& il_jmin0, il_jmax0, &
399	& il_rho(:) )
400
401	IF( .NOT. ASSOCIATED(tl_multi%t_mpp) )THEN
402	CALL logger_error("CREATE BATHY: no mpp file to work on. "//&
403	& "check cn_varfile in namelist.")
404	ELSE
405
406	ALLOCATE( tl_var( tl_multi%i_nvar ) )
407	jk=0
408	DO ji=1,tl_multi%i_nmpp
409
410	WRITE(cl_data,'(a,i2.2)') 'data-',jk+1
411	IF( .NOT. ASSOCIATED(tl_multi%t_mpp(ji)%t_proc(1)%t_var) )THEN
412
413	CALL logger_fatal("CREATE BATHY: no variable to work on for "//&
414	& "mpp file"//TRIM(tl_multi%t_mpp(ji)%c_name)//&
415	& ". check cn_varfile in namelist.")
416
417	ELSEIF( TRIM(tl_multi%t_mpp(ji)%c_name) == TRIM(cl_data) )THEN
418
419	!- use input matrix to initialise variable
420	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
421	jk=jk+1
422	tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
423
424	tl_var(jk)=create_bathy_matrix(tl_tmp, tl_coord1)
425	ENDDO
426	! clean
427	CALL var_clean(tl_tmp)
428
429	ELSE
430
431	tl_mpp=mpp_init( file_init(TRIM(tl_multi%t_mpp(ji)%c_name)) )
432	CALL grid_get_info(tl_mpp)
433
434	! open mpp file
435	CALL iom_mpp_open(tl_mpp)
436
437	! get or check depth value
438	CALL create_bathy_check_depth( tl_mpp, tl_depth )
439
440	! get or check time value
441	CALL create_bathy_check_time( tl_mpp, tl_time )
442
443	! close mpp file
444	CALL iom_mpp_close(tl_mpp)
445
446	IF( ANY(tl_mpp%t_dim(1:2)%i_len /= tl_coord0%t_dim(1:2)%i_len).OR.&
447	& ALL(il_rho(:)==1) )THEN
448	!- extract bathymetry from fine grid bathymetry
449	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
450	jk=jk+1
451	tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
452
453	tl_var(jk)=create_bathy_extract( tl_tmp, tl_mpp, &
454	& tl_coord1 )
455	ENDDO
456	! clean
457	CALL var_clean(tl_tmp)
458	ELSE
459	!- get bathymetry from coarse grid bathymetry
460	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
461	jk=jk+1
462	tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
463
464	il_offset(:,:)= grid_get_fine_offset(tl_coord0, &
465	& il_imin0, il_jmin0, &
466	& il_imax0, il_jmax0, &
467	& tl_coord1, &
468	& il_rho(:), &
469	& TRIM(tl_tmp%c_point) )
470
471	tl_var(jk)=create_bathy_get_var( tl_tmp, tl_mpp, &
472	& il_imin0, il_jmin0, &
473	& il_imax0, il_jmax0, &
474	& il_offset(:,:), &
475	& il_rho(:) )
476	ENDDO
477	! clean
478	CALL var_clean(tl_tmp)
479	ENDIF
480
481	! clean structure
482	CALL mpp_clean(tl_mpp)
483
484	ENDIF
485	ENDDO
486	ENDIF
487
488	! use additional request
489	DO jk=1,tl_multi%i_nvar
490
491	! change unit and apply factor
492	CALL var_chg_unit(tl_var(jk))
493
494	! forced min and max value
495	CALL var_limit_value(tl_var(jk))
496
497	! fill closed sea
498	IF( ll_fillclosed )THEN
499	ALLOCATE( il_mask(tl_var(jk)%t_dim(1)%i_len, &
500	& tl_var(jk)%t_dim(2)%i_len) )
501
502	! split domain in N sea subdomain
503	il_mask(:,:)=grid_split_domain(tl_var(jk))
504
505	! fill smallest domain
506	CALL grid_fill_small_dom( tl_var(jk), il_mask(:,:) )
507
508	DEALLOCATE( il_mask )
509	ENDIF
510
511	! filter
512	CALL filter_fill_value(tl_var(jk))
513
514	! check bathymetry
515	dl_minbat=MINVAL(tl_var(jk)%d_value(:,:,:,:))
516	IF( TRIM(tl_var(jk)%c_stdname) == 'bathymetry' .AND. &
517	& dl_minbat <= 0._dp )THEN
518	CALL logger_debug("CREATE BATHY: min value "//TRIM(fct_str(dl_minbat)))
519	CALL logger_fatal("CREATE BATHY: Bathymetry has value <= 0")
520	ENDIF
521
522	ENDDO
523
524	! create file
525	tl_fileout=file_init(TRIM(cn_fileout))
526
527	! add dimension
528	tl_dim(:)=var_max_dim(tl_var(:))
529
530	DO ji=1,ip_maxdim
531	IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
532	ENDDO
533
534	! add variables
535	IF( ALL( tl_dim(1:2)%l_use ) )THEN
536
537	! open mpp files
538	CALL iom_mpp_open(tl_coord1)
539
540	! add longitude
541	tl_lon=iom_mpp_read_var(tl_coord1,'longitude')
542	CALL file_add_var(tl_fileout, tl_lon)
543	CALL var_clean(tl_lon)
544
545	! add latitude
546	tl_lat=iom_mpp_read_var(tl_coord1,'latitude')
547	CALL file_add_var(tl_fileout, tl_lat)
548	CALL var_clean(tl_lat)
549
550	! close mpp files
551	CALL iom_mpp_close(tl_coord1)
552
553	ENDIF
554
555	IF( tl_dim(3)%l_use )THEN
556	! add depth
557	CALL file_add_var(tl_fileout, tl_depth)
558	CALL var_clean(tl_depth)
559	ENDIF
560
561	IF( tl_dim(4)%l_use )THEN
562	! add time
563	CALL file_add_var(tl_fileout, tl_time)
564	CALL var_clean(tl_time)
565	ENDIF
566
567	! add other variables
568	DO jk=tl_multi%i_nvar,1,-1
569	CALL file_add_var(tl_fileout, tl_var(jk))
570	CALL var_clean(tl_var(jk))
571	ENDDO
572	DEALLOCATE(tl_var)
573
574	! add some attribute
575	tl_att=att_init("Created_by","SIREN create_bathy")
576	CALL file_add_att(tl_fileout, tl_att)
577
578	cl_date=date_print(date_now())
579	tl_att=att_init("Creation_date",cl_date)
580	CALL file_add_att(tl_fileout, tl_att)
581
582	! add attribute periodicity
583	il_attid=0
584	IF( ASSOCIATED(tl_fileout%t_att) )THEN
585	il_attid=att_get_index(tl_fileout%t_att(:),'periodicity')
586	ENDIF
587	IF( tl_coord1%i_perio >= 0 .AND. il_attid == 0 )THEN
588	tl_att=att_init('periodicity',tl_coord1%i_perio)
589	CALL file_add_att(tl_fileout,tl_att)
590	ENDIF
591
592	! add attribute east west overlap
593	il_attid=0
594	IF( ASSOCIATED(tl_fileout%t_att) )THEN
595	il_attid=att_get_index(tl_fileout%t_att(:),'ew_overlap')
596	ENDIF
597	IF( tl_coord1%i_ew >= 0 .AND. il_attid == 0 )THEN
598	tl_att=att_init('ew_overlap',tl_coord1%i_ew)
599	CALL file_add_att(tl_fileout,tl_att)
600	ENDIF
601
602	! create file
603	CALL iom_create(tl_fileout)
604
605	! write file
606	CALL iom_write_file(tl_fileout)
607
608	! close file
609	CALL iom_close(tl_fileout)
610
611	! clean
612	CALL att_clean(tl_att)
613
614	CALL file_clean(tl_fileout)
615	CALL mpp_clean(tl_coord1)
616	CALL mpp_clean(tl_coord0)
617
618	! close log file
619	CALL logger_footer()
620	CALL logger_close()
621
622	CONTAINS
623	!-------------------------------------------------------------------
624	!> @brief
625	!> This function create variable, filled with matrix value
626	!>
627	!> @details
628	!> A variable is create with the same name that the input variable,
629	!> and with dimension of the coordinate file.<br/>
630	!> Then the variable array of value is split into equal subdomain.
631	!> Each subdomain is filled with the corresponding value of the matrix.
632	!>
633	!> @author J.Paul
634	!> @date November, 2013 - Initial Version
635	!>
636	!> @param[in] td_var variable structure
637	!> @param[in] td_coord coordinate file structure
638	!> @return variable structure
639	!-------------------------------------------------------------------
640	FUNCTION create_bathy_matrix(td_var, td_coord)
641	IMPLICIT NONE
642	! Argument
643	TYPE(TVAR), INTENT(IN) :: td_var
644	TYPE(TMPP), INTENT(IN) :: td_coord
645
646	! function
647	TYPE(TVAR) :: create_bathy_matrix
648
649	! local variable
650	INTEGER(i4) , DIMENSION(2,2) :: il_xghost
651	INTEGER(i4) , DIMENSION(3) :: il_dim
652	INTEGER(i4) , DIMENSION(3) :: il_size
653	INTEGER(i4) , DIMENSION(3) :: il_rest
654
655	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ishape
656	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_jshape
657	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_kshape
658
659	REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
660
661	TYPE(TVAR) :: tl_lon
662	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
663
664	TYPE(TMPP) :: tl_coord
665
666	! loop indices
667	INTEGER(i4) :: ji
668	INTEGER(i4) :: jj
669	INTEGER(i4) :: jk
670	!----------------------------------------------------------------
671
672	! copy structure
673	tl_coord=mpp_copy(td_coord)
674
675	! use only edge processor
676	CALL mpp_get_contour(tl_coord)
677
678	! open useful processor
679	CALL iom_mpp_open(tl_coord)
680
681	! read output grid
682	tl_lon=iom_mpp_read_var(tl_coord,'longitude')
683
684	! look for ghost cell
685	il_xghost(:,:)=grid_get_ghost( tl_coord )
686
687	! close processor
688	CALL iom_mpp_close(tl_coord)
689	! clean
690	CALL mpp_clean(tl_coord)
691
692	! remove ghost cell
693	CALL grid_del_ghost(tl_lon, il_xghost(:,:))
694
695	! write value on grid
696	! get matrix dimension
697	il_dim(:)=td_var%t_dim(1:3)%i_len
698	! output dimension
699	tl_dim(:)=dim_copy(tl_lon%t_dim(:))
700	! clean
701	CALL var_clean(tl_lon)
702
703	! split output domain in N subdomain depending of matrix dimension
704	il_size(:) = tl_dim(1:3)%i_len / il_dim(:)
705	il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:))
706
707	ALLOCATE( il_ishape(il_dim(1)+1) )
708	il_ishape(:)=0
709	DO ji=2,il_dim(1)+1
710	il_ishape(ji)=il_ishape(ji-1)+il_size(1)
711	ENDDO
712	! add rest to last cell
713	il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1)
714
715	ALLOCATE( il_jshape(il_dim(2)+1) )
716	il_jshape(:)=0
717	DO jj=2,il_dim(2)+1
718	il_jshape(jj)=il_jshape(jj-1)+il_size(2)
719	ENDDO
720	! add rest to last cell
721	il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2)
722
723	ALLOCATE( il_kshape(il_dim(3)+1) )
724	il_kshape(:)=0
725	DO jk=2,il_dim(3)+1
726	il_kshape(jk)=il_kshape(jk-1)+il_size(3)
727	ENDDO
728	! add rest to last cell
729	il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3)
730
731	! write ouput array of value
732	ALLOCATE(dl_value( tl_dim(1)%i_len, &
733	& tl_dim(2)%i_len, &
734	& tl_dim(3)%i_len, &
735	& tl_dim(4)%i_len) )
736
737	dl_value(:,:,:,:)=0
738
739	DO jk=2,il_dim(3)+1
740	DO jj=2,il_dim(2)+1
741	DO ji=2,il_dim(1)+1
742
743	dl_value( 1+il_ishape(ji-1):il_ishape(ji), &
744	& 1+il_jshape(jj-1):il_jshape(jj), &
745	& 1+il_kshape(jk-1):il_kshape(jk), &
746	& 1 ) = td_var%d_value(ji-1,jj-1,jk-1,1)
747
748	ENDDO
749	ENDDO
750	ENDDO
751
752	! initialise variable with value
753	create_bathy_matrix=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:))
754
755	DEALLOCATE(dl_value)
756
757	! add ghost cell
758	CALL grid_add_ghost(create_bathy_matrix, il_xghost(:,:))
759
760	! clean
761	CALL dim_clean(tl_dim(:))
762
763	END FUNCTION create_bathy_matrix
764	!-------------------------------------------------------------------
765	!> @brief
766	!> This function extract variable from file over coordinate domain and
767	!> return variable structure
768	!>
769	!> @author J.Paul
770	!> @date November, 2013 - Initial Version
771	!>
772	!> @param[in] td_var variable structure
773	!> @param[in] td_mpp mpp file structure
774	!> @param[in] td_coord coordinate file structure
775	!> @return variable structure
776	!-------------------------------------------------------------------
777	FUNCTION create_bathy_extract(td_var, td_mpp, &
778	& td_coord)
779	IMPLICIT NONE
780	! Argument
781	TYPE(TVAR), INTENT(IN) :: td_var
782	TYPE(TMPP), INTENT(IN) :: td_mpp
783	TYPE(TMPP), INTENT(IN) :: td_coord
784
785	! function
786	TYPE(TVAR) :: create_bathy_extract
787
788	! local variable
789	INTEGER(i4), DIMENSION(2,2) :: il_ind
790
791	INTEGER(i4) :: il_imin
792	INTEGER(i4) :: il_jmin
793	INTEGER(i4) :: il_imax
794	INTEGER(i4) :: il_jmax
795
796	TYPE(TMPP) :: tl_mpp
797
798	TYPE(TATT) :: tl_att
799
800	TYPE(TVAR) :: tl_var
801
802	TYPE(TDOM) :: tl_dom
803	! loop indices
804	!----------------------------------------------------------------
805
806	IF( .NOT. ASSOCIATED(td_mpp%t_proc) )THEN
807	CALL logger_error("CREATE BATHY EXTRACT: no processor associated "//&
808	& "to mpp "//TRIM(td_mpp%c_name))
809	ELSE
810
811	!init
812	tl_mpp=mpp_copy(td_mpp)
813
814	! compute file grid indices around coord grid
815	il_ind(:,:)=grid_get_coarse_index(tl_mpp, td_coord )
816
817	il_imin=il_ind(1,1) ; il_imax=il_ind(1,2)
818	il_jmin=il_ind(2,1) ; il_jmax=il_ind(2,2)
819
820	! check grid coincidence
821	CALL grid_check_coincidence( tl_mpp, td_coord, &
822	& il_imin, il_imax, &
823	& il_jmin, il_jmax, &
824	& (/1, 1, 1/) )
825
826	! compute domain
827	tl_dom=dom_init(tl_mpp, &
828	& il_imin, il_imax, &
829	& il_jmin, il_jmax)
830
831	! open mpp files over domain
832	CALL iom_dom_open(tl_mpp, tl_dom)
833
834	! read variable on domain
835	tl_var=iom_dom_read_var(tl_mpp,TRIM(td_var%c_name),tl_dom)
836
837	! close mpp file
838	CALL iom_dom_close(tl_mpp)
839
840	! add ghost cell
841	CALL grid_add_ghost(tl_var,tl_dom%i_ghost(:,:))
842
843	! check result
844	IF( ANY( tl_var%t_dim(:)%l_use .AND. &
845	& tl_var%t_dim(:)%i_len /= td_coord%t_dim(:)%i_len) )THEN
846	CALL logger_debug("CREATE BATHY EXTRACT: "//&
847	& "dimensoin of variable "//TRIM(td_var%c_name)//" "//&
848	& TRIM(fct_str(tl_var%t_dim(1)%i_len))//","//&
849	& TRIM(fct_str(tl_var%t_dim(2)%i_len))//","//&
850	& TRIM(fct_str(tl_var%t_dim(3)%i_len))//","//&
851	& TRIM(fct_str(tl_var%t_dim(4)%i_len)) )
852	CALL logger_debug("CREATE BATHY EXTRACT: "//&
853	& "dimensoin of coordinate file "//&
854	& TRIM(fct_str(td_coord%t_dim(1)%i_len))//","//&
855	& TRIM(fct_str(td_coord%t_dim(2)%i_len))//","//&
856	& TRIM(fct_str(td_coord%t_dim(3)%i_len))//","//&
857	& TRIM(fct_str(td_coord%t_dim(4)%i_len)) )
858	CALL logger_fatal("CREATE BATHY EXTRACT: "//&
859	& "dimensoin of extracted "//&
860	& "variable and coordinate file dimension differ")
861	ENDIF
862
863	! add attribute to variable
864	tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
865	CALL var_move_att(tl_var, tl_att)
866
867	tl_att=att_init('src_i_indices',(/tl_dom%i_imin, tl_dom%i_imax/))
868	CALL var_move_att(tl_var, tl_att)
869
870	tl_att=att_init('src_j_indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
871	CALL var_move_att(tl_var, tl_att)
872
873	! save result
874	create_bathy_extract=var_copy(tl_var)
875
876	! clean structure
877	CALL att_clean(tl_att)
878	CALL var_clean(tl_var)
879	CALL mpp_clean(tl_mpp)
880	ENDIF
881
882	END FUNCTION create_bathy_extract
883	!-------------------------------------------------------------------
884	!> @brief
885	!> This function get coarse grid variable, interpolate variable, and return
886	!> variable structure over fine grid
887	!>
888	!> @author J.Paul
889	!> @date November, 2013 - Initial Version
890	!>
891	!> @param[in] td_var variable structure
892	!> @param[in] td_mpp mpp file structure
893	!> @param[in] id_imin i-direction lower left corner indice
894	!> @param[in] id_imax i-direction upper right corner indice
895	!> @param[in] id_jmin j-direction lower left corner indice
896	!> @param[in] id_jmax j-direction upper right corner indice
897	!> @param[in] id_offset offset between fine grid and coarse grid
898	!> @param[in] id_rho array of refinement factor
899	!> @return variable structure
900	!-------------------------------------------------------------------
901	FUNCTION create_bathy_get_var(td_var, td_mpp, &
902	& id_imin, id_jmin, &
903	& id_imax, id_jmax, &
904	& id_offset, &
905	& id_rho )
906	IMPLICIT NONE
907	! Argument
908	TYPE(TVAR) , INTENT(IN) :: td_var
909	TYPE(TMPP) , INTENT(IN) :: td_mpp
910	INTEGER(i4) , INTENT(IN) :: id_imin
911	INTEGER(i4) , INTENT(IN) :: id_imax
912	INTEGER(i4) , INTENT(IN) :: id_jmin
913	INTEGER(i4) , INTENT(IN) :: id_jmax
914	INTEGER(i4), DIMENSION(:,:), INTENT(IN) :: id_offset
915	INTEGER(i4), DIMENSION(:) , INTENT(IN) :: id_rho
916
917	! function
918	TYPE(TVAR) :: create_bathy_get_var
919
920	! local variable
921	TYPE(TMPP) :: tl_mpp
922	TYPE(TATT) :: tl_att
923	TYPE(TVAR) :: tl_var
924	TYPE(TDOM) :: tl_dom
925
926	INTEGER(i4) :: il_size
927	INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
928
929	! loop indices
930	!----------------------------------------------------------------
931	IF( ANY(SHAPE(id_offset(:,:)) /= 2) )THEN
932	CALL logger_error("CREATE BATHY GET VAR: invalid dimension of "//&
933	& "offset array")
934	ENDIF
935
936	! copy structure
937	tl_mpp=mpp_copy(td_mpp)
938
939	!- compute domain
940	tl_dom=dom_init(tl_mpp, &
941	& id_imin, id_imax, &
942	& id_jmin, id_jmax)
943
944	!- add extra band (if possible) to compute interpolation
945	CALL dom_add_extra(tl_dom)
946
947	!- open mpp files over domain
948	CALL iom_dom_open(tl_mpp, tl_dom)
949
950	!- read variable value on domain
951	tl_var=iom_dom_read_var(tl_mpp,TRIM(td_var%c_name),tl_dom)
952
953	!- close mpp files
954	CALL iom_dom_close(tl_mpp)
955
956	il_size=SIZE(id_rho(:))
957	ALLOCATE( il_rho(il_size) )
958	il_rho(:)=id_rho(:)
959
960	!- interpolate variable
961	CALL create_bathy_interp(tl_var, il_rho(:), id_offset(:,:))
962
963	!- remove extraband added to domain
964	CALL dom_del_extra( tl_var, tl_dom, il_rho(:) )
965
966	CALL dom_clean_extra( tl_dom )
967
968	!- add ghost cell
969	CALL grid_add_ghost(tl_var,tl_dom%i_ghost(:,:))
970
971	!- add attribute to variable
972	tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
973	CALL var_move_att(tl_var, tl_att)
974
975	tl_att=att_init('src_i_indices',(/tl_dom%i_imin, tl_dom%i_imax/))
976	CALL var_move_att(tl_var, tl_att)
977
978	tl_att=att_init('src_j_indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
979	CALL var_move_att(tl_var, tl_att)
980
981	IF( .NOT. ALL(id_rho(:)==1) )THEN
982	tl_att=att_init("refinment_factor",(/id_rho(jp_I),id_rho(jp_J)/))
983	CALL var_move_att(tl_var, tl_att)
984	ENDIF
985
986	DEALLOCATE( il_rho )
987
988	!- save result
989	create_bathy_get_var=var_copy(tl_var)
990
991	!- clean structure
992	CALL att_clean(tl_att)
993	CALL var_clean(tl_var)
994	CALL mpp_clean(tl_mpp)
995
996	END FUNCTION create_bathy_get_var
997	!-------------------------------------------------------------------
998	!> @brief
999	!> This subroutine interpolate variable
1000	!>
1001	!> @author J.Paul
1002	!> @date November, 2013 - Initial Version
1003	!>
1004	!> @param[inout] td_var variable structure
1005	!> @param[in] id_rho array of refinment factor
1006	!> @param[in] id_offset array of offset between fine and coarse grid
1007	!> @param[in] id_iext i-direction size of extra bands (default=im_minext)
1008	!> @param[in] id_jext j-direction size of extra bands (default=im_minext)
1009	!-------------------------------------------------------------------
1010	SUBROUTINE create_bathy_interp( td_var, &
1011	& id_rho, &
1012	& id_offset, &
1013	& id_iext, id_jext)
1014
1015	IMPLICIT NONE
1016
1017	! Argument
1018	TYPE(TVAR) , INTENT(INOUT) :: td_var
1019	INTEGER(i4), DIMENSION(:) , INTENT(IN ) :: id_rho
1020	INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset
1021	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext
1022	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext
1023
1024	! local variable
1025	TYPE(TVAR) :: tl_mask
1026
1027	INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
1028
1029	INTEGER(i4) :: il_iext
1030	INTEGER(i4) :: il_jext
1031
1032	! loop indices
1033	!----------------------------------------------------------------
1034
1035	!WARNING: two extrabands are required for cubic interpolation
1036	il_iext=3
1037	IF( PRESENT(id_iext) ) il_iext=id_iext
1038
1039	il_jext=3
1040	IF( PRESENT(id_jext) ) il_jext=id_jext
1041
1042	IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1043	CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
1044	& "on two points are required with cubic interpolation ")
1045	il_iext=2
1046	ENDIF
1047
1048	IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1049	CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
1050	& "on two points are required with cubic interpolation ")
1051	il_jext=2
1052	ENDIF
1053
1054	! work on mask
1055	! create mask
1056	ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
1057	& td_var%t_dim(2)%i_len, &
1058	& td_var%t_dim(3)%i_len, &
1059	& td_var%t_dim(4)%i_len) )
1060
1061	bl_mask(:,:,:,:)=1
1062	WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0
1063
1064	SELECT CASE(TRIM(td_var%c_point))
1065	CASE DEFAULT ! 'T'
1066	tl_mask=var_init('tmask', bl_mask(:,:,:,:), td_dim=td_var%t_dim(:), &
1067	& id_ew=td_var%i_ew )
1068	CASE('U','V','F')
1069	CALL logger_fatal("CREATE BATHY INTERP: can not computed "//&
1070	& "interpolation on "//TRIM(td_var%c_point)//&
1071	& " grid point (variable "//TRIM(td_var%c_name)//&
1072	& "). check namelist.")
1073	END SELECT
1074
1075	DEALLOCATE(bl_mask)
1076
1077	! interpolate mask
1078	CALL interp_fill_value( tl_mask, id_rho(:), &
1079	& id_offset=id_offset(:,:) )
1080
1081	! work on variable
1082	! add extraband
1083	CALL extrap_add_extrabands(td_var, il_iext, il_jext)
1084
1085	! extrapolate variable
1086	CALL extrap_fill_value( td_var )
1087
1088	! interpolate Bathymetry
1089	CALL interp_fill_value( td_var, id_rho(:), &
1090	& id_offset=id_offset(:,:) )
1091
1092	! remove extraband
1093	CALL extrap_del_extrabands(td_var, il_iextid_rho(jp_I), il_jextid_rho(jp_J))
1094
1095	! keep original mask
1096	WHERE( tl_mask%d_value(:,:,:,:) == 0 )
1097	td_var%d_value(:,:,:,:)=td_var%d_fill
1098	END WHERE
1099
1100	! clean variable structure
1101	CALL var_clean(tl_mask)
1102
1103	END SUBROUTINE create_bathy_interp
1104	!-------------------------------------------------------------------
1105	!> @brief
1106	!> This subroutine get depth variable value in an open mpp structure
1107	!> and check if agree with already input depth variable.
1108	!>
1109	!> @details
1110	!>
1111	!> @author J.Paul
1112	!> @date January, 2016 - Initial Version
1113	!>
1114	!> @param[in] td_mpp mpp structure
1115	!> @param[inout] td_depth depth variable structure
1116	!-------------------------------------------------------------------
1117	SUBROUTINE create_bathy_check_depth( td_mpp, td_depth )
1118
1119	IMPLICIT NONE
1120
1121	! Argument
1122	TYPE(TMPP) , INTENT(IN ) :: td_mpp
1123	TYPE(TVAR) , INTENT(INOUT) :: td_depth
1124
1125	! local variable
1126	INTEGER(i4) :: il_varid
1127	TYPE(TVAR) :: tl_depth
1128	! loop indices
1129	!----------------------------------------------------------------
1130
1131	! get or check depth value
1132	IF( td_mpp%t_proc(1)%i_depthid /= 0 )THEN
1133
1134	il_varid=td_mpp%t_proc(1)%i_depthid
1135	IF( ASSOCIATED(td_depth%d_value) )THEN
1136
1137	tl_depth=iom_mpp_read_var(td_mpp, il_varid)
1138
1139	IF( ANY( td_depth%d_value(:,:,:,:) /= &
1140	& tl_depth%d_value(:,:,:,:) ) )THEN
1141
1142	CALL logger_warn("CREATE BATHY: depth value from "//&
1143	& TRIM(td_mpp%c_name)//" not conform "//&
1144	& " to those from former file(s).")
1145
1146	ENDIF
1147	CALL var_clean(tl_depth)
1148
1149	ELSE
1150	td_depth=iom_mpp_read_var(td_mpp,il_varid)
1151	ENDIF
1152
1153	ENDIF
1154
1155	END SUBROUTINE create_bathy_check_depth
1156	!-------------------------------------------------------------------
1157	!> @brief
1158	!> This subroutine get date and time in an open mpp structure
1159	!> and check if agree with date and time already read.
1160	!>
1161	!> @details
1162	!>
1163	!> @author J.Paul
1164	!> @date January, 2016 - Initial Version
1165	!>
1166	!> @param[in] td_mpp mpp structure
1167	!> @param[inout] td_time time variable structure
1168	!-------------------------------------------------------------------
1169	SUBROUTINE create_bathy_check_time( td_mpp, td_time )
1170
1171	IMPLICIT NONE
1172
1173	! Argument
1174	TYPE(TMPP), INTENT(IN ) :: td_mpp
1175	TYPE(TVAR), INTENT(INOUT) :: td_time
1176
1177	! local variable
1178	INTEGER(i4) :: il_varid
1179	TYPE(TVAR) :: tl_time
1180
1181	TYPE(TDATE) :: tl_date1
1182	TYPE(TDATE) :: tl_date2
1183	! loop indices
1184	!----------------------------------------------------------------
1185
1186	! get or check depth value
1187	IF( td_mpp%t_proc(1)%i_timeid /= 0 )THEN
1188
1189	il_varid=td_mpp%t_proc(1)%i_timeid
1190	IF( ASSOCIATED(td_time%d_value) )THEN
1191
1192	tl_time=iom_mpp_read_var(td_mpp, il_varid)
1193
1194	tl_date1=var_to_date(td_time)
1195	tl_date2=var_to_date(tl_time)
1196	IF( tl_date1 - tl_date2 /= 0 )THEN
1197
1198	CALL logger_warn("CREATE BATHY: date from "//&
1199	& TRIM(td_mpp%c_name)//" not conform "//&
1200	& " to those from former file(s).")
1201
1202	ENDIF
1203	CALL var_clean(tl_time)
1204
1205	ELSE
1206	td_time=iom_mpp_read_var(td_mpp,il_varid)
1207	ENDIF
1208
1209	ENDIF
1210
1211	END SUBROUTINE create_bathy_check_time
1212	END PROGRAM create_bathy

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

Download in other formats: