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

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

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/create_coord.f90 @ 5214

Last change on this file since 5214 was 5214, checked in by davestorkey, 9 years ago

Merge in changes from the trunk up to rev 5107.

File size: 22.2 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5!
6! PROGRAM: create_coord
7!
8! DESCRIPTION:
9!> @file
10!> @brief
11!> This program create fine grid coordinate file.
12!>
13!> @details
14!> @section sec1 method
15!>    All variables from the input coordinates coarse grid file, are extracted
16!>    and interpolated to create fine grid coordinates files.<br/>
17!>    @note
18!>       interpolation method could be different for each variable.
19!>
20!> @section sec2 how to
21!>    to create fine grid coordinates files:<br/>
22!> @code{.sh}
23!>    ./SIREN/bin/create_coord create_coord.nam
24!> @endcode
25!>   
26!>    create_coord.nam comprise 6 namelists:<br/>
27!>       - logger namelist (namlog)
28!>       - config namelist (namcfg)
29!>       - coarse grid namelist (namcrs)
30!>       - variable namelist (namvar)
31!>       - nesting namelist (namnst)
32!>       - output namelist (namout)
33!>   
34!>    @note
35!>       All namelists have to be in file create_coord.nam,
36!>       however variables of those namelists are all optional.
37!>
38!>    * _logger namelist (namlog)_:<br/>
39!>       - cn_logfile   : log filename
40!>       - cn_verbosity : verbosity ('trace','debug','info',
41!> 'warning','error','fatal')
42!>       - in_maxerror  : maximum number of error allowed
43!>
44!>    * _config namelist (namcfg)_:<br/>
45!>       - cn_varcfg : variable configuration file
46!> (see ./SIREN/cfg/variable.cfg)
47!>
48!>    * _coarse grid namelist (namcrs)_:<br/>
49!>       - cn_coord0 : coordinate file
50!>       - in_perio0 : NEMO periodicity index (see Model Boundary Condition in
51!> [NEMO documentation](http://www.nemo-ocean.eu/About-NEMO/Reference-manuals))
52!>
53!>    * _variable namelist (namvar)_:<br/>
54!>       - cn_varinfo : list of variable and extra information about request(s)
55!> to be used.<br/>
56!>          each elements of *cn_varinfo* is a string character.<br/>
57!>          it is composed of the variable name follow by ':',
58!>          then request(s) to be used on this variable.<br/>
59!>          request could be:
60!>             - interpolation method
61!>             - extrapolation method
62!>             - filter method
63!>
64!>                requests must be separated by ';' .<br/>
65!>                order of requests does not matter.<br/>
66!>
67!>          informations about available method could be find in @ref interp,
68!>          @ref extrap and @ref filter modules.<br/>
69!>
70!>          Example: 'votemper: linear; hann(2,3); dist_weight',
71!>          'vosaline: cubic'<br/>
72!>          @note
73!>             If you do not specify a method which is required,
74!>             default one is applied.
75!>
76!>    * _nesting namelist (namnst)_:<br/>
77!>       - in_imin0 : i-direction lower left  point indice
78!> of coarse grid subdomain to be used
79!>       - in_imax0 : i-direction upper right point indice
80!> of coarse grid subdomain to be used
81!>       - in_jmin0 : j-direction lower left  point indice
82!> of coarse grid subdomain to be used
83!>       - in_jmax0 : j-direction upper right point indice
84!> of coarse grid subdomain to be used
85!>       - in_rhoi  : refinement factor in i-direction
86!>       - in_rhoj  : refinement factor in j-direction<br/>
87!>
88!>       \image html  grid_zoom_40.png
89!>       \image latex grid_zoom_40.png
90!>
91!>    * _output namelist (namout)_:
92!>       - cn_fileout : output coordinate file
93!>
94!> @author J.Paul
95! REVISION HISTORY:
96!> @date November, 2013 - Initial Version
97!> @date September, 2014
98!> - add header for user
99!> - compute offset considering grid point
100!> - add global attributes in output file
101!>
102!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
103!----------------------------------------------------------------------
104PROGRAM create_coord
105
106   USE global                          ! global variable
107   USE kind                            ! F90 kind parameter
108   USE logger                          ! log file manager
109   USE fct                             ! basic useful function
110   USE date                            ! date manager
111   USE att                             ! attribute manager
112   USE dim                             ! dimension manager
113   USE var                             ! variable manager
114   USE file                            ! file manager
115   USE iom                             ! I/O manager
116   USE grid                            ! grid manager
117   USE extrap                          ! extrapolation manager
118   USE interp                          ! interpolation manager
119   USE filter                          ! filter manager
120   USE mpp                             ! MPP manager
121   USE dom                             ! domain manager
122   USE iom_mpp                         ! MPP I/O manager
123   USE iom_dom                         ! DOM I/O manager
124
125   IMPLICIT NONE
126
127   ! local variable
128   CHARACTER(LEN=lc)                                    :: cl_namelist
129   CHARACTER(LEN=lc)                                    :: cl_date
130
131   INTEGER(i4)                                          :: il_narg
132   INTEGER(i4)                                          :: il_status
133   INTEGER(i4)                                          :: il_fileid
134   INTEGER(i4)                                          :: il_attid
135   INTEGER(i4)                                          :: il_ind
136   INTEGER(i4)                                          :: il_nvar
137   INTEGER(i4)                                          :: il_ew
138   INTEGER(i4)      , DIMENSION(ip_maxdim)              :: il_rho
139   INTEGER(i4)      , DIMENSION(2,2,ip_npoint)          :: il_offset
140
141   LOGICAL                                              :: ll_exist
142
143   TYPE(TATT)                                           :: tl_att
144
145   TYPE(TDOM)                                           :: tl_dom
146
147   TYPE(TVAR)       , DIMENSION(:)        , ALLOCATABLE :: tl_var
148
149   TYPE(TDIM)       , DIMENSION(ip_maxdim)              :: tl_dim
150
151   TYPE(TMPP)                                           :: tl_coord0
152   TYPE(TFILE)                                          :: tl_fileout
153
154   ! check
155!   INTEGER(i4)                                          :: il_imin0
156!   INTEGER(i4)                                          :: il_imax0
157!   INTEGER(i4)                                          :: il_jmin0
158!   INTEGER(i4)                                          :: il_jmax0
159!   INTEGER(i4)      , DIMENSION(2,2)                    :: il_ind2
160!   TYPE(TMPP)                                           :: tl_mppout
161
162   ! loop indices
163   INTEGER(i4) :: ji
164   INTEGER(i4) :: jj
165
166   ! namelist variable
167   CHARACTER(LEN=lc) :: cn_logfile = 'create_coord.log' 
168   CHARACTER(LEN=lc) :: cn_verbosity = 'warning' 
169   INTEGER(i4)       :: in_maxerror = 5
170
171   CHARACTER(LEN=lc) :: cn_coord0 = '' 
172   INTEGER(i4)       :: in_perio0 = -1
173
174   CHARACTER(LEN=lc) :: cn_varcfg = '../cfg/variable.cfg' 
175
176   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''
177
178   INTEGER(i4)       :: in_imin0 = 0
179   INTEGER(i4)       :: in_imax0 = 0
180   INTEGER(i4)       :: in_jmin0 = 0
181   INTEGER(i4)       :: in_jmax0 = 0
182   INTEGER(i4)       :: in_rhoi  = 1
183   INTEGER(i4)       :: in_rhoj  = 1
184
185   CHARACTER(LEN=lc) :: cn_fileout= 'coord_fine.nc'
186   !-------------------------------------------------------------------
187
188   NAMELIST /namlog/ &  !  logger namelist
189   &  cn_logfile,    &  !< logger file name
190   &  cn_verbosity,  &  !< logger verbosity
191   &  in_maxerror       !< logger maximum error
192
193   NAMELIST /namcfg/ &  !  config namelist
194   &  cn_varcfg         !< variable configuration file
195
196   NAMELIST /namcrs/ &  !  coarse grid namelist
197   &  cn_coord0 , &     !< coordinate file
198   &  in_perio0         !< periodicity index
199
200   NAMELIST /namvar/ &  !  variable namelist
201   &  cn_varinfo        !< list of variable and extra information about
202                        !< interpolation, extrapolation or filter method to be used.
203                        !< (ex: 'votemper:linear,hann,dist_weight','vosaline:cubic' )
204   
205   NAMELIST /namnst/ &  !  nesting namelist
206   &  in_imin0,   &     !< i-direction lower left  point indice
207   &  in_imax0,   &     !< i-direction upper right point indice
208   &  in_jmin0,   &     !< j-direction lower left  point indice
209   &  in_jmax0,   &     !< j-direction upper right point indice
210   &  in_rhoi,    &     !< refinement factor in i-direction
211   &  in_rhoj           !< refinement factor in j-direction
212
213   NAMELIST /namout/ &  !  output namelist
214   &  cn_fileout        !< fine grid coordinate file   
215   !-------------------------------------------------------------------
216
217   ! namelist
218   ! get namelist
219   il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
220   IF( il_narg/=1 )THEN
221      PRINT *,"ERROR in create_coord: need a namelist"
222      STOP
223   ELSE
224      CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec
225   ENDIF
226   
227   ! read namelist
228   INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
229   IF( ll_exist )THEN
230 
231      il_fileid=fct_getunit()
232
233      OPEN( il_fileid, FILE=TRIM(cl_namelist),   &
234      &                FORM='FORMATTED',         &
235      &                ACCESS='SEQUENTIAL',      &
236      &                STATUS='OLD',             &
237      &                ACTION='READ',            &
238      &                IOSTAT=il_status)
239      CALL fct_err(il_status)
240      IF( il_status /= 0 )THEN
241         PRINT *,"ERROR in create_coord: error opening "//TRIM(cl_namelist)
242         STOP
243      ENDIF
244
245      READ( il_fileid, NML = namlog )
246      ! define logger file
247      CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
248      CALL logger_header()
249
250      READ( il_fileid, NML = namcfg )
251      ! get variable extra information on configuration file
252      CALL var_def_extra(TRIM(cn_varcfg))
253
254      READ( il_fileid, NML = namcrs )
255      READ( il_fileid, NML = namvar )
256      ! add user change in extra information
257      CALL var_chg_extra( cn_varinfo )
258
259      READ( il_fileid, NML = namnst )
260      READ( il_fileid, NML = namout )
261
262      CLOSE( il_fileid, IOSTAT=il_status )
263      CALL fct_err(il_status)
264      IF( il_status /= 0 )THEN
265         CALL logger_error("CREATE COORD: closing "//TRIM(cl_namelist))
266      ENDIF
267
268   ELSE
269
270      PRINT *,"ERROR in create_coord: can't find "//TRIM(cl_namelist)
271      STOP
272
273   ENDIF
274
275   ! open files
276   IF( cn_coord0 /= '' )THEN
277      tl_coord0=mpp_init( file_init(TRIM(cn_coord0)), id_perio=in_perio0)
278      CALL grid_get_info(tl_coord0)
279   ELSE
280      CALL logger_fatal("CREATE COORD: no coarse grid coordinate found. "//&
281      &     "check namelist")     
282   ENDIF
283
284   ! check
285   ! check output file do not already exist
286   INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
287   IF( ll_exist )THEN
288      CALL logger_fatal("CREATE COORD: output file "//TRIM(cn_fileout)//&
289      &  " already exist.")
290   ENDIF
291
292   ! check nesting parameters
293   IF( in_imin0 < 0 .OR. in_imax0 < 0 .OR. in_jmin0 < 0 .OR. in_jmax0 < 0)THEN
294      CALL logger_fatal("CREATE COORD: invalid points indices."//&
295      &  " check namelist "//TRIM(cl_namelist))
296   ENDIF
297
298   il_rho(:)=1
299   IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
300      CALL logger_error("CREATE COORD: invalid refinement factor."//&
301      &  " check namelist "//TRIM(cl_namelist))
302   ELSE
303      il_rho(jp_I)=in_rhoi
304      il_rho(jp_J)=in_rhoj     
305
306      il_offset(:,:,:)=create_coord_get_offset(il_rho(:))
307
308   ENDIF
309
310   ! check domain validity
311   CALL grid_check_dom(tl_coord0, in_imin0, in_imax0, in_jmin0, in_jmax0 )
312
313   ! compute domain
314   tl_dom=dom_init( tl_coord0,         &
315   &                in_imin0, in_imax0,&
316   &                in_jmin0, in_jmax0 )
317
318   ! add extra band (if need be) to compute interpolation
319   CALL dom_add_extra(tl_dom)
320
321   ! open mpp files
322   CALL iom_dom_open(tl_coord0, tl_dom)
323
324   il_nvar=tl_coord0%t_proc(1)%i_nvar
325   ALLOCATE( tl_var(il_nvar) )
326   DO ji=1,il_nvar
327
328      tl_var(ji)=iom_dom_read_var(tl_coord0, &
329      &                          TRIM(tl_coord0%t_proc(1)%t_var(ji)%c_name),&
330      &                          tl_dom)
331
332      SELECT CASE(TRIM(tl_var(ji)%c_point))
333         CASE('T')
334            jj=jp_T
335         CASE('U')
336            jj=jp_U
337         CASE('V')
338            jj=jp_V
339         CASE('F')
340            jj=jp_F
341      END SELECT
342
343      ! interpolate variables
344      CALL create_coord_interp( tl_var(ji), il_rho(:), &
345      &                         il_offset(:,:,jj) )
346
347      ! remove extraband added to domain
348      CALL dom_del_extra( tl_var(ji), tl_dom, il_rho(:), .true. )
349
350      ! do not add ghost cell.
351      ! ghost cell already replace by value for coordinates
352      ! CALL grid_add_ghost(tl_var(ji),tl_dom%i_ghost(:,:))
353
354      ! filter
355      CALL filter_fill_value(tl_var(ji))     
356
357   ENDDO
358
359   ! close mpp files
360   CALL iom_dom_close(tl_coord0)
361
362   ! clean
363   CALL mpp_clean(tl_coord0)
364
365   ! create file
366   tl_fileout=file_init(TRIM(cn_fileout))
367
368   ! add dimension
369   ! save biggest dimension
370   tl_dim(:)=var_max_dim(tl_var(:))
371
372   DO ji=1,ip_maxdim
373      IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
374   ENDDO
375
376   ! add variables
377   DO ji=1,il_nvar
378      CALL file_add_var(tl_fileout, tl_var(ji))
379   ENDDO
380
381   ! recompute some attribute
382
383   ! add some attribute
384   tl_att=att_init("Created_by","SIREN create_coord")
385   CALL file_add_att(tl_fileout, tl_att)
386
387   cl_date=date_print(date_now())
388   tl_att=att_init("Creation_date",cl_date)
389   CALL file_add_att(tl_fileout, tl_att)
390
391   tl_att=att_init("src_file",TRIM(fct_basename(cn_coord0)))
392   CALL file_add_att(tl_fileout, tl_att)   
393
394   tl_att=att_init("src_i_indices",(/in_imin0,in_imax0/))
395   CALL file_add_att(tl_fileout, tl_att)   
396   tl_att=att_init("src_j_indices",(/in_jmin0,in_jmax0/))
397   CALL file_add_att(tl_fileout, tl_att)
398   IF( .NOT. ALL(il_rho(:)==1) )THEN
399      tl_att=att_init("refinment_factor",(/il_rho(jp_I),il_rho(jp_J)/))
400      CALL file_add_att(tl_fileout, tl_att)
401   ENDIF
402
403   ! add attribute periodicity
404   il_attid=0
405   IF( ASSOCIATED(tl_fileout%t_att) )THEN
406      il_attid=att_get_index(tl_fileout%t_att(:),'periodicity')
407   ENDIF
408   IF( tl_dom%i_perio >= 0 .AND. il_attid == 0 )THEN
409      tl_att=att_init('periodicity',tl_dom%i_perio)
410      CALL file_add_att(tl_fileout,tl_att)
411   ENDIF
412
413   ! add attribute east west overlap
414   il_attid=0
415   IF( ASSOCIATED(tl_fileout%t_att) )THEN
416      il_attid=att_get_index(tl_fileout%t_att(:),'ew_overlap')
417   ENDIF
418   IF( il_attid == 0 )THEN
419      il_ind=var_get_index(tl_fileout%t_var(:),'longitude')
420      il_ew=grid_get_ew_overlap(tl_fileout%t_var(il_ind))
421      IF( il_ew >= 0 )THEN
422         tl_att=att_init('ew_overlap',il_ew)
423         CALL file_add_att(tl_fileout,tl_att)
424      ENDIF
425   ENDIF
426
427   ! create file
428   CALL iom_create(tl_fileout)
429
430   ! write file
431   CALL iom_write_file(tl_fileout)
432
433   ! close file
434   CALL iom_close(tl_fileout)
435
436   ! clean
437   CALL att_clean(tl_att)
438   CALL var_clean(tl_var(:))
439   DEALLOCATE( tl_var) 
440
441   CALL file_clean(tl_fileout)
442
443!   ! check domain
444!   tl_coord0=mpp_init( file_init(TRIM(cn_coord0)), id_perio=in_perio0)
445!   tl_mppout=mpp_init( file_init(TRIM(cn_fileout)) )
446!   CALL grid_get_info(tl_coord0)
447!   CALL iom_mpp_open(tl_mppout)
448!
449!   il_ind2(:,:)=grid_get_coarse_index( tl_coord0, tl_mppout, &
450!   &                                   id_rho=il_rho(:) )
451!
452!   il_imin0=il_ind2(1,1) ; il_imax0=il_ind2(1,2)
453!   il_jmin0=il_ind2(2,1) ; il_jmax0=il_ind2(2,2)
454!
455!   IF( il_imin0 /= in_imin0 .OR. &
456!   &   il_imax0 /= in_imax0 .OR. &
457!   &   il_jmin0 /= in_jmin0 .OR. &
458!   &   il_jmax0 /= in_jmax0 )THEN
459!      CALL logger_debug("CREATE COORD: output indices ("//&
460!      &                 TRIM(fct_str(il_imin0))//","//&
461!      &                 TRIM(fct_str(il_imax0))//") ("//&
462!      &                 TRIM(fct_str(il_jmin0))//","//&
463!      &                 TRIM(fct_str(il_jmax0))//")" )
464!      CALL logger_debug("CREATE COORD: input indices ("//&
465!      &                 TRIM(fct_str(in_imin0))//","//&
466!      &                 TRIM(fct_str(in_imax0))//") ("//&
467!      &                 TRIM(fct_str(in_jmin0))//","//&
468!      &                 TRIM(fct_str(in_jmax0))//")" )
469!      CALL logger_fatal("CREATE COORD: output domain not confrom "//&
470!      &                 "with input indices")
471!   ENDIF
472!
473!   CALL iom_mpp_close(tl_coord0)
474!   CALL iom_mpp_close(tl_mppout)
475
476   ! close log file
477   CALL logger_footer()
478   CALL logger_close() 
479
480CONTAINS
481   !-------------------------------------------------------------------
482   !> @brief
483   !> This function compute offset over Arakawa grid points,
484   !> given refinement factor.
485   !>
486   !> @author J.Paul
487   !> @date August, 2014 - Initial Version
488   !>
489   !> @param[in] id_rho array of refinement factor
490   !> @return array of offset
491   !-------------------------------------------------------------------
492   FUNCTION create_coord_get_offset( id_rho )
493      IMPLICIT NONE
494      ! Argument     
495      INTEGER(i4), DIMENSION(:), INTENT(IN) :: id_rho
496
497      ! function
498      INTEGER(i4), DIMENSION(2,2,ip_npoint) :: create_coord_get_offset
499      ! local variable
500      ! loop indices
501      !----------------------------------------------------------------
502
503      ! case 'T'
504      create_coord_get_offset(jp_I,:,jp_T)=FLOOR(REAL(id_rho(jp_I)-1,dp)*0.5)
505      create_coord_get_offset(jp_J,:,jp_T)=FLOOR(REAL(id_rho(jp_J)-1,dp)*0.5)
506      ! case 'U'
507      create_coord_get_offset(jp_I,1,jp_U)=0
508      create_coord_get_offset(jp_I,2,jp_U)=id_rho(jp_I)-1
509      create_coord_get_offset(jp_J,:,jp_U)=FLOOR(REAL(id_rho(jp_J)-1,dp)*0.5)
510      ! case 'V'
511      create_coord_get_offset(jp_I,:,jp_V)=FLOOR(REAL(id_rho(jp_I)-1,dp)*0.5)
512      create_coord_get_offset(jp_J,1,jp_V)=0
513      create_coord_get_offset(jp_J,2,jp_V)=id_rho(jp_J)-1
514      ! case 'F'
515      create_coord_get_offset(jp_I,1,jp_F)=0
516      create_coord_get_offset(jp_I,2,jp_F)=id_rho(jp_I)-1
517      create_coord_get_offset(jp_J,1,jp_F)=0
518      create_coord_get_offset(jp_J,2,jp_F)=id_rho(jp_J)-1
519
520
521   END FUNCTION create_coord_get_offset
522   !-------------------------------------------------------------------
523   !> @brief
524   !> This subroutine interpolate variable, given refinment factor.
525   !>
526   !> @details
527   !>  Optionaly, you could specify number of points
528   !>    to be extrapolated in i- and j-direction.<br/>
529   !>  variable mask is first computed (using _FillValue) and interpolated.<br/>
530   !>  variable is then extrapolated, and interpolated.<br/>
531   !>  Finally interpolated mask is applied on refined variable.
532   !>
533   !> @author J.Paul
534   !> @date November, 2013 - Initial Version
535   !>
536   !> @param[inout] td_var variable strcuture
537   !> @param[in] id_rho    array of refinement factor
538   !> @param[in] id_offset offset between fine grid and coarse grid
539   !> @param[in] id_iext   number of points to be extrapolated in i-direction
540   !> @param[in] id_jext   number of points to be extrapolated in j-direction
541   !-------------------------------------------------------------------
542   SUBROUTINE create_coord_interp( td_var,          &
543   &                               id_rho,          &
544   &                               id_offset,       &
545   &                               id_iext, id_jext)
546
547      IMPLICIT NONE
548
549      ! Argument
550      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
551      INTEGER(i4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
552      INTEGER(i4), DIMENSION(:,:), INTENT(IN)    :: id_offset
553      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_iext
554      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_jext
555
556      ! local variable
557      TYPE(TVAR)  :: tl_mask
558
559      INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
560
561      INTEGER(i4) :: il_iext
562      INTEGER(i4) :: il_jext
563
564      ! loop indices
565      !----------------------------------------------------------------
566
567      IF( ANY(SHAPE(id_offset(:,:)) /= 2) )THEN
568         CALL logger_error("CREATE COORD INTERP: invalid dimension of "//&
569         &                 "offset array")
570      ENDIF
571
572      !WARNING: two extrabands are required for cubic interpolation
573      il_iext=2
574      IF( PRESENT(id_iext) ) il_iext=id_iext
575
576      il_jext=2
577      IF( PRESENT(id_jext) ) il_jext=id_jext
578     
579      IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
580         CALL logger_warn("CREATE COORD INTERP: at least extrapolation "//&
581         &  "on two points are required with cubic interpolation ")
582         il_iext=2
583      ENDIF
584
585      IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
586         CALL logger_warn("CREATE COORD INTERP: at least extrapolation "//&
587         &  "on two points are required with cubic interpolation ")
588         il_jext=2
589      ENDIF
590
591      IF( ANY(id_rho(:)>1) )THEN
592         ! work on mask
593         ! create mask
594         ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
595         &                td_var%t_dim(2)%i_len, &
596         &                td_var%t_dim(3)%i_len, &
597         &                td_var%t_dim(4)%i_len) )
598
599         bl_mask(:,:,:,:)=1
600         WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0     
601
602         SELECT CASE(TRIM(td_var%c_point))
603         CASE DEFAULT ! 'T'
604            tl_mask=var_init('tmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
605            &                id_ew=td_var%i_ew )
606         CASE('U')
607            tl_mask=var_init('umask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
608            &                id_ew=td_var%i_ew )
609         CASE('V')
610            tl_mask=var_init('vmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
611            &                id_ew=td_var%i_ew )
612         CASE('F')
613            tl_mask=var_init('fmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
614            &                id_ew=td_var%i_ew )
615         END SELECT         
616
617         DEALLOCATE(bl_mask)
618
619         ! interpolate mask
620         CALL interp_fill_value( tl_mask, id_rho(:), &
621         &                       id_offset=id_offset(:,:) )
622
623         ! work on variable
624         ! add extraband
625         CALL extrap_add_extrabands(td_var, il_iext, il_jext)
626
627         ! extrapolate variable
628         CALL extrap_fill_value( td_var, id_iext=il_iext, id_jext=il_jext )
629
630         ! interpolate variable
631         CALL interp_fill_value( td_var, id_rho(:), &
632         &                       id_offset=id_offset(:,:))
633
634         ! remove extraband
635         CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J))
636
637         ! keep original mask
638         WHERE( tl_mask%d_value(:,:,:,:) == 0 )
639            td_var%d_value(:,:,:,:)=td_var%d_fill
640         END WHERE
641      ENDIF
642
643      ! clean variable structure
644      CALL var_clean(tl_mask)
645
646   END SUBROUTINE create_coord_interp
647END PROGRAM create_coord
Note: See TracBrowser for help on using the repository browser.