source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/TOOLS/SIREN/src/create_bathy.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: 37.2 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! PROGRAM: create_bathy
6!
7! DESCRIPTION:
8!> @file
9!> @brief
10!> This program create fine grid bathymetry file.
11!>
12!> @details
13!> @section sec1 method
14!> Bathymetry could be extracted from fine grid Bathymetry file, interpolated
15!> from coarse grid Bathymetry file, or manually written.
16!>
17!> @section sec2 how to
18!>    to create fine grid bathymetry file:<br/>
19!> @code{.sh}
20!>    ./SIREN/bin/create_bathy create_bathy.nam
21!> @endcode
22!>   
23!>    create_bathy.nam comprise 7 namelists:<br/>
24!>       - logger namelist (namlog)
25!>       - config namelist (namcfg)
26!>       - coarse grid namelist (namcrs)
27!>       - fine grid namelist (namfin)
28!>       - variable namelist (namvar)
29!>       - nesting namelist (namnst)
30!>       - output namelist (namout)
31!>   
32!>    @note
33!>       All namelists have to be in file create_bathy.nam, however variables of
34!>       those namelists are all optional.
35!>
36!>    * _logger namelist (namlog)_:<br/>
37!>       - cn_logfile   : log filename
38!>       - cn_verbosity : verbosity ('trace','debug','info',
39!> 'warning','error','fatal')
40!>       - in_maxerror  : maximum number of error allowed
41!>
42!>    * _config namelist (namcfg)_:<br/>
43!>       - cn_varcfg : variable configuration file
44!> (see ./SIREN/cfg/variable.cfg)
45!>
46!>    * _coarse grid namelist (namcrs)_:<br/>
47!>       - cn_coord0 : coordinate file
48!>       - in_perio0 : NEMO periodicity index (see Model Boundary Condition in
49!> [NEMO documentation](http://www.nemo-ocean.eu/About-NEMO/Reference-manuals))
50!>
51!>    * _fine grid namelist (namfin)_:<br/>
52!>       - cn_coord1 : coordinate file
53!>       - in_perio1 : periodicity index
54!>       - ln_fillclosed : fill closed sea or not
55!>
56!>    * _variable namelist (namvar)_:<br/>
57!>       - cn_varinfo : list of variable and extra information about request(s)
58!>       to be used.<br/>
59!>          each elements of *cn_varinfo* is a string character.<br/>
60!>          it is composed of the variable name follow by ':',
61!>          then request(s) to be used on this variable.<br/>
62!>          request could be:
63!>             - interpolation method
64!>             - extrapolation method
65!>             - filter method
66!>             - > minimum value
67!>             - < maximum value
68!>
69!>                requests must be separated by ';'.<br/>
70!>                order of requests does not matter.<br/>
71!>
72!>          informations about available method could be find in @ref interp,
73!>          @ref extrap and @ref filter modules.<br/>
74!>          Example: 'Bathymetry: 2*hamming(2,3); > 0'
75!>          @note
76!>             If you do not specify a method which is required,
77!>             default one is apply.
78!>          @warning
79!>             variable name must be __Bathymetry__ here.
80!>       - cn_varfile : list of variable, and corresponding file.<br/>
81!>          *cn_varfile* is the path and filename of the file where find
82!>          variable.
83!>          @note
84!>             *cn_varfile* could be a matrix of value, if you want to filled
85!>             manually variable value.<br/>
86!>             the variable array of value is split into equal subdomain.<br/>
87!>             Each subdomain is filled with the corresponding value
88!>             of the matrix.<br/>         
89!>             separators used to defined matrix are:
90!>                - ',' for line
91!>                - '/' for row
92!>                - '\' for level<br/>
93!>                Example:<br/>
94!>                   3,2,3/1,4,5  =>  @f$ \left( \begin{array}{ccc}
95!>                                         3 & 2 & 3 \\
96!>                                         1 & 4 & 5 \end{array} \right) @f$
97!>
98!>          Examples:
99!>             - 'Bathymetry:gridT.nc'
100!>             - 'Bathymetry:5000,5000,5000/5000,3000,5000/5000,5000,5000'
101!>
102!>          \image html  bathy_40.png
103!>          \image latex bathy_30.png
104!>
105!>    * _nesting namelist (namnst)_:<br/>
106!>       - in_rhoi  : refinement factor in i-direction
107!>       - in_rhoj  : refinement factor in j-direction
108!>       @note
109!>          coarse grid indices will be deduced from fine grid
110!>          coordinate file.
111!>
112!>    * _output namelist (namout)_:<br/>
113!>       - cn_fileout : output bathymetry file
114!>
115!> @author J.Paul
116! REVISION HISTORY:
117!> @date November, 2013 - Initial Version
118!> @date Sepember, 2014
119!> - add header for user
120!> - Bug fix, compute offset depending of grid point
121!
122!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
123!----------------------------------------------------------------------
124PROGRAM create_bathy
125
126   USE global                          ! global variable
127   USE kind                            ! F90 kind parameter
128   USE logger                          ! log file manager
129   USE fct                             ! basic useful function
130   USE date                            ! date manager
131   USE att                             ! attribute manager
132   USE dim                             ! dimension manager
133   USE var                             ! variable manager
134   USE file                            ! file manager
135   USE multi                           ! multi file manager
136   USE iom                             ! I/O manager
137   USE grid                            ! grid manager
138   USE extrap                          ! extrapolation manager
139   USE interp                          ! interpolation manager
140   USE filter                          ! filter manager
141   USE mpp                             ! MPP manager
142   USE dom                             ! domain manager
143   USE iom_mpp                         ! MPP I/O manager
144   USE iom_dom                         ! DOM I/O manager
145
146   IMPLICIT NONE
147
148   ! local variable
149   CHARACTER(LEN=lc)                                  :: cl_namelist
150   CHARACTER(LEN=lc)                                  :: cl_date
151   CHARACTER(LEN=lc)                                  :: cl_data
152
153   INTEGER(i4)                                        :: il_narg
154   INTEGER(i4)                                        :: il_status
155   INTEGER(i4)                                        :: il_fileid
156   INTEGER(i4)                                        :: il_varid
157   INTEGER(i4)                                        :: il_attid
158   INTEGER(i4)                                        :: il_imin0
159   INTEGER(i4)                                        :: il_imax0
160   INTEGER(i4)                                        :: il_jmin0
161   INTEGER(i4)                                        :: il_jmax0
162   INTEGER(i4)      , DIMENSION(ip_maxdim)            :: il_rho
163   INTEGER(i4)      , DIMENSION(2,2)                  :: il_offset
164   INTEGER(i4)      , DIMENSION(2,2)                  :: il_ind
165   INTEGER(i4)      , DIMENSION(:,:)    , ALLOCATABLE :: il_mask
166
167   LOGICAL                                            :: ll_exist
168
169   TYPE(TMPP)                                         :: tl_coord0
170   TYPE(TMPP)                                         :: tl_coord1
171   TYPE(TMPP)                                         :: tl_mpp
172   TYPE(TFILE)                                        :: tl_fileout
173
174   TYPE(TATT)                                         :: tl_att
175   
176   TYPE(TVAR)                                         :: tl_lon
177   TYPE(TVAR)                                         :: tl_lat
178   TYPE(TVAR)                                         :: tl_depth
179   TYPE(TVAR)                                         :: tl_time
180
181   TYPE(TVAR)                                         :: tl_tmp
182   TYPE(TVAR)       , DIMENSION(:), ALLOCATABLE       :: tl_var
183   
184   TYPE(TDIM)       , DIMENSION(ip_maxdim)            :: tl_dim
185
186   TYPE(TMULTI)                                       :: tl_multi
187
188   REAL(dp)                                           :: dl_minbat
189
190   ! loop indices
191   INTEGER(i4) :: ji
192   INTEGER(i4) :: jj
193   INTEGER(i4) :: jk
194
195   ! namelist variable
196   ! namlog
197   CHARACTER(LEN=lc) :: cn_logfile = 'create_bathy.log' 
198   CHARACTER(LEN=lc) :: cn_verbosity = 'warning' 
199   INTEGER(i4)       :: in_maxerror = 5
200
201   ! namcfg
202   CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg' 
203
204   ! namcrs
205   CHARACTER(LEN=lc) :: cn_coord0 = '' 
206   INTEGER(i4)       :: in_perio0 = -1
207
208   ! namfin
209   CHARACTER(LEN=lc) :: cn_coord1 = ''
210   INTEGER(i4)       :: in_perio1 = -1
211   LOGICAL           :: ln_fillclosed = .TRUE.
212
213   ! namvar
214   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
215   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = ''
216
217   ! namnst
218   INTEGER(i4)       :: in_rhoi  = 1
219   INTEGER(i4)       :: in_rhoj  = 1
220
221   ! namout
222   CHARACTER(LEN=lc) :: cn_fileout = 'bathy_fine.nc' 
223   !-------------------------------------------------------------------
224
225   NAMELIST /namlog/ &   !< logger namelist
226   &  cn_logfile,    &   !< log file
227   &  cn_verbosity,  &   !< log verbosity
228   &  in_maxerror        !< logger maximum error
229
230   NAMELIST /namcfg/ &   !< configuration namelist
231   &  cn_varcfg          !< variable configuration file
232
233   NAMELIST /namcrs/ &   !< coarse grid namelist
234   &  cn_coord0,  &      !< coordinate file
235   &  in_perio0          !< periodicity index
236   
237   NAMELIST /namfin/ &   !< fine grid namelist
238   &  cn_coord1,   &     !< coordinate file
239   &  in_perio1,   &     !< periodicity index
240   &  ln_fillclosed      !< fill closed sea
241 
242   NAMELIST /namvar/ &   !< variable namelist
243   &  cn_varinfo, &      !< list of variable and interpolation method to be used. (ex: 'votemper:linear','vosaline:cubic' )
244   &  cn_varfile         !< list of variable file
245   
246   NAMELIST /namnst/ &   !< nesting namelist
247   &  in_rhoi,    &      !< refinement factor in i-direction
248   &  in_rhoj            !< refinement factor in j-direction
249
250   NAMELIST /namout/ &   !< output namlist
251   &  cn_fileout         !< fine grid bathymetry file
252   !-------------------------------------------------------------------
253
254   ! namelist
255   ! get namelist
256   il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
257   IF( il_narg/=1 )THEN
258      PRINT *,"ERROR in create_bathy: need a namelist"
259      STOP
260   ELSE
261      CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
262   ENDIF
263 
264   ! read namelist
265   INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
266   IF( ll_exist )THEN
267 
268      il_fileid=fct_getunit()
269
270      OPEN( il_fileid, FILE=TRIM(cl_namelist), &
271      &                FORM='FORMATTED',       &
272      &                ACCESS='SEQUENTIAL',    &
273      &                STATUS='OLD',           &
274      &                ACTION='READ',          &
275      &                IOSTAT=il_status)
276      CALL fct_err(il_status)
277      IF( il_status /= 0 )THEN
278         PRINT *,"ERROR in create_bathy: error opening "//TRIM(cl_namelist)
279         STOP
280      ENDIF
281
282      READ( il_fileid, NML = namlog )
283      ! define log file
284      CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
285      CALL logger_header()
286
287      READ( il_fileid, NML = namcfg )
288      ! get variable extra information
289      CALL var_def_extra(TRIM(cn_varcfg))
290
291      READ( il_fileid, NML = namcrs )
292      READ( il_fileid, NML = namfin )
293      READ( il_fileid, NML = namvar )
294      ! add user change in extra information
295      CALL var_chg_extra( cn_varinfo )
296      ! match variable with file
297      tl_multi=multi_init(cn_varfile)
298     
299      READ( il_fileid, NML = namnst )
300      READ( il_fileid, NML = namout )
301
302      CLOSE( il_fileid, IOSTAT=il_status )
303      CALL fct_err(il_status)
304      IF( il_status /= 0 )THEN
305         CALL logger_error("CREATE BATHY: closing "//TRIM(cl_namelist))
306      ENDIF
307
308   ELSE
309
310      PRINT *,"ERROR in create_bathy: can't find "//TRIM(cl_namelist)
311
312   ENDIF
313
314   CALL multi_print(tl_multi)
315
316   ! open files
317   IF( cn_coord0 /= '' )THEN
318      tl_coord0=mpp_init( file_init(TRIM(cn_coord0)), id_perio=in_perio0)
319      CALL grid_get_info(tl_coord0)
320   ELSE
321      CALL logger_fatal("CREATE BATHY: no coarse grid coordinate found. "//&
322      &     "check namelist")     
323   ENDIF
324
325   IF( TRIM(cn_coord1) /= '' )THEN
326      tl_coord1=mpp_init( file_init(TRIM(cn_coord1)),id_perio=in_perio1)
327      CALL grid_get_info(tl_coord1)
328   ELSE
329      CALL logger_fatal("CREATE BATHY: no fine grid coordinate found. "//&
330      &     "check namelist")
331   ENDIF
332
333   ! check
334   ! check output file do not already exist
335   INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
336   IF( ll_exist )THEN
337      CALL logger_fatal("CREATE BATHY: output file "//TRIM(cn_fileout)//&
338      &  " already exist.")
339   ENDIF
340
341   ! check namelist
342   ! check refinement factor
343   il_rho(:)=1
344   IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
345      CALL logger_error("CREATE BATHY: invalid refinement factor."//&
346      &  " check namelist "//TRIM(cl_namelist))
347   ELSE
348      il_rho(jp_I)=in_rhoi
349      il_rho(jp_J)=in_rhoj
350   ENDIF
351
352   ! check domain indices
353   ! compute coarse grid indices around fine grid
354   il_ind(:,:)=grid_get_coarse_index( tl_coord0, tl_coord1, &
355   &                                  id_rho=il_rho(:) )
356
357   il_imin0=il_ind(jp_I,1) ; il_imax0=il_ind(jp_I,2)
358   il_jmin0=il_ind(jp_J,1) ; il_jmax0=il_ind(jp_J,2)
359
360   ! check domain validity
361   CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0)
362
363   ! check coincidence between coarse and fine grid
364   CALL grid_check_coincidence( tl_coord0, tl_coord1, &
365   &                            il_imin0, il_imax0, &
366   &                            il_jmin0, il_jmax0, &
367   &                            il_rho(:) )
368
369   IF( .NOT. ASSOCIATED(tl_multi%t_mpp) )THEN
370      CALL logger_error("CREATE BATHY: no mpp file to work on. "//&
371      &                 "check cn_varfile in namelist.")
372   ELSE
373
374      ALLOCATE( tl_var( tl_multi%i_nvar ) )
375      jk=0
376      DO ji=1,tl_multi%i_nmpp
377     
378         WRITE(cl_data,'(a,i2.2)') 'data-',jk+1
379         IF( .NOT. ASSOCIATED(tl_multi%t_mpp(ji)%t_proc(1)%t_var) )THEN
380
381            CALL logger_fatal("CREATE BATHY: no variable to work on for "//&
382            &                 "mpp file"//TRIM(tl_multi%t_mpp(ji)%c_name)//&
383            &                 ". check cn_varfile in namelist.")
384
385         ELSEIF( TRIM(tl_multi%t_mpp(ji)%c_name) == TRIM(cl_data) )THEN
386
387            !- use input matrix to initialise variable
388            DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
389               jk=jk+1
390               tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
391           
392               tl_var(jk)=create_bathy_matrix(tl_tmp, tl_coord1)
393            ENDDO
394            ! clean
395            CALL var_clean(tl_tmp)
396
397         ELSE
398
399            tl_mpp=mpp_init( file_init(TRIM(tl_multi%t_mpp(ji)%c_name)) )
400            CALL grid_get_info(tl_mpp)
401
402            ! open mpp file
403            CALL iom_mpp_open(tl_mpp)
404
405            ! get or check depth value
406            IF( tl_multi%t_mpp(ji)%t_proc(1)%i_depthid /= 0 )THEN
407               il_varid=tl_multi%t_mpp(ji)%t_proc(1)%i_depthid
408               IF( ASSOCIATED(tl_depth%d_value) )THEN
409                  tl_tmp=iom_mpp_read_var(tl_mpp,il_varid)
410                  IF( ANY( tl_depth%d_value(:,:,:,:) /= &
411                  &        tl_tmp%d_value(:,:,:,:) ) )THEN
412                     CALL logger_fatal("CREATE BATHY: depth value from "//&
413                     &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//&
414                     &  " to those from former file(s).")
415                  ENDIF
416                  CALL var_clean(tl_tmp)
417               ELSE
418                  tl_depth=iom_mpp_read_var(tl_mpp,il_varid)
419               ENDIF
420            ENDIF
421
422            ! get or check time value
423            IF( tl_multi%t_mpp(ji)%t_proc(1)%i_timeid /= 0 )THEN
424               il_varid=tl_multi%t_mpp(ji)%t_proc(1)%i_timeid
425               IF( ASSOCIATED(tl_time%d_value) )THEN
426                  tl_tmp=iom_mpp_read_var(tl_mpp,il_varid)
427                  IF( ANY( tl_time%d_value(:,:,:,:) /= &
428                  &        tl_tmp%d_value(:,:,:,:) ) )THEN
429                     CALL logger_fatal("CREATE BATHY: time value from "//&
430                     &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//&
431                     &  " to those from former file(s).")
432                  ENDIF
433                  CALL var_clean(tl_tmp)
434               ELSE
435                  tl_time=iom_mpp_read_var(tl_mpp,il_varid)
436               ENDIF
437            ENDIF
438
439            ! close mpp file
440            CALL iom_mpp_close(tl_mpp)
441
442            IF( ANY( tl_mpp%t_dim(1:2)%i_len /= &
443            &        tl_coord0%t_dim(1:2)%i_len) )THEN
444               !- extract bathymetry from fine grid bathymetry
445               DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
446                  jk=jk+1
447                  tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
448               
449                  tl_var(jk)=create_bathy_extract( tl_tmp, tl_mpp, &
450                  &                                tl_coord1 )
451               ENDDO
452               ! clean
453               CALL var_clean(tl_tmp)
454            ELSE
455               !- get bathymetry from coarse grid bathymetry
456               DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar
457                  jk=jk+1
458                  tl_tmp=var_copy(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj))
459
460                  il_offset(:,:)= grid_get_fine_offset(tl_coord0,    &
461                  &                                    il_imin0, il_jmin0, &
462                  &                                    il_imax0, il_jmax0, &
463                  &                                    tl_coord1,          &
464                  &                                    il_rho(:),          &
465                  &                                    TRIM(tl_tmp%c_point) )
466
467                  tl_var(jk)=create_bathy_get_var( tl_tmp, tl_mpp,     &
468                  &                                il_imin0, il_jmin0, &
469                  &                                il_imax0, il_jmax0, &
470                  &                                il_offset(:,:),  &
471                  &                                il_rho(:) )
472               ENDDO
473               ! clean
474               CALL var_clean(tl_tmp)
475            ENDIF
476
477            ! clean structure
478            CALL mpp_clean(tl_mpp)
479
480         ENDIF
481      ENDDO
482   ENDIF
483
484   DO jk=1,tl_multi%i_nvar
485         ! forced min and max value
486         CALL var_limit_value(tl_var(jk))
487
488         ! fill closed sea
489         IF( ln_fillclosed )THEN
490            ALLOCATE( il_mask(tl_var(jk)%t_dim(1)%i_len, &
491            &                 tl_var(jk)%t_dim(2)%i_len) )
492
493            ! split domain in N sea subdomain
494            il_mask(:,:)=grid_split_domain(tl_var(jk))
495
496            !  fill smallest domain
497            CALL grid_fill_small_dom( tl_var(jk), il_mask(:,:) )
498
499            DEALLOCATE( il_mask )
500         ENDIF
501
502         ! filter
503         CALL filter_fill_value(tl_var(jk))
504
505         ! check bathymetry
506         dl_minbat=MINVAL(tl_var(jk)%d_value(:,:,:,:))
507         IF( TRIM(tl_var(jk)%c_stdname) == 'bathymetry' .AND. &
508         &   dl_minbat <= 0._dp  )THEN
509            CALL logger_debug("CREATE BATHY: min value "//TRIM(fct_str(dl_minbat)))
510            CALL logger_error("CREATE BATHY: Bathymetry has value <= 0")
511         ENDIF
512
513   ENDDO
514
515   ! create file
516   tl_fileout=file_init(TRIM(cn_fileout))
517
518   ! add dimension
519   tl_dim(:)=var_max_dim(tl_var(:))
520
521   DO ji=1,ip_maxdim
522      IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
523   ENDDO
524
525   ! add variables
526   IF( ALL( tl_dim(1:2)%l_use ) )THEN
527
528      ! open mpp files
529      CALL iom_mpp_open(tl_coord1)
530
531      ! add longitude
532      tl_lon=iom_mpp_read_var(tl_coord1,'longitude')
533      CALL file_add_var(tl_fileout, tl_lon)
534      CALL var_clean(tl_lon)
535
536      ! add latitude
537      tl_lat=iom_mpp_read_var(tl_coord1,'latitude')
538      CALL file_add_var(tl_fileout, tl_lat)
539      CALL var_clean(tl_lat)
540
541      ! close mpp files
542      CALL iom_mpp_close(tl_coord1)
543
544   ENDIF
545
546   IF( tl_dim(3)%l_use )THEN
547      ! add depth
548      CALL file_add_var(tl_fileout, tl_depth)
549      CALL var_clean(tl_depth)
550   ENDIF
551
552   IF( tl_dim(4)%l_use )THEN
553      ! add time
554      CALL file_add_var(tl_fileout, tl_time)
555      CALL var_clean(tl_time)
556   ENDIF
557
558   ! add other variables
559   DO jk=1,tl_multi%i_nvar
560      CALL file_add_var(tl_fileout, tl_var(jk))
561      CALL var_clean(tl_var(jk))
562   ENDDO
563   DEALLOCATE(tl_var)
564
565   ! add some attribute
566   tl_att=att_init("Created_by","SIREN create_bathy")
567   CALL file_add_att(tl_fileout, tl_att)
568
569   cl_date=date_print(date_now())
570   tl_att=att_init("Creation_date",cl_date)
571   CALL file_add_att(tl_fileout, tl_att)
572
573   ! add attribute periodicity
574   il_attid=0
575   IF( ASSOCIATED(tl_fileout%t_att) )THEN
576      il_attid=att_get_index(tl_fileout%t_att(:),'periodicity')
577   ENDIF
578   IF( tl_coord1%i_perio >= 0 .AND. il_attid == 0 )THEN
579      tl_att=att_init('periodicity',tl_coord1%i_perio)
580      CALL file_add_att(tl_fileout,tl_att)
581   ENDIF
582
583   ! add attribute east west overlap
584   il_attid=0
585   IF( ASSOCIATED(tl_fileout%t_att) )THEN
586      il_attid=att_get_index(tl_fileout%t_att(:),'ew_overlap')
587   ENDIF
588   IF( tl_coord1%i_ew >= 0 .AND. il_attid == 0 )THEN
589      tl_att=att_init('ew_overlap',tl_coord1%i_ew)
590      CALL file_add_att(tl_fileout,tl_att)
591   ENDIF
592   
593   ! create file
594   CALL iom_create(tl_fileout)
595
596   ! write file
597   CALL iom_write_file(tl_fileout)
598
599   ! close file
600   CALL iom_close(tl_fileout)
601
602   ! clean
603   CALL att_clean(tl_att)
604
605   CALL file_clean(tl_fileout)
606   CALL mpp_clean(tl_coord1)
607   CALL mpp_clean(tl_coord0)
608
609   ! close log file
610   CALL logger_footer()
611   CALL logger_close()
612
613CONTAINS
614   !-------------------------------------------------------------------
615   !> @brief
616   !> This function create variable, filled with matrix value
617   !>
618   !> @details
619   !> A variable is create with the same name that the input variable,
620   !> and with dimension of the coordinate file.<br/>
621   !> Then the variable array of value is split into equal subdomain.
622   !> Each subdomain is filled with the corresponding value of the matrix.
623   !>
624   !> @author J.Paul
625   !> - November, 2013- Initial Version
626   !>
627   !> @param[in] td_var    variable structure
628   !> @param[in] td_coord  coordinate file structure
629   !> @return variable structure
630   !-------------------------------------------------------------------
631   FUNCTION create_bathy_matrix(td_var, td_coord)
632      IMPLICIT NONE
633      ! Argument
634      TYPE(TVAR), INTENT(IN) :: td_var
635      TYPE(TMPP), INTENT(IN) :: td_coord
636
637      ! function
638      TYPE(TVAR) :: create_bathy_matrix
639
640      ! local variable
641      INTEGER(i4)      , DIMENSION(2,2)                  :: il_xghost
642      INTEGER(i4)      , DIMENSION(3)                    :: il_dim
643      INTEGER(i4)      , DIMENSION(3)                    :: il_size
644      INTEGER(i4)      , DIMENSION(3)                    :: il_rest
645
646      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_ishape
647      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_jshape
648      INTEGER(i4)      , DIMENSION(:)      , ALLOCATABLE :: il_kshape
649
650      REAL(dp)         , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value
651
652      TYPE(TVAR)                                         :: tl_lon
653      TYPE(TDIM)       , DIMENSION(ip_maxdim)            :: tl_dim
654
655      TYPE(TMPP)                                         :: tl_coord
656
657      ! loop indices
658      INTEGER(i4) :: ji
659      INTEGER(i4) :: jj
660      INTEGER(i4) :: jk
661      !----------------------------------------------------------------
662
663      ! copy structure
664      tl_coord=mpp_copy(td_coord)
665
666      ! use only edge processor
667      CALL mpp_get_contour(tl_coord)
668
669      ! open useful processor
670      CALL iom_mpp_open(tl_coord)
671
672      ! read output grid
673      tl_lon=iom_mpp_read_var(tl_coord,'longitude')
674
675      ! look for ghost cell
676      il_xghost(:,:)=grid_get_ghost( tl_coord )
677
678      ! close processor
679      CALL iom_mpp_close(tl_coord)
680      ! clean
681      CALL mpp_clean(tl_coord)
682
683      ! remove ghost cell
684      CALL grid_del_ghost(tl_lon, il_xghost(:,:))
685
686      ! write value on grid
687      ! get matrix dimension
688      il_dim(:)=td_var%t_dim(1:3)%i_len
689      ! output dimension
690      tl_dim(:)=dim_copy(tl_lon%t_dim(:))
691      ! clean
692      CALL var_clean(tl_lon)
693
694      ! split output domain in N subdomain depending of matrix dimension
695      il_size(:) = tl_dim(1:3)%i_len / il_dim(:)
696      il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:))
697
698      ALLOCATE( il_ishape(il_dim(1)+1) )
699      il_ishape(:)=0
700      DO ji=2,il_dim(1)+1
701         il_ishape(ji)=il_ishape(ji-1)+il_size(1)
702      ENDDO
703      ! add rest to last cell
704      il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1)
705
706      ALLOCATE( il_jshape(il_dim(2)+1) )
707      il_jshape(:)=0
708      DO jj=2,il_dim(2)+1
709         il_jshape(jj)=il_jshape(jj-1)+il_size(2)
710      ENDDO
711      ! add rest to last cell
712      il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2)
713
714      ALLOCATE( il_kshape(il_dim(3)+1) )
715      il_kshape(:)=0
716      DO jk=2,il_dim(3)+1
717         il_kshape(jk)=il_kshape(jk-1)+il_size(3)
718      ENDDO
719      ! add rest to last cell
720      il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3)
721
722      ! write ouput array of value
723      ALLOCATE(dl_value( tl_dim(1)%i_len, &
724      &                  tl_dim(2)%i_len, &
725      &                  tl_dim(3)%i_len, &
726      &                  tl_dim(4)%i_len) )
727
728      dl_value(:,:,:,:)=0
729
730      DO jk=2,il_dim(3)+1
731         DO jj=2,il_dim(2)+1
732            DO ji=2,il_dim(1)+1
733               
734               dl_value( 1+il_ishape(ji-1):il_ishape(ji), &
735               &         1+il_jshape(jj-1):il_jshape(jj), &
736               &         1+il_kshape(jk-1):il_kshape(jk), &
737               &         1 ) = td_var%d_value(ji-1,jj-1,jk-1,1)
738
739            ENDDO
740         ENDDO
741      ENDDO
742
743      ! initialise variable with value
744      create_bathy_matrix=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:))
745
746      DEALLOCATE(dl_value)
747
748      ! add ghost cell
749      CALL grid_add_ghost(create_bathy_matrix, il_xghost(:,:))
750
751      ! clean
752      CALL dim_clean(tl_dim(:))
753
754   END FUNCTION create_bathy_matrix
755   !-------------------------------------------------------------------
756   !> @brief
757   !> This function extract variable from file over coordinate domain and
758   !> return variable structure
759   !>
760   !> @author J.Paul
761   !> - November, 2013- Initial Version
762   !>
763   !> @param[in] td_var    variable structure
764   !> @param[in] td_mpp    mpp file structure
765   !> @param[in] td_coord  coordinate file structure
766   !> @return variable structure
767   !-------------------------------------------------------------------
768   FUNCTION create_bathy_extract(td_var, td_mpp, &
769   &                             td_coord)
770      IMPLICIT NONE
771      ! Argument
772      TYPE(TVAR), INTENT(IN) :: td_var 
773      TYPE(TMPP), INTENT(IN) :: td_mpp
774      TYPE(TMPP), INTENT(IN) :: td_coord
775
776      ! function
777      TYPE(TVAR) :: create_bathy_extract
778
779      ! local variable
780      INTEGER(i4), DIMENSION(2,2) :: il_ind
781
782      INTEGER(i4) :: il_imin
783      INTEGER(i4) :: il_jmin
784      INTEGER(i4) :: il_imax
785      INTEGER(i4) :: il_jmax
786
787      TYPE(TMPP)  :: tl_mpp
788
789      TYPE(TATT)  :: tl_att
790
791      TYPE(TVAR)  :: tl_var
792     
793      TYPE(TDOM)  :: tl_dom
794      ! loop indices
795      !----------------------------------------------------------------
796
797      IF( .NOT. ASSOCIATED(td_mpp%t_proc) )THEN
798         CALL logger_error("CREATE BATHY EXTRACT: no processor associated "//&
799         &  "to mpp "//TRIM(td_mpp%c_name))
800      ELSE
801
802         !init
803         tl_mpp=mpp_copy(td_mpp)
804
805         ! compute file grid indices around coord grid
806         il_ind(:,:)=grid_get_coarse_index(tl_mpp, td_coord )
807
808         il_imin=il_ind(1,1) ; il_imax=il_ind(1,2)
809         il_jmin=il_ind(2,1) ; il_jmax=il_ind(2,2)
810
811         ! check grid coincidence
812         CALL grid_check_coincidence( tl_mpp, td_coord, &
813         &                            il_imin, il_imax, &
814         &                            il_jmin, il_jmax, &
815         &                            (/1, 1, 1/) )
816
817         ! compute domain
818         tl_dom=dom_init(tl_mpp,           &
819         &               il_imin, il_imax, &
820         &               il_jmin, il_jmax)
821
822         ! open mpp files over domain
823         CALL iom_dom_open(tl_mpp, tl_dom)
824
825         ! read variable on domain
826         tl_var=iom_dom_read_var(tl_mpp,TRIM(td_var%c_name),tl_dom)
827
828         ! close mpp file
829         CALL iom_dom_close(tl_mpp)
830
831         ! add ghost cell
832         CALL grid_add_ghost(tl_var,tl_dom%i_ghost(:,:))
833
834         ! check result
835         IF( ANY( tl_var%t_dim(:)%l_use .AND. &
836         &        tl_var%t_dim(:)%i_len /= td_coord%t_dim(:)%i_len) )THEN
837            CALL logger_debug("CREATE BATHY EXTRACT: "//&
838            &        "dimensoin of variable "//TRIM(td_var%c_name)//" "//&
839            &        TRIM(fct_str(tl_var%t_dim(1)%i_len))//","//&
840            &        TRIM(fct_str(tl_var%t_dim(2)%i_len))//","//&
841            &        TRIM(fct_str(tl_var%t_dim(3)%i_len))//","//&
842            &        TRIM(fct_str(tl_var%t_dim(4)%i_len)) )
843            CALL logger_debug("CREATE BATHY EXTRACT: "//&
844            &        "dimensoin of coordinate file "//&
845            &        TRIM(fct_str(td_coord%t_dim(1)%i_len))//","//&
846            &        TRIM(fct_str(td_coord%t_dim(2)%i_len))//","//&
847            &        TRIM(fct_str(td_coord%t_dim(3)%i_len))//","//&
848            &        TRIM(fct_str(td_coord%t_dim(4)%i_len)) )
849            CALL logger_fatal("CREATE BATHY EXTRACT: "//&
850            &  "dimensoin of extracted "//&
851            &  "variable and coordinate file dimension differ")
852         ENDIF
853
854         ! add attribute to variable
855         tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
856         CALL var_move_att(tl_var, tl_att)         
857
858         tl_att=att_init('src_i_indices',(/tl_dom%i_imin, tl_dom%i_imax/))
859         CALL var_move_att(tl_var, tl_att)
860
861         tl_att=att_init('src_j_indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
862         CALL var_move_att(tl_var, tl_att)
863
864         ! save result
865         create_bathy_extract=var_copy(tl_var)
866
867         ! clean structure
868         CALL att_clean(tl_att)
869         CALL var_clean(tl_var)
870         CALL mpp_clean(tl_mpp)
871      ENDIF
872
873   END FUNCTION create_bathy_extract
874   !-------------------------------------------------------------------
875   !> @brief
876   !> This function get coarse grid variable, interpolate variable, and return
877   !> variable structure over fine grid
878   !>
879   !> @author J.Paul
880   !> - November, 2013- Initial Version
881   !>
882   !> @param[in] td_var    variable structure
883   !> @param[in] td_mpp    mpp file structure
884   !> @param[in] id_imin   i-direction lower left  corner indice
885   !> @param[in] id_imax   i-direction upper right corner indice
886   !> @param[in] id_jmin   j-direction lower left  corner indice
887   !> @param[in] id_jmax   j-direction upper right corner indice
888   !> @param[in] id_offset offset between fine grid and coarse grid
889   !> @param[in] id_rho    array of refinement factor
890   !> @return variable structure
891   !-------------------------------------------------------------------
892   FUNCTION create_bathy_get_var(td_var, td_mpp,   &
893            &                    id_imin, id_jmin, &
894            &                    id_imax, id_jmax, &
895            &                    id_offset,        &
896            &                    id_rho )
897      IMPLICIT NONE
898      ! Argument
899      TYPE(TVAR) , INTENT(IN) :: td_var 
900      TYPE(TMPP) , INTENT(IN) :: td_mpp 
901      INTEGER(i4), INTENT(IN) :: id_imin
902      INTEGER(i4), INTENT(IN) :: id_imax
903      INTEGER(i4), INTENT(IN) :: id_jmin
904      INTEGER(i4), INTENT(IN) :: id_jmax
905      INTEGER(i4), DIMENSION(:,:), INTENT(IN) :: id_offset
906      INTEGER(i4), DIMENSION(:)  , INTENT(IN) :: id_rho
907
908      ! function
909      TYPE(TVAR) :: create_bathy_get_var
910
911      ! local variable
912      TYPE(TMPP)  :: tl_mpp
913      TYPE(TATT)  :: tl_att
914      TYPE(TVAR)  :: tl_var
915      TYPE(TDOM)  :: tl_dom
916
917      INTEGER(i4) :: il_size
918      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho
919
920      ! loop indices
921      !----------------------------------------------------------------
922      IF( ANY(SHAPE(id_offset(:,:)) /= 2) )THEN
923         CALL logger_error("CREATE BATHY GET VAR: invalid dimension of "//&
924         &                 "offset array")
925      ENDIF
926
927      ! copy structure
928      tl_mpp=mpp_copy(td_mpp)
929
930      !- compute domain
931      tl_dom=dom_init(tl_mpp,           &
932      &               id_imin, id_imax, &
933      &               id_jmin, id_jmax)
934
935      !- add extra band (if possible) to compute interpolation
936      CALL dom_add_extra(tl_dom)
937
938      !- open mpp files over domain
939      CALL iom_dom_open(tl_mpp, tl_dom)
940
941      !- read variable value on domain
942      tl_var=iom_dom_read_var(tl_mpp,TRIM(td_var%c_name),tl_dom)
943
944      !- close mpp files
945      CALL iom_dom_close(tl_mpp)
946
947      il_size=SIZE(id_rho(:))
948      ALLOCATE( il_rho(il_size) )
949      il_rho(:)=id_rho(:)
950     
951      !- interpolate variable
952      CALL create_bathy_interp(tl_var, il_rho(:), id_offset(:,:))
953
954      !- remove extraband added to domain
955      CALL dom_del_extra( tl_var, tl_dom, il_rho(:) )
956
957      !- add ghost cell
958      CALL grid_add_ghost(tl_var,tl_dom%i_ghost(:,:))
959 
960      !- add attribute to variable
961      tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name)))
962      CALL var_move_att(tl_var, tl_att)
963
964      tl_att=att_init('src_i_indices',(/tl_dom%i_imin, tl_dom%i_imax/))
965      CALL var_move_att(tl_var, tl_att)
966
967      tl_att=att_init('src_j_indices',(/tl_dom%i_jmin, tl_dom%i_jmax/))
968      CALL var_move_att(tl_var, tl_att)
969
970      IF( .NOT. ALL(id_rho(:)==1) )THEN
971         tl_att=att_init("refinment_factor",(/id_rho(jp_I),id_rho(jp_J)/))
972         CALL var_move_att(tl_var, tl_att)
973      ENDIF
974
975      DEALLOCATE( il_rho )
976
977      !- save result
978      create_bathy_get_var=var_copy(tl_var)
979
980      !- clean structure
981      CALL att_clean(tl_att)
982      CALL var_clean(tl_var)
983      CALL mpp_clean(tl_mpp)
984 
985   END FUNCTION create_bathy_get_var
986   !-------------------------------------------------------------------
987   !> @brief
988   !> This subroutine interpolate variable
989   !>
990   !> @author J.Paul
991   !> - November, 2013- Initial Version
992   !>
993   !> @param[inout] td_var variable structure
994   !> @param[in] id_rho    array of refinment factor
995   !> @param[in] id_offset array of offset between fine and coarse grid
996   !> @param[in] id_iext   i-direction size of extra bands (default=im_minext)
997   !> @param[in] id_jext   j-direction size of extra bands (default=im_minext)
998   !-------------------------------------------------------------------
999   SUBROUTINE create_bathy_interp( td_var,         &
1000   &                               id_rho,          &
1001   &                               id_offset,       &
1002   &                               id_iext, id_jext)
1003
1004      IMPLICIT NONE
1005
1006      ! Argument
1007      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
1008      INTEGER(i4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
1009      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_offset
1010      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_iext
1011      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_jext
1012
1013      ! local variable
1014      TYPE(TVAR)  :: tl_mask
1015
1016      INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
1017
1018      INTEGER(i4) :: il_iext
1019      INTEGER(i4) :: il_jext
1020
1021      ! loop indices
1022      !----------------------------------------------------------------
1023
1024      !WARNING: two extrabands are required for cubic interpolation
1025      il_iext=3
1026      IF( PRESENT(id_iext) ) il_iext=id_iext
1027
1028      il_jext=3
1029      IF( PRESENT(id_jext) ) il_jext=id_jext
1030
1031      IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1032         CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
1033         &  "on two points are required with cubic interpolation ")
1034         il_iext=2
1035      ENDIF
1036
1037      IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
1038         CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//&
1039         &  "on two points are required with cubic interpolation ")
1040         il_jext=2
1041      ENDIF
1042
1043      ! work on mask
1044      ! create mask
1045      ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
1046      &                td_var%t_dim(2)%i_len, &
1047      &                td_var%t_dim(3)%i_len, &
1048      &                td_var%t_dim(4)%i_len) )
1049
1050      bl_mask(:,:,:,:)=1
1051      WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0     
1052
1053      SELECT CASE(TRIM(td_var%c_point))
1054      CASE DEFAULT ! 'T'
1055         tl_mask=var_init('tmask', bl_mask(:,:,:,:), td_dim=td_var%t_dim(:), &
1056         &                id_ew=td_var%i_ew )
1057      CASE('U','V','F')
1058         CALL logger_fatal("CREATE BATHY INTERP: can not computed "//&
1059         &                 "interpolation on "//TRIM(td_var%c_point)//&
1060         &                 " grid point (variable "//TRIM(td_var%c_name)//&
1061         &                 "). check namelist.")
1062      END SELECT
1063
1064      DEALLOCATE(bl_mask)
1065
1066      ! interpolate mask
1067      CALL interp_fill_value( tl_mask, id_rho(:), &
1068      &                       id_offset=id_offset(:,:) )
1069
1070      ! work on variable
1071      ! add extraband
1072      CALL extrap_add_extrabands(td_var, il_iext, il_jext)
1073
1074      ! extrapolate variable
1075      CALL extrap_fill_value( td_var, id_offset=id_offset(:,:), &
1076      &                               id_rho=id_rho(:),         &
1077      &                               id_iext=il_iext, id_jext=il_jext )
1078
1079      ! interpolate Bathymetry
1080      CALL interp_fill_value( td_var, id_rho(:), &
1081      &                       id_offset=id_offset(:,:) )
1082
1083      ! remove extraband
1084      CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J))
1085
1086      ! keep original mask
1087      WHERE( tl_mask%d_value(:,:,:,:) == 0 )
1088         td_var%d_value(:,:,:,:)=td_var%d_fill
1089      END WHERE
1090
1091      ! clean variable structure
1092      CALL var_clean(tl_mask)
1093
1094   END SUBROUTINE create_bathy_interp
1095END PROGRAM create_bathy
Note: See TracBrowser for help on using the repository browser.