Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

create_boundary.f90 @ 9987

Last change on this file since 9987 was 9987, checked in by emmafiedler, 6 years ago
Merge with GO6 FOAMv14 package branch r9288
File size: 61.5 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	!
6	! PROGRAM: create_boundary
7	!
8	! DESCRIPTION:
9	!> @file
10	!> @brief
11	!> This program create boundary files.
12	!>
13	!> @details
14	!> @section sec1 method
15	!> Variables are read from coarse grid standard output
16	!> and interpolated on fine grid or manually written.<br/>
17	!> @note
18	!> method could be different for each variable.
19	!>
20	!> @section sec2 how to
21	!> to create boundaries files:<br/>
22	!> @code{.sh}
23	!> ./SIREN/bin/create_boundary create_boundary.nam
24	!> @endcode
25	!> <br/>
26	!> \image html boundary_NEATL36_70.png
27	!> \image latex boundary_NEATL36_70.png
28	!>
29	!> @note
30	!> you could find a template of the namelist in templates directory.
31	!>
32	!> create_boundary.nam comprise 9 namelists:<br/>
33	!> - logger namelist (namlog)
34	!> - config namelist (namcfg)
35	!> - coarse grid namelist (namcrs)
36	!> - fine grid namelist (namfin)
37	!> - variable namelist (namvar)
38	!> - nesting namelist (namnst)
39	!> - boundary namelist (nambdy)
40	!> - vertical grid namelist (namzgr)
41	!> - output namelist (namout)
42	!>
43	!> @note
44	!> All namelists have to be in file create_boundary.nam,
45	!> however variables of those namelists are all optional.
46	!>
47	!> * _logger namelist (namlog)_:<br/>
48	!> - cn_logfile : log filename
49	!> - cn_verbosity : verbosity ('trace','debug','info',
50	!> 'warning','error','fatal','none')
51	!> - in_maxerror : maximum number of error allowed
52	!>
53	!> * _config namelist (namcfg)_:<br/>
54	!> - cn_varcfg : variable configuration file
55	!> (see ./SIREN/cfg/variable.cfg)
56	!>
57	!> * _coarse grid namelist (namcrs)_:<br/>
58	!> - cn_coord0 : coordinate file
59	!> - in_perio0 : NEMO periodicity index (see Model Boundary Condition in
60	!> [NEMO documentation](http://www.nemo-ocean.eu/About-NEMO/Reference-manuals))
61	!>
62	!> * _fine grid namelist (namfin)_:<br/>
63	!> - cn_coord1 : coordinate file
64	!> - cn_bathy1 : bathymetry file
65	!> - in_perio1 : periodicity index
66	!>
67	!> * _vertical grid namelist (namzgr)_:<br/>
68	!> - dn_pp_to_be_computed :
69	!> - dn_ppsur :
70	!> - dn_ppa0 :
71	!> - dn_ppa1 :
72	!> - dn_ppa2 :
73	!> - dn_ppkth :
74	!> - dn_ppkth2 :
75	!> - dn_ppacr :
76	!> - dn_ppacr2 :
77	!> - dn_ppdzmin :
78	!> - dn_pphmax :
79	!> - in_nlevel : number of vertical level
80	!>
81	!> * _partial step namelist (namzps)_:<br/>
82	!> - dn_e3zps_mi :
83	!> - dn_e3zps_rat :
84	!>
85	!> * _variable namelist (namvar)_:<br/>
86	!> - cn_varinfo : list of variable and extra information about request(s)
87	!> to be used (separated by ',').<br/>
88	!> each elements of cn_varinfo is a string character.<br/>
89	!> it is composed of the variable name follow by ':',
90	!> then request(s) to be used on this variable.<br/>
91	!> request could be:
92	!> - int = interpolation method
93	!> - ext = extrapolation method
94	!> - flt = filter method
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.
100	!>
101	!> informations about available method could be find in @ref interp,
102	!> @ref extrap and @ref filter.<br/>
103	!>
104	!> Example: 'votemper:int=linear;flt=hann;ext=dist_weight', 'vosaline:int=cubic'
105	!> @note
106	!> If you do not specify a method which is required,
107	!> default one is apply.
108	!> - cn_varfile : list of variable, and corresponding file<br/>
109	!> cn_varfile is the path and filename of the file where find
110	!> variable.<br/>
111	!> @note
112	!> cn_varfile could be a matrix of value, if you want to filled
113	!> manually variable value.<br/>
114	!> the variable array of value is split into equal subdomain.<br/>
115	!> Each subdomain is filled with the corresponding value
116	!> of the matrix.<br/>
117	!> separators used to defined matrix are:
118	!> - ',' for line
119	!> - '/' for row
120	!> - '\' for level<br/>
121	!> Example:<br/>
122	!> 3,2,3/1,4,5 => @f$ \left( \begin{array}{ccc}
123	!> 3 & 2 & 3 \\
124	!> 1 & 4 & 5 \end{array} \right) @f$
125	!> @warning
126	!> the same matrix is used for all boundaries.
127	!>
128	!> Examples:
129	!> - 'votemper:gridT.nc', 'vozocrtx:gridU.nc'
130	!> - 'votemper:10\25', 'vozocrtx:gridU.nc'
131	!>
132	!> * _nesting namelist (namnst)_:<br/>
133	!> - in_rhoi : refinement factor in i-direction
134	!> - in_rhoj : refinement factor in j-direction
135	!>
136	!> * _boundary namelist (nambdy)_:<br/>
137	!> - ln_north : use north boundary
138	!> - ln_south : use south boundary
139	!> - ln_east : use east boundary
140	!> - ln_west : use west boundary
141	!> - cn_north : north boundary indices on fine grid
142	!> cn_north is a string character defining boundary
143	!> segmentation.<br/>
144	!> segments are separated by '\|'.<br/>
145	!> each segments of the boundary is composed of:
146	!> - indice of velocity (orthogonal to boundary .ie.
147	!> for north boundary, J-indice).
148	!> - indice of segemnt start (I-indice for north boundary)
149	!> - indice of segment end (I-indice for north boundary)<br/>
150	!> indices must be separated by ':' .<br/>
151	!> - optionally, boundary size could be added between '(' and ')'
152	!> in the first segment defined.
153	!> @note
154	!> boundary width is the same for all segments of one boundary.
155	!>
156	!> Examples:
157	!> - cn_north='index1,first1:last1(width)'
158	!> - cn_north='index1(width),first1:last1\|index2,first2:last2'
159	!> \image html boundary_50.png
160	!> \image latex boundary_50.png
161	!> - cn_south : south boundary indices on fine grid
162	!> - cn_east : east boundary indices on fine grid
163	!> - cn_west : west boundary indices on fine grid
164	!> - ln_oneseg : use only one segment for each boundary or not
165	!>
166	!> * _output namelist (namout)_:<br/>
167	!> - cn_fileout : fine grid boundary basename
168	!> (cardinal and segment number will be automatically added)
169	!> - dn_dayofs : date offset in day (change only ouput file name)
170	!> - ln_extrap : extrapolate land point or not
171	!>
172	!> Examples:
173	!> - cn_fileout=boundary.nc<br/>
174	!> if time_counter (16/07/2015 00h) is read on input file (see varfile),
175	!> west boundary will be named boundary_west_y2015m07d16
176	!> - dn_dayofs=-2.<br/>
177	!> if you use day offset you get boundary_west_y2015m07d14
178	!>
179	!>
180	!> @author J.Paul
181	! REVISION HISTORY:
182	!> @date November, 2013 - Initial Version
183	!> @date September, 2014
184	!> - add header for user
185	!> - take into account grid point to compue boundaries
186	!> - reorder output dimension for north and south boundaries
187	!> @date June, 2015
188	!> - extrapolate all land points, and add ln_extrap in namelist.
189	!> - allow to change unit.
190	!> @date July, 2015
191	!> - add namelist parameter to shift date of output file name.
192	!>
193	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
194	!----------------------------------------------------------------------
195	PROGRAM create_boundary
196
197	USE netcdf ! nf90 library
198	USE global ! global variable
199	USE phycst ! physical constant
200	USE kind ! F90 kind parameter
201	USE logger ! log file manager
202	USE fct ! basic useful function
203	USE date ! date manager
204	USE att ! attribute manager
205	USE dim ! dimension manager
206	USE var ! variable manager
207	USE file ! file manager
208	USE multi ! multi file manager
209	USE boundary ! boundary manager
210	USE iom ! I/O manager
211	USE dom ! domain manager
212	USE grid ! grid manager
213	USE vgrid ! vertical grid manager
214	USE extrap ! extrapolation manager
215	USE interp ! interpolation manager
216	USE filter ! filter manager
217	USE mpp ! MPP manager
218	USE iom_mpp ! MPP I/O manager
219
220	IMPLICIT NONE
221
222	! local variable
223	CHARACTER(LEN=lc) :: cl_namelist
224	CHARACTER(LEN=lc) :: cl_date
225	CHARACTER(LEN=lc) :: cl_name
226	CHARACTER(LEN=lc) :: cl_bdyout
227	CHARACTER(LEN=lc) :: cl_data
228	CHARACTER(LEN=lc) :: cl_dimorder
229	CHARACTER(LEN=lc) :: cl_point
230	CHARACTER(LEN=lc) :: cl_fmt
231
232	INTEGER(i4) :: il_narg
233	INTEGER(i4) :: il_status
234	INTEGER(i4) :: il_fileid
235	INTEGER(i4) :: il_imin0
236	INTEGER(i4) :: il_imax0
237	INTEGER(i4) :: il_jmin0
238	INTEGER(i4) :: il_jmax0
239	INTEGER(i4) :: il_shift
240	INTEGER(i4) , DIMENSION(ip_maxdim) :: il_rho
241	INTEGER(i4) , DIMENSION(2,2) :: il_offset
242	INTEGER(i4) , DIMENSION(2,2) :: il_ind
243
244	LOGICAL :: ll_exist
245
246	TYPE(TATT) :: tl_att
247
248	TYPE(TVAR) :: tl_depth
249	TYPE(TVAR) :: tl_time
250	TYPE(TVAR) :: tl_var1
251	TYPE(TVAR) :: tl_var0
252	TYPE(TVAR) :: tl_lon1
253	TYPE(TVAR) :: tl_lat1
254	TYPE(TVAR) :: tl_lvl1
255	TYPE(TVAR) , DIMENSION(:) , ALLOCATABLE :: tl_level
256	TYPE(TVAR) , DIMENSION(:,:,:) , ALLOCATABLE :: tl_seglvl1
257	TYPE(TVAR) , DIMENSION(:,:,:) , ALLOCATABLE :: tl_segvar1
258
259	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
260
261	TYPE(TDATE) :: tl_date
262
263	TYPE(TBDY) , DIMENSION(ip_ncard) :: tl_bdy
264
265	TYPE(TDOM) :: tl_dom0
266	TYPE(TDOM) :: tl_dom1
267	TYPE(TDOM) , DIMENSION(:,:,:) , ALLOCATABLE :: tl_segdom1
268
269	TYPE(TFILE) :: tl_fileout
270
271	TYPE(TMPP) :: tl_coord0
272	TYPE(TMPP) :: tl_coord1
273	TYPE(TMPP) :: tl_bathy1
274	TYPE(TMPP) :: tl_mpp
275
276	TYPE(TMULTI) :: tl_multi
277
278	! loop indices
279	INTEGER(i4) :: jvar
280	INTEGER(i4) :: jpoint
281	INTEGER(i4) :: ji
282	INTEGER(i4) :: jj
283	INTEGER(i4) :: jk
284	INTEGER(i4) :: jl
285
286	! namelist variable
287	! namlog
288	CHARACTER(LEN=lc) :: cn_logfile = 'create_boundary.log'
289	CHARACTER(LEN=lc) :: cn_verbosity = 'warning'
290	INTEGER(i4) :: in_maxerror = 5
291
292	! namcfg
293	CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg'
294
295	! namcrs
296	CHARACTER(LEN=lc) :: cn_coord0 = ''
297	INTEGER(i4) :: in_perio0 = -1
298
299	! namfin
300	CHARACTER(LEN=lc) :: cn_coord1 = ''
301	CHARACTER(LEN=lc) :: cn_bathy1 = ''
302	INTEGER(i4) :: in_perio1 = -1
303
304	!namzgr
305	REAL(dp) :: dn_pp_to_be_computed = 0._dp
306	REAL(dp) :: dn_ppsur = -3958.951371276829_dp
307	REAL(dp) :: dn_ppa0 = 103.9530096000000_dp
308	REAL(dp) :: dn_ppa1 = 2.4159512690000_dp
309	REAL(dp) :: dn_ppa2 = 100.7609285000000_dp
310	REAL(dp) :: dn_ppkth = 15.3510137000000_dp
311	REAL(dp) :: dn_ppkth2 = 48.0298937200000_dp
312	REAL(dp) :: dn_ppacr = 7.0000000000000_dp
313	REAL(dp) :: dn_ppacr2 = 13.000000000000_dp
314	REAL(dp) :: dn_ppdzmin = 6._dp
315	REAL(dp) :: dn_pphmax = 5750._dp
316	INTEGER(i4) :: in_nlevel = 75
317
318	!namzps
319	REAL(dp) :: dn_e3zps_min = 25._dp
320	REAL(dp) :: dn_e3zps_rat = 0.2_dp
321
322	! namvar
323	CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
324	CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = ''
325
326	! namnst
327	INTEGER(i4) :: in_rhoi = 0
328	INTEGER(i4) :: in_rhoj = 0
329
330	! nambdy
331	LOGICAL :: ln_north = .TRUE.
332	LOGICAL :: ln_south = .TRUE.
333	LOGICAL :: ln_east = .TRUE.
334	LOGICAL :: ln_west = .TRUE.
335	CHARACTER(LEN=lc) :: cn_north = ''
336	CHARACTER(LEN=lc) :: cn_south = ''
337	CHARACTER(LEN=lc) :: cn_east = ''
338	CHARACTER(LEN=lc) :: cn_west = ''
339	LOGICAL :: ln_oneseg = .TRUE.
340
341	! namout
342	CHARACTER(LEN=lc) :: cn_fileout = 'boundary.nc'
343	REAL(dp) :: dn_dayofs = 0._dp
344	LOGICAL :: ln_extrap = .FALSE.
345	!-------------------------------------------------------------------
346
347	NAMELIST /namlog/ & !< logger namelist
348	& cn_logfile, & !< log file
349	& cn_verbosity, & !< log verbosity
350	& in_maxerror
351
352	NAMELIST /namcfg/ & !< config namelist
353	& cn_varcfg !< variable configuration file
354
355	NAMELIST /namcrs/ & !< coarse grid namelist
356	& cn_coord0, & !< coordinate file
357	& in_perio0 !< periodicity index
358
359	NAMELIST /namfin/ & !< fine grid namelist
360	& cn_coord1, & !< coordinate file
361	& cn_bathy1, & !< bathymetry file
362	& in_perio1 !< periodicity index
363
364	NAMELIST /namzgr/ &
365	& dn_pp_to_be_computed, &
366	& dn_ppsur, &
367	& dn_ppa0, &
368	& dn_ppa1, &
369	& dn_ppa2, &
370	& dn_ppkth, &
371	& dn_ppkth2, &
372	& dn_ppacr, &
373	& dn_ppacr2, &
374	& dn_ppdzmin, &
375	& dn_pphmax, &
376	& in_nlevel !< number of vertical level
377
378	NAMELIST /namzps/ &
379	& dn_e3zps_min, &
380	& dn_e3zps_rat
381
382	NAMELIST /namvar/ & !< variable namelist
383	& cn_varinfo, & !< list of variable and method to apply on. (ex: 'votemper:linear','vosaline:cubic' )
384	& cn_varfile !< list of variable and file where find it. (ex: 'votemper:GLORYS_gridT.nc' )
385
386	NAMELIST /namnst/ & !< nesting namelist
387	& in_rhoi, & !< refinement factor in i-direction
388	& in_rhoj !< refinement factor in j-direction
389
390	NAMELIST /nambdy/ & !< boundary namelist
391	& ln_north, & !< use north boundary
392	& ln_south, & !< use south boundary
393	& ln_east , & !< use east boundary
394	& ln_west , & !< use west boundary
395	& cn_north, & !< north boundary indices on fine grid
396	& cn_south, & !< south boundary indices on fine grid
397	& cn_east , & !< east boundary indices on fine grid
398	& cn_west , & !< west boundary indices on fine grid
399	& ln_oneseg !< use only one segment for each boundary or not
400
401	NAMELIST /namout/ & !< output namelist
402	& cn_fileout, & !< fine grid boundary file basename
403	& dn_dayofs, & !< date offset in day (change only ouput file name)
404	& ln_extrap !< extrapolate or not
405	!-------------------------------------------------------------------
406
407	! namelist
408	! get namelist
409	il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
410	IF( il_narg/=1 )THEN
411	PRINT *,"CREATE BOUNDARY: ERROR. need a namelist"
412	STOP
413	ELSE
414	CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
415	ENDIF
416
417	! read namelist
418	INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
419	IF( ll_exist )THEN
420
421	il_fileid=fct_getunit()
422
423	OPEN( il_fileid, FILE=TRIM(cl_namelist), &
424	& FORM='FORMATTED', &
425	& ACCESS='SEQUENTIAL', &
426	& STATUS='OLD', &
427	& ACTION='READ', &
428	& IOSTAT=il_status)
429	CALL fct_err(il_status)
430	IF( il_status /= 0 )THEN
431	PRINT *,"CREATE BOUNDARY: ERROR opening "//TRIM(cl_namelist)
432	STOP
433	ENDIF
434
435	READ( il_fileid, NML = namlog )
436	! define log file
437	CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
438	CALL logger_header()
439
440	READ( il_fileid, NML = namcfg )
441	! get variable extra information
442	CALL var_def_extra(TRIM(cn_varcfg))
443
444	READ( il_fileid, NML = namcrs )
445	READ( il_fileid, NML = namfin )
446	READ( il_fileid, NML = namzgr )
447	READ( il_fileid, NML = namvar )
448	! add user change in extra information
449	CALL var_chg_extra(cn_varinfo)
450	! match variable with file
451	tl_multi=multi_init(cn_varfile)
452
453	READ( il_fileid, NML = namnst )
454	READ( il_fileid, NML = nambdy )
455	READ( il_fileid, NML = namout )
456
457	CLOSE( il_fileid, IOSTAT=il_status )
458	CALL fct_err(il_status)
459	IF( il_status /= 0 )THEN
460	CALL logger_error("CREATE BOUNDARY: ERROR closing "//TRIM(cl_namelist))
461	ENDIF
462
463	ELSE
464
465	PRINT *,"CREATE BOUNDARY: ERROR. can not find "//TRIM(cl_namelist)
466	STOP
467
468	ENDIF
469
470	CALL multi_print(tl_multi)
471	IF( tl_multi%i_nvar <= 0 )THEN
472	CALL logger_fatal("CREATE BOUNDARY: no variable to be used."//&
473	& " check namelist.")
474	ENDIF
475
476	! open files
477	IF( TRIM(cn_coord0) /= '' )THEN
478	tl_coord0=mpp_init( file_init(TRIM(cn_coord0)), id_perio=in_perio0)
479	CALL grid_get_info(tl_coord0)
480	ELSE
481	CALL logger_fatal("CREATE BOUNDARY: can not find coarse grid "//&
482	& "coordinate file. check namelist")
483	ENDIF
484
485	IF( TRIM(cn_coord1) /= '' )THEN
486	tl_coord1=mpp_init( file_init(TRIM(cn_coord1)), id_perio=in_perio1)
487	CALL grid_get_info(tl_coord1)
488	ELSE
489	CALL logger_fatal("CREATE BOUNDARY: can not find fine grid coordinate "//&
490	& "file. check namelist")
491	ENDIF
492
493	IF( TRIM(cn_bathy1) /= '' )THEN
494	tl_bathy1=mpp_init( file_init(TRIM(cn_bathy1)), id_perio=in_perio1)
495	CALL grid_get_info(tl_bathy1)
496	ELSE
497	CALL logger_fatal("CREATE BOUNDARY: can not find fine grid bathymetry "//&
498	& "file. check namelist")
499	ENDIF
500
501	! check
502	! check output file do not already exist
503	! WARNING: do not work when use time to create output file name
504	DO jk=1,ip_ncard
505	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
506	& TRIM(cp_card(jk)), 1 )
507	INQUIRE(FILE=TRIM(cl_bdyout), EXIST=ll_exist)
508	IF( ll_exist )THEN
509	CALL logger_fatal("CREATE BOUNDARY: output file "//TRIM(cl_bdyout)//&
510	& " already exist.")
511	ENDIF
512
513	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
514	& TRIM(cp_card(jk)) )
515	INQUIRE(FILE=TRIM(cl_bdyout), EXIST=ll_exist)
516	IF( ll_exist )THEN
517	CALL logger_fatal("CREATE BOUNDARY: output file "//TRIM(cl_bdyout)//&
518	& " already exist.")
519	ENDIF
520	ENDDO
521
522	! check namelist
523	! check refinement factor
524	il_rho(:)=1
525	IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
526	CALL logger_error("CREATE BOUNDARY: invalid refinement factor."//&
527	& " check namelist "//TRIM(cl_namelist))
528	ELSE
529	il_rho(jp_I)=in_rhoi
530	il_rho(jp_J)=in_rhoj
531	ENDIF
532
533	!
534	! compute coarse grid indices around fine grid
535	il_ind(:,:)=grid_get_coarse_index(tl_coord0, tl_coord1, &
536	& id_rho=il_rho(:))
537
538	il_imin0=il_ind(1,1) ; il_imax0=il_ind(1,2)
539	il_jmin0=il_ind(2,1) ; il_jmax0=il_ind(2,2)
540
541	! check domain validity
542	CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0)
543
544	! check coordinate file
545	CALL grid_check_coincidence( tl_coord0, tl_coord1, &
546	& il_imin0, il_imax0, &
547	& il_jmin0, il_jmax0, &
548	& il_rho(:) )
549
550	! read or compute boundary
551	CALL mpp_get_contour(tl_bathy1)
552
553	CALL iom_mpp_open(tl_bathy1)
554
555	tl_var1=iom_mpp_read_var(tl_bathy1,'Bathymetry')
556
557	CALL iom_mpp_close(tl_bathy1)
558
559	! get boundaries indices
560	tl_bdy(:)=boundary_init(tl_var1, ln_north, ln_south, ln_east, ln_west, &
561	& cn_north, cn_south, cn_east, cn_west, &
562	& ln_oneseg )
563
564	CALL var_clean(tl_var1)
565
566	! compute level
567	ALLOCATE(tl_level(ip_npoint))
568	tl_level(:)=vgrid_get_level(tl_bathy1, cl_namelist )
569
570	! get coordinate for each segment of each boundary
571	ALLOCATE( tl_segdom1(ip_npoint,ip_maxseg,ip_ncard) )
572	ALLOCATE( tl_seglvl1(ip_npoint,ip_maxseg,ip_ncard) )
573
574	DO jl=1,ip_ncard
575	IF( tl_bdy(jl)%l_use )THEN
576	DO jk=1,tl_bdy(jl)%i_nseg
577
578	! get fine grid segment domain
579	tl_segdom1(:,jk,jl)=create_boundary_get_dom( tl_bathy1, &
580	& tl_bdy(jl), jk )
581
582	IF( .NOT. ln_extrap )THEN
583	! get fine grid level
584	tl_seglvl1(:,jk,jl)= &
585	& create_boundary_get_level( tl_level(:), &
586	& tl_segdom1(:,jk,jl))
587	ENDIF
588
589	! add extra band to fine grid domain (if possible)
590	! to avoid dimension of one and so be able to compute offset
591	DO jj=1,ip_npoint
592	CALL dom_add_extra(tl_segdom1(jj,jk,jl), &
593	& il_rho(jp_I), il_rho(jp_J))
594	ENDDO
595
596	ENDDO
597	ENDIF
598	ENDDO
599
600	! clean
601	CALL var_clean(tl_level(:))
602	DEALLOCATE(tl_level)
603
604	! clean bathy
605	CALL mpp_clean(tl_bathy1)
606
607	ALLOCATE( tl_segvar1(tl_multi%i_nvar,ip_maxseg,ip_ncard) )
608	! compute boundary for variable to be used (see namelist)
609	IF( .NOT. ASSOCIATED(tl_multi%t_mpp) )THEN
610	CALL logger_error("CREATE BOUNDARY: no file to work on. "//&
611	& "check cn_varfile in namelist.")
612	ELSE
613
614	jvar=0
615	! for each file
616	DO ji=1,tl_multi%i_nmpp
617
618	WRITE(cl_data,'(a,i2.2)') 'data-',jvar+1
619
620	IF( .NOT. ASSOCIATED(tl_multi%t_mpp(ji)%t_proc(1)%t_var) )THEN
621
622	CALL logger_error("CREATE BOUNDARY: no variable to work on for "//&
623	& "mpp "//TRIM(tl_multi%t_mpp(ji)%c_name)//&
624	& ". check cn_varfile in namelist.")
625
626	ELSEIF( TRIM(tl_multi%t_mpp(ji)%c_name) == TRIM(cl_data) )THEN
627	!- use input matrix to fill variable
628
629	WRITE(*,'(a)') "work on data"
630	! for each variable initialise from matrix
631	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
632
633	jvar=jvar+1
634	WRITE(*,'(2x,a,a)') "work on variable "//&
635	& TRIM(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name)
636
637	tl_var1=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
638
639	SELECT CASE(TRIM(tl_var1%c_point))
640	CASE DEFAULT !'T'
641	jpoint=jp_T
642	CASE('U')
643	jpoint=jp_U
644	CASE('V')
645	jpoint=jp_V
646	CASE('F')
647	jpoint=jp_F
648	END SELECT
649
650	WRITE(*,'(4x,a,a)') 'work on '//TRIM(tl_var1%c_name)
651	DO jl=1,ip_ncard
652	IF( tl_bdy(jl)%l_use )THEN
653
654	DO jk=1,tl_bdy(jl)%i_nseg
655
656	! fill value with matrix data
657	tl_segvar1(jvar,jk,jl)=create_boundary_matrix( &
658	& tl_var1, &
659	& tl_segdom1(jpoint,jk,jl), &
660	& in_nlevel )
661
662	!del extra
663	CALL dom_del_extra( tl_segvar1(jvar,jk,jl), &
664	& tl_segdom1(jpoint,jk,jl) )
665
666	ENDDO
667
668	ENDIF
669	ENDDO
670
671	! clean
672	CALL var_clean(tl_var1)
673
674	ENDDO
675
676	!- end of use input matrix to fill variable
677	ELSE
678	!- use file to fill variable
679
680	WRITE(*,'(a)') "work on file "//TRIM(tl_multi%t_mpp(ji)%c_name)
681	!
682	tl_mpp=mpp_init(file_init(TRIM(tl_multi%t_mpp(ji)%t_proc(1)%c_name)))
683	CALL grid_get_info(tl_mpp)
684
685	! check vertical dimension
686	IF( tl_mpp%t_dim(jp_K)%l_use .AND. &
687	& tl_mpp%t_dim(jp_K)%i_len /= in_nlevel )THEN
688	CALL logger_error("CREATE BOUNDARY: dimension in file "//&
689	& TRIM(tl_mpp%c_name)//" not agree with namelist in_nlevel ")
690	ENDIF
691
692	! open mpp file
693	CALL iom_mpp_open(tl_mpp)
694
695	! get or check depth value
696	CALL create_boundary_check_depth( tl_mpp, tl_depth )
697
698	! get or check time value
699	CALL create_boundary_check_time( tl_mpp, tl_time )
700
701	! close mpp file
702	CALL iom_mpp_close(tl_mpp)
703
704	IF( ANY( tl_mpp%t_dim(1:2)%i_len /= &
705	& tl_coord0%t_dim(1:2)%i_len) )THEN
706	!- extract value from fine grid
707
708	IF( ANY( tl_mpp%t_dim(1:2)%i_len <= &
709	& tl_coord1%t_dim(1:2)%i_len) )THEN
710	CALL logger_fatal("CREATE BOUNDARY: dimension in file "//&
711	& TRIM(tl_mpp%c_name)//" smaller than those in fine"//&
712	& " grid coordinates.")
713	ENDIF
714
715	DO jl=1,ip_ncard
716	IF( tl_bdy(jl)%l_use )THEN
717
718	WRITE(*,'(2x,a,a)') 'work on '//TRIM(tl_bdy(jl)%c_card)//&
719	& ' boundary'
720	DO jk=1,tl_bdy(jl)%i_nseg
721	! compute domain on fine grid
722
723	! for each variable of this file
724	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
725
726	cl_name=tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name
727	WRITE(*,'(4x,a,a)') "work on (extract) variable "//&
728	& TRIM(cl_name)
729
730	cl_point=tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_point
731	! open mpp file on domain
732	SELECT CASE(TRIM(cl_point))
733	CASE DEFAULT !'T'
734	jpoint=jp_T
735	CASE('U')
736	jpoint=jp_U
737	CASE('V')
738	jpoint=jp_V
739	CASE('F')
740	jpoint=jp_F
741	END SELECT
742
743	tl_dom1=dom_copy(tl_segdom1(jpoint,jk,jl))
744
745	! open mpp files
746	CALL iom_dom_open(tl_mpp, tl_dom1)
747
748	!7-5 read variable over domain
749	tl_segvar1(jvar+jj,jk,jl)=iom_dom_read_var( &
750	& tl_mpp, TRIM(cl_name), tl_dom1)
751
752	! del extra point
753	CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), &
754	& tl_dom1 )
755
756	! clean extra point information on fine grid domain
757	CALL dom_clean_extra( tl_dom1 )
758
759	! add attribute to variable
760	tl_att=att_init('src_file', &
761	& TRIM(fct_basename(tl_mpp%c_name)))
762	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att)
763
764	tl_att=att_init('src_i_indices', &
765	& (/tl_dom1%i_imin, tl_dom1%i_imax/))
766	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att)
767
768	tl_att=att_init('src_j_indices', &
769	& (/tl_dom1%i_jmin, tl_dom1%i_jmax/))
770	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att)
771
772	! clean structure
773	CALL att_clean(tl_att)
774	CALL dom_clean(tl_dom1)
775
776	! close mpp files
777	CALL iom_dom_close(tl_mpp)
778
779	! clean
780	CALL var_clean(tl_lvl1)
781
782	ENDDO ! jj
783	ENDDO ! jk
784
785	ENDIF
786	ENDDO ! jl
787
788	! clean
789	CALL mpp_clean(tl_mpp)
790
791	jvar=jvar+tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
792
793	!- end of extract value from fine grid
794	ELSE
795	!- interpolate value from coarse grid
796
797	DO jl=1,ip_ncard
798	IF( tl_bdy(jl)%l_use )THEN
799
800	WRITE(*,'(2x,a,a)') 'work on '//TRIM(tl_bdy(jl)%c_card)//&
801	& ' boundary'
802	DO jk=1,tl_bdy(jl)%i_nseg
803
804	! for each variable of this file
805	DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
806
807	WRITE(*,'(4x,a,a)') "work on (interp) variable "//&
808	& TRIM(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name)
809
810	tl_var0=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
811	! open mpp file on domain
812	SELECT CASE(TRIM(tl_var0%c_point))
813	CASE DEFAULT !'T'
814	jpoint=jp_T
815	CASE('U')
816	jpoint=jp_U
817	CASE('V')
818	jpoint=jp_V
819	CASE('F')
820	jpoint=jp_F
821	END SELECT
822
823	tl_dom1=dom_copy(tl_segdom1(jpoint,jk,jl))
824
825	CALL create_boundary_get_coord( tl_coord1, tl_dom1, &
826	& tl_var0%c_point, &
827	& tl_lon1, tl_lat1 )
828
829	! get coarse grid indices of this segment
830	il_ind(:,:)=grid_get_coarse_index(tl_coord0, &
831	& tl_lon1, tl_lat1, &
832	& id_rho=il_rho(:) )
833
834	IF( ANY(il_ind(:,:)==0) )THEN
835	CALL logger_error("CREATE BOUNDARY: error "//&
836	& "computing coarse grid indices")
837	ELSE
838	il_imin0=il_ind(1,1)
839	il_imax0=il_ind(1,2)
840
841	il_jmin0=il_ind(2,1)
842	il_jmax0=il_ind(2,2)
843	ENDIF
844
845	il_offset(:,:)= grid_get_fine_offset( &
846	& tl_coord0, &
847	& il_imin0, il_jmin0,&
848	& il_imax0, il_jmax0,&
849	& tl_lon1%d_value(:,:,1,1),&
850	& tl_lat1%d_value(:,:,1,1),&
851	& il_rho(:),&
852	& TRIM(tl_var0%c_point) )
853
854	! compute coarse grid segment domain
855	tl_dom0=dom_init( tl_coord0, &
856	& il_imin0, il_imax0,&
857	& il_jmin0, il_jmax0 )
858
859	! add extra band (if possible) to compute
860	! interpolation
861	CALL dom_add_extra(tl_dom0)
862
863	! read variables on domain
864	! open mpp files
865	CALL iom_dom_open(tl_mpp, tl_dom0)
866
867	cl_name=tl_var0%c_name
868	! read variable value on domain
869	tl_segvar1(jvar+jj,jk,jl)= &
870	& iom_dom_read_var(tl_mpp, TRIM(cl_name), tl_dom0)
871
872	! work on variable
873	CALL create_boundary_interp( &
874	& tl_segvar1(jvar+jj,jk,jl),&
875	& il_rho(:), il_offset(:,:) )
876
877	! remove extraband added to domain
878	CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), &
879	& tl_dom0, il_rho(:) )
880
881	! del extra point on fine grid
882	CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), &
883	& tl_dom1 )
884	! clean extra point information on coarse grid domain
885	CALL dom_clean_extra( tl_dom0 )
886
887	! add attribute to variable
888	tl_att=att_init('src_file',&
889	& TRIM(fct_basename(tl_mpp%c_name)))
890	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
891	& tl_att)
892
893	! use clean extra avt creer attribut
894	tl_att=att_init('src_i-indices',&
895	& (/tl_dom0%i_imin, tl_dom0%i_imax/))
896	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
897	& tl_att)
898
899	tl_att=att_init('src_j-indices', &
900	& (/tl_dom0%i_jmin, tl_dom0%i_jmax/))
901	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
902	& tl_att)
903
904	IF( ANY(il_rho(:)/=1) )THEN
905	tl_att=att_init("refinment_factor", &
906	& (/il_rho(jp_I),il_rho(jp_J)/))
907	CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
908	& tl_att)
909	ENDIF
910
911	! clean structure
912	CALL att_clean(tl_att)
913
914	! clean
915	CALL dom_clean(tl_dom0)
916	CALL dom_clean(tl_dom1)
917
918	! close mpp files
919	CALL iom_dom_close(tl_mpp)
920
921	! clean structure
922	CALL var_clean(tl_lon1)
923	CALL var_clean(tl_lat1)
924	CALL var_clean(tl_lvl1)
925
926	ENDDO ! jj
927
928	! clean
929	CALL var_clean(tl_var0)
930
931	ENDDO ! jk
932
933	ENDIF
934	ENDDO ! jl
935
936	jvar=jvar+tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
937
938	!- end of interpolate value from coarse grid
939	ENDIF
940
941	! clean
942	CALL mpp_clean(tl_mpp)
943
944	!- end of use file to fill variable
945	ENDIF
946	ENDDO
947	ENDIF
948
949	IF( jvar /= tl_multi%i_nvar )THEN
950	CALL logger_error("CREATE BOUNDARY: it seems some variable "//&
951	& "can not be read")
952	ENDIF
953
954	! write file for each segment of each boundary
955	DO jl=1,ip_ncard
956	IF( tl_bdy(jl)%l_use )THEN
957
958	DO jk=1,tl_bdy(jl)%i_nseg
959	!-
960	CALL create_boundary_get_coord( tl_coord1, tl_segdom1(jp_T,jk,jl),&
961	& 'T', tl_lon1, tl_lat1 )
962
963	! force to use nav_lon, nav_lat as variable name
964	tl_lon1%c_name='nav_lon'
965	tl_lat1%c_name='nav_lat'
966
967	! del extra point on fine grid
968	CALL dom_del_extra( tl_lon1, tl_segdom1(jp_T,jk,jl) )
969	CALL dom_del_extra( tl_lat1, tl_segdom1(jp_T,jk,jl) )
970
971	! clean
972	DO jpoint=1,ip_npoint
973	CALL dom_clean(tl_segdom1(jpoint,jk,jl))
974	ENDDO
975
976	! swap array
977	CALL boundary_swap(tl_lon1, tl_bdy(jl))
978	CALL boundary_swap(tl_lat1, tl_bdy(jl))
979	DO jvar=1,tl_multi%i_nvar
980
981	! use additional request
982	! change unit and apply factor
983	CALL var_chg_unit(tl_segvar1(jvar,jk,jl))
984
985	! forced min and max value
986	CALL var_limit_value(tl_segvar1(jvar,jk,jl))
987
988	! filter
989	CALL filter_fill_value(tl_segvar1(jvar,jk,jl))
990
991	IF( .NOT. ln_extrap )THEN
992	! use mask
993	SELECT CASE(TRIM(tl_segvar1(jvar,jk,jl)%c_point))
994	CASE DEFAULT !'T'
995	jpoint=jp_T
996	CASE('U')
997	jpoint=jp_U
998	CASE('V')
999	jpoint=jp_V
1000	CASE('F')
1001	jpoint=jp_F
1002	END SELECT
1003
1004	CALL create_boundary_use_mask(tl_segvar1(jvar,jk,jl), &
1005	& tl_seglvl1(jpoint,jk,jl))
1006	ENDIF
1007
1008	! swap dimension order
1009	CALL boundary_swap(tl_segvar1(jvar,jk,jl), tl_bdy(jl))
1010
1011	ENDDO
1012
1013	! create file
1014	! create file structure
1015	! set file namearray of level variable structure
1016	IF( tl_bdy(jl)%i_nseg > 1 )THEN
1017	IF( ASSOCIATED(tl_time%d_value) )THEN
1018	cl_fmt="('y',i0.4,'m',i0.2,'d',i0.2)"
1019	tl_date=var_to_date(tl_time)
1020	tl_date=tl_date+dn_dayofs
1021	cl_date=date_print( tl_date, cl_fmt )
1022
1023	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1024	& TRIM(tl_bdy(jl)%c_card), jk,&
1025	& cd_date=TRIM(cl_date) )
1026	ELSE
1027	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1028	& TRIM(tl_bdy(jl)%c_card), jk )
1029	ENDIF
1030	ELSE
1031	IF( ASSOCIATED(tl_time%d_value) )THEN
1032	cl_fmt="('y',i0.4,'m',i0.2,'d',i0.2)"
1033	tl_date=var_to_date(tl_time)
1034	tl_date=tl_date+dn_dayofs
1035	cl_date=date_print( tl_date, cl_fmt )
1036
1037	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1038	& TRIM(tl_bdy(jl)%c_card), &
1039	& cd_date=TRIM(cl_date) )
1040	ELSE
1041	cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1042	& TRIM(tl_bdy(jl)%c_card) )
1043	ENDIF
1044	ENDIF
1045	!
1046	tl_fileout=file_init(TRIM(cl_bdyout),id_perio=in_perio1)
1047
1048	! add dimension
1049	tl_dim(:)=var_max_dim(tl_segvar1(:,jk,jl))
1050
1051	SELECT CASE(TRIM(tl_bdy(jl)%c_card))
1052	CASE DEFAULT ! 'north','south'
1053	cl_dimorder='xyzt'
1054	CASE('east','west')
1055	cl_dimorder='yxzt'
1056	END SELECT
1057
1058	DO ji=1,ip_maxdim
1059	IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
1060	ENDDO
1061
1062	! add variables
1063	IF( ALL( tl_dim(1:2)%l_use ) )THEN
1064	! add longitude
1065	CALL file_add_var(tl_fileout, tl_lon1)
1066	CALL var_clean(tl_lon1)
1067
1068	! add latitude
1069	CALL file_add_var(tl_fileout, tl_lat1)
1070	CALL var_clean(tl_lat1)
1071	ENDIF
1072
1073
1074
1075	IF( tl_dim(3)%l_use )THEN
1076	IF( ASSOCIATED(tl_depth%d_value) )THEN
1077	! add depth
1078	CALL file_add_var(tl_fileout, tl_depth)
1079	ENDIF
1080	ENDIF
1081
1082	IF( tl_dim(4)%l_use )THEN
1083	IF( ASSOCIATED(tl_time%d_value) )THEN
1084	! add time
1085	CALL file_add_var(tl_fileout, tl_time)
1086	ENDIF
1087	ENDIF
1088
1089	! add other variable
1090	DO jvar=tl_multi%i_nvar,1,-1
1091	CALL file_add_var(tl_fileout, tl_segvar1(jvar,jk,jl))
1092	CALL var_clean(tl_segvar1(jvar,jk,jl))
1093	ENDDO
1094
1095	! add some attribute
1096	tl_att=att_init("Created_by","SIREN create_boundary")
1097	CALL file_add_att(tl_fileout, tl_att)
1098
1099	cl_date=date_print(date_now())
1100	tl_att=att_init("Creation_date",cl_date)
1101	CALL file_add_att(tl_fileout, tl_att)
1102
1103	! add shift on north and east boundary
1104	! boundary compute on T point but express on U or V point
1105	SELECT CASE(TRIM(tl_bdy(jl)%c_card))
1106	CASE DEFAULT ! 'south','west'
1107	il_shift=0
1108	CASE('north','east')
1109	il_shift=1
1110	END SELECT
1111
1112	! add indice of velocity row or column
1113	tl_att=att_init('bdy_ind',tl_bdy(jl)%t_seg(jk)%i_index-il_shift)
1114	CALL file_move_att(tl_fileout, tl_att)
1115
1116	! add width of the relaxation zone
1117	tl_att=att_init('bdy_width',tl_bdy(jl)%t_seg(jk)%i_width)
1118	CALL file_move_att(tl_fileout, tl_att)
1119
1120	! add indice of segment start
1121	tl_att=att_init('bdy_deb',tl_bdy(jl)%t_seg(jk)%i_first)
1122	CALL file_move_att(tl_fileout, tl_att)
1123
1124	! add indice of segment end
1125	tl_att=att_init('bdy_end',tl_bdy(jl)%t_seg(jk)%i_last)
1126	CALL file_move_att(tl_fileout, tl_att)
1127
1128	! clean
1129	CALL att_clean(tl_att)
1130
1131	! create file
1132	CALL iom_create(tl_fileout)
1133
1134	! write file
1135	CALL iom_write_file(tl_fileout, cl_dimorder)
1136
1137	! close file
1138	CALL iom_close(tl_fileout)
1139	CALL file_clean(tl_fileout)
1140
1141	ENDDO ! jk
1142
1143	ENDIF
1144	! clean
1145	CALL boundary_clean(tl_bdy(jl))
1146	ENDDO !jl
1147
1148	! clean
1149	IF( ASSOCIATED(tl_depth%d_value) ) CALL var_clean(tl_depth)
1150	IF( ASSOCIATED(tl_time%d_value) ) CALL var_clean(tl_time)
1151	DEALLOCATE( tl_segdom1 )
1152	DEALLOCATE( tl_segvar1 )
1153	CALL var_clean(tl_seglvl1(:,:,:))
1154	DEALLOCATE( tl_seglvl1 )
1155
1156
1157	CALL mpp_clean(tl_coord1)
1158	CALL mpp_clean(tl_coord0)
1159
1160	CALL multi_clean(tl_multi)
1161
1162	! close log file
1163	CALL logger_footer()
1164	CALL logger_close()
1165
1166	CONTAINS
1167	!-------------------------------------------------------------------
1168	!> @brief
1169	!> This subroutine compute boundary domain for each grid point (T,U,V,F)
1170	!>
1171	!> @author J.Paul
1172	!> @date November, 2013 - Initial Version
1173	!> @date September, 2014
1174	!> - take into account grid point to compute boundary indices
1175	!>
1176	!> @param[in] td_bathy1 file structure
1177	!> @param[in] td_bdy boundary structure
1178	!> @param[in] id_seg segment indice
1179	!> @return array of domain structure
1180	!-------------------------------------------------------------------
1181	FUNCTION create_boundary_get_dom( td_bathy1, td_bdy, id_seg )
1182
1183	IMPLICIT NONE
1184
1185	! Argument
1186	TYPE(TMPP) , INTENT(IN ) :: td_bathy1
1187	TYPE(TBDY) , INTENT(IN ) :: td_bdy
1188	INTEGER(i4), INTENT(IN ) :: id_seg
1189
1190	! function
1191	TYPE(TDOM), DIMENSION(ip_npoint) :: create_boundary_get_dom
1192
1193	! local variable
1194	INTEGER(i4) :: il_imin1
1195	INTEGER(i4) :: il_imax1
1196	INTEGER(i4) :: il_jmin1
1197	INTEGER(i4) :: il_jmax1
1198
1199	INTEGER(i4) :: il_imin
1200	INTEGER(i4) :: il_imax
1201	INTEGER(i4) :: il_jmin
1202	INTEGER(i4) :: il_jmax
1203
1204	INTEGER(i4), DIMENSION(ip_npoint) :: il_ishift
1205	INTEGER(i4), DIMENSION(ip_npoint) :: il_jshift
1206
1207	! loop indices
1208	INTEGER(i4) :: ji
1209	INTEGER(i4) :: jk
1210	!----------------------------------------------------------------
1211	! init
1212	jk=id_seg
1213
1214	il_ishift(:)=0
1215	il_jshift(:)=0
1216
1217	! get boundary definition
1218	SELECT CASE(TRIM(td_bdy%c_card))
1219	CASE('north')
1220
1221	il_imin1=td_bdy%t_seg(jk)%i_first
1222	il_imax1=td_bdy%t_seg(jk)%i_last
1223	il_jmin1=td_bdy%t_seg(jk)%i_index-(td_bdy%t_seg(jk)%i_width-1)
1224	il_jmax1=td_bdy%t_seg(jk)%i_index
1225
1226	il_jshift(jp_V)=-1
1227	il_jshift(jp_F)=-1
1228
1229	CASE('south')
1230
1231	il_imin1=td_bdy%t_seg(jk)%i_first
1232	il_imax1=td_bdy%t_seg(jk)%i_last
1233	il_jmin1=td_bdy%t_seg(jk)%i_index
1234	il_jmax1=td_bdy%t_seg(jk)%i_index+(td_bdy%t_seg(jk)%i_width-1)
1235
1236	CASE('east')
1237
1238	il_imin1=td_bdy%t_seg(jk)%i_index-(td_bdy%t_seg(jk)%i_width-1)
1239	il_imax1=td_bdy%t_seg(jk)%i_index
1240	il_jmin1=td_bdy%t_seg(jk)%i_first
1241	il_jmax1=td_bdy%t_seg(jk)%i_last
1242
1243	il_ishift(jp_U)=-1
1244	il_ishift(jp_F)=-1
1245
1246	CASE('west')
1247
1248	il_imin1=td_bdy%t_seg(jk)%i_index
1249	il_imax1=td_bdy%t_seg(jk)%i_index+(td_bdy%t_seg(jk)%i_width-1)
1250	il_jmin1=td_bdy%t_seg(jk)%i_first
1251	il_jmax1=td_bdy%t_seg(jk)%i_last
1252
1253	END SELECT
1254
1255	!-read fine grid domain
1256	DO ji=1,ip_npoint
1257
1258	! shift domain
1259	il_imin=il_imin1+il_ishift(ji)
1260	il_imax=il_imax1+il_ishift(ji)
1261
1262	il_jmin=il_jmin1+il_jshift(ji)
1263	il_jmax=il_jmax1+il_jshift(ji)
1264
1265	! compute domain
1266	create_boundary_get_dom(ji)=dom_init( td_bathy1, &
1267	& il_imin, il_imax,&
1268	& il_jmin, il_jmax,&
1269	& TRIM(td_bdy%c_card) )
1270
1271	ENDDO
1272
1273	END FUNCTION create_boundary_get_dom
1274	!-------------------------------------------------------------------
1275	!> @brief
1276	!> This subroutine get coordinates over boundary domain
1277	!>
1278	!> @author J.Paul
1279	!> @date November, 2013 - Initial Version
1280	!> @date September, 2014
1281	!> - take into account grid point
1282	!>
1283	!> @param[in] td_coord1 coordinates file structure
1284	!> @param[in] td_dom1 boundary domain structure
1285	!> @param[in] cd_point grid point
1286	!> @param[out] td_lon1 longitude variable structure
1287	!> @param[out] td_lat1 latitude variable structure
1288	!-------------------------------------------------------------------
1289	SUBROUTINE create_boundary_get_coord( td_coord1, td_dom1, cd_point, &
1290	& td_lon1, td_lat1 )
1291
1292	IMPLICIT NONE
1293	! Argument
1294	TYPE(TMPP) , INTENT(IN ) :: td_coord1
1295	TYPE(TDOM) , INTENT(IN ) :: td_dom1
1296	TYPE(TVAR) , INTENT( OUT) :: td_lon1
1297	TYPE(TVAR) , INTENT( OUT) :: td_lat1
1298	CHARACTER(LEN=*), INTENT(IN ), OPTIONAL :: cd_point
1299
1300	! local variable
1301	TYPE(TMPP) :: tl_coord1
1302
1303	CHARACTER(LEN=lc) :: cl_name
1304	! loop indices
1305	!----------------------------------------------------------------
1306	!read variables on domain (ugly way to do it, have to work on it)
1307	! init mpp structure
1308	tl_coord1=mpp_copy(td_coord1)
1309
1310	! open mpp files
1311	CALL iom_dom_open(tl_coord1, td_dom1)
1312
1313	! read variable value on domain
1314	WRITE(cl_name,*) 'longitude_'//TRIM(cd_point)
1315	td_lon1=iom_dom_read_var( tl_coord1, TRIM(cl_name), td_dom1)
1316	WRITE(cl_name,*) 'latitude_'//TRIM(cd_point)
1317	td_lat1=iom_dom_read_var( tl_coord1, TRIM(cl_name), td_dom1)
1318
1319	! close mpp files
1320	CALL iom_dom_close(tl_coord1)
1321
1322	! clean structure
1323	CALL mpp_clean(tl_coord1)
1324
1325	END SUBROUTINE create_boundary_get_coord
1326	!-------------------------------------------------------------------
1327	!> @brief
1328	!> This subroutine interpolate variable on boundary
1329	!>
1330	!> @details
1331	!>
1332	!> @author J.Paul
1333	!> @date November, 2013 - Initial Version
1334	!>
1335	!> @param[inout] td_var variable structure
1336	!> @param[in] id_rho array of refinment factor
1337	!> @param[in] id_offset array of offset between fine and coarse grid
1338	!> @param[in] id_iext i-direction size of extra bands (default=im_minext)
1339	!> @param[in] id_jext j-direction size of extra bands (default=im_minext)
1340	!-------------------------------------------------------------------
1341	SUBROUTINE create_boundary_interp( td_var, &
1342	& id_rho, &
1343	& id_offset, &
1344	& id_iext, id_jext )
1345
1346	IMPLICIT NONE
1347
1348	! Argument
1349	TYPE(TVAR) , INTENT(INOUT) :: td_var
1350	INTEGER(I4), DIMENSION(:) , INTENT(IN ) :: id_rho
1351	INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset
1352
1353	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext
1354	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext
1355
1356
1357	! local variable
1358	INTEGER(i4) :: il_iext
1359	INTEGER(i4) :: il_jext
1360	! loop indices
1361	!----------------------------------------------------------------
1362
1363	!WARNING: at least two extrabands are required for cubic interpolation
1364	il_iext=2
1365	IF( PRESENT(id_iext) ) il_iext=id_iext
1366
1367	il_jext=2
1368	IF( PRESENT(id_jext) ) il_jext=id_jext
1369
1370	IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1371	CALL logger_warn("CREATE BOUNDARY INTERP: at least extrapolation "//&
1372	& "on two points are required with cubic interpolation ")
1373	il_iext=2
1374	ENDIF
1375
1376	IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1377	CALL logger_warn("CREATE BOUNDARY INTERP: at least extrapolation "//&
1378	& "on two points are required with cubic interpolation ")
1379	il_jext=2
1380	ENDIF
1381
1382	! work on variable
1383	! add extraband
1384	CALL extrap_add_extrabands(td_var, il_iext, il_jext)
1385
1386	! extrapolate variable
1387	CALL extrap_fill_value( td_var )
1388
1389	! interpolate Bathymetry
1390	CALL interp_fill_value( td_var, id_rho(:), &
1391	& id_offset=id_offset(:,:) )
1392
1393	! remove extraband
1394	CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), &
1395	& il_jext*id_rho(jp_J))
1396
1397	END SUBROUTINE create_boundary_interp
1398	!-------------------------------------------------------------------
1399	!> @brief
1400	!> This function create variable, filled with matrix value
1401	!>
1402	!> @details
1403	!> A variable is create with the same name that the input variable,
1404	!> and with dimension of the coordinate file.
1405	!> Then the variable array of value is split into equal subdomain.
1406	!> Each subdomain is fill with the linked value of the matrix.
1407	!>
1408	!> @author J.Paul
1409	!> @date November, 2013 - Initial Version
1410	!>
1411	!> @param[in] td_var variable structure
1412	!> @param[in] td_dom domain structure
1413	!> @param[in] id_nlevel number of levels
1414	!> @return variable structure
1415	!-------------------------------------------------------------------
1416	FUNCTION create_boundary_matrix(td_var, td_dom, id_nlevel)
1417	IMPLICIT NONE
1418	! Argument
1419	TYPE(TVAR) , INTENT(IN) :: td_var
1420	TYPE(TDOM) , INTENT(IN) :: td_dom
1421	INTEGER(i4), INTENT(IN) :: id_nlevel
1422
1423	! function
1424	TYPE(TVAR) :: create_boundary_matrix
1425
1426	! local variable
1427	INTEGER(i4) , DIMENSION(3) :: il_dim
1428	INTEGER(i4) , DIMENSION(3) :: il_size
1429	INTEGER(i4) , DIMENSION(3) :: il_rest
1430
1431	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ishape
1432	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_jshape
1433	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_kshape
1434
1435	REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
1436
1437	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
1438
1439	! loop indices
1440	INTEGER(i4) :: ji
1441	INTEGER(i4) :: jj
1442	INTEGER(i4) :: jk
1443	!----------------------------------------------------------------
1444
1445	! write value on grid
1446	! get matrix dimension
1447	il_dim(:)=td_var%t_dim(1:3)%i_len
1448
1449	tl_dim(jp_I:jp_J)=dim_copy(td_dom%t_dim(jp_I:jp_J))
1450	tl_dim(jp_K)%i_len=id_nlevel
1451
1452	! split output domain in N subdomain depending of matrix dimension
1453	il_size(:) = tl_dim(1:3)%i_len / il_dim(:)
1454	il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:))
1455
1456	ALLOCATE( il_ishape(il_dim(1)+1) )
1457	il_ishape(:)=0
1458	DO ji=2,il_dim(1)+1
1459	il_ishape(ji)=il_ishape(ji-1)+il_size(1)
1460	ENDDO
1461	! add rest to last cell
1462	il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1)
1463
1464	ALLOCATE( il_jshape(il_dim(2)+1) )
1465	il_jshape(:)=0
1466	DO jj=2,il_dim(2)+1
1467	il_jshape(jj)=il_jshape(jj-1)+il_size(2)
1468	ENDDO
1469	! add rest to last cell
1470	il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2)
1471
1472	ALLOCATE( il_kshape(il_dim(3)+1) )
1473	il_kshape(:)=0
1474	DO jk=2,il_dim(3)+1
1475	il_kshape(jk)=il_kshape(jk-1)+il_size(3)
1476	ENDDO
1477	! add rest to last cell
1478	il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3)
1479
1480	! write ouput array of value
1481	ALLOCATE(dl_value( tl_dim(1)%i_len, &
1482	& tl_dim(2)%i_len, &
1483	& tl_dim(3)%i_len, &
1484	& tl_dim(4)%i_len) )
1485
1486	dl_value(:,:,:,:)=0
1487
1488	DO jk=2,il_dim(3)+1
1489	DO jj=2,il_dim(2)+1
1490	DO ji=2,il_dim(1)+1
1491
1492	dl_value( 1+il_ishape(ji-1):il_ishape(ji), &
1493	& 1+il_jshape(jj-1):il_jshape(jj), &
1494	& 1+il_kshape(jk-1):il_kshape(jk), &
1495	& 1 ) = td_var%d_value(ji-1,jj-1,jk-1,1)
1496
1497	ENDDO
1498	ENDDO
1499	ENDDO
1500
1501	! initialise variable with value
1502	create_boundary_matrix=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:))
1503
1504	DEALLOCATE(dl_value)
1505
1506	END FUNCTION create_boundary_matrix
1507	!-------------------------------------------------------------------
1508	!> @brief
1509	!> This subroutine use mask to filled land point with _FillValue
1510	!>
1511	!> @details
1512	!>
1513	!> @author J.Paul
1514	!> @date November, 2013 - Initial Version
1515	!>
1516	!> @param[inout] td_var variable structure
1517	!> @param[in] td_mask mask variable structure
1518	!-------------------------------------------------------------------
1519	SUBROUTINE create_boundary_use_mask( td_var, td_mask )
1520
1521	IMPLICIT NONE
1522
1523	! Argument
1524	TYPE(TVAR), INTENT(INOUT) :: td_var
1525	TYPE(TVAR), INTENT(IN ) :: td_mask
1526
1527	! local variable
1528	INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask
1529
1530	! loop indices
1531	INTEGER(i4) :: jk
1532	INTEGER(i4) :: jl
1533	!----------------------------------------------------------------
1534
1535	IF( ANY(td_var%t_dim(1:2)%i_len /= &
1536	& td_mask%t_dim(1:2)%i_len) )THEN
1537	CALL logger_debug(" mask dimension ( "//&
1538	& TRIM(fct_str(td_mask%t_dim(1)%i_len))//","//&
1539	& TRIM(fct_str(td_mask%t_dim(2)%i_len))//")" )
1540	CALL logger_debug(" variable dimension ( "//&
1541	& TRIM(fct_str(td_var%t_dim(1)%i_len))//","//&
1542	& TRIM(fct_str(td_var%t_dim(2)%i_len))//")" )
1543	CALL logger_fatal("CREATE BOUNDARY USE MASK: mask and "//&
1544	& "variable dimension differ." )
1545	ENDIF
1546
1547	ALLOCATE( il_mask(td_var%t_dim(1)%i_len, &
1548	& td_var%t_dim(2)%i_len) )
1549
1550	il_mask(:,:)=INT(td_mask%d_value(:,:,1,1))
1551
1552	DO jl=1,td_var%t_dim(4)%i_len
1553	DO jk=1,td_var%t_dim(3)%i_len
1554	WHERE( il_mask(:,:) < jk ) td_var%d_value(:,:,jk,jl)=td_var%d_fill
1555	ENDDO
1556	ENDDO
1557
1558	DEALLOCATE( il_mask )
1559
1560	END SUBROUTINE create_boundary_use_mask
1561	!-------------------------------------------------------------------
1562	!> @brief
1563	!> This function extract level over domain on each grid point, and return
1564	!> array of variable structure
1565	!>
1566	!> @author J.Paul
1567	!> @date November, 2013 - Initial Version
1568	!>
1569	!> @param[in] td_level array of level variable structure
1570	!> @param[in] td_dom array of domain structure
1571	!> @return array of variable structure
1572	!-------------------------------------------------------------------
1573	FUNCTION create_boundary_get_level(td_level, td_dom)
1574	IMPLICIT NONE
1575	! Argument
1576	TYPE(TVAR), DIMENSION(:), INTENT(IN) :: td_level
1577	TYPE(TDOM), DIMENSION(:), INTENT(IN) :: td_dom
1578
1579	! function
1580	TYPE(TVAR), DIMENSION(ip_npoint) :: create_boundary_get_level
1581
1582	! local variable
1583	TYPE(TVAR), DIMENSION(ip_npoint) :: tl_var
1584
1585	! loop indices
1586	INTEGER(i4) :: ji
1587	!----------------------------------------------------------------
1588
1589	IF( SIZE(td_level(:)) /= ip_npoint .OR. &
1590	& SIZE(td_dom(:)) /= ip_npoint )THEN
1591	CALL logger_error("CREATE BDY GET LEVEL: invalid dimension. "//&
1592	& "check input array of level and domain.")
1593	ELSE
1594
1595	DO ji=1,ip_npoint
1596
1597	tl_var(ji)=var_copy(td_level(ji))
1598
1599	IF( ASSOCIATED(tl_var(ji)%d_value) ) DEALLOCATE(tl_var(ji)%d_value)
1600
1601	tl_var(ji)%t_dim(1)%i_len=td_dom(ji)%t_dim(1)%i_len
1602	tl_var(ji)%t_dim(2)%i_len=td_dom(ji)%t_dim(2)%i_len
1603	ALLOCATE(tl_var(ji)%d_value(tl_var(ji)%t_dim(1)%i_len, &
1604	& tl_var(ji)%t_dim(2)%i_len, &
1605	& tl_var(ji)%t_dim(3)%i_len, &
1606	& tl_var(ji)%t_dim(4)%i_len) )
1607
1608	tl_var(ji)%d_value(:,:,:,:) = &
1609	& td_level(ji)%d_value( td_dom(ji)%i_imin:td_dom(ji)%i_imax, &
1610	& td_dom(ji)%i_jmin:td_dom(ji)%i_jmax, :, : )
1611
1612	ENDDO
1613	! save result
1614	create_boundary_get_level(:)=var_copy(tl_var(:))
1615
1616	! clean
1617	CALL var_clean(tl_var(:))
1618
1619	ENDIF
1620	END FUNCTION create_boundary_get_level
1621	!-------------------------------------------------------------------
1622	!> @brief
1623	!> This subroutine get depth variable value in an open mpp structure
1624	!> and check if agree with already input depth variable.
1625	!>
1626	!> @details
1627	!>
1628	!> @author J.Paul
1629	!> @date November, 2014 - Initial Version
1630	!>
1631	!> @param[in] td_mpp mpp structure
1632	!> @param[inout] td_depth depth variable structure
1633	!-------------------------------------------------------------------
1634	SUBROUTINE create_boundary_check_depth( td_mpp, td_depth )
1635
1636	IMPLICIT NONE
1637
1638	! Argument
1639	TYPE(TMPP) , INTENT(IN ) :: td_mpp
1640	TYPE(TVAR) , INTENT(INOUT) :: td_depth
1641
1642	! local variable
1643	INTEGER(i4) :: il_varid
1644	TYPE(TVAR) :: tl_depth
1645	! loop indices
1646	!----------------------------------------------------------------
1647
1648	! get or check depth value
1649	IF( td_mpp%t_proc(1)%i_depthid /= 0 )THEN
1650
1651	il_varid=td_mpp%t_proc(1)%i_depthid
1652	IF( ASSOCIATED(td_depth%d_value) )THEN
1653
1654	tl_depth=iom_mpp_read_var(td_mpp, il_varid)
1655	IF( ANY( td_depth%d_value(:,:,:,:) /= &
1656	& tl_depth%d_value(:,:,:,:) ) )THEN
1657
1658	CALL logger_fatal("CREATE BOUNDARY: depth value from "//&
1659	& TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//&
1660	& " to those from former file(s).")
1661
1662	ENDIF
1663	CALL var_clean(tl_depth)
1664
1665	ELSE
1666	td_depth=iom_mpp_read_var(td_mpp,il_varid)
1667	ENDIF
1668
1669	ENDIF
1670
1671	END SUBROUTINE create_boundary_check_depth
1672	!-------------------------------------------------------------------
1673	!> @brief
1674	!> This subroutine get date and time in an open mpp structure
1675	!> and check if agree with date and time already read.
1676	!>
1677	!> @details
1678	!>
1679	!> @author J.Paul
1680	!> @date November, 2014 - Initial Version
1681	!>
1682	!> @param[in] td_mpp mpp structure
1683	!> @param[inout] td_time time variable structure
1684	!-------------------------------------------------------------------
1685	SUBROUTINE create_boundary_check_time( td_mpp, td_time )
1686
1687	IMPLICIT NONE
1688
1689	! Argument
1690	TYPE(TMPP), INTENT(IN ) :: td_mpp
1691	TYPE(TVAR), INTENT(INOUT) :: td_time
1692
1693	! local variable
1694	INTEGER(i4) :: il_varid
1695	TYPE(TVAR) :: tl_time
1696
1697	TYPE(TDATE) :: tl_date1
1698	TYPE(TDATE) :: tl_date2
1699	! loop indices
1700	!----------------------------------------------------------------
1701
1702	! get or check depth value
1703	IF( td_mpp%t_proc(1)%i_timeid /= 0 )THEN
1704
1705	il_varid=td_mpp%t_proc(1)%i_timeid
1706	IF( ASSOCIATED(td_time%d_value) )THEN
1707
1708	tl_time=iom_mpp_read_var(td_mpp, il_varid)
1709
1710	tl_date1=var_to_date(td_time)
1711	tl_date2=var_to_date(tl_time)
1712	IF( tl_date1 - tl_date2 /= 0 )THEN
1713
1714	CALL logger_fatal("CREATE BOUNDARY: date from "//&
1715	& TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//&
1716	& " to those from former file(s).")
1717
1718	ENDIF
1719	CALL var_clean(tl_time)
1720
1721	ELSE
1722	td_time=iom_mpp_read_var(td_mpp,il_varid)
1723	ENDIF
1724
1725	ENDIF
1726
1727	END SUBROUTINE create_boundary_check_time
1728	END PROGRAM create_boundary

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/TOOLS/SIREN/src/create_boundary.f90 @ 9987

Download in other formats: