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/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/SIREN/src/create_coord.f90 @ 8862

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

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

Download in other formats: