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 utils/tools/SIREN/src – NEMO

source: utils/tools/SIREN/src/create_boundary.F90

Last change on this file was 13369, checked in by jpaul, 4 years ago

update: cf changelog inside documentation

File size: 78.4 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! DESCRIPTION:
6!> @file
7!> @brief
8!> This program creates boundary files.
9!>
10!> @details
11!> @section sec1 method
12!> Variables are read from source/coarse grid standard output,
13!> extracted or interpolated on target/fine grid.
14!> Variables could also be manually written.<br/>
15!> @note
16!>    method could be different for each variable.
17!>
18!>  <br/>
19!> @image html  boundary_NEATL36_70.png
20!> <center>@image latex boundary_NEATL36_70.png
21!> </center>
22!>
23!> @section sec2 how to
24!> USAGE: create_boundary create_bounary.nam [-v] [-h]<br/>
25!>    - positional arguments:<br/>
26!>       - create_boundary.nam<br/>
27!>          namelist of create_boundary
28!>          @note
29!>             a template of the namelist could be created running (in templates directory):
30!>             @code{.sh}
31!>                python create_templates.py create_boundary
32!>             @endcode
33!>
34!>    - optional arguments:<br/>
35!>       - -h, --help<br/>
36!>          show this help message (and exit)<br/>
37!>       - -v, --version<br/>
38!>          show Siren's version   (and exit)<br/>
39!>    @note
40!>       compiled with @a key_mpp_mpi, could be run on multi processor :<br/>
41!>       USAGE: create_boundary create_bounary.nam create_bounary2.nam ... [-v] [-h]<br/>
42!>
43!> @section sec_boundary create_boundary.nam
44!>    create_boundary.nam contains 9 namelists:<br/>
45!>       - **namlog** to set logger parameters
46!>       - **namcfg** to set configuration file parameters
47!>       - **namsrc** to set source/coarse grid parameters
48!>       - **namtgt** to set target/fine grid parameters
49!>       - **namvar** to set variable parameters
50!>       - **namnst** to set sub domain and nesting paramters
51!>       - **nambdy** to set boundary parameters
52!>       - **namzgr** to set vertical grid parameters
53!>       - **namout** to set output parameters
54!>
55!>    here after, each sub-namelist parameters is detailed.
56!>    @note
57!>       default values are specified between brackets
58!>
59!> @subsection sublog namlog
60!>    the logger sub-namelist parameters are :
61!>
62!>    - **cn_logfile** [@a create_boundary.log]<br/>
63!>       logger filename
64!>
65!>    - **cn_verbosity** [@a warning]<br/>
66!>       verbosity level, choose between :
67!>          - trace
68!>          - debug
69!>          - info
70!>          - warning
71!>          - error
72!>          - fatal
73!>          - none
74!>
75!>    - **in_maxerror** [@a 5]<br/>
76!>       maximum number of error allowed
77!>
78!> @subsection subcfg namcfg
79!>    the configuration sub-namelist parameters are :
80!>
81!>    - **cn_varcfg** [@a ./cfg/variable.cfg]<br/>
82!>       path to the variable configuration file.<br/>
83!>       the variable configuration file defines standard name,
84!>       default interpolation method, axis,...
85!>       to be used for some known variables.<br/>
86!>
87!>    - **cn_dimcfg** [@a ./cfg/dimension.cfg]<br/>
88!>       path to the dimension configuration file.<br/>
89!>       the dimension configuration file defines dimensions allowed.<br/>
90!>
91!>    - **cn_dumcfg** [@a ./cfg/dummy.cfg]<br/>
92!>       path to the useless (dummy) configuration file.<br/>
93!>       the dummy configuration file defines useless
94!>       dimension or variable. these dimension(s) or variable(s) will not be
95!>       processed.<br/>
96!>
97!> @subsection subcrs namcrs
98!>    the source/coarse grid sub-namelist parameters are :
99!>
100!>    - **cn_coord0** [@a ]<br/>
101!>       path to the coordinate file
102!>
103!>    - **in_perio0** [@a ]<br/>
104!>       NEMO periodicity index<br/>
105!>       the NEMO periodicity could be choose between 0 to 6:
106!>       <dl>
107!>          <dt>in_perio=0</dt>
108!>          <dd>standard regional model</dd>
109!>          <dt>in_perio=1</dt>
110!>          <dd>east-west cyclic model</dd>
111!>          <dt>in_perio=2</dt>
112!>          <dd>model with symmetric boundary condition across the equator</dd>
113!>          <dt>in_perio=3</dt>
114!>          <dd>regional model with North fold boundary and T-point pivot</dd>
115!>          <dt>in_perio=4</dt>
116!>          <dd>global model with a T-point pivot.<br/>
117!>          example: ORCA2, ORCA025, ORCA12</dd>
118!>          <dt>in_perio=5</dt>
119!>          <dd>regional model with North fold boundary and F-point pivot</dd>
120!>          <dt>in_perio=6</dt>
121!>          <dd>global model with a F-point pivot<br/>
122!>          example: ORCA05</dd>
123!>          </dd>
124!>       </dl>
125!>       @sa For more information see @ref md_src_docsrc_6_perio
126!>       and Model Boundary Condition paragraph in the
127!>       [NEMO documentation](https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/manual/pdf/NEMO_manual.pdf)
128!>
129!> @subsection subfin namfin
130!>    the target/fine grid sub-namelist parameters are :
131!>
132!>    - **cn_coord1** [@a ]<br/>
133!>       path to coordinate file
134!>
135!>    - **cn_bathy1** [@a ]<br/>
136!>       path to bathymetry file
137!>       @warning
138!>
139!>    - **in_perio1** [@a ]<br/>
140!>       NEMO periodicity index (see above)
141!>    @note if the fine/target coordinates file (cn_coord1) was created by SIREN, you do
142!>    not need to fill this parameter. SIREN will read it on the global attributes of
143!>    the coordinates file.
144!>
145!> @subsection subzgr namzgr
146!>    the vertical grid sub-namelist parameters are :
147!>
148!>    - **dn_pp_to_be_computed** [@a 0]<br/>
149!>
150!>    - **dn_ppsur** [@a -3958.951371276829]<br/>
151!>       coefficient to compute vertical grid
152!>
153!>    - **dn_ppa0** [@a 103.953009600000]<br/>
154!>       coefficient to compute vertical grid
155!>
156!>    - **dn_ppa1** [@a 2.415951269000]<br/>
157!>       coefficient to compute vertical grid
158!>
159!>    - **dn_ppa2** [@a 100.760928500000]<br/>
160!>       double tanh function parameter
161!>
162!>    - **dn_ppkth** [@a 15.351013700000]<br/>
163!>       coefficient to compute vertical grid
164!>
165!>    - **dn_ppkth2** [@a 48.029893720000]<br/>
166!>       double tanh function parameter
167!>
168!>    - **dn_ppacr** [@a 7.000000000000]<br/>
169!>       coefficient to compute vertical grid
170!>
171!>    - **dn_ppacr2** [@a 13.000000000000]<br/>
172!>       double tanh function parameter
173!>
174!>    - **dn_ppdzmin** [@a 6.]<br/>
175!>       minimum vertical spacing
176!>
177!>    - **dn_pphmax** [@a 5750.]<br/>
178!>       maximum depth
179!>
180!>    - **in_nlevel** [@a 75]<br/>
181!>       number of vertical level
182!>
183!>     @note
184!>       If *dn_ppa1*, *dn_ppa0* and *dn_ppsur* are undefined,
185!>       NEMO will compute them from *dn_ppdzmin, dn_pphmax, dn_ppkth, dn_ppacr*
186!>
187!> @subsection subzps namzps
188!>    the partial step sub-namelist parameters are :
189!>
190!>    - **dn_e3zps_min** [@a 25.]<br/>
191!>       minimum thickness of partial step level (meters)
192!>    - **dn_e3zps_rat** [@a 0.2]<br/>
193!>       minimum thickness ratio of partial step level
194!>
195!> @subsection subvar namvar
196!>    the variable sub-namelist parameters are :
197!>
198!>    - **cn_varfile** [@a ]<br/>
199!>       list of variable, and associated file
200!>
201!>       *cn_varfile* is the path and filename of the file where find
202!>       variable.
203!>       @note
204!>          *cn_varfile* could be a matrix of value, if you want to handwrite
205!>          variable value.<br/>
206!>          the variable array of value is split into equal subdomain.<br/>
207!>          each subdomain is filled with the corresponding value
208!>          of the matrix.<br/>
209!>          separators used to defined matrix are:
210!>             - ',' for line
211!>             - '/' for row
212!>             - '\' for level<br/>
213!>             Example:<br/>
214!>                3,2,3/1,4,5  =>  @f$ \left( \begin{array}{ccc}
215!>                                      3 & 2 & 3 \\
216!>                                      1 & 4 & 5 \end{array} \right) @f$
217!>
218!>          @warning
219!>             the same matrix is used for all boundaries.
220!>
221!>       Examples:
222!>          - 'votemper:gridT.nc', 'vozocrtx:gridU.nc'
223!>          - 'votemper:10\25', 'vozocrtx:gridU.nc'<br/>
224!>
225!>       @note
226!>          Optionnaly, NEMO periodicity could be added following the filename.
227!>          the periodicity must be separated by ';'
228!>
229!>       Example:
230!>          - 'votemper:gridT.nc ; perio=4'
231!>
232!>    - **cn_varinfo** [@a ]<br/>
233!>       list of variable and extra information about request(s) to be used<br/>
234!>
235!>       each elements of *cn_varinfo* is a string character (separated by ',').<br/>
236!>       it is composed of the variable name follow by ':',
237!>       then request(s) to be used on this variable.<br/>
238!>       request could be:
239!>          - int = interpolation method
240!>          - ext = extrapolation method
241!>          - flt = filter method
242!>          - min = minimum value
243!>          - max = maximum value
244!>          - unt = new units
245!>          - unf = unit scale factor (linked to new units)
246!>
247!>             requests must be separated by ';'.<br/>
248!>             order of requests does not matter.<br/>
249!>
250!>       informations about available method could be find in @ref interp,
251!>       @ref extrap and @ref filter modules.<br/>
252!>       Example:
253!>          - 'votemper: int=linear; flt=hann; ext=dist_weight',
254!>            'vosaline: int=cubic'
255!>
256!>       @note
257!>          If you do not specify a method which is required,
258!>          default one is apply.
259!>
260!> @subsection subnst namnst
261!>    the nesting sub-namelist parameters are :
262!>
263!>    - **in_rhoi**  [@a 1]<br/>
264!>       refinement factor in i-direction
265!>
266!>    - **in_rhoj**  [@a 1]<br/>
267!>       refinement factor in j-direction
268!>
269!>    @note
270!>       source/coarse grid indices will be deduced from target/fine grid
271!>       coordinate file.
272!>
273!> @subsection subbdy nambdy
274!>    the boundary sub-namelist parameters are :
275!>
276!>    - **ln_north** [@a .TRUE.]<br/>
277!>       logical to use north boundary or not
278!>    - **ln_south** [@a .TRUE.]<br/>
279!>       logical to use south boundary or not
280!>    - **ln_east**  [@a .TRUE.]<br/>
281!>       logical to use east boundary or not
282!>    - **ln_west**  [@a .TRUE.]<br/>
283!>       logical to use west  boundary or not
284!>    <br/> <br/>
285!>    - **cn_north** [@a ]<br/>
286!>       north boundary indices on target/fine grid<br/>
287!>    - **cn_south** [@a ]<br/>
288!>       south boundary indices on target/fine grid<br/>
289!>    - **cn_east**  [@a ]<br/>
290!>       east  boundary indices on target/fine grid<br/>
291!>    - **cn_west**  [@a ]<br/>
292!>       west  boundary indices on target/fine grid<br/>
293!>
294!>       *cn_north* is a string character defining boundary
295!>       segmentation.<br/>
296!>       segments are separated by '|'.<br/>
297!>       each segments of the boundary is composed of:
298!>          - indice of velocity (orthogonal to boundary .ie.
299!>             for north boundary, J-indice).
300!>          - indice of segment start (I-indice for north boundary)
301!>          - indice of segment end   (I-indice for north boundary)<br/>
302!>             indices must be separated by ':' .<br/>
303!>          - optionally, boundary size could be added between '(' and ')'
304!>          in the first segment defined.
305!>             @note
306!>                boundary size is the same for all segments of one boundary.
307!>
308!>       Examples:
309!>          - cn_north='index1,first1:last1(width)'
310!>          - cn_north='index1(width),first1:last1|index2,first2:last2'
311!>
312!>       @image html  boundary_50.png
313!>       <center>@image latex boundary_50.png
314!>       </center>
315!>
316!>    - **ln_oneseg** [@a .TRUE.]<br/>
317!>       logical to use only one segment for each boundary or not
318!>
319!>    @note
320!>       the number of point(s) with source/coarse value save at boundaries is
321!>       defined with the *weight* variable (see @ref merge_bathy)
322!>
323!> @subsection subout namout
324!>    the output sub-namelist parameter is :
325!>
326!>    - **cn_fileout** [@a boundary.nc]<br/>
327!>       output bathymetry filename
328!>
329!>       @note
330!>          cardinal point and segment number will be automatically added
331!>
332!>    - **ln_extrap** [@a .FALSE.]<br/>
333!>       extrapolate on land point
334!>
335!>    - **dn_dayofs** [@a 0]<br/>
336!>       date offset in day (change only ouput file name)
337!>
338!>       Examples:
339!>          - cn_fileout='boundary.nc'<br/>
340!>             if time_counter (16/07/2015 00h) is read on input file (see varfile),
341!>             west boundary will be named boundary_west_y2015m07d16
342!>          - dn_dayofs=-2.<br/>
343!>             if you use day offset you get boundary_west_y2015m07d14
344!>
345!> @subsection sub_nambdy How to fill Lateral Boundary Condition in NEMO namelist
346!>    To use boundary condition within NEMO, you need to fill the NEMO namelist.<br/>
347!>    As this is a little bit messy for lateral boundary condition, here after
348!>    is an explanation of how to do it.
349!>
350!>    This will be done in 3 steps.
351!>
352!>    @subsubsection ss_nambdy nambdy
353!>       The *nambdy* NEMO sub-namelist defines open boundaries.<br/>
354!>       Here we indicate the number of open boundary (**nb_bdy**).
355!>
356!>       @note
357!>          we have to fill most of the parameters with as many elements as there are open boundaries
358!>
359!>       Regarding the width of the relaxation zone **nn_rimwidth**,
360!>       this information is available as a global attribute (**bdy_width**)
361!>       in the metadata of boundary files created with SIREN
362!>
363!> @code{.sh}
364!>    ncdump -h boundary_east.nc
365!> @endcode
366!>       @warning
367!>          The order of the boundaries must stay unchanged, in parameters list as well as
368!>          in the next sub-namelsits
369!>
370!>    Example:<br/>
371!>       here is an example for a domain with two boundaries East and North
372!>
373!> @code{.sh}
374!> !-----------------------------------------------------------------------
375!> &nambdy        !  unstructured open boundaries                          ("key_bdy")
376!> !-----------------------------------------------------------------------
377!>   nb_bdy         = 2                    !  number of open boundary sets
378!>   ln_coords_file = .false.,.false.      !  =T : read bdy coordinates from file
379!>   cn_coords_file = '',''                !  bdy coordinates files
380!>   ln_mask_file   = .false.              !  =T : read mask from file
381!>   cn_mask_file   = ''                   !  name of mask file (if ln_mask_file=.TRUE.)
382!>   cn_dyn2d       = 'flather','flather'  !
383!>   nn_dyn2d_dta   = 1,1                  !  = 0, bdy data are equal to the initial state
384!>                                         !  = 1, bdy data are read in 'bdydata   .nc' files
385!>                                         !  = 2, use tidal harmonic forcing data from files
386!>                                         !  = 3, use external data AND tidal harmonic forcing
387!>   cn_dyn3d       = 'specified','specified' !
388!>   nn_dyn3d_dta   = 1,1                  !  = 0, bdy data are equal to the initial state
389!>                                         !  = 1, bdy data are read in 'bdydata   .nc' files
390!>   cn_tra         = 'specified','specified' !
391!>   nn_tra_dta     = 1,1                  !  = 0, bdy data are equal to the initial state
392!>                                         !  = 1, bdy data are read in 'bdydata   .nc' files
393!>                                         !
394!>   ln_tra_dmp    =.true.,.true.          !  open boudaries conditions for tracers
395!>   ln_dyn3d_dmp  =.true.,.true.          !  open boundary condition for baroclinic velocities
396!>   rn_time_dmp     =  1.,1.              ! Damping time scale in days
397!>   rn_time_dmp_out =  1.,1.              ! Outflow damping time scale
398!>   nn_rimwidth   = 10,10                 !  width of the relaxation zone
399!>   ln_vol        = .false.               !  total volume correction (see nn_volctl parameter)
400!>   nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero
401!> /
402!> @endcode
403!>
404!>    @subsubsection ss_nambdy_index nambdy_index
405!>       The *nambdy_index* NEMO sub-namelist describes the boundaries we will use.
406!>
407!>       @warning
408!>          We have to add as many as sub namelist *nambdy_index* than open boundaries (nb_bdy),
409!>          and keep them in the same order as above
410!>
411!>          Here we indicate if the open boundary is North, South, East, or West (**ctypebdy**).<br/>
412!>          We also indicate indice of segment start and end (respectively **nbdybeg**  and **nbdyend**)
413!>          as well as indice of velocity row or column (**nbdyind**).<br/>
414!>
415!>          Those informations are available as global attributes
416!>          (respectively **bdy_deb, bdy_end, bdy_ind**) in the metadata of our boundary files
417!>          created with SIREN.
418!>
419!>    Example:<br/>
420!>       here is an example for a domain with two boundaries East and North
421!>
422!> @code{.sh}
423!>    !-----------------------------------------------------------------------
424!>    &nambdy_index  !  structured open boundaries definition     ("key_bdy")
425!>    !-----------------------------------------------------------------------
426!>      ctypebdy ='E'                   ! Open boundary type (W,E,S or N)
427!>      nbdyind  = 407                  ! indice of velocity row or column
428!>                                      ! if ==-1, set obc at the domain boundary
429!>                                      !        , discard start and end indices
430!>      nbdybeg  = 32                   ! indice of segment start
431!>      nbdyend  = 300                  ! indice of segment end
432!>    /
433!>    !-----------------------------------------------------------------------
434!>    &nambdy_index  !  structured open boundaries definition     ("key_bdy")
435!>    !-----------------------------------------------------------------------
436!>      ctypebdy ='N'                   ! Open boundary type (W,E,S or N)
437!>      nbdyind  = 299                  ! indice of velocity row or column
438!>                                      ! if ==-1, set obc at the domain boundary
439!>                                      !        , discard start and end indices
440!>      nbdybeg  = 200                  ! indice of segment start
441!>      nbdyend  = 408                  ! indice of segment end
442!>    /
443!> @endcode
444!>
445!>    @subsubsection ss_nambdy_dat nambdy_dta
446!>       The *nambdy_dta* NEMO sub-namelists describes the boundary data and files to be used.<br/>
447!>       @warning
448!>          We have to add as many as sub namelist *nambdy_dta* than open boundaries (nb_bdy),
449!>          and keep them in the same order as above
450!>
451!>    Example:<br/>
452!>       here is an example for a domain with two boundaries East and North
453!>
454!> @code{.sh}
455!> !-----------------------------------------------------------------------
456!> &nambdy_dta      !  open boundaries - external data           ("key_bdy")
457!> !-----------------------------------------------------------------------
458!> !          ! file name !  freq (hours)   ! variable ! time interp. ! clim  ! 'yearly'/ ! weights...
459!> !          !           ! (if < 0 months) !   name   !   (logical)  ! (T/F) ! 'monthly' ! filename
460!>   bn_ssh = 'boundary_east' , -12   , 'sossheig' , .false. , .true. , 'yearly' , '', '', ''
461!>   bn_u2d = 'boundary_east' , -12   , 'vobtcrtx' , .false. , .true. , 'yearly' , '', '', ''
462!>   bn_v2d = 'boundary_east' , -12   , 'vobtcrty' , .false. , .true. , 'yearly' , '', '', ''
463!>   bn_u3d = 'boundary_east' , -12   , 'vozocrtx' , .false. , .true. , 'yearly' , '', '', ''
464!>   bn_v3d = 'boundary_east' , -12   , 'vomecrty' , .false. , .true. , 'yearly' , '', '', ''
465!>   bn_tem = 'boundary_east' , -12   , 'votemper' , .false. , .true. , 'yearly' , '', '', ''
466!>   bn_sal = 'boundary_east' , -12   , 'vosaline' , .false. , .true. , 'yearly' , '', '', ''
467!>   cn_dir = './'
468!>   ln_full_vel = .true.
469!> /
470!> !-----------------------------------------------------------------------
471!> &nambdy_dta      !  open boundaries - external data           ("key_bdy")
472!> !-----------------------------------------------------------------------
473!> !          ! file name !  freq (hours)   ! variable ! time interp. ! clim  ! 'yearly'/ ! weights...
474!> !          !           ! (if < 0 months) !   name   !   (logical)  ! (T/F) ! 'monthly' ! filename
475!>   bn_ssh = 'boundary_north' , -12   , 'sossheig' ,  .false. , .true. , 'yearly' , '', '', ''
476!>   bn_u2d = 'boundary_north' , -12   , 'vobtcrtx' ,  .false. , .true. , 'yearly' , '', '', ''
477!>   bn_v2d = 'boundary_north' , -12   , 'vobtcrty' ,  .false. , .true. , 'yearly' , '', '', ''
478!>   bn_u3d = 'boundary_north' , -12   , 'vozocrtx' ,  .false. , .true. , 'yearly' , '', '', ''
479!>   bn_v3d = 'boundary_north' , -12   , 'vomecrty' ,  .false. , .true. , 'yearly' , '', '', ''
480!>   bn_tem = 'boundary_north' , -12   , 'votemper' ,  .false. , .true. , 'yearly' , '', '', ''
481!>   bn_sal = 'boundary_north' , -12   , 'vosaline' ,  .false. , .true. , 'yearly' , '', '', ''
482!>   cn_dir = './'
483!>   ln_full_vel = .true.
484!> /
485!> @endcode
486!>
487!> <hr>
488!> @author J.Paul
489!>
490!> @date November, 2013 - Initial Version
491!> @date September, 2014
492!> - add header for user
493!> - take into account grid point to compue boundaries
494!> - reorder output dimension for north and south boundaries
495!> @date June, 2015
496!> - extrapolate all land points, and add ln_extrap in namelist.
497!> - allow to change unit.
498!> @date July, 2015
499!> - add namelist parameter to shift date of output file name.
500!> @date September, 2015
501!> - manage useless (dummy) variable, attributes, and dimension
502!> - allow to run on multi processors with key_mpp_mpi
503!> @date January, 2016
504!> - same process use for variable extracted or interpolated from input file.
505!> @date October, 2016
506!> - dimension to be used select from configuration file
507!> @date January, 2019
508!> - add url path to global attributes of output file(s)
509!> - create and clean file structure to avoid memory leaks
510!> - explain how to fill Lateral Boundary Condition in NEMO namelist
511!> @date February, 2019
512!> - rename sub namelist namcrs to namsrc
513!> - rename sub namelist namfin to namtgt
514!> @date August, 2019
515!> - use periodicity read from namelist, and store in multi structure
516!> @date Ocober, 2019
517!> - add help and version optional arguments
518!>
519!> @todo
520!> - rewitre using meshmask instead of bathymetry and coordinates files.
521!>
522!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
523!----------------------------------------------------------------------
524PROGRAM create_boundary
525
526   USE netcdf                          ! nf90 library
527   USE global                          ! global variable
528   USE phycst                          ! physical constant
529   USE kind                            ! F90 kind parameter
530   USE fct                             ! basic useful function
531   USE date                            ! date manager
532   USE att                             ! attribute manager
533   USE dim                             ! dimension manager
534   USE var                             ! variable manager
535   USE file                            ! file manager
536   USE multi                           ! multi file manager
537   USE boundary                        ! boundary manager
538   USE iom                             ! I/O manager
539   USE dom                             ! domain manager
540   USE grid                            ! grid manager
541   USE vgrid                           ! vertical grid manager
542   USE extrap                          ! extrapolation manager
543   USE interp                          ! interpolation manager
544   USE filter                          ! filter manager
545   USE mpp                             ! MPP manager
546   USE iom_mpp                         ! MPP I/O manager
547
548   IMPLICIT NONE
549
550   ! parameters
551   CHARACTER(LEN=lc), PARAMETER  :: cp_myname = "create_boundary"
552
553   ! local variable
554   CHARACTER(LEN=lc)                                  :: cl_arg
555   CHARACTER(LEN=lc)                                  :: cl_errormsg
556
557   INTEGER(i4)                                        :: il_narg
558
559#if defined key_mpp_mpi
560   ! mpp variable
561   CHARACTER(LEN=lc), DIMENSION(:)      , ALLOCATABLE :: cl_args
562
563   INTEGER(i4)                                        :: ierror
564   INTEGER(i4)                                        :: iproc
565   INTEGER(i4)                                        :: nproc
566   INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_nprog
567
568   ! loop indices
569   INTEGER(i4) :: jm
570#else
571   CHARACTER(LEN=lc)                                  :: cl_namelist
572#endif
573   !-------------------------------------------------------------------
574#if defined key_mpp_mpi
575   INCLUDE 'mpif.h'
576#endif
577   !-------------------------------------------------------------------
578
579   !
580   ! Initialisation
581   ! --------------
582   !
583   il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
584
585#if ! defined key_mpp_mpi
586
587   ! Traitement des arguments fournis
588   ! --------------------------------
589   IF( il_narg /= 1 )THEN
590      WRITE(cl_errormsg,*) ' ERROR : one argument is needed '
591      CALL fct_help(cp_myname,cl_errormsg)
592      CALL EXIT(1)
593   ELSE
594
595      CALL GET_COMMAND_ARGUMENT(1,cl_arg) !f03 intrinsec
596      SELECT CASE (cl_arg)
597         CASE ('-v', '--version')
598
599            CALL fct_version(cp_myname)
600            CALL EXIT(0)
601
602         CASE ('-h', '--help')
603
604            CALL fct_help(cp_myname)
605            CALL EXIT(0)
606
607         CASE DEFAULT
608
609            cl_namelist=cl_arg
610
611            CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
612            CALL create_boundary__mono(cl_namelist)
613
614      END SELECT
615   ENDIF
616#else
617
618   ! Initialize MPI
619   CALL mpi_init(ierror)
620   CALL mpi_comm_rank(mpi_comm_world,iproc,ierror)
621   CALL mpi_comm_size(mpi_comm_world,nproc,ierror)
622
623   ! Traitement des arguments fournis
624   ! --------------------------------
625   IF( il_narg == 0 )THEN
626      WRITE(cl_errormsg,*) ' ERROR : at least one argument is needed '
627      CALL fct_help(cp_myname,cl_errormsg)
628      CALL EXIT(1)
629   ELSE
630
631      ALLOCATE(cl_args(il_narg))
632      DO jm=1,il_narg
633
634         CALL GET_COMMAND_ARGUMENT(jm,cl_arg) !f03 intrinsec
635         SELECT CASE (cl_arg)
636            CASE ('-v', '--version')
637
638               CALL fct_version(cp_myname)
639               CALL EXIT(0)
640
641            CASE ('-h', '--help')
642
643               CALL fct_help(cp_myname)
644               CALL EXIT(0)
645
646            CASE DEFAULT
647
648               cl_args(jm)=TRIM(cl_arg)
649
650         END SELECT
651      ENDDO
652   ENDIF
653
654   ALLOCATE(il_nprog(il_narg))
655   DO jm=1, il_narg
656      il_nprog(jm)= MOD(jm,nproc)
657   ENDDO
658
659   DO jm=1, il_narg
660      IF ( il_nprog(jm) .eq. iproc ) THEN
661         CALL create_boundary__mono(cl_args(jm))
662      ENDIF
663   ENDDO
664
665   CALL mpi_finalize(ierror)
666
667   DEALLOCATE(cl_args)
668   DEALLOCATE(il_nprog)
669#endif
670
671CONTAINS
672   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
673   SUBROUTINE create_boundary__mono(cd_namelist)
674   !-------------------------------------------------------------------
675   !> @brief
676   !> This subroutine create boundary files.
677   !>
678   !> @details
679   !>
680   !> @author J.Paul
681   !> @date January, 2016 - Initial Version
682   !>
683   !> @param[in] cd_namelist namelist file
684   !-------------------------------------------------------------------
685
686   USE logger                          ! log file manager
687
688   IMPLICIT NONE
689   ! Argument
690   CHARACTER(LEN=lc), INTENT(IN) :: cd_namelist
691
692   ! local variable
693   CHARACTER(LEN=lc)                                  :: cl_date
694   CHARACTER(LEN=lc)                                  :: cl_name
695   CHARACTER(LEN=lc)                                  :: cl_bdyout
696   CHARACTER(LEN=lc)                                  :: cl_data
697   CHARACTER(LEN=lc)                                  :: cl_dimorder
698   CHARACTER(LEN=lc)                                  :: cl_fmt
699   CHARACTER(LEN=lc)                                  :: cl_url
700
701   INTEGER(i4)                                        :: il_status
702   INTEGER(i4)                                        :: il_fileid
703   INTEGER(i4)                                        :: il_imin0
704   INTEGER(i4)                                        :: il_imax0
705   INTEGER(i4)                                        :: il_jmin0
706   INTEGER(i4)                                        :: il_jmax0
707   INTEGER(i4)                                        :: il_shift
708   INTEGER(i4)      , DIMENSION(ip_maxdim)            :: il_rho
709   INTEGER(i4)      , DIMENSION(2,2)                  :: il_offset
710   INTEGER(i4)      , DIMENSION(2,2)                  :: il_ind
711
712   LOGICAL                                            :: ll_exist
713
714   TYPE(TATT)                                         :: tl_att
715
716   TYPE(TVAR)                                         :: tl_depth
717   TYPE(TVAR)                                         :: tl_time
718   TYPE(TVAR)                                         :: tl_var1
719   TYPE(TVAR)                                         :: tl_var0
720   TYPE(TVAR)                                         :: tl_lon1
721   TYPE(TVAR)                                         :: tl_lat1
722   TYPE(TVAR)                                         :: tl_lvl1
723   TYPE(TVAR)       , DIMENSION(:)      , ALLOCATABLE :: tl_level
724   TYPE(TVAR)       , DIMENSION(:,:,:)  , ALLOCATABLE :: tl_seglvl1
725   TYPE(TVAR)       , DIMENSION(:,:,:)  , ALLOCATABLE :: tl_segvar1
726
727   TYPE(TDIM)       , DIMENSION(ip_maxdim)            :: tl_dim
728
729   TYPE(TDATE)                                        :: tl_date
730
731   TYPE(TBDY)       , DIMENSION(ip_ncard)             :: tl_bdy
732
733   TYPE(TDOM)                                         :: tl_dom0
734   TYPE(TDOM)                                         :: tl_dom1
735   TYPE(TDOM)       , DIMENSION(:,:,:)  , ALLOCATABLE :: tl_segdom1
736
737   TYPE(TFILE)                                        :: tl_file
738   TYPE(TFILE)                                        :: tl_fileout
739
740   TYPE(TMPP)                                         :: tl_coord0
741   TYPE(TMPP)                                         :: tl_coord1
742   TYPE(TMPP)                                         :: tl_bathy1
743   TYPE(TMPP)                                         :: tl_mpp
744
745   TYPE(TMULTI)                                       :: tl_multi
746
747   ! loop indices
748   INTEGER(i4) :: jvar
749   INTEGER(i4) :: jpoint
750   INTEGER(i4) :: ji
751   INTEGER(i4) :: jj
752   INTEGER(i4) :: jk
753   INTEGER(i4) :: jl
754
755   ! namelist variable
756   ! namlog
757   CHARACTER(LEN=lc)                       :: cn_logfile    = 'create_boundary.log'
758   CHARACTER(LEN=lc)                       :: cn_verbosity  = 'warning'
759   INTEGER(i4)                             :: in_maxerror   = 5
760
761   ! namcfg
762   CHARACTER(LEN=lc)                       :: cn_varcfg  = './cfg/variable.cfg'
763   CHARACTER(LEN=lc)                       :: cn_dimcfg  = './cfg/dimension.cfg'
764   CHARACTER(LEN=lc)                       :: cn_dumcfg  = './cfg/dummy.cfg'
765
766   ! namsrc
767   CHARACTER(LEN=lc)                       :: cn_coord0  = ''
768   INTEGER(i4)                             :: in_perio0  = -1
769
770   ! namtgt
771   CHARACTER(LEN=lc)                       :: cn_coord1  = ''
772   CHARACTER(LEN=lc)                       :: cn_bathy1  = ''
773   INTEGER(i4)                             :: in_perio1  = -1
774
775   !namzgr
776   REAL(dp)                                :: dn_pp_to_be_computed = 0._dp
777   REAL(dp)                                :: dn_ppsur   = -3958.951371276829_dp
778   REAL(dp)                                :: dn_ppa0    =   103.953009600000_dp
779   REAL(dp)                                :: dn_ppa1    =     2.415951269000_dp
780   REAL(dp)                                :: dn_ppa2    =   100.760928500000_dp
781   REAL(dp)                                :: dn_ppkth   =    15.351013700000_dp
782   REAL(dp)                                :: dn_ppkth2  =    48.029893720000_dp
783   REAL(dp)                                :: dn_ppacr   =     7.000000000000_dp
784   REAL(dp)                                :: dn_ppacr2  =    13.000000000000_dp
785   REAL(dp)                                :: dn_ppdzmin = 6._dp
786   REAL(dp)                                :: dn_pphmax  = 5750._dp
787   INTEGER(i4)                             :: in_nlevel  = 75
788
789   !namzps
790   REAL(dp)                                :: dn_e3zps_min  = 25._dp
791   REAL(dp)                                :: dn_e3zps_rat  = 0.2_dp
792
793   ! namvar
794   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = ''
795   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
796
797   ! namnst
798   INTEGER(i4)                             :: in_rhoi    = 1
799   INTEGER(i4)                             :: in_rhoj    = 1
800
801   ! nambdy
802   LOGICAL                                 :: ln_north   = .TRUE.
803   LOGICAL                                 :: ln_south   = .TRUE.
804   LOGICAL                                 :: ln_east    = .TRUE.
805   LOGICAL                                 :: ln_west    = .TRUE.
806   LOGICAL                                 :: ln_oneseg  = .TRUE.
807   CHARACTER(LEN=lc)                       :: cn_north   = ''
808   CHARACTER(LEN=lc)                       :: cn_south   = ''
809   CHARACTER(LEN=lc)                       :: cn_east    = ''
810   CHARACTER(LEN=lc)                       :: cn_west    = ''
811
812   ! namout
813   CHARACTER(LEN=lc)                       :: cn_fileout = 'boundary.nc'
814   REAL(dp)                                :: dn_dayofs  = 0._dp
815   LOGICAL                                 :: ln_extrap  = .FALSE.
816   !-------------------------------------------------------------------
817
818   NAMELIST /namlog/ &  !< logger namelist
819   &  cn_logfile,    &  !< log file
820   &  cn_verbosity,  &  !< log verbosity
821   &  in_maxerror       !< logger maximum error
822
823   NAMELIST /namcfg/ &  !< configuration namelist
824   &  cn_varcfg,     &  !< variable configuration file
825   &  cn_dimcfg,     &  !< dimension configuration file
826   &  cn_dumcfg         !< dummy configuration file
827
828   NAMELIST /namsrc/ &  !< source/coarse grid namelist
829   &  cn_coord0,     &  !< coordinate file
830   &  in_perio0         !< periodicity index
831
832   NAMELIST /namtgt/ &  !< target/fine grid namelist
833   &  cn_coord1,     &  !< coordinate file
834   &  cn_bathy1,     &  !< bathymetry file
835   &  in_perio1         !< periodicity index
836
837   NAMELIST /namzgr/ &
838   &  dn_pp_to_be_computed, &
839   &  dn_ppsur,      &
840   &  dn_ppa0,       &
841   &  dn_ppa1,       &
842   &  dn_ppa2,       &
843   &  dn_ppkth,      &
844   &  dn_ppkth2,     &
845   &  dn_ppacr,      &
846   &  dn_ppacr2,     &
847   &  dn_ppdzmin,    &
848   &  dn_pphmax,     &
849   &  in_nlevel         !< number of vertical level
850
851   NAMELIST /namzps/ &
852   &  dn_e3zps_min,  &
853   &  dn_e3zps_rat
854
855   NAMELIST /namvar/ &  !< variable namelist
856   &  cn_varfile,    &  !< list of variable and file where find it. (ex: 'votemper:GLORYS_gridT.nc' )
857   &  cn_varinfo        !< list of variable and method to apply on. (ex: 'votemper:linear','vosaline:cubic' )
858
859   NAMELIST /namnst/ &  !< nesting namelist
860   &  in_rhoi,       &  !< refinement factor in i-direction
861   &  in_rhoj           !< refinement factor in j-direction
862
863   NAMELIST /nambdy/ &  !< boundary namelist
864   &  ln_north,      &  !< use north boundary
865   &  ln_south,      &  !< use south boundary
866   &  ln_east ,      &  !< use east  boundary
867   &  ln_west ,      &  !< use west  boundary
868   &  cn_north,      &  !< north boundary indices on fine grid
869   &  cn_south,      &  !< south boundary indices on fine grid
870   &  cn_east ,      &  !< east  boundary indices on fine grid
871   &  cn_west ,      &  !< west  boundary indices on fine grid
872   &  ln_oneseg         !< use only one segment for each boundary or not
873
874   NAMELIST /namout/ &  !< output namelist
875   &  cn_fileout,    &  !< fine grid boundary file basename
876   &  dn_dayofs,     &  !< date offset in day (change only ouput file name)
877   &  ln_extrap         !< extrapolate or not
878   !-------------------------------------------------------------------
879
880   ! read namelist
881   INQUIRE(FILE=TRIM(cd_namelist), EXIST=ll_exist)
882
883   IF( ll_exist )THEN
884
885      il_fileid=fct_getunit()
886
887      OPEN( il_fileid, FILE=TRIM(cd_namelist), &
888      &                FORM='FORMATTED',       &
889      &                ACCESS='SEQUENTIAL',    &
890      &                STATUS='OLD',           &
891      &                ACTION='READ',          &
892      &                IOSTAT=il_status)
893      CALL fct_err(il_status)
894      IF( il_status /= 0 )THEN
895         PRINT *,"CREATE BOUNDARY: ERROR opening "//TRIM(cd_namelist)
896         STOP
897      ENDIF
898
899      READ( il_fileid, NML = namlog )
900      ! define log file
901      CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
902      CALL logger_header()
903
904      READ( il_fileid, NML = namcfg )
905      ! get variable extra information
906      CALL var_def_extra(TRIM(cn_varcfg))
907
908      ! get dimension allowed
909      CALL dim_def_extra(TRIM(cn_dimcfg))
910
911      ! get dummy variable
912      CALL var_get_dummy(TRIM(cn_dumcfg))
913      ! get dummy dimension
914      CALL dim_get_dummy(TRIM(cn_dumcfg))
915      ! get dummy attribute
916      CALL att_get_dummy(TRIM(cn_dumcfg))
917
918      READ( il_fileid, NML = namsrc )
919      READ( il_fileid, NML = namtgt )
920      READ( il_fileid, NML = namzgr )
921      READ( il_fileid, NML = namvar )
922      ! add user change in extra information
923      CALL var_chg_extra(cn_varinfo)
924      ! match variable with file
925      tl_multi=multi_init(cn_varfile)
926
927      READ( il_fileid, NML = namnst )
928      READ( il_fileid, NML = nambdy )
929      READ( il_fileid, NML = namout )
930
931      CLOSE( il_fileid, IOSTAT=il_status )
932      CALL fct_err(il_status)
933      IF( il_status /= 0 )THEN
934         CALL logger_error("CREATE BOUNDARY: ERROR closing "//TRIM(cd_namelist))
935      ENDIF
936
937   ELSE
938
939      WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cd_namelist)
940      CALL fct_help(cp_myname,cl_errormsg)
941      CALL EXIT(1)
942
943   ENDIF
944
945   CALL multi_print(tl_multi)
946
947   IF( tl_multi%i_nvar <= 0 )THEN
948      CALL logger_fatal("CREATE BOUNDARY: no variable to be used."//&
949      &  " check namelist.")
950   ENDIF
951
952   ! open files
953   IF( TRIM(cn_coord0) /= '' )THEN
954      tl_file= file_init(TRIM(cn_coord0))
955      tl_coord0=mpp_init( tl_file, id_perio=in_perio0)
956      ! clean
957      CALL file_clean(tl_file)
958      CALL grid_get_info(tl_coord0)
959   ELSE
960      CALL logger_fatal("CREATE BOUNDARY: can not find source/coarse grid "//&
961      &  "coordinate file. check namelist")
962   ENDIF
963
964   IF( TRIM(cn_coord1) /= '' )THEN
965      tl_file=file_init(TRIM(cn_coord1))
966      tl_coord1=mpp_init( tl_file, id_perio=in_perio1)
967      ! clean
968      CALL file_clean(tl_file)
969      CALL grid_get_info(tl_coord1)
970   ELSE
971      CALL logger_fatal("CREATE BOUNDARY: can not find fine grid coordinate "//&
972      &  "file. check namelist")
973   ENDIF
974
975   IF( TRIM(cn_bathy1) /= '' )THEN
976      tl_file=file_init(TRIM(cn_bathy1))
977      tl_bathy1=mpp_init( tl_file, id_perio=in_perio1)
978      ! clean
979      CALL file_clean(tl_file)
980      CALL grid_get_info(tl_bathy1)
981   ELSE
982      CALL logger_fatal("CREATE BOUNDARY: can not find fine grid bathymetry "//&
983      &  "file. check namelist")
984   ENDIF
985
986   ! check
987   ! check output file do not already exist
988   ! WARNING: do not work when use time to create output file name
989   DO jk=1,ip_ncard
990      cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
991      &                                TRIM(cp_card(jk)), 1 )
992      INQUIRE(FILE=TRIM(cl_bdyout), EXIST=ll_exist)
993      IF( ll_exist )THEN
994         CALL logger_fatal("CREATE BOUNDARY: output file "//TRIM(cl_bdyout)//&
995         &  " already exist.")
996      ENDIF
997
998      cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
999         &                             TRIM(cp_card(jk)) )
1000      INQUIRE(FILE=TRIM(cl_bdyout), EXIST=ll_exist)
1001      IF( ll_exist )THEN
1002         CALL logger_fatal("CREATE BOUNDARY: output file "//TRIM(cl_bdyout)//&
1003         &  " already exist.")
1004      ENDIF
1005   ENDDO
1006
1007   ! check namelist
1008   ! check refinement factor
1009   il_rho(:)=1
1010   IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
1011      CALL logger_error("CREATE BOUNDARY: invalid refinement factor."//&
1012         &  " check namelist "//TRIM(cd_namelist))
1013   ELSE
1014      il_rho(jp_I)=in_rhoi
1015      il_rho(jp_J)=in_rhoj
1016   ENDIF
1017
1018   !
1019   ! compute coarse grid indices around fine grid
1020   il_ind(:,:)=grid_get_coarse_index(tl_coord0, tl_coord1, &
1021      &                              id_rho=il_rho(:))
1022
1023   il_imin0=il_ind(1,1) ; il_imax0=il_ind(1,2)
1024   il_jmin0=il_ind(2,1) ; il_jmax0=il_ind(2,2)
1025
1026   ! check domain validity
1027   CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0)
1028
1029   ! check coordinate file
1030   CALL grid_check_coincidence( tl_coord0, tl_coord1, &
1031      &                         il_imin0, il_imax0, &
1032      &                         il_jmin0, il_jmax0, &
1033      &                         il_rho(:) )
1034
1035   ! read or compute boundary
1036   CALL mpp_get_contour(tl_bathy1)
1037
1038   CALL iom_mpp_open(tl_bathy1)
1039
1040   tl_var1=iom_mpp_read_var(tl_bathy1,'Bathymetry')
1041
1042   CALL iom_mpp_close(tl_bathy1)
1043
1044   ! get boundaries indices
1045   tl_bdy(:)=boundary_init(tl_var1, ln_north, ln_south, ln_east, ln_west, &
1046      &                             cn_north, cn_south, cn_east, cn_west, &
1047      &                             ln_oneseg )
1048
1049
1050   CALL var_clean(tl_var1)
1051
1052   ! compute level
1053   ALLOCATE(tl_level(ip_npoint))
1054   tl_level(:)=vgrid_get_level(tl_bathy1, cd_namelist )
1055
1056   ! get coordinate for each segment of each boundary
1057   ALLOCATE( tl_segdom1(ip_npoint,ip_maxseg,ip_ncard) )
1058   ALLOCATE( tl_seglvl1(ip_npoint,ip_maxseg,ip_ncard) )
1059
1060   DO jl=1,ip_ncard
1061      IF( tl_bdy(jl)%l_use )THEN
1062         DO jk=1,tl_bdy(jl)%i_nseg
1063
1064            ! get fine grid segment domain
1065            tl_segdom1(:,jk,jl)=create_boundary_get_dom( tl_bathy1, &
1066               &                                         tl_bdy(jl), jk )
1067
1068            IF( .NOT. ln_extrap )THEN
1069               ! get fine grid level
1070               tl_seglvl1(:,jk,jl)= &
1071                  & create_boundary_get_level( tl_level(:), &
1072                  &                            tl_segdom1(:,jk,jl))
1073            ENDIF
1074
1075            ! add extra band to fine grid domain (if possible)
1076            ! to avoid dimension of one and so be able to compute offset
1077            DO jj=1,ip_npoint
1078               CALL dom_add_extra(tl_segdom1(jj,jk,jl), &
1079                  &               il_rho(jp_I), il_rho(jp_J))
1080            ENDDO
1081
1082         ENDDO
1083      ENDIF
1084   ENDDO
1085
1086   ! clean
1087   CALL var_clean(tl_level(:))
1088   DEALLOCATE(tl_level)
1089
1090   ! clean bathy
1091   CALL mpp_clean(tl_bathy1)
1092
1093   ALLOCATE( tl_segvar1(tl_multi%i_nvar,ip_maxseg,ip_ncard) )
1094   ! compute boundary for variable to be used (see namelist)
1095   IF( .NOT. ASSOCIATED(tl_multi%t_mpp) )THEN
1096      CALL logger_error("CREATE BOUNDARY: no file to work on. "//&
1097         &              "check cn_varfile in namelist.")
1098   ELSE
1099
1100      jvar=0
1101      ! for each file
1102      DO ji=1,tl_multi%i_nmpp
1103
1104         WRITE(cl_data,'(a,i2.2)') 'data-',jvar+1
1105         IF( .NOT. ASSOCIATED(tl_multi%t_mpp(ji)%t_proc(1)%t_var) )THEN
1106
1107            CALL logger_error("CREATE BOUNDARY: no variable to work on for "//&
1108               &              "mpp "//TRIM(tl_multi%t_mpp(ji)%c_name)//&
1109               &              ". check cn_varfile in namelist.")
1110
1111         ELSEIF( TRIM(tl_multi%t_mpp(ji)%c_name) == TRIM(cl_data) )THEN
1112         !- use input matrix to fill variable
1113
1114            WRITE(*,'(a)') "work on data"
1115            ! for each variable initialise from matrix
1116            DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
1117
1118               jvar=jvar+1
1119               WRITE(*,'(2x,a,a)') "work on variable "//&
1120               &  TRIM(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name)
1121
1122               tl_var1=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
1123
1124               SELECT CASE(TRIM(tl_var1%c_point))
1125               CASE DEFAULT !'T'
1126                  jpoint=jp_T
1127               CASE('U')
1128                  jpoint=jp_U
1129               CASE('V')
1130                  jpoint=jp_V
1131               CASE('F')
1132                  jpoint=jp_F
1133               END SELECT
1134
1135               WRITE(*,'(4x,a,a)') 'work on '//TRIM(tl_var1%c_name)
1136               DO jl=1,ip_ncard
1137                  IF( tl_bdy(jl)%l_use )THEN
1138
1139                     DO jk=1,tl_bdy(jl)%i_nseg
1140
1141                        ! fill value with matrix data
1142                        tl_segvar1(jvar,jk,jl)=create_boundary_matrix( &
1143                           &                          tl_var1, &
1144                           &                          tl_segdom1(jpoint,jk,jl), &
1145                           &                          in_nlevel )
1146
1147                        !del extra
1148                        CALL dom_del_extra( tl_segvar1(jvar,jk,jl), &
1149                           &                tl_segdom1(jpoint,jk,jl) )
1150
1151                     ENDDO
1152
1153                  ENDIF
1154               ENDDO
1155
1156               ! clean
1157               CALL var_clean(tl_var1)
1158
1159            ENDDO
1160
1161         !- end of use input matrix to fill variable
1162         ELSE
1163         !- use mpp file to fill variable
1164
1165            WRITE(*,'(a)') "work on file "//TRIM(tl_multi%t_mpp(ji)%c_name)
1166            !
1167            tl_file=file_init(TRIM(tl_multi%t_mpp(ji)%t_proc(1)%c_name), &
1168               &              id_perio=tl_multi%t_mpp(ji)%i_perio)
1169            tl_mpp=mpp_init( tl_file )
1170            !tl_mpp=mpp_init( tl_file, id_perio=tl_multi%t_mpp(ji)%t_proc(1)%i_perio)
1171            ! clean
1172            CALL file_clean(tl_file)
1173            CALL grid_get_info(tl_mpp)
1174
1175            DO jl=1,ip_ncard
1176               IF( tl_bdy(jl)%l_use )THEN
1177
1178                  WRITE(*,'(2x,a,a)') 'work on '//TRIM(tl_bdy(jl)%c_card)//&
1179                     &  ' boundary'
1180                  DO jk=1,tl_bdy(jl)%i_nseg
1181
1182                     ! for each variable of this file
1183                     DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
1184
1185                        WRITE(*,'(4x,a,a)') "work on variable "//&
1186                        &  TRIM(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name)
1187
1188                        tl_var0=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
1189
1190                        ! open mpp file
1191                        CALL iom_mpp_open(tl_mpp)
1192
1193                        ! get or check depth value
1194                        CALL create_boundary_check_depth( tl_var0, tl_mpp, &
1195                        &                                 in_nlevel, tl_depth )
1196
1197                        ! get or check time value
1198                        CALL create_boundary_check_time( tl_var0, tl_mpp, &
1199                        &                                tl_time )
1200
1201                        ! close mpp file
1202                        CALL iom_mpp_close(tl_mpp)
1203
1204                        ! open mpp file on domain
1205                        SELECT CASE(TRIM(tl_var0%c_point))
1206                           CASE DEFAULT !'T'
1207                              jpoint=jp_T
1208                           CASE('U')
1209                              jpoint=jp_U
1210                           CASE('V')
1211                              jpoint=jp_V
1212                           CASE('F')
1213                              jpoint=jp_F
1214                        END SELECT
1215
1216                        tl_dom1=dom_copy(tl_segdom1(jpoint,jk,jl))
1217
1218                        CALL create_boundary_get_coord( tl_coord1, tl_dom1, &
1219                        &                               tl_var0%c_point,    &
1220                        &                               tl_lon1, tl_lat1 )
1221
1222                        ! get source/coarse grid indices of this segment
1223                        il_ind(:,:)=grid_get_coarse_index(tl_coord0, &
1224                        &                                 tl_lon1, tl_lat1, &
1225                        &                                 id_rho=il_rho(:) )
1226
1227                        IF( ANY(il_ind(:,:)==0) )THEN
1228                           CALL logger_error("CREATE BOUNDARY: error "//&
1229                           &  "computing source/coarse grid indices")
1230                        ELSE
1231                           il_imin0=il_ind(1,1)
1232                           il_imax0=il_ind(1,2)
1233
1234                           il_jmin0=il_ind(2,1)
1235                           il_jmax0=il_ind(2,2)
1236                        ENDIF
1237
1238                        il_offset(:,:)= grid_get_fine_offset( &
1239                        &                    tl_coord0, &
1240                        &                    il_imin0, il_jmin0,&
1241                        &                    il_imax0, il_jmax0,&
1242                        &                    tl_lon1%d_value(:,:,1,1),&
1243                        &                    tl_lat1%d_value(:,:,1,1),&
1244                        &                    il_rho(:),&
1245                        &                    TRIM(tl_var0%c_point) )
1246
1247                        ! compute source/coarse grid segment domain
1248                        tl_dom0=dom_init( tl_coord0,         &
1249                        &                 il_imin0, il_imax0,&
1250                        &                 il_jmin0, il_jmax0 )
1251
1252                        ! add extra band (if possible) to compute interpolation
1253                        CALL dom_add_extra(tl_dom0)
1254
1255                        ! open mpp files
1256                        CALL iom_dom_open(tl_mpp, tl_dom0)
1257
1258                        cl_name=tl_var0%c_name
1259                        ! read variable value on domain
1260                        tl_segvar1(jvar+jj,jk,jl)= &
1261                        &    iom_dom_read_var(tl_mpp, TRIM(cl_name), tl_dom0)
1262
1263                        IF( ANY(il_rho(:)/=1) )THEN
1264                           WRITE(*,'(4x,a,a)') "interp variable "//TRIM(cl_name)
1265                           ! work on variable
1266                           CALL create_boundary_interp( &
1267                           &                 tl_segvar1(jvar+jj,jk,jl),&
1268                           &                 il_rho(:), il_offset(:,:) )
1269                        ENDIF
1270
1271                        ! remove extraband added to domain
1272                        CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), &
1273                        &                   tl_dom0, il_rho(:) )
1274
1275                        ! del extra point on fine grid
1276                        CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), &
1277                        &                   tl_dom1 )
1278
1279                        ! clean extra point information on source/coarse grid domain
1280                        CALL dom_clean_extra( tl_dom0 )
1281
1282                        ! add attribute to variable
1283                        tl_att=att_init('src_file',&
1284                        &  TRIM(fct_basename(tl_mpp%c_name)))
1285                        CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
1286                        &                 tl_att)
1287
1288                        !
1289                        tl_att=att_init('src_i_indices',&
1290                        &  (/tl_dom0%i_imin, tl_dom0%i_imax/))
1291                        CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
1292                        &                 tl_att)
1293
1294                        tl_att=att_init('src_j_indices', &
1295                        &  (/tl_dom0%i_jmin, tl_dom0%i_jmax/))
1296                        CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
1297                        &                 tl_att)
1298
1299                        IF( ANY(il_rho(:)/=1) )THEN
1300                           tl_att=att_init("refinment_factor", &
1301                           &               (/il_rho(jp_I),il_rho(jp_J)/))
1302                           CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), &
1303                           &                 tl_att)
1304                        ENDIF
1305
1306                        ! clean structure
1307                        CALL att_clean(tl_att)
1308
1309                        ! clean
1310                        CALL dom_clean(tl_dom0)
1311                        CALL dom_clean(tl_dom1)
1312
1313                        ! close mpp files
1314                        CALL iom_dom_close(tl_mpp)
1315
1316                        ! clean structure
1317                        CALL var_clean(tl_lon1)
1318                        CALL var_clean(tl_lat1)
1319                        CALL var_clean(tl_lvl1)
1320
1321                     ENDDO ! jj
1322
1323                     ! clean
1324                     CALL var_clean(tl_var0)
1325
1326                  ENDDO ! jk
1327
1328               ENDIF
1329            ENDDO ! jl
1330
1331            jvar=jvar+tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
1332
1333            ! clean
1334            CALL mpp_clean(tl_mpp)
1335
1336         !- end of use file to fill variable
1337         ENDIF
1338      ENDDO ! ji
1339   ENDIF
1340
1341   IF( jvar /= tl_multi%i_nvar )THEN
1342      CALL logger_error("CREATE BOUNDARY: it seems some variable "//&
1343         &  "can not be read")
1344   ENDIF
1345
1346   ! write file for each segment of each boundary
1347   DO jl=1,ip_ncard
1348      IF( tl_bdy(jl)%l_use )THEN
1349
1350         DO jk=1,tl_bdy(jl)%i_nseg
1351            !-
1352            CALL create_boundary_get_coord( tl_coord1, tl_segdom1(jp_T,jk,jl),&
1353               &                           'T', tl_lon1, tl_lat1 )
1354
1355            ! force to use nav_lon, nav_lat as variable name
1356            tl_lon1%c_name='nav_lon'
1357            tl_lat1%c_name='nav_lat'
1358
1359            ! del extra point on fine grid
1360            CALL dom_del_extra( tl_lon1, tl_segdom1(jp_T,jk,jl) )
1361            CALL dom_del_extra( tl_lat1, tl_segdom1(jp_T,jk,jl) )
1362
1363            ! clean
1364            DO jpoint=1,ip_npoint
1365               CALL dom_clean(tl_segdom1(jpoint,jk,jl))
1366            ENDDO
1367
1368            ! swap array
1369            CALL boundary_swap(tl_lon1, tl_bdy(jl))
1370            CALL boundary_swap(tl_lat1, tl_bdy(jl))
1371            DO jvar=1,tl_multi%i_nvar
1372
1373               ! use additional request
1374               ! change unit and apply factor
1375               CALL var_chg_unit(tl_segvar1(jvar,jk,jl))
1376
1377               ! forced min and max value
1378               CALL var_limit_value(tl_segvar1(jvar,jk,jl))
1379
1380               ! filter
1381               CALL filter_fill_value(tl_segvar1(jvar,jk,jl))
1382
1383               IF( .NOT. ln_extrap )THEN
1384                  ! use mask
1385                  SELECT CASE(TRIM(tl_segvar1(jvar,jk,jl)%c_point))
1386                  CASE DEFAULT !'T'
1387                     jpoint=jp_T
1388                  CASE('U')
1389                     jpoint=jp_U
1390                  CASE('V')
1391                     jpoint=jp_V
1392                  CASE('F')
1393                     jpoint=jp_F
1394                  END SELECT
1395
1396                  CALL create_boundary_use_mask(tl_segvar1(jvar,jk,jl), &
1397                  &                             tl_seglvl1(jpoint,jk,jl))
1398               ENDIF
1399
1400               ! swap dimension order
1401               CALL boundary_swap(tl_segvar1(jvar,jk,jl), tl_bdy(jl))
1402
1403            ENDDO
1404
1405            ! create file
1406            ! create file structure
1407            ! set file namearray of level variable structure
1408            IF( tl_bdy(jl)%i_nseg > 1 )THEN
1409               IF( ASSOCIATED(tl_time%d_value) )THEN
1410                  cl_fmt="('y',i0.4,'m',i0.2,'d',i0.2)"
1411                  tl_date=var_to_date(tl_time)
1412                  tl_date=tl_date+dn_dayofs
1413                  cl_date=date_print( tl_date, cl_fmt )
1414
1415                  cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1416                  &                                TRIM(tl_bdy(jl)%c_card), jk,&
1417                  &                                cd_date=TRIM(cl_date) )
1418               ELSE
1419                  cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1420                  &                                TRIM(tl_bdy(jl)%c_card), jk )
1421               ENDIF
1422            ELSE
1423               IF( ASSOCIATED(tl_time%d_value) )THEN
1424                  cl_fmt="('y',i0.4,'m',i0.2,'d',i0.2)"
1425                  tl_date=var_to_date(tl_time)
1426                  tl_date=tl_date+dn_dayofs
1427                  cl_date=date_print( tl_date, cl_fmt )
1428
1429                  cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1430                  &                                TRIM(tl_bdy(jl)%c_card), &
1431                  &                                cd_date=TRIM(cl_date) )
1432               ELSE
1433                  cl_bdyout=boundary_set_filename( TRIM(cn_fileout), &
1434                  &                                TRIM(tl_bdy(jl)%c_card) )
1435               ENDIF
1436            ENDIF
1437            !
1438            tl_fileout=file_init(TRIM(cl_bdyout),id_perio=in_perio1)
1439
1440            ! add dimension
1441            tl_dim(:)=var_max_dim(tl_segvar1(:,jk,jl))
1442
1443            SELECT CASE(TRIM(tl_bdy(jl)%c_card))
1444               CASE DEFAULT ! 'north','south'
1445                  cl_dimorder='xyzt'
1446               CASE('east','west')
1447                  cl_dimorder='yxzt'
1448            END SELECT
1449
1450            DO ji=1,ip_maxdim
1451               IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
1452            ENDDO
1453
1454            ! add variables
1455            IF( ALL( tl_dim(1:2)%l_use ) )THEN
1456               ! add longitude
1457               CALL file_add_var(tl_fileout, tl_lon1)
1458               CALL var_clean(tl_lon1)
1459
1460               ! add latitude
1461               CALL file_add_var(tl_fileout, tl_lat1)
1462               CALL var_clean(tl_lat1)
1463            ENDIF
1464
1465
1466
1467            IF( tl_dim(3)%l_use )THEN
1468               IF( ASSOCIATED(tl_depth%d_value) )THEN
1469                  ! add depth
1470                  CALL file_add_var(tl_fileout, tl_depth)
1471               ENDIF
1472            ENDIF
1473
1474            IF( tl_dim(4)%l_use )THEN
1475               IF( ASSOCIATED(tl_time%d_value) )THEN
1476                  ! add time
1477                  CALL file_add_var(tl_fileout, tl_time)
1478               ENDIF
1479            ENDIF
1480
1481            ! add other variable
1482            DO jvar=tl_multi%i_nvar,1,-1
1483               CALL file_add_var(tl_fileout, tl_segvar1(jvar,jk,jl))
1484               CALL var_clean(tl_segvar1(jvar,jk,jl))
1485            ENDDO
1486
1487            ! add some attribute
1488            tl_att=att_init("Created_by","SIREN create_boundary")
1489            CALL file_add_att(tl_fileout, tl_att)
1490
1491            !add source url
1492            cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
1493            tl_att=att_init("SIREN_url",cl_url)
1494            CALL file_add_att(tl_fileout, tl_att)
1495
1496            ! add date of creation
1497            cl_date=date_print(date_now())
1498            tl_att=att_init("Creation_date",cl_date)
1499            CALL file_add_att(tl_fileout, tl_att)
1500
1501            ! add shift on north and east boundary
1502            ! boundary compute on T point but express on U or V point
1503            SELECT CASE(TRIM(tl_bdy(jl)%c_card))
1504            CASE DEFAULT ! 'south','west'
1505               il_shift=0
1506            CASE('north','east')
1507               il_shift=1
1508            END SELECT
1509
1510            ! add indice of velocity row or column
1511            tl_att=att_init('bdy_ind',tl_bdy(jl)%t_seg(jk)%i_index-il_shift)
1512            CALL file_move_att(tl_fileout, tl_att)
1513
1514            ! add width of the relaxation zone
1515            tl_att=att_init('bdy_width',tl_bdy(jl)%t_seg(jk)%i_width)
1516            CALL file_move_att(tl_fileout, tl_att)
1517
1518            ! add indice of segment start
1519            tl_att=att_init('bdy_deb',tl_bdy(jl)%t_seg(jk)%i_first)
1520            CALL file_move_att(tl_fileout, tl_att)
1521
1522            ! add indice of segment end
1523            tl_att=att_init('bdy_end',tl_bdy(jl)%t_seg(jk)%i_last)
1524            CALL file_move_att(tl_fileout, tl_att)
1525
1526            ! clean
1527            CALL att_clean(tl_att)
1528
1529            ! create file
1530            CALL iom_create(tl_fileout)
1531
1532            ! write file
1533            CALL iom_write_file(tl_fileout, cl_dimorder)
1534
1535            ! close file
1536            CALL iom_close(tl_fileout)
1537            CALL file_clean(tl_fileout)
1538
1539         ENDDO ! jk
1540
1541      ENDIF
1542      ! clean
1543      CALL boundary_clean(tl_bdy(jl))
1544   ENDDO !jl
1545
1546   ! clean
1547   IF( ASSOCIATED(tl_depth%d_value) ) CALL var_clean(tl_depth)
1548   IF( ASSOCIATED(tl_time%d_value) ) CALL var_clean(tl_time)
1549   DEALLOCATE( tl_segdom1 )
1550   DEALLOCATE( tl_segvar1 )
1551   CALL var_clean(tl_seglvl1(:,:,:))
1552   DEALLOCATE( tl_seglvl1 )
1553
1554
1555   CALL mpp_clean(tl_coord1)
1556   CALL mpp_clean(tl_coord0)
1557   CALL var_clean_extra()
1558
1559   CALL multi_clean(tl_multi)
1560
1561   ! close log file
1562   CALL logger_footer()
1563   CALL logger_close()
1564   CALL logger_clean()
1565
1566   END SUBROUTINE create_boundary__mono
1567   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1568   FUNCTION create_boundary_get_dom(td_bathy1, td_bdy, id_seg) &
1569         & RESULT (tf_dom)
1570   !-------------------------------------------------------------------
1571   !> @brief
1572   !> This subroutine compute boundary domain for each grid point (T,U,V,F)
1573   !>
1574   !> @author J.Paul
1575   !> @date November, 2013 - Initial Version
1576   !> @date September, 2014
1577   !> - take into account grid point to compute boundary indices
1578   !>
1579   !> @param[in] td_bathy1 file structure
1580   !> @param[in] td_bdy    boundary structure
1581   !> @param[in] id_seg    segment indice
1582   !> @return array of domain structure
1583   !-------------------------------------------------------------------
1584
1585      IMPLICIT NONE
1586
1587      ! Argument
1588      TYPE(TMPP) , INTENT(IN   ) :: td_bathy1
1589      TYPE(TBDY) , INTENT(IN   ) :: td_bdy
1590      INTEGER(i4), INTENT(IN   ) :: id_seg
1591
1592      ! function
1593      TYPE(TDOM), DIMENSION(ip_npoint) :: tf_dom
1594
1595      ! local variable
1596      INTEGER(i4) :: il_imin1
1597      INTEGER(i4) :: il_imax1
1598      INTEGER(i4) :: il_jmin1
1599      INTEGER(i4) :: il_jmax1
1600
1601      INTEGER(i4) :: il_imin
1602      INTEGER(i4) :: il_imax
1603      INTEGER(i4) :: il_jmin
1604      INTEGER(i4) :: il_jmax
1605
1606      INTEGER(i4), DIMENSION(ip_npoint) :: il_ishift
1607      INTEGER(i4), DIMENSION(ip_npoint) :: il_jshift
1608
1609      ! loop indices
1610      INTEGER(i4) :: ji
1611      INTEGER(i4) :: jk
1612      !----------------------------------------------------------------
1613      ! init
1614      jk=id_seg
1615
1616      il_ishift(:)=0
1617      il_jshift(:)=0
1618
1619      ! get boundary definition
1620      SELECT CASE(TRIM(td_bdy%c_card))
1621         CASE('north')
1622
1623            il_imin1=td_bdy%t_seg(jk)%i_first
1624            il_imax1=td_bdy%t_seg(jk)%i_last
1625            il_jmin1=td_bdy%t_seg(jk)%i_index-(td_bdy%t_seg(jk)%i_width-1)
1626            il_jmax1=td_bdy%t_seg(jk)%i_index
1627
1628            il_jshift(jp_V)=-1
1629            il_jshift(jp_F)=-1
1630
1631         CASE('south')
1632
1633            il_imin1=td_bdy%t_seg(jk)%i_first
1634            il_imax1=td_bdy%t_seg(jk)%i_last
1635            il_jmin1=td_bdy%t_seg(jk)%i_index
1636            il_jmax1=td_bdy%t_seg(jk)%i_index+(td_bdy%t_seg(jk)%i_width-1)
1637
1638         CASE('east')
1639
1640            il_imin1=td_bdy%t_seg(jk)%i_index-(td_bdy%t_seg(jk)%i_width-1)
1641            il_imax1=td_bdy%t_seg(jk)%i_index
1642            il_jmin1=td_bdy%t_seg(jk)%i_first
1643            il_jmax1=td_bdy%t_seg(jk)%i_last
1644
1645            il_ishift(jp_U)=-1
1646            il_ishift(jp_F)=-1
1647
1648         CASE('west')
1649
1650            il_imin1=td_bdy%t_seg(jk)%i_index
1651            il_imax1=td_bdy%t_seg(jk)%i_index+(td_bdy%t_seg(jk)%i_width-1)
1652            il_jmin1=td_bdy%t_seg(jk)%i_first
1653            il_jmax1=td_bdy%t_seg(jk)%i_last
1654
1655      END SELECT
1656
1657      !-read fine grid domain
1658      DO ji=1,ip_npoint
1659
1660         ! shift domain
1661         il_imin=il_imin1+il_ishift(ji)
1662         il_imax=il_imax1+il_ishift(ji)
1663
1664         il_jmin=il_jmin1+il_jshift(ji)
1665         il_jmax=il_jmax1+il_jshift(ji)
1666
1667         ! compute domain
1668         tf_dom(ji)=dom_init(td_bathy1,         &
1669            &                il_imin, il_imax,  &
1670            &                il_jmin, il_jmax,  &
1671            &                TRIM(td_bdy%c_card))
1672
1673      ENDDO
1674
1675   END FUNCTION create_boundary_get_dom
1676   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1677   SUBROUTINE create_boundary_get_coord(td_coord1, td_dom1, cd_point, &
1678      &                                 td_lon1, td_lat1)
1679   !-------------------------------------------------------------------
1680   !> @brief
1681   !> This subroutine get coordinates over boundary domain
1682   !>
1683   !> @author J.Paul
1684   !> @date November, 2013 - Initial Version
1685   !> @date September, 2014
1686   !> - take into account grid point
1687   !>
1688   !> @param[in] td_coord1 coordinates file structure
1689   !> @param[in] td_dom1   boundary domain structure
1690   !> @param[in] cd_point  grid point
1691   !> @param[out] td_lon1  longitude variable structure
1692   !> @param[out] td_lat1  latitude variable structure
1693   !-------------------------------------------------------------------
1694
1695      IMPLICIT NONE
1696
1697      ! Argument
1698      TYPE(TMPP)      , INTENT(IN   ) :: td_coord1
1699      TYPE(TDOM)      , INTENT(IN   ) :: td_dom1
1700      CHARACTER(LEN=*), INTENT(IN   ) :: cd_point
1701      TYPE(TVAR)      , INTENT(  OUT) :: td_lon1
1702      TYPE(TVAR)      , INTENT(  OUT) :: td_lat1
1703
1704      ! local variable
1705      TYPE(TMPP)        :: tl_coord1
1706
1707      CHARACTER(LEN=lc) :: cl_name
1708      ! loop indices
1709      !----------------------------------------------------------------
1710      !read variables on domain (ugly way to do it, have to work on it)
1711      ! init mpp structure
1712      tl_coord1=mpp_copy(td_coord1)
1713
1714      ! open mpp files
1715      CALL iom_dom_open(tl_coord1, td_dom1)
1716
1717      ! read variable value on domain
1718      WRITE(cl_name,*) 'longitude_'//TRIM(cd_point)
1719      td_lon1=iom_dom_read_var( tl_coord1, TRIM(cl_name), td_dom1)
1720      WRITE(cl_name,*) 'latitude_'//TRIM(cd_point)
1721      td_lat1=iom_dom_read_var( tl_coord1, TRIM(cl_name), td_dom1)
1722
1723      ! close mpp files
1724      CALL iom_dom_close(tl_coord1)
1725
1726      ! clean structure
1727      CALL mpp_clean(tl_coord1)
1728
1729   END SUBROUTINE create_boundary_get_coord
1730   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1731   SUBROUTINE create_boundary_interp(td_var, id_rho, id_offset, &
1732         &                           id_iext, id_jext)
1733   !-------------------------------------------------------------------
1734   !> @brief
1735   !> This subroutine interpolate variable on boundary
1736   !>
1737   !> @details
1738   !>
1739   !> @author J.Paul
1740   !> @date November, 2013 - Initial Version
1741   !>
1742   !> @param[inout] td_var variable structure
1743   !> @param[in] id_rho    array of refinment factor
1744   !> @param[in] id_offset array of offset between fine and source/coarse grid
1745   !> @param[in] id_iext   i-direction size of extra bands (default=im_minext)
1746   !> @param[in] id_jext   j-direction size of extra bands (default=im_minext)
1747   !-------------------------------------------------------------------
1748
1749      IMPLICIT NONE
1750
1751      ! Argument
1752      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
1753      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
1754      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_offset
1755
1756      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_iext
1757      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_jext
1758
1759
1760      ! local variable
1761      INTEGER(i4) :: il_iext
1762      INTEGER(i4) :: il_jext
1763      ! loop indices
1764      !----------------------------------------------------------------
1765
1766      !WARNING: at least two extrabands are required for cubic interpolation
1767      il_iext=2
1768      IF( PRESENT(id_iext) ) il_iext=id_iext
1769
1770      il_jext=2
1771      IF( PRESENT(id_jext) ) il_jext=id_jext
1772
1773      IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1774         CALL logger_warn("CREATE BOUNDARY INTERP: at least extrapolation "//&
1775         &  "on two points are required with cubic interpolation ")
1776         il_iext=2
1777      ENDIF
1778
1779      IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1780         CALL logger_warn("CREATE BOUNDARY INTERP: at least extrapolation "//&
1781         &  "on two points are required with cubic interpolation ")
1782         il_jext=2
1783      ENDIF
1784
1785      ! work on variable
1786      ! add extraband
1787      CALL extrap_add_extrabands(td_var, il_iext, il_jext)
1788
1789      ! extrapolate variable
1790      CALL extrap_fill_value( td_var )
1791
1792      ! interpolate variable
1793      CALL interp_fill_value( td_var, id_rho(:), &
1794      &                       id_offset=id_offset(:,:) )
1795
1796      ! remove extraband
1797      CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), &
1798         &                               il_jext*id_rho(jp_J))
1799
1800   END SUBROUTINE create_boundary_interp
1801   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1802   FUNCTION create_boundary_matrix(td_var, td_dom, id_nlevel) &
1803         & RESULT (tf_var)
1804   !-------------------------------------------------------------------
1805   !> @brief
1806   !> This function create variable, filled with matrix value
1807   !>
1808   !> @details
1809   !> A variable is create with the same name that the input variable,
1810   !> and with dimension of the coordinate file.
1811   !> Then the variable array of value is split into equal subdomain.
1812   !> Each subdomain is fill with the associated value of the matrix.
1813   !>
1814   !> @author J.Paul
1815   !> @date November, 2013 - Initial Version
1816   !>
1817   !> @param[in] td_var    variable structure
1818   !> @param[in] td_dom    domain structure
1819   !> @param[in] id_nlevel number of levels
1820   !> @return variable structure
1821   !-------------------------------------------------------------------
1822
1823      IMPLICIT NONE
1824
1825      ! Argument
1826      TYPE(TVAR) , INTENT(IN) :: td_var
1827      TYPE(TDOM) , INTENT(IN) :: td_dom
1828      INTEGER(i4), INTENT(IN) :: id_nlevel
1829
1830      ! function
1831      TYPE(TVAR)              :: tf_var
1832
1833      ! local variable
1834      INTEGER(i4)      , DIMENSION(3)                    :: il_dim
1835      INTEGER(i4)      , DIMENSION(3)                    :: il_size
1836      INTEGER(i4)      , DIMENSION(3)                    :: il_rest
1837
1838      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_ishape
1839      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_jshape
1840      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_kshape
1841
1842      REAL(dp)         , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
1843
1844      TYPE(TDIM)       , DIMENSION(ip_maxdim)            :: tl_dim
1845
1846      ! loop indices
1847      INTEGER(i4) :: ji
1848      INTEGER(i4) :: jj
1849      INTEGER(i4) :: jk
1850      !----------------------------------------------------------------
1851
1852      ! write value on grid
1853      ! get matrix dimension
1854      il_dim(:)=td_var%t_dim(1:3)%i_len
1855
1856      tl_dim(jp_I:jp_J)=dim_copy(td_dom%t_dim(jp_I:jp_J))
1857      tl_dim(jp_K)%i_len=id_nlevel
1858
1859      ! split output domain in N subdomain depending of matrix dimension
1860      il_size(:) = tl_dim(1:3)%i_len / il_dim(:)
1861      il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:))
1862
1863      ALLOCATE( il_ishape(il_dim(1)+1) )
1864      il_ishape(:)=0
1865      DO ji=2,il_dim(1)+1
1866         il_ishape(ji)=il_ishape(ji-1)+il_size(1)
1867      ENDDO
1868      ! add rest to last cell
1869      il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1)
1870
1871      ALLOCATE( il_jshape(il_dim(2)+1) )
1872      il_jshape(:)=0
1873      DO jj=2,il_dim(2)+1
1874         il_jshape(jj)=il_jshape(jj-1)+il_size(2)
1875      ENDDO
1876      ! add rest to last cell
1877      il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2)
1878
1879      ALLOCATE( il_kshape(il_dim(3)+1) )
1880      il_kshape(:)=0
1881      DO jk=2,il_dim(3)+1
1882         il_kshape(jk)=il_kshape(jk-1)+il_size(3)
1883      ENDDO
1884      ! add rest to last cell
1885      il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3)
1886
1887      ! write ouput array of value
1888      ALLOCATE(dl_value( tl_dim(1)%i_len, &
1889      &                  tl_dim(2)%i_len, &
1890      &                  tl_dim(3)%i_len, &
1891      &                  tl_dim(4)%i_len) )
1892
1893      dl_value(:,:,:,:)=0
1894
1895      DO jk=2,il_dim(3)+1
1896         DO jj=2,il_dim(2)+1
1897            DO ji=2,il_dim(1)+1
1898
1899               dl_value( 1+il_ishape(ji-1):il_ishape(ji), &
1900               &         1+il_jshape(jj-1):il_jshape(jj), &
1901               &         1+il_kshape(jk-1):il_kshape(jk), &
1902               &         1 ) = td_var%d_value(ji-1,jj-1,jk-1,1)
1903
1904            ENDDO
1905         ENDDO
1906      ENDDO
1907
1908      ! initialise variable with value
1909      tf_var=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:))
1910
1911      DEALLOCATE(dl_value)
1912
1913   END FUNCTION create_boundary_matrix
1914   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1915   SUBROUTINE create_boundary_use_mask(td_var, td_mask)
1916   !-------------------------------------------------------------------
1917   !> @brief
1918   !> This subroutine use mask to filled land point with _FillValue
1919   !>
1920   !> @details
1921   !>
1922   !> @author J.Paul
1923   !> @date November, 2013 - Initial Version
1924   !>
1925   !> @param[inout] td_var variable structure
1926   !> @param[in] td_mask   mask variable structure
1927   !-------------------------------------------------------------------
1928
1929      IMPLICIT NONE
1930
1931      ! Argument
1932      TYPE(TVAR), INTENT(INOUT) :: td_var
1933      TYPE(TVAR), INTENT(IN   ) :: td_mask
1934
1935      ! local variable
1936      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask
1937
1938      ! loop indices
1939      INTEGER(i4) :: jk
1940      INTEGER(i4) :: jl
1941      !----------------------------------------------------------------
1942
1943      IF( ANY(td_var%t_dim(1:2)%i_len /= &
1944      &       td_mask%t_dim(1:2)%i_len) )THEN
1945         CALL logger_debug("     mask dimension ( "//&
1946         &              TRIM(fct_str(td_mask%t_dim(1)%i_len))//","//&
1947         &              TRIM(fct_str(td_mask%t_dim(2)%i_len))//")" )
1948         CALL logger_debug(" variable dimension ( "//&
1949         &              TRIM(fct_str(td_var%t_dim(1)%i_len))//","//&
1950         &              TRIM(fct_str(td_var%t_dim(2)%i_len))//")" )
1951         CALL logger_fatal("CREATE BOUNDARY USE MASK: mask and "//&
1952         &                 "variable dimension differ."   )
1953      ENDIF
1954
1955      ALLOCATE( il_mask(td_var%t_dim(1)%i_len, &
1956      &                 td_var%t_dim(2)%i_len) )
1957
1958      il_mask(:,:)=INT(td_mask%d_value(:,:,1,1))
1959
1960      DO jl=1,td_var%t_dim(4)%i_len
1961         DO jk=1,td_var%t_dim(3)%i_len
1962            WHERE( il_mask(:,:) < jk ) td_var%d_value(:,:,jk,jl)=td_var%d_fill
1963         ENDDO
1964      ENDDO
1965
1966      DEALLOCATE( il_mask )
1967
1968   END SUBROUTINE create_boundary_use_mask
1969   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1970   FUNCTION create_boundary_get_level(td_level, td_dom) &
1971         & RESULT (tf_var)
1972   !-------------------------------------------------------------------
1973   !> @brief
1974   !> This function extract level over domain on each grid point, and return
1975   !> array of variable structure
1976   !>
1977   !> @author J.Paul
1978   !> @date November, 2013 - Initial Version
1979   !>
1980   !> @param[in] td_level  array of level variable structure
1981   !> @param[in] td_dom    array of domain structure
1982   !> @return array of variable structure
1983   !-------------------------------------------------------------------
1984
1985      IMPLICIT NONE
1986
1987      ! Argument
1988      TYPE(TVAR), DIMENSION(:), INTENT(IN) :: td_level
1989      TYPE(TDOM), DIMENSION(:), INTENT(IN) :: td_dom
1990
1991      ! function
1992      TYPE(TVAR), DIMENSION(ip_npoint)     :: tf_var
1993
1994      ! local variable
1995      ! loop indices
1996      INTEGER(i4) :: ji
1997      !----------------------------------------------------------------
1998
1999      IF( SIZE(td_level(:)) /= ip_npoint .OR. &
2000      &   SIZE(td_dom(:)) /= ip_npoint )THEN
2001         CALL logger_error("CREATE BDY GET LEVEL: invalid dimension. "//&
2002         &  "check input array of level and domain.")
2003      ELSE
2004
2005         DO ji=1,ip_npoint
2006
2007            tf_var(ji)=var_copy(td_level(ji))
2008
2009            IF( ASSOCIATED(tf_var(ji)%d_value) ) DEALLOCATE(tf_var(ji)%d_value)
2010
2011            tf_var(ji)%t_dim(1)%i_len=td_dom(ji)%t_dim(1)%i_len
2012            tf_var(ji)%t_dim(2)%i_len=td_dom(ji)%t_dim(2)%i_len
2013            ALLOCATE(tf_var(ji)%d_value(tf_var(ji)%t_dim(1)%i_len, &
2014            &                           tf_var(ji)%t_dim(2)%i_len, &
2015            &                           tf_var(ji)%t_dim(3)%i_len, &
2016            &                           tf_var(ji)%t_dim(4)%i_len) )
2017
2018            tf_var(ji)%d_value(:,:,:,:) = &
2019            &  td_level(ji)%d_value( td_dom(ji)%i_imin:td_dom(ji)%i_imax, &
2020            &                        td_dom(ji)%i_jmin:td_dom(ji)%i_jmax, :, : )
2021
2022         ENDDO
2023
2024      ENDIF
2025
2026   END FUNCTION create_boundary_get_level
2027   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2028   SUBROUTINE create_boundary_check_depth(td_var, td_mpp, id_nlevel, td_depth)
2029   !-------------------------------------------------------------------
2030   !> @brief
2031   !> This subroutine check if variable need depth dimension,
2032   !> get depth variable value in an open mpp structure
2033   !> and check if agree with already input depth variable.
2034   !>
2035   !> @details
2036   !>
2037   !> @author J.Paul
2038   !> @date November, 2014 - Initial Version
2039   !> @date January, 2016
2040   !> - check if variable need/use depth dimension
2041   !>
2042   !> @param[in] td_var       variable structure
2043   !> @param[in] td_mpp       mpp structure
2044   !> @param[in] id_nlevel    mpp structure
2045   !> @param[inout] td_depth  depth variable structure
2046   !-------------------------------------------------------------------
2047
2048      IMPLICIT NONE
2049
2050      ! Argument
2051      TYPE(TVAR) , INTENT(IN   ) :: td_var
2052      TYPE(TMPP) , INTENT(IN   ) :: td_mpp
2053      INTEGER(i4), INTENT(IN   ) :: id_nlevel
2054      TYPE(TVAR) , INTENT(INOUT) :: td_depth
2055
2056      ! local variable
2057      INTEGER(i4) :: il_varid
2058      TYPE(TVAR)  :: tl_depth
2059      ! loop indices
2060      !----------------------------------------------------------------
2061
2062      IF( td_var%t_dim(jp_K)%l_use .AND. &
2063      &   ( TRIM(td_var%c_axis) == '' .OR. &
2064      &     INDEX(TRIM(td_var%c_axis),'Z') /= 0 )&
2065      & )THEN
2066
2067         ! check vertical dimension
2068         IF( td_mpp%t_dim(jp_K)%l_use )THEN
2069            IF( td_mpp%t_dim(jp_K)%i_len /= id_nlevel .AND. &
2070              & td_mpp%t_dim(jp_K)%i_len /= 1 )THEN
2071               CALL logger_error("CREATE BOUNDARY: dimension in file "//&
2072               &  TRIM(td_mpp%c_name)//" not agree with namelist in_nlevel ")
2073            ENDIF
2074         ENDIF
2075
2076         ! get or check depth value
2077         IF( td_mpp%t_proc(1)%i_depthid /= 0 )THEN
2078
2079            il_varid=td_mpp%t_proc(1)%i_depthid
2080            IF( ASSOCIATED(td_depth%d_value) )THEN
2081
2082               tl_depth=iom_mpp_read_var(td_mpp, il_varid)
2083
2084               IF( ANY( td_depth%d_value(:,:,:,:) /= &
2085               &        tl_depth%d_value(:,:,:,:) ) )THEN
2086
2087                  CALL logger_error("CREATE BOUNDARY: depth value "//&
2088                  &  "for variable "//TRIM(td_var%c_name)//&
2089                  &  "from "//TRIM(td_mpp%c_name)//" not conform "//&
2090                  &  " to those from former file(s).")
2091
2092               ENDIF
2093               CALL var_clean(tl_depth)
2094
2095            ELSE
2096               td_depth=iom_mpp_read_var(td_mpp,il_varid)
2097            ENDIF
2098
2099         ENDIF
2100      ELSE
2101         CALL logger_debug("CREATE BOUNDARY: no depth dimension use"//&
2102         &                 " for variable "//TRIM(td_var%c_name))
2103      ENDIF
2104
2105   END SUBROUTINE create_boundary_check_depth
2106   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2107   SUBROUTINE create_boundary_check_time(td_var, td_mpp, td_time)
2108   !-------------------------------------------------------------------
2109   !> @brief
2110   !> This subroutine check if variable need time dimension,
2111   !> get date and time in an open mpp structure
2112   !> and check if agree with date and time already read.
2113   !>
2114   !> @details
2115   !>
2116   !> @author J.Paul
2117   !> @date November, 2014 - Initial Version
2118   !> @date January, 2016
2119   !> - check if variable need/use time dimension
2120   !>
2121   !> @param[in] td_var       variable structure
2122   !> @param[in] td_mpp      mpp structure
2123   !> @param[inout] td_time  time variable structure
2124   !-------------------------------------------------------------------
2125
2126      IMPLICIT NONE
2127
2128      ! Argument
2129      TYPE(TVAR), INTENT(IN   ) :: td_var
2130      TYPE(TMPP), INTENT(IN   ) :: td_mpp
2131      TYPE(TVAR), INTENT(INOUT) :: td_time
2132
2133      ! local variable
2134      INTEGER(i4) :: il_varid
2135      TYPE(TVAR)  :: tl_time
2136
2137      TYPE(TDATE) :: tl_date1
2138      TYPE(TDATE) :: tl_date2
2139      ! loop indices
2140      !----------------------------------------------------------------
2141      IF( td_var%t_dim(jp_L)%l_use .AND. &
2142      &   ( TRIM(td_var%c_axis) == '' .OR. &
2143      &     INDEX(TRIM(td_var%c_axis),'T') /= 0 )&
2144      & )THEN
2145
2146         ! get or check time value
2147         IF( td_mpp%t_proc(1)%i_timeid /= 0 )THEN
2148
2149            il_varid=td_mpp%t_proc(1)%i_timeid
2150            IF( ASSOCIATED(td_time%d_value) )THEN
2151
2152               tl_time=iom_mpp_read_var(td_mpp, il_varid)
2153
2154               tl_date1=var_to_date(td_time)
2155               tl_date2=var_to_date(tl_time)
2156               IF( tl_date1 - tl_date2 /= 0 )THEN
2157
2158                  CALL logger_warn("CREATE BOUNDARY: date from "//&
2159                  &  TRIM(td_mpp%c_name)//" not conform "//&
2160                  &  " to those from former file(s).")
2161
2162               ENDIF
2163               CALL var_clean(tl_time)
2164
2165            ELSE
2166               td_time=iom_mpp_read_var(td_mpp,il_varid)
2167            ENDIF
2168
2169         ENDIF
2170
2171      ELSE
2172         CALL logger_debug("CREATE BOUNDARY: no time dimension use"//&
2173         &                 " for variable "//TRIM(td_var%c_name))
2174      ENDIF
2175
2176   END SUBROUTINE create_boundary_check_time
2177   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2178END PROGRAM create_boundary
Note: See TracBrowser for help on using the repository browser.