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_boundary.F90 in branches/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/SIREN/src/create_boundary.F90 @ 8862

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

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

Download in other formats: