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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
create_bathy.f90 in utils/tools/SIREN/src – NEMO

source: utils/tools/SIREN/src/create_bathy.f90 @ 12080

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

update nemo trunk

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