Context Navigation

interp.f90 @ 14712

Last change on this file since 14712 was 13369, checked in by jpaul, 4 years ago
update: cf changelog inside documentation
File size: 34.2 KB

Rev	Line
[4213]	1	!----------------------------------------------------------------------
	2	! NEMO system team, System and Interface for oceanic RElocable Nesting
	3	!----------------------------------------------------------------------
	4	!
	5	! DESCRIPTION:
[13369]	6	!> @brief
[4213]	7	!> This module manage interpolation on regular grid.
	8	!>
[13369]	9	!> @details Interpolation method to be used is specify inside variable
	10	!> strcuture, as array of string character.<br/>
[5037]	11	!> - td_var\%c_interp(1) string character is the interpolation name choose between:
	12	!> - 'nearest'
	13	!> - 'cubic '
	14	!> - 'linear '
[13369]	15	!> - td_var\%c_interp(2) string character is an operation to be used
[5037]	16	!> on interpolated value.<br/>
	17	!> operation have to be mulitplication '*' or division '/'.<br/>
[13369]	18	!> coefficient have to be refinement factor following i-direction 'rhoi',
[5037]	19	!> j-direction 'rhoj', or k-direction 'rhok'.<br/>
	20	!>
	21	!> Examples: '*rhoi', '/rhoj'.
[13369]	22	!>
[5037]	23	!> @note Those informations are read from namelist or variable configuration file (default).<br/>
	24	!> Interplation method could be specify for each variable in namelist _namvar_,
	25	!> defining string character _cn\_varinfo_.<br/>
	26	!> Example:
[13369]	27	!> - cn_varinfo='varname1:int=cubic/rhoi', 'varname2:int=linear'
[5037]	28	!>
	29	!> to create mixed grid (with coarse grid point needed to compute
	30	!> interpolation):<br/>
	31	!> @code
	32	!> CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] )
	33	!> @endcode
	34	!> - td_var is coarse grid variable (should be extrapolated)
	35	!> - td_mix is mixed grid variable structure [output]
	36	!> - id_rho is array of refinment factor [optional]
	37	!>
	38	!> to detected point to be interpolated:<br/>
	39	!> @code
	40	!> il_detect(:,:,:)=interp_detect( td_mix [,id_rho] )
	41	!> @endcode
	42	!> - il_detect(:,:,:) is 3D array of detected point to be interpolated
	43	!> - td_mix is mixed grid variable
	44	!> - id_rho is array of refinement factor [optional]
	45	!>
	46	!> to interpolate variable value:<br/>
	47	!> @code
	48	!> CALL interp_fill_value( td_var [,id_rho] [,id_offset] )
	49	!> @endcode
	50	!> - td_var is variable structure
	51	!> - id_rho is array of refinement factor [optional]
	52	!> - id_offset is array of offset between fine and coarse grid [optional]
	53	!>
	54	!> to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/>
	55	!> @code
	56	!> CALL interp_clean_mixed_grid( td_mix, td_var, id_rho )
	57	!> @endcode
	58	!> - td_mix is mixed grid variable structure
	59	!> - td_var is variable structure [output]
	60	!> - id_rho is array of refinement factor [optional]
	61	!> - id_offset is array of offset between fine and coarse grid [optional]
	62	!>
	63	!> @note It use to work on ORCA grid, as we work only with grid indices.
	64	!>
[4213]	65	!> @warning due to the use of second derivative when using cubic interpolation
[5037]	66	!> you should add at least 2 extrabands.
[4213]	67	!>
	68	!> @author
	69	!> J.Paul
[13369]	70	!>
[5037]	71	!> @date November, 2013 - Initial Version
	72	!> @date September, 2014
	73	!> - add header
	74	!> - use interpolation method modules
	75	!>
[12080]	76	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[4213]	77	!----------------------------------------------------------------------
	78	MODULE interp
	79
	80	USE netcdf ! nf90 library
	81	USE global ! global variable
	82	USE kind ! F90 kind parameter
[5037]	83	USE logger ! log file manager
[4213]	84	USE fct ! basic useful function
[5037]	85	USE date ! date manager
[4213]	86	USE att ! attribute manager
	87	USE dim ! dimension manager
	88	USE var ! variable manager
	89	USE grid ! grid manager
	90	USE extrap ! extrapolation manager
[5037]	91	USE interp_cubic ! cubic interpolation manager
	92	USE interp_linear ! linear interpolation manager
	93	USE interp_nearest ! nearest interpolation manager
[4213]	94
	95	IMPLICIT NONE
	96	! NOTE_avoid_public_variables_if_possible
	97
	98	! type and variable
	99
	100	! function and subroutine
	101	PUBLIC :: interp_detect !< detected point to be interpolated
[5037]	102	PUBLIC :: interp_fill_value !< interpolate value
[4213]	103	PUBLIC :: interp_create_mixed_grid !< create mixed grid
	104	PUBLIC :: interp_clean_mixed_grid !< clean mixed grid
	105
[5037]	106	PRIVATE :: interp__detect ! detected point to be interpolated
	107	PRIVATE :: interp__detect_wrapper ! detected point to be interpolated
	108	PRIVATE :: interp__fill_value_wrapper ! interpolate value over detectected point
	109	PRIVATE :: interp__fill_value ! interpolate value over detectected point
	110	PRIVATE :: interp__clean_even_grid ! clean even mixed grid
	111	PRIVATE :: interp__check_method ! check if interpolation method available
[13369]	112
[4213]	113	TYPE TINTERP
	114	CHARACTER(LEN=lc) :: c_name = '' !< interpolation method name
	115	CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor
	116	CHARACTER(LEN=lc) :: c_divisor= '' !< interpolation divisor
	117	END TYPE TINTERP
	118
	119	INTERFACE interp_detect
[13369]	120	MODULE PROCEDURE interp__detect_wrapper
[4213]	121	END INTERFACE interp_detect
	122
	123	INTERFACE interp_fill_value
[13369]	124	MODULE PROCEDURE interp__fill_value_wrapper
[4213]	125	END INTERFACE interp_fill_value
	126
	127	CONTAINS
[12080]	128	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	129	FUNCTION interp__check_method(cd_method) &
	130	& RESULT (lf_avail)
[4213]	131	!-------------------------------------------------------------------
	132	!> @brief
[5037]	133	!> This function check if interpolation method is available.
[13369]	134	!>
	135	!> @details
[5037]	136	!> check if name of interpolation method is present in global list of string
	137	!> character cp_interp_list (see global.f90).
[4213]	138	!>
	139	!> @author J.Paul
[5617]	140	!> @date November, 2013 - Initial Version
[4213]	141	!
[5037]	142	!> @param[in] cd_method interpolation method
[12080]	143	!> @return true if interpolation method is available
[4213]	144	!-------------------------------------------------------------------
[12080]	145
[4213]	146	IMPLICIT NONE
[12080]	147
[4213]	148	! Argument
	149	CHARACTER(LEN=lc) :: cd_method
	150
	151	! function
[12080]	152	LOGICAL :: lf_avail
[4213]	153
	154	! local variable
	155	CHARACTER(LEN=lc) :: cl_interp
	156	CHARACTER(LEN=lc) :: cl_method
	157
	158	! loop indices
	159	INTEGER(I4) :: ji
	160	!----------------------------------------------------------------
	161
	162	cl_method=fct_lower(cd_method)
	163
[12080]	164	lf_avail=.FALSE.
[5037]	165	DO ji=1,ip_ninterp
	166	cl_interp=fct_lower(cp_interp_list(ji))
[4213]	167	IF( TRIM(cl_interp) == TRIM(cl_method) )THEN
[12080]	168	lf_avail=.TRUE.
[4213]	169	EXIT
	170	ENDIF
	171	ENDDO
	172
	173	END FUNCTION interp__check_method
[12080]	174	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	175	FUNCTION interp__detect_wrapper(td_mix, id_rho) &
	176	& RESULT (if_detect)
[4213]	177	!-------------------------------------------------------------------
	178	!> @brief
	179	!> This function detected point to be interpolated.
[13369]	180	!>
	181	!> @details
[4213]	182	!> Actually it checks, the number of dimension used for this variable
	183	!> and launch interp__detect which detected point to be interpolated.
	184	!>
	185	!> @author J.Paul
[5617]	186	!> @date November, 2013 - Initial Version
[5037]	187	!>
	188	!> @param[in] td_mix mixed grid variable (to interpolate)
	189	!> @param[in] id_rho array of refinement factor
[13369]	190	!> @return 3D array of detected point to be interpolated
[4213]	191	!-------------------------------------------------------------------
[12080]	192
[4213]	193	IMPLICIT NONE
[12080]	194
[4213]	195	! Argument
	196	TYPE(TVAR) , INTENT(IN) :: td_mix
	197	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
	198
	199	! function
	200	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
[12080]	201	& td_mix%t_dim(2)%i_len,&
	202	& td_mix%t_dim(3)%i_len ) :: if_detect
[4213]	203
	204	! local variable
	205	! loop indices
	206	!----------------------------------------------------------------
	207
	208	IF( .NOT. ANY(td_mix%t_dim(1:3)%l_use) )THEN
	209	! no dimension I-J-K used
	210	CALL logger_debug(" INTERP DETECT: nothing done for variable"//&
[12080]	211	& TRIM(td_mix%c_name) )
[13369]	212
[12080]	213	if_detect(:,:,:)=0
[4213]	214
	215	ELSE IF( ALL(td_mix%t_dim(1:3)%l_use) )THEN
[13369]	216
[4213]	217	! detect point to be interpolated on I-J-K
	218	CALL logger_debug(" INTERP DETECT: detect point "//&
[12080]	219	& TRIM(td_mix%c_point)//" for variable "//&
	220	& TRIM(td_mix%c_name) )
[13369]	221
[12080]	222	if_detect(:,:,:)=interp__detect( td_mix, id_rho(:) )
[4213]	223
	224	ELSE IF( ALL(td_mix%t_dim(1:2)%l_use) )THEN
	225
	226	! detect point to be interpolated on I-J
	227	CALL logger_debug(" INTERP DETECT: detect point "//&
[12080]	228	& TRIM(td_mix%c_point)//" for variable "//&
	229	& TRIM(td_mix%c_name) )
[13369]	230
[12080]	231	if_detect(:,:,1:1)=interp__detect( td_mix, id_rho(:))
[4213]	232
	233	ELSE IF( td_mix%t_dim(3)%l_use )THEN
[13369]	234
[4213]	235	! detect point to be interpolated on K
	236	CALL logger_debug(" INTERP DETECT: detect vertical point "//&
[12080]	237	& " for variable "//TRIM(td_mix%c_name) )
[13369]	238
[12080]	239	if_detect(1:1,1:1,:)=interp__detect( td_mix, id_rho(:) )
[4213]	240
[13369]	241	ENDIF
[4213]	242
	243	END FUNCTION interp__detect_wrapper
[12080]	244	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	245	FUNCTION interp__detect(td_mix, id_rho) &
	246	& RESULT (if_detect)
[4213]	247	!-------------------------------------------------------------------
	248	!> @brief
	249	!> This function detected point to be interpolated.
[13369]	250	!>
	251	!> @details
[4213]	252	!> A special case is done for even refinement on ARAKAWA-C grid.
	253	!>
	254	!> @author J.Paul
[5617]	255	!> @date November, 2013 - Initial Version
[4213]	256	!
[5037]	257	!> @param[in] td_mix mixed grid variable (to interpolate)
	258	!> @param[in] id_rho array of refinement factor
[13369]	259	!> @return 3D array of detected point to be interpolated
[4213]	260	!-------------------------------------------------------------------
[12080]	261
[4213]	262	IMPLICIT NONE
[12080]	263
[4213]	264	! Argument
	265	TYPE(TVAR) , INTENT(IN) :: td_mix
	266	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
	267
	268	! function
	269	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
[12080]	270	& td_mix%t_dim(2)%i_len,&
	271	& td_mix%t_dim(3)%i_len ) :: if_detect
[4213]	272
	273	! local variable
	274	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
	275
	276	INTEGER(I4) :: il_xextra
	277	INTEGER(I4) :: il_yextra
	278	INTEGER(I4) :: il_zextra
	279
	280	INTEGER(i4), DIMENSION(3) :: il_dim
	281
	282	LOGICAL , DIMENSION(3) :: ll_even
	283
	284	! loop indices
	285	INTEGER(I4) :: ji
	286	INTEGER(I4) :: jj
	287	INTEGER(I4) :: jk
	288	!----------------------------------------------------------------
	289	ALLOCATE( il_rho(ip_maxdim) )
	290	il_rho(:)=1
[5037]	291	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
[4213]	292
	293	! special case for even refinement on ARAKAWA-C grid
	294	ll_even(:)=.FALSE.
	295	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
	296	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
[13369]	297	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
[4213]	298
	299	SELECT CASE(TRIM(td_mix%c_point))
	300	CASE('U')
	301	ll_even(1)=.FALSE.
	302	CASE('V')
	303	ll_even(2)=.FALSE.
	304	CASE('F')
	305	ll_even(:)=.FALSE.
	306	END SELECT
	307
	308	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
	309	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
	310	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
	311
	312	! special case for cubic interpolation
	313	il_xextra=0
	314	il_yextra=0
	315	il_zextra=0
	316	SELECT CASE(TRIM(td_mix%c_interp(1)))
	317	CASE('cubic')
	318	! those points can not be compute cause cubic interpolation
	319	! need second derivative.
	320	IF( il_rho(jp_I) /= 1 ) il_xextra=3*il_rho(jp_I)
	321	IF( il_rho(jp_J) /= 1 ) il_yextra=3*il_rho(jp_J)
	322	IF( il_rho(jp_K) /= 1 ) il_zextra=3*il_rho(jp_K)
	323	END SELECT
	324
	325	il_dim(:)=td_mix%t_dim(1:3)%i_len
	326
[5037]	327	! init
[12080]	328	if_detect(:,:,:)=1
[4213]	329
	330	! do not compute coarse grid point
[12080]	331	if_detect(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
	332	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
	333	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K) ) = 0
[4213]	334
	335	! do not compute point near fill value
[5037]	336	DO jk=1,il_dim(3),il_rho(jp_K)
	337	DO jj=1,il_dim(2),il_rho(jp_J)
	338	DO ji=1,il_dim(1),il_rho(jp_I)
[4213]	339
[5037]	340	IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN
[4213]	341
[5037]	342	! i-direction
[12080]	343	if_detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
	344	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
[13369]	345	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
[5037]	346	! j-direction
[12080]	347	if_detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
	348	& MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
	349	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
[5037]	350	! k-direction
[12080]	351	if_detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
	352	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
[13369]	353	& MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0
[4213]	354
[5037]	355	ENDIF
	356
	357	ENDDO
	358	ENDDO
	359	ENDDO
	360
	361	DEALLOCATE( il_rho )
	362
[4213]	363	END FUNCTION interp__detect
[12080]	364	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	365	SUBROUTINE interp_create_mixed_grid(td_var, td_mix, id_rho)
[4213]	366	!-------------------------------------------------------------------
	367	!> @brief
	368	!> This subroutine create mixed grid.
[13369]	369	!>
	370	!> @details
[4213]	371	!> Created grid is fine resolution grid.
[13369]	372	!> First and last point are coasre grid point.
[4213]	373	!>
	374	!> A special case is done for even refinement on ARAKAWA-C grid.
	375	!>
	376	!> @author J.Paul
[5617]	377	!> @date November, 2013 - Initial Version
[4213]	378	!>
[5037]	379	!> @param[in] td_var coarse grid variable (should be extrapolated)
	380	!> @param[out] td_mix mixed grid variable
[13369]	381	!> @param[in] id_rho array of refinment factor (default 1)
[4213]	382	!-------------------------------------------------------------------
[12080]	383
[4213]	384	IMPLICIT NONE
[12080]	385
[4213]	386	! Argument
[13369]	387	TYPE(TVAR) , INTENT(IN ) :: td_var
	388	TYPE(TVAR) , INTENT( OUT) :: td_mix
[4213]	389	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
	390
	391	! local variable
	392	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
	393
	394	INTEGER(i4) :: il_xextra
	395	INTEGER(i4) :: il_yextra
	396	INTEGER(i4) :: il_zextra
	397
	398	LOGICAL, DIMENSION(3) :: ll_even
	399
	400	! loop indices
	401	!----------------------------------------------------------------
	402	ALLOCATE(il_rho(ip_maxdim))
	403	il_rho(:)=1
[5037]	404	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
[4213]	405
	406	! special case for even refinement on ARAKAWA-C grid
	407	ll_even(:)=.FALSE.
	408	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
	409	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
[13369]	410	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
[4213]	411
	412	SELECT CASE(TRIM(td_var%c_point))
	413	CASE('U')
	414	ll_even(1)=.FALSE.
	415	CASE('V')
	416	ll_even(2)=.FALSE.
	417	CASE('F')
	418	ll_even(:)=.FALSE.
	419	END SELECT
	420
	421	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
	422	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
	423	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
	424
	425	! copy variable
[5037]	426	td_mix=var_copy(td_var)
[4213]	427
	428	! compute new dimension length
	429	il_xextra=il_rho(jp_I)-1
	430	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
	431
	432	il_yextra=il_rho(jp_J)-1
	433	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
	434
	435	il_zextra=il_rho(jp_K)-1
	436	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
	437
	438	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
	439	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
	440	& td_mix%t_dim(2)%i_len, &
	441	& td_mix%t_dim(3)%i_len, &
	442	& td_mix%t_dim(4)%i_len) )
	443
	444	! initialise to FillValue
	445	td_mix%d_value(:,:,:,:)=td_mix%d_fill
	446
	447	! quid qd coord ou bathy fourni par user ?? (offset ??)
	448	td_mix%d_value(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
	449	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
	450	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K), :) = &
	451	& td_var%d_value(:,:,:,:)
	452
[5037]	453	DEALLOCATE(il_rho)
	454
[4213]	455	END SUBROUTINE interp_create_mixed_grid
[12080]	456	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	457	SUBROUTINE interp__clean_even_grid(td_mix, id_rho)
[4213]	458	!-------------------------------------------------------------------
	459	!> @brief
	460	!> This subroutine remove points added to mixed grid to compute
	461	!> interpolation in the special case of even refinement on ARAKAWA-C grid.
	462	!>
[13369]	463	!> @details
	464	!>
[4213]	465	!> @author J.Paul
[5617]	466	!> @date November, 2013 - Initial Version
[4213]	467	!>
[5037]	468	!> @param[inout] td_mix mixed grid variable
[13369]	469	!> @param[in] id_rho array of refinment factor
[4213]	470	!-------------------------------------------------------------------
[12080]	471
[4213]	472	IMPLICIT NONE
[12080]	473
[4213]	474	! Argument
[13369]	475	TYPE(TVAR) , INTENT(INOUT) :: td_mix
[4213]	476	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
	477
	478	! local variable
	479	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
	480
	481	INTEGER(i4) :: il_xextra
	482	INTEGER(i4) :: il_yextra
	483	INTEGER(i4) :: il_zextra
	484
	485	LOGICAL, DIMENSION(3) :: ll_even
	486
	487	LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ll_mask
	488
	489	REAL(dp), DIMENSION(:) , ALLOCATABLE :: dl_vect
	490
	491	TYPE(TVAR) :: tl_mix
	492
	493	! loop indices
	494	!----------------------------------------------------------------
	495	ALLOCATE(il_rho(ip_maxdim))
	496	il_rho(:)=1
	497	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
	498
	499	! special case for even refinement on ARAKAWA-C grid
	500	ll_even(:)=.FALSE.
	501	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
	502	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
[13369]	503	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
[4213]	504
	505	SELECT CASE(TRIM(td_mix%c_point))
	506	CASE('U')
	507	ll_even(1)=.FALSE.
	508	CASE('V')
	509	ll_even(2)=.FALSE.
	510	CASE('F')
	511	ll_even(:)=.FALSE.
	512	END SELECT
	513
	514	! remove some point only if refinement in some direction is even
	515	IF( ANY(ll_even(:)) )THEN
	516
[5037]	517	! copy variable
	518	tl_mix=var_copy(td_mix)
	519
[4213]	520	ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, &
	521	& tl_mix%t_dim(2)%i_len, &
	522	& tl_mix%t_dim(3)%i_len, &
	523	& tl_mix%t_dim(4)%i_len) )
	524
	525	ll_mask(:,:,:,:)=.TRUE.
	526
	527	IF( tl_mix%t_dim(1)%l_use .AND. ll_even(1) )THEN
	528
	529	il_rho(jp_I)=il_rho(jp_I)+1
	530
	531	! locate wrong point on mixed grid
	532	ll_mask(1:td_mix%t_dim(1)%i_len:il_rho(jp_I),:,:,:)=.FALSE.
	533
	534	! compute coasre grid dimension length
	535	il_xextra=il_rho(jp_I)-1
	536	td_mix%t_dim(1)%i_len=(tl_mix%t_dim(1)%i_len+il_xextra)/il_rho(jp_I)
	537
	538	il_rho(jp_I)=il_rho(jp_I)-1
	539	! compute right fine grid dimension length
	540	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
	541
	542	ENDIF
	543
	544	IF( tl_mix%t_dim(2)%l_use .AND. ll_even(2) )THEN
	545
	546	il_rho(jp_J)=il_rho(jp_J)+1
	547
	548	! locate wrong point on mixed grid
	549	ll_mask(:,1:tl_mix%t_dim(2)%i_len:il_rho(jp_J),:,:)=.FALSE.
	550
	551	! compute coasre grid dimension length
	552	il_yextra=il_rho(jp_J)-1
	553	td_mix%t_dim(2)%i_len=(tl_mix%t_dim(2)%i_len+il_yextra)/il_rho(jp_J)
	554
	555	il_rho(jp_J)=il_rho(jp_J)-1
	556	! compute right fine grid dimension length
	557	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
[13369]	558
[4213]	559	ENDIF
	560
	561	IF( tl_mix%t_dim(3)%l_use .AND. ll_even(3) )THEN
	562
	563	il_rho(jp_K)=il_rho(jp_K)+1
	564
	565	! locate wrong point on mixed grid
	566	ll_mask(:,:,1:tl_mix%t_dim(3)%i_len:il_rho(jp_K),:)=.FALSE.
	567
	568	! compute coasre grid dimension length
	569	il_zextra=il_rho(jp_K)-1
	570	td_mix%t_dim(3)%i_len=(tl_mix%t_dim(3)%i_len+il_zextra)/il_rho(jp_K)
	571
	572	il_rho(jp_K)=il_rho(jp_K)-1
	573	! compute right fine grid dimension length
	574	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
	575
[13369]	576	ENDIF
[4213]	577
	578	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
	579	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
	580	& td_mix%t_dim(2)%i_len, &
	581	& td_mix%t_dim(3)%i_len, &
	582	& td_mix%t_dim(4)%i_len) )
	583
	584	! initialise to FillValue
	585	td_mix%d_value(:,:,:,:)=td_mix%d_fill
	586
	587	IF( COUNT(ll_mask(:,:,:,:)) /= SIZE(td_mix%d_value(:,:,:,:)) )THEN
	588
	589	CALL logger_error("INTERP CLEAN EVEN GRID: output value size "//&
	590	& " and mask count differ ")
	591
	592	ELSE
	593
	594	ALLOCATE( dl_vect(COUNT(ll_mask(:,:,:,:))) )
	595
	596	dl_vect(:)= PACK( tl_mix%d_value(:,:,:,:), &
	597	& MASK=ll_mask(:,:,:,:) )
	598
	599	td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), &
[5037]	600	& SHAPE=(/td_mix%t_dim(1)%i_len, &
	601	& td_mix%t_dim(2)%i_len, &
	602	& td_mix%t_dim(3)%i_len, &
	603	& td_mix%t_dim(4)%i_len/) )
[4213]	604
	605	DEALLOCATE( dl_vect )
	606
	607	ENDIF
	608
	609	DEALLOCATE( ll_mask )
	610
[5037]	611	! clean
	612	CALL var_clean(tl_mix)
	613
[4213]	614	ENDIF
	615
[5037]	616	! clean
	617	DEALLOCATE(il_rho)
[4213]	618
	619	END SUBROUTINE interp__clean_even_grid
[12080]	620	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	621	SUBROUTINE interp_clean_mixed_grid(td_mix, td_var, &
	622	& id_rho, id_offset)
[4213]	623	!-------------------------------------------------------------------
	624	!> @brief
[13369]	625	!> This subroutine remove points added on mixed grid
	626	!> to compute interpolation. And save interpolated value over domain.
[4213]	627	!>
[13369]	628	!> @details
	629	!>
[4213]	630	!> @author J.Paul
[5617]	631	!> @date November, 2013 - Initial Version
[5037]	632	!> @date September, 2014
	633	!> - use offset to save useful domain
	634	!>
	635	!> @param[in] td_mix mixed grid variable structure
	636	!> @param[out] td_var variable structure
	637	!> @param[in] id_rho array of refinement factor (default 1)
	638	!> @param[in] id_offset 2D array of offset between fine and coarse grid
[4213]	639	!-------------------------------------------------------------------
[12080]	640
[4213]	641	IMPLICIT NONE
[12080]	642
[4213]	643	! Argument
	644	TYPE(TVAR) , INTENT(IN ) :: td_mix
	645	TYPE(TVAR) , INTENT( OUT) :: td_var
[5037]	646	INTEGER(I4), DIMENSION(:) , INTENT(IN ) :: id_rho
	647	INTEGER(I4), DIMENSION(2,2), INTENT(IN ) :: id_offset
[4213]	648
	649	! local variable
	650	INTEGER(i4) :: il_imin0
	651	INTEGER(i4) :: il_imax0
	652	INTEGER(i4) :: il_jmin0
	653	INTEGER(i4) :: il_jmax0
	654
	655	INTEGER(i4) :: il_imin1
	656	INTEGER(i4) :: il_jmin1
	657	INTEGER(i4) :: il_imax1
	658	INTEGER(i4) :: il_jmax1
	659
	660	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
	661
	662	TYPE(TVAR) :: tl_mix
	663
	664	! loop indices
	665	!----------------------------------------------------------------
	666	! copy mixed variable in temporary structure
[5037]	667	tl_mix=var_copy(td_mix)
[4213]	668
[12080]	669	! remove useless points over mixed grid for even refinement
[4213]	670	CALL interp__clean_even_grid(tl_mix, id_rho(:))
	671
	672	! copy cleaned mixed variable
[5037]	673	td_var=var_copy(tl_mix)
[4213]	674
[5037]	675	! delete array of value
	676	CALL var_del_value(td_var)
[4213]	677
[5037]	678	! compute domain indices in i-direction
	679	il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len
[13369]	680
[5037]	681	IF( td_var%t_dim(1)%l_use )THEN
	682	il_imin1=il_imin0+id_offset(Jp_I,1)
	683	il_imax1=il_imax0-id_offset(Jp_I,2)
	684	ELSE
[4213]	685
[5037]	686	il_imin1=il_imin0
	687	il_imax1=il_imax0
[4213]	688
[5037]	689	ENDIF
[4213]	690
[5037]	691	! compute domain indices in j-direction
	692	il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len
[13369]	693
[5037]	694	IF( td_var%t_dim(2)%l_use )THEN
	695	il_jmin1=il_jmin0+id_offset(Jp_J,1)
	696	il_jmax1=il_jmax0-id_offset(Jp_J,2)
	697	ELSE
[4213]	698
[5037]	699	il_jmin1=il_jmin0
	700	il_jmax1=il_jmax0
[13369]	701
[5037]	702	ENDIF
[4213]	703
[5037]	704	! compute new dimension
	705	td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
	706	td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
[13369]	707
[5037]	708	ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
	709	& td_var%t_dim(2)%i_len, &
	710	& td_var%t_dim(3)%i_len, &
	711	& td_var%t_dim(4)%i_len) )
[4213]	712
[5037]	713	dl_value( 1:td_var%t_dim(1)%i_len, &
	714	& 1:td_var%t_dim(2)%i_len, &
	715	& :,:) = tl_mix%d_value( il_imin1:il_imax1, &
	716	& il_jmin1:il_jmax1, &
	717	& :, : )
[4213]	718
[5037]	719	! add variable value
	720	CALL var_add_value(td_var,dl_value(:,:,:,:))
[4213]	721
[5037]	722	DEALLOCATE(dl_value)
[13369]	723
[5037]	724	! clean
[4213]	725	CALL var_clean(tl_mix)
	726
	727	END SUBROUTINE interp_clean_mixed_grid
[12080]	728	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[13369]	729	SUBROUTINE interp__fill_value_wrapper(td_var, &
[12080]	730	& id_rho, &
	731	& id_offset)
[4213]	732	!-------------------------------------------------------------------
	733	!> @brief
[5037]	734	!> This subroutine interpolate variable value.
[13369]	735	!>
	736	!> @details
[4213]	737	!> Actually it checks, the number of dimension used for this variable
	738	!> and launch interp__fill_value.
	739	!>
	740	!> @author J.Paul
[5617]	741	!> @date November, 2013 - Initial Version
[5037]	742	!>
	743	!> @param[inout] td_var variable structure
	744	!> @param[in] id_rho array of refinement factor
	745	!> @param[in] id_offset 2D array of offset between fine and coarse grid
[4213]	746	!-------------------------------------------------------------------
[12080]	747
[4213]	748	IMPLICIT NONE
[12080]	749
[4213]	750	! Argument
	751	TYPE(TVAR) , INTENT(INOUT) :: td_var
[5037]	752	INTEGER(I4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
[4213]	753	INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
	754
	755	! local variable
[5037]	756	CHARACTER(LEN=lc) :: cl_method
	757	INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho
	758	INTEGER(i4) , DIMENSION(2,2) :: il_offset
[4213]	759
	760	! loop indices
	761	!----------------------------------------------------------------
	762
	763	ALLOCATE( il_rho(ip_maxdim) )
	764	il_rho(:)=1
[5037]	765	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
[4213]	766	IF( ANY(il_rho(:) < 0) )THEN
	767	CALL logger_error("INTERP FILL VALUE: invalid "//&
	768	& " refinement factor ")
	769	ENDIF
	770
	771	il_offset(:,:)=0
	772	IF( PRESENT(id_offset) )THEN
	773	IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN
[5037]	774	CALL logger_error("INTERP FILL VALUE: invalid array of offset")
[4213]	775	ELSE
	776	il_offset(:,:)=id_offset(:,:)
	777	ENDIF
	778	ENDIF
	779
	780	IF( (il_rho(jp_I) /= 1 .AND. td_var%t_dim(1)%l_use) .OR. &
	781	& (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. &
	782	& (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN
	783
[5037]	784	SELECT CASE(TRIM(td_var%c_interp(1)))
	785	CASE('cubic','linear','nearest')
	786	cl_method=TRIM(td_var%c_interp(1))
	787	CASE DEFAULT
	788	CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//&
	789	& " use linear interpolation")
	790	cl_method='linear'
	791	! update variable structure value
	792	td_var%c_interp(1)='linear'
	793	END SELECT
	794
[4213]	795	CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//&
	796	& " using "//TRIM(cl_method)//" method." )
	797	CALL logger_info("INTERP FILL: refinement factor "//&
	798	& TRIM(fct_str(il_rho(jp_I)))//&
	799	& " "//TRIM(fct_str(il_rho(jp_J)))//&
	800	& " "//TRIM(fct_str(il_rho(jp_K))) )
[13369]	801
	802	CALL interp__fill_value( td_var, cl_method, &
[5037]	803	& il_rho(:), il_offset(:,:) )
[4213]	804
	805	SELECT CASE(TRIM(td_var%c_interp(2)))
	806	CASE('/rhoi')
	807	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
	808	td_var%d_value(:,:,:,:) = &
	809	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_I),dp)
	810	END WHERE
	811	CASE('/rhoj')
	812	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
	813	td_var%d_value(:,:,:,:) = &
	814	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_J),dp)
	815	END WHERE
	816	CASE('/rhok')
	817	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
	818	td_var%d_value(:,:,:,:) = &
	819	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_K),dp)
	820	END WHERE
	821	CASE('*rhoi')
	822	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
	823	td_var%d_value(:,:,:,:) = &
	824	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_I),dp)
	825	END WHERE
	826	CASE('*rhoj')
[13369]	827	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
[4213]	828	td_var%d_value(:,:,:,:) = &
	829	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_J),dp)
	830	END WHERE
	831	CASE('*rhok')
[13369]	832	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
[4213]	833	td_var%d_value(:,:,:,:) = &
	834	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_K),dp)
	835	END WHERE
	836	CASE DEFAULT
	837	td_var%c_interp(2)=''
	838	END SELECT
	839
[5037]	840	ELSE
	841	td_var%c_interp(:)=''
[4213]	842	ENDIF
	843
[5037]	844	DEALLOCATE(il_rho)
	845
[4213]	846	END SUBROUTINE interp__fill_value_wrapper
[12080]	847	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	848	SUBROUTINE interp__fill_value(td_var, cd_method, &
	849	& id_rho, id_offset)
[4213]	850	!-------------------------------------------------------------------
	851	!> @brief
	852	!> This subroutine interpolate value over mixed grid.
[13369]	853	!>
[5037]	854	!> @author J.Paul
[5617]	855	!> @date November, 2013 - Initial Version
[5037]	856	!> @date September, 2014
	857	!> - use interpolation method modules
[4213]	858	!>
[5037]	859	!> @param[inout] td_var variable structure
[13369]	860	!> @param[in] cd_method interpolation method
	861	!> @param[in] id_rho array of refinment factor
[5037]	862	!> @param[in] id_offset 2D array of offset between fine and coarse grid
[4213]	863	!-------------------------------------------------------------------
[12080]	864
[4213]	865	IMPLICIT NONE
[12080]	866
[4213]	867	! Argument
[13369]	868	TYPE(TVAR) , INTENT(INOUT) :: td_var
[4213]	869	CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
	870	INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho
[5037]	871	INTEGER(I4) , DIMENSION(2,2), INTENT(IN ) :: id_offset
[4213]	872
	873	! local variable
	874	CHARACTER(LEN=lc) :: cl_interp
	875
	876	INTEGER(I4) , DIMENSION(:) , ALLOCATABLE :: il_rho
	877	INTEGER(i4) , DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
[13369]	878
[4213]	879	REAL(dp) :: dl_min
	880	REAL(dp) :: dl_max
	881
	882	LOGICAL , DIMENSION(3) :: ll_even
	883	LOGICAL :: ll_discont
[13369]	884
[4213]	885	TYPE(TVAR) :: tl_mix
	886
	887	TYPE(TATT) :: tl_att
[13369]	888
[4213]	889	! loop indices
	890	!----------------------------------------------------------------
	891
	892	!1- create mixed grid
	893	CALL interp_create_mixed_grid( td_var, tl_mix, id_rho(:) )
	894
	895	! clean variable structure
	896	CALL var_clean(td_var)
	897
	898	!2- detect point to be interpolated
	899	ALLOCATE( il_detect( tl_mix%t_dim(1)%i_len, &
	900	& tl_mix%t_dim(2)%i_len, &
	901	& tl_mix%t_dim(3)%i_len) )
	902
	903	il_detect(:,:,:)=0
	904
	905	il_detect(:,:,:)=interp_detect(tl_mix, id_rho(:) )
	906
	907	! add attribute to variable
	908	cl_interp=fct_concat(tl_mix%c_interp(:))
	909	tl_att=att_init('interpolation',cl_interp)
	910	CALL var_move_att(tl_mix, tl_att)
	911
[5037]	912	! clean
	913	CALL att_clean(tl_att)
	914
[4213]	915	! special case for even refinement on ARAKAWA-C grid
	916	ll_even(:)=.FALSE.
	917	IF( MOD(id_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
	918	IF( MOD(id_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
[13369]	919	IF( MOD(id_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
[4213]	920
	921	SELECT CASE(TRIM(tl_mix%c_point))
	922	CASE('U')
	923	ll_even(1)=.FALSE.
	924	CASE('V')
	925	ll_even(2)=.FALSE.
	926	CASE('F')
	927	ll_even(:)=.FALSE.
	928	END SELECT
	929
	930	ALLOCATE(il_rho(ip_maxdim))
	931	il_rho(:)=id_rho(:)
	932
	933	IF( ll_even(1) ) il_rho(jp_I)=id_rho(jp_I)+1
	934	IF( ll_even(2) ) il_rho(jp_J)=id_rho(jp_J)+1
	935	IF( ll_even(3) ) il_rho(jp_K)=id_rho(jp_K)+1
	936
	937	! special case for longitude
	938	ll_discont=.FALSE.
	939	IF( TRIM(tl_mix%c_units) == 'degrees_east' )THEN
	940	dl_min=MINVAL( tl_mix%d_value(:,:,:,:), &
	941	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
	942	dl_max=MAXVAL( tl_mix%d_value(:,:,:,:), &
[13369]	943	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
[4213]	944	IF( dl_min < -170_dp .AND. dl_max > 170_dp .OR. &
	945	& dl_min < 10_dp .AND. dl_max > 350_dp )THEN
	946	ll_discont=.TRUE.
	947	ENDIF
	948	ENDIF
	949
	950	!3- interpolate
[5037]	951	CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
	952	& " discont "//TRIM(fct_str(ll_discont)) )
	953	SELECT CASE(TRIM(cd_method))
	954	CASE('cubic')
	955	CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
	956	& il_detect(:,:,:), &
	957	& il_rho(:), ll_even(:), ll_discont )
	958	CASE('nearest')
	959	CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), &
	960	& il_detect(:,:,:), &
	961	& il_rho(:) )
	962	CASE DEFAULT ! linear
	963	CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
	964	& il_detect(:,:,:), &
	965	& il_rho(:), ll_even(:), ll_discont )
[13369]	966	END SELECT
[4213]	967
	968	IF( ANY(il_detect(:,:,:)==1) )THEN
	969	CALL logger_warn("INTERP FILL: some points can not be interpolated "//&
	970	& "for variable "//TRIM(tl_mix%c_name) )
	971	ENDIF
	972
	973	DEALLOCATE(il_detect)
[5609]	974
[4213]	975	!4- save useful domain (remove offset)
	976	CALL interp_clean_mixed_grid( tl_mix, td_var, &
[13369]	977	& id_rho(:), id_offset(:,:) )
[4213]	978
	979	! clean variable structure
[5037]	980	DEALLOCATE(il_rho)
[4213]	981	CALL var_clean(tl_mix)
	982
	983	END SUBROUTINE interp__fill_value
[12080]	984	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[4213]	985	END MODULE interp

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

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

Context Navigation

source: utils/tools/SIREN/src/interp.f90 @ 14712

Download in other formats: