source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/TOOLS/SIREN/src/create_boundary.f90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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