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/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_coord.f90 @ 7025

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: