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 @ 7153

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

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

Download in other formats: