Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

interp_linear.f90

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

Rev	Line
[5037]	1	!----------------------------------------------------------------------
	2	! NEMO system team, System and Interface for oceanic RElocable Nesting
	3	!----------------------------------------------------------------------
	4	!
	5	! DESCRIPTION:
[13369]	6	!> @brief
[5037]	7	!> This module manage linear interpolation on regular grid.
	8	!>
	9	!> @details
	10	!> to compute linear interpolation:<br/>
	11	!> @code
	12	!> CALL interp_linear_fill(dd_value, dd_fill, id_detect, id_rho, ld_even [,ld_discont] )
	13	!> @endcode
[13369]	14	!> - dd_value is 2D array of variable value
[5037]	15	!> - dd_fill is the FillValue of variable
[13369]	16	!> - id_detect is 2D array of point to be interpolated (see interp module)
[5037]	17	!> - id_rho is array of refinment factor
[13369]	18	!> - ld_even indicates even refinment or not
[5037]	19	!> - ld_discont indicates longitudinal discontinuity (-180°/180°, 0°/360°) or not
	20	!>
	21	!> @author
	22	!> J.Paul
[12080]	23	!>
[5617]	24	!> @date September, 2014 - Initial version
[5037]	25	!>
[12080]	26	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[5037]	27	!----------------------------------------------------------------------
	28	MODULE interp_linear
	29
	30	USE netcdf ! nf90 library
	31	USE global ! global variable
	32	USE kind ! F90 kind parameter
	33	USE logger ! log file manager
	34	USE fct ! basic useful function
	35	USE extrap ! extrapolation manager
	36
	37	IMPLICIT NONE
	38	! NOTE_avoid_public_variables_if_possible
	39
	40	! type and variable
	41
	42	! function and subroutine
	43	PUBLIC :: interp_linear_fill !< compute interpolation using linear method
	44
	45	PRIVATE :: interp_linear__2D !< compute bilinear interpolation on 2D gid
	46	PRIVATE :: interp_linear__1D !< compute linear interpolation on 1D gid
	47	PRIVATE :: interp_linear__2D_coef !< compute coefficient for bilinear interpolation
	48	PRIVATE :: interp_linear__2D_fill !< fill value using bilinear interpolation
	49	PRIVATE :: interp_linear__1D_coef !< compute coefficient for linear interpolation
	50	PRIVATE :: interp_linear__1D_fill !< fill value using linear interpolation
	51	PRIVATE :: interp_linear__get_weight2D !< compute interpoaltion weight for 2D array.
	52	PRIVATE :: interp_linear__get_weight1D !< compute interpoaltion weight for 1D array.
	53
[12080]	54	CONTAINS
	55	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	56	SUBROUTINE interp_linear_fill(dd_value, dd_fill, id_detect, &
	57	& id_rho, ld_even, ld_discont)
[5037]	58	!-------------------------------------------------------------------
	59	!> @brief
[13369]	60	!> This subroutine compute horizontal linear interpolation on 4D array of value.
[5037]	61	!>
[13369]	62	!> @details
	63	!>
[5037]	64	!> @author J.Paul
[5617]	65	!> @date September, 2014 - Initial Version
[5609]	66	!> @date July, 2015 - reinitialise detect array for each level
[5037]	67	!>
[13369]	68	!> @param[inout] dd_value 2D array of variable value
[5037]	69	!> @param[in] dd_fill FillValue of variable
[13369]	70	!> @param[inout] id_detect 2D array of point to be interpolated
[5037]	71	!> @param[in] id_rho array of refinment factor
[13369]	72	!> @param[in] ld_even even refinment or not
[5037]	73	!> @param[in] ld_discont longitudinal discontinuity (-180°/180°, 0°/360°) or not
	74	!-------------------------------------------------------------------
[12080]	75
[5037]	76	IMPLICIT NONE
[12080]	77
[5037]	78	! Argument
[13369]	79	REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value
	80	REAL(dp) , INTENT(IN ) :: dd_fill
[5037]	81	INTEGER(I4) , DIMENSION(:,:,:) , INTENT(INOUT) :: id_detect
	82	INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho
	83	LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even
	84	LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont
	85
	86	! local variable
[5609]	87	INTEGER(i4), DIMENSION(4) :: il_shape
[5037]	88
[5609]	89	INTEGER(I4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect
	90
	91	LOGICAL :: ll_discont
[13369]	92
[5609]	93	REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_weight_IJ
	94	REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_weight_I
	95	REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_weight_J
[13369]	96
[5037]	97	! loop indices
	98	INTEGER(i4) :: ji
	99	INTEGER(i4) :: jj
	100	INTEGER(i4) :: jk
	101	INTEGER(i4) :: jl
	102	!----------------------------------------------------------------
	103	ll_discont=.FALSE.
	104	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
	105
	106	il_shape(:)=SHAPE(dd_value)
	107
	108	! compute vect2D
[5609]	109	ALLOCATE(dl_weight_IJ(4,((id_rho(jp_I)+1)*(id_rho(jp_J)+1))) )
[5037]	110	CALL interp_linear__get_weight2D(dl_weight_IJ(:,:), &
	111	& id_rho(:), ld_even(:))
	112
[5609]	113	ALLOCATE( dl_weight_I( 2,((id_rho(jp_I)+1) )) )
	114	ALLOCATE( dl_weight_J( 2,( (id_rho(jp_J)+1))) )
[5037]	115	CALL interp_linear__get_weight1D(dl_weight_I(:,:), &
	116	& id_rho(jp_I), ld_even(jp_I))
	117	CALL interp_linear__get_weight1D(dl_weight_J(:,:), &
	118	& id_rho(jp_J), ld_even(jp_J))
	119
[5609]	120	ALLOCATE(il_detect(il_shape(1),il_shape(2),il_shape(3)))
	121
[5037]	122	DO jl=1,il_shape(4)
[5609]	123	il_detect(:,:,:)=id_detect(:,:,:)
[5037]	124	! loop on vertical level
	125	DO jk=1,il_shape(3)
	126
	127	! I-J plan
	128	CALL interp_linear__2D(dd_value(:,:,jk,jl), dd_fill,&
[5609]	129	& il_detect(:,:,jk), &
[5037]	130	& dl_weight_IJ(:,:), &
	131	& id_rho(jp_I), id_rho(jp_J), &
[13369]	132	& ll_discont)
[5609]	133	IF( ANY(il_detect(:,:,jk)==1) )THEN
[5037]	134	! I direction
	135	DO jj=1,il_shape(2)
	136	CALL interp_linear__1D( dd_value(:,jj,jk,jl), dd_fill,&
[5609]	137	& il_detect(:,jj,jk), &
[5037]	138	& dl_weight_I(:,:), &
	139	& id_rho(jp_I), ll_discont )
	140	ENDDO
[5609]	141	IF( ALL(il_detect(:,:,jk)==0) )THEN
[5037]	142	CYCLE
	143	ELSE
	144	! J direction
	145	DO ji=1,il_shape(1)
	146	CALL interp_linear__1D( dd_value(ji,:,jk,jl), dd_fill,&
[5609]	147	& il_detect(ji,:,jk), &
[5037]	148	& dl_weight_J(:,:), &
	149	& id_rho(jp_J), ll_discont )
	150	ENDDO
	151	ENDIF
	152	ENDIF
	153
	154	ENDDO
	155	ENDDO
	156
[5609]	157	id_detect(:,:,:)=il_detect(:,:,:)
	158	DEALLOCATE(il_detect)
	159
[5037]	160	DEALLOCATE(dl_weight_IJ)
	161	DEALLOCATE(dl_weight_I)
	162	DEALLOCATE(dl_weight_J)
[13369]	163
[5037]	164	END SUBROUTINE interp_linear_fill
[12080]	165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	166	SUBROUTINE interp_linear__2D(dd_value, dd_fill, &
	167	& id_detect, &
	168	& dd_weight, &
	169	& id_rhoi, id_rhoj, &
	170	& ld_discont)
[5037]	171	!-------------------------------------------------------------------
	172	!> @brief
[13369]	173	!> This subroutine compute linear interpolation on 2D array of value.
[5037]	174	!>
[13369]	175	!> @details
	176	!>
[5037]	177	!> @author J.Paul
[5617]	178	!> @date September, 2014 - Initial Version
[5037]	179	!>
[13369]	180	!> @param[inout] dd_value 2D array of variable value
[5037]	181	!> @param[in] dd_fill FillValue of variable
[13369]	182	!> @param[inout] id_detect 2D array of point to be interpolated
[5037]	183	!> @param[in] id_rhoi refinment factor in i-direction
	184	!> @param[in] id_rhoj refinment factor in j-direction
	185	!> @param[in] id_rhok refinment factor in k-direction
[13369]	186	!> @param[in] ld_even even refinment or not
[5037]	187	!> @param[in] ld_discont longitudinal discontinuity (-180°/180°, 0°/360°) or not
	188	!-------------------------------------------------------------------
	189
	190	IMPLICIT NONE
[12080]	191
[5037]	192	! Argument
[13369]	193	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value
	194	REAL(dp) , INTENT(IN ) :: dd_fill
[5037]	195	INTEGER(I4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect
	196	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_weight
	197	INTEGER(I4) , INTENT(IN ) :: id_rhoi
	198	INTEGER(I4) , INTENT(IN ) :: id_rhoj
	199	LOGICAL , INTENT(IN ) :: ld_discont
	200
	201	! local variable
	202	INTEGER(I4) :: il_xextra
	203	INTEGER(I4) :: il_yextra
	204	INTEGER(i4), DIMENSION(2) :: il_shape
	205	INTEGER(i4), DIMENSION(2) :: il_dim
	206
[13369]	207	REAL(dp) :: dl_min
	208	REAL(dp) :: dl_max
[5037]	209	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_coef
	210	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_coarse
	211	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_tmp
	212
	213	! loop indices
	214	INTEGER(i4) :: ji
	215	INTEGER(i4) :: jj
	216	INTEGER(i4) :: ii
	217	INTEGER(i4) :: ij
	218
	219	!----------------------------------------------------------------
	220
	221	IF( ANY(id_detect(:,:)==1) )THEN
	222	il_shape(:)=SHAPE(dd_value)
	223
	224	! compute coarse grid dimension
	225	il_xextra=id_rhoi-1
	226	il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
	227
	228	il_yextra=id_rhoj-1
	229	il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj
	230
	231	ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )
	232
	233	! value on coarse grid
	234	dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &
	235	& 1:il_shape(2):id_rhoj )
	236
	237	ALLOCATE( dl_tmp(2,2) )
	238	ALLOCATE( dl_coef(4) )
	239
	240	DO jj=1,il_shape(2)-1,id_rhoj
	241	ij=((jj-1)/id_rhoj)+1
	242	DO ji=1,il_shape(1)-1,id_rhoi
	243	ii=((ji-1)/id_rhoi)+1
[13369]	244
[5037]	245	! check if point to be interpolated
	246	IF( ALL(id_detect(ji:ji+id_rhoi, &
	247	& jj:jj+id_rhoj)==0) ) CYCLE
[5609]	248	! check data needed to interpolate
[5037]	249	IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE
	250	! check longitude discontinuity
	251	dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
	252	IF( ld_discont )THEN
	253
	254	dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
	255	dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
	256	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
[13369]	257	WHERE( dl_tmp(:,:) < 0_dp )
[5037]	258	dl_tmp(:,:) = dl_tmp(:,:)+360._dp
	259	END WHERE
	260	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
[13369]	261	WHERE( dl_tmp(:,:) > 180_dp )
[5037]	262	dl_tmp(:,:) = dl_tmp(:,:)-180._dp
	263	END WHERE
	264	ENDIF
	265
	266	ENDIF
	267
	268	! compute bilinear coefficient
	269	dl_coef(:)=interp_linear__2D_coef(dl_tmp(:,:),&
	270	& dd_fill )
	271
	272	! compute value on detetected point
	273	CALL interp_linear__2D_fill(dd_value( ji:ji+id_rhoi, &
	274	& jj:jj+id_rhoj ), &
	275	& id_detect(ji:ji+id_rhoi, &
	276	& jj:jj+id_rhoj ), &
	277	& dd_weight(:,:), dl_coef(:),&
	278	& dd_fill, id_rhoi, id_rhoj )
	279
	280	IF( ld_discont )THEN
	281	WHERE( dd_value( ji:ji+id_rhoi, &
	282	& jj:jj+id_rhoj ) >= 180._dp .AND. &
	283	& dd_value( ji:ji+id_rhoi, &
	284	& jj:jj+id_rhoj ) /= dd_fill )
	285	dd_value( ji:ji+id_rhoi, &
	286	& jj:jj+id_rhoj ) = &
	287	& dd_value( ji:ji+id_rhoi, &
	288	& jj:jj+id_rhoj ) - 360._dp
	289	END WHERE
	290	ENDIF
	291
	292	ENDDO
	293	ENDDO
	294
	295	DEALLOCATE(dl_coef)
	296	DEALLOCATE(dl_tmp )
	297
	298	DEALLOCATE( dl_coarse )
	299	ENDIF
	300
	301	END SUBROUTINE interp_linear__2D
[12080]	302	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	303	SUBROUTINE interp_linear__1D(dd_value, dd_fill, &
	304	& id_detect, &
	305	& dd_weight, &
	306	& id_rhoi, &
	307	& ld_discont)
[5037]	308	!-------------------------------------------------------------------
	309	!> @brief
[13369]	310	!> This subroutine compute linear interpolation on 1D array of value.
[5037]	311	!>
[13369]	312	!> @details
	313	!>
[5037]	314	!> @author J.Paul
[5617]	315	!> @date September, 2014 - Initial Version
[5037]	316	!>
[13369]	317	!> @param[inout] dd_value 1D array of variable value
[5037]	318	!> @param[in] dd_fill FillValue of variable
[13369]	319	!> @param[inout] id_detect 1D array of point to be interpolated
[5037]	320	!> @param[in] id_rhoi refinment factor
[13369]	321	!> @param[in] ld_even even refinment or not
[5037]	322	!> @param[in] ld_discont longitudinal discontinuity (-180°/180°, 0°/360°) or not
	323	!-------------------------------------------------------------------
	324
	325	IMPLICIT NONE
[12080]	326
[5037]	327	! Argument
[13369]	328	REAL(dp) , DIMENSION(:) , INTENT(INOUT) :: dd_value
	329	REAL(dp) , INTENT(IN ) :: dd_fill
[5037]	330	INTEGER(I4) , DIMENSION(:) , INTENT(INOUT) :: id_detect
	331	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_weight
	332	INTEGER(I4) , INTENT(IN ) :: id_rhoi
	333	LOGICAL , INTENT(IN ) :: ld_discont
	334
	335	! local variable
	336	INTEGER(I4) :: il_xextra
	337	INTEGER(i4), DIMENSION(1) :: il_shape
	338	INTEGER(i4), DIMENSION(1) :: il_dim
	339
[13369]	340	REAL(dp) :: dl_min
	341	REAL(dp) :: dl_max
[5037]	342	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coef
	343	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coarse
	344	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_tmp
	345
	346	! loop indices
	347	INTEGER(i4) :: ji
	348	INTEGER(i4) :: ii
	349
	350	!----------------------------------------------------------------
	351
	352	IF( ANY(id_detect(:)==1) )THEN
	353	il_shape(:)=SHAPE(dd_value)
	354
	355	! compute coarse grid dimension
	356	il_xextra=id_rhoi-1
	357	il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
	358
	359	ALLOCATE( dl_coarse(il_dim(1)) )
	360
	361	! value on coarse grid
	362	dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )
	363
	364	ALLOCATE( dl_tmp(2) )
	365	ALLOCATE( dl_coef(4) )
	366
	367	DO ji=1,il_shape(1)-1,id_rhoi
	368	ii=((ji-1)/id_rhoi)+1
[13369]	369
[5037]	370	! check if point to be interpolated
	371	IF( ALL(id_detect(ji:ji+id_rhoi)==0) ) CYCLE
	372	! check data needed to interpolate
	373	IF( ANY(dl_coarse(ii:ii+1)==dd_fill) ) CYCLE
	374	! check longitude discontinuity
	375	dl_tmp(:)=dl_coarse(ii:ii+1)
	376	IF( ld_discont )THEN
	377
	378	dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
	379	dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
	380	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
[13369]	381	WHERE( dl_tmp(:) < 0_dp )
[5037]	382	dl_tmp(:) = dl_tmp(:)+360._dp
	383	END WHERE
	384	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
[13369]	385	WHERE( dl_tmp(:) > 180_dp )
[5037]	386	dl_tmp(:) = dl_tmp(:)-180._dp
	387	END WHERE
	388	ENDIF
	389
	390	ENDIF
	391
	392	! compute bilinear coefficient
	393	dl_coef(:)=interp_linear__1D_coef(dl_tmp(:), &
	394	& dd_fill )
	395
	396	! compute value on detetected point
	397	CALL interp_linear__1D_fill( dd_value( ji:ji+id_rhoi ), &
	398	& id_detect(ji:ji+id_rhoi ), &
	399	& dd_weight(:,:), dl_coef(:),&
	400	& dd_fill, id_rhoi )
	401
	402	IF( ld_discont )THEN
	403	WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
	404	& dd_value( ji:ji+id_rhoi ) /= dd_fill )
	405	dd_value(ji:ji+id_rhoi) = dd_value(ji:ji+id_rhoi) - 360._dp
	406	END WHERE
	407	ENDIF
	408
	409	ENDDO
	410
	411	DEALLOCATE(dl_coef)
	412	DEALLOCATE(dl_tmp )
	413
	414	DEALLOCATE( dl_coarse )
	415	ENDIF
	416
	417	END SUBROUTINE interp_linear__1D
[12080]	418	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	419	FUNCTION interp_linear__2D_coef(dd_value, dd_fill) &
	420	& RESULT (df_coef)
[5037]	421	!-------------------------------------------------------------------
	422	!> @brief
[13369]	423	!> This subroutine compute 2D array of coefficient for linear interpolation.
	424	!>
[5037]	425	!> @author J.Paul
[5617]	426	!> @date September, 2014 - Initial Version
[5037]	427	!>
	428	!> @param[in] dd_value 2D array of value
	429	!> @param[in] dd_fill FillValue of variable
	430	!-------------------------------------------------------------------
[12080]	431
[5037]	432	IMPLICIT NONE
[12080]	433
[5037]	434	! Argument
[13369]	435	REAL(dp), DIMENSION(:,:) , INTENT(IN) :: dd_value
	436	REAL(dp) , INTENT(IN) :: dd_fill
[5037]	437
	438	! function
[12080]	439	REAL(dp), DIMENSION(4) :: df_coef
[5037]	440
	441	! local variable
	442	REAL(dp), DIMENSION(4,4), PARAMETER :: dl_matrix = RESHAPE( &
[13369]	443	& (/ 1 ,-1 ,-1 , 1 ,&
	444	0 , 1 , 0 ,-1 ,&
	445	0 , 0 , 1 ,-1 ,&
[5037]	446	0 , 0 , 0 , 1 /), &
	447	& (/ 4, 4 /) )
	448
	449	REAL(dp), DIMENSION(4) :: dl_vect
	450
	451	!----------------------------------------------------------------
	452	! init
[12080]	453	df_coef(:)=dd_fill
[5037]	454
	455	dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
[12080]	456	df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
[5037]	457
	458	END FUNCTION interp_linear__2D_coef
[12080]	459	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	460	SUBROUTINE interp_linear__2D_fill(dd_value, id_detect, &
	461	& dd_weight, dd_coef, &
	462	& dd_fill, id_rhoi, id_rhoj)
[5037]	463	!-------------------------------------------------------------------
	464	!> @brief
[13369]	465	!> This subroutine compute linear interpolation of a 2D array of value.
	466	!>
[5037]	467	!> @author J.Paul
[5617]	468	!> @date September, 2014 - Initial Version
[5609]	469	!>
[5037]	470	!> @param[inout] dd_value 2D array of mixed grid value
	471	!> @param[inout] id_detect 2D array of point to be interpolated
	472	!> @param[in] dd_coef 2D array of coefficient
	473	!> @param[in] dd_fill FillValue of variable
	474	!> @param[in] ld_even even refinment or not
[13369]	475	!> @param[in] id_rhoi refinement factor in i-direction
[5037]	476	!> @param[in] id_rhoj refinement factor in j-direction
	477	!-------------------------------------------------------------------
[12080]	478
[5037]	479	IMPLICIT NONE
[12080]	480
[5037]	481	! Argument
[13369]	482	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value
[5037]	483	INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect
	484	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_weight
[13369]	485	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef
	486	REAL(dp) , INTENT(IN ) :: dd_fill
[5037]	487	INTEGER(I4) , INTENT(IN ) :: id_rhoi
	488	INTEGER(I4) , INTENT(IN ) :: id_rhoj
	489
	490	! local variable
	491
	492	! loop indices
	493	INTEGER(i4) :: ii
	494
	495	INTEGER(i4) :: ji
	496	INTEGER(i4) :: jj
	497	!----------------------------------------------------------------
	498
[5609]	499	IF( ANY( dd_coef(:)==dd_fill ) )THEN
	500	CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//&
	501	& "can not compute interpolation.")
	502	ELSE
[5037]	503
[5609]	504	ii=0
	505	DO jj=1,id_rhoj+1
	506	DO ji=1,id_rhoi+1
[5037]	507
[5609]	508	ii=ii+1
	509	IF(id_detect(ji,jj)==1)THEN
[5037]	510
[5609]	511	dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ii))
	512	id_detect(ji,jj)=0
[5037]	513
[5609]	514	ENDIF
[5037]	515
	516	ENDDO
[5609]	517	ENDDO
[5037]	518
[5609]	519	ENDIF
[5037]	520
	521	END SUBROUTINE interp_linear__2D_fill
[12080]	522	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	523	FUNCTION interp_linear__1D_coef(dd_value, dd_fill) &
	524	& RESULT (df_coef)
[5037]	525	!-------------------------------------------------------------------
	526	!> @brief
[13369]	527	!> This subroutine compute 1D array of coefficient for linear interpolation.
	528	!>
[5037]	529	!> @author J.Paul
[5617]	530	!> @date September, 2014 - Initial Version
[5037]	531	!>
	532	!> @param[in] dd_value 1D array of value
	533	!> @param[in] dd_fill FillValue of variable
	534	!-------------------------------------------------------------------
[12080]	535
[5037]	536	IMPLICIT NONE
[12080]	537
[5037]	538	! Argument
[13369]	539	REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value
	540	REAL(dp) , INTENT(IN) :: dd_fill
[5037]	541
	542	! function
[12080]	543	REAL(dp), DIMENSION(2) :: df_coef
[5037]	544
	545	! local variable
	546	REAL(dp), DIMENSION(2,2), PARAMETER :: dl_matrix = RESHAPE( &
[13369]	547	& (/ 1 ,-1 ,&
[5037]	548	0 , 1 /), &
	549	& (/ 2, 2 /) )
	550
	551	REAL(dp), DIMENSION(2) :: dl_vect
	552
	553	!----------------------------------------------------------------
	554	! init
[12080]	555	df_coef(:)=dd_fill
[5037]	556
	557	dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
[12080]	558	df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
[5037]	559
	560	END FUNCTION interp_linear__1D_coef
[12080]	561	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	562	SUBROUTINE interp_linear__1D_fill(dd_value, id_detect, &
	563	& dd_weight, dd_coef, &
	564	& dd_fill, id_rho)
[5037]	565	!-------------------------------------------------------------------
	566	!> @brief
[13369]	567	!> This subroutine compute linear interpolation of a 1D array of value.
	568	!>
[5037]	569	!> @author J.Paul
[5617]	570	!> @date September, 2014 - Initial Version
[13369]	571	!>
[5037]	572	!> @param[inout] dd_value 1D array of mixed grid value
	573	!> @param[inout] id_detect 1D array of point to be interpolated
	574	!> @param[in] dd_coef 1D array of coefficient
	575	!> @param[in] dd_fill FillValue of variable
	576	!> @param[in] ld_even even refinment or not
	577	!> @param[in] id_rho refinement factor
	578	!-------------------------------------------------------------------
[12080]	579
[5037]	580	IMPLICIT NONE
[12080]	581
[5037]	582	! Argument
[13369]	583	REAL(dp) , DIMENSION(:) , INTENT(INOUT) :: dd_value
[5037]	584	INTEGER(i4) , DIMENSION(:) , INTENT(INOUT) :: id_detect
	585	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_weight
[13369]	586	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef
	587	REAL(dp) , INTENT(IN ) :: dd_fill
	588	INTEGER(I4) , INTENT(IN ) :: id_rho
[5037]	589
	590	! local variable
	591
	592	! loop indices
	593	INTEGER(i4) :: ji
	594	!----------------------------------------------------------------
	595
	596	IF( ANY( dd_coef(:)==dd_fill ) )THEN
	597	CALL logger_error("INTERP LINEAR FILL: fill value detected. "//&
	598	& "can not compute interpolation")
	599	ELSE
	600
	601	DO ji=1,id_rho+1
[13369]	602
[5037]	603	IF(id_detect(ji)==1)THEN
	604
	605	dd_value(ji)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ji))
	606	id_detect(ji)=0
	607
	608	ENDIF
	609
	610	ENDDO
	611
	612	ENDIF
	613
	614	END SUBROUTINE interp_linear__1D_fill
[12080]	615	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	616	SUBROUTINE interp_linear__get_weight2D(dd_weight, id_rho, ld_even)
[5037]	617	!-------------------------------------------------------------------
	618	!> @brief
[13369]	619	!> This subroutine compute interpoaltion weight for 2D array.
	620	!>
[5037]	621	!> @author J.Paul
[5617]	622	!> @date September, 2014 - Initial Version
[5037]	623	!>
	624	!> @param[in] dd_weight interpolation weight of 2D array
	625	!> @param[in] ld_even even refinment or not
[13369]	626	!> @param[in] id_rho refinement factor
[5037]	627	!-------------------------------------------------------------------
[12080]	628
[5037]	629	IMPLICIT NONE
[12080]	630
[5037]	631	! Argument
[12080]	632
[5037]	633	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
	634	INTEGER(I4), DIMENSION(:) , INTENT(IN ) :: id_rho
	635	LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even
[12080]	636
[5037]	637	! local variable
	638	REAL(dp) :: dl_dx
[13369]	639	REAL(dp) :: dl_x
[5037]	640	REAL(dp) :: dl_dy
[13369]	641	REAL(dp) :: dl_y
	642
[5037]	643	! loop indices
[13369]	644	INTEGER(i4) :: ii
[5037]	645
[13369]	646	INTEGER(i4) :: ji
[5037]	647	INTEGER(i4) :: jj
	648	!----------------------------------------------------------------
[13369]	649
[5037]	650	IF( ld_even(jp_I) )THEN
[6393]	651	dl_dx=1._dp/REAL(id_rho(jp_I)-1,dp)
[5037]	652	ELSE ! odd refinement
[6393]	653	dl_dx=1._dp/REAL(id_rho(jp_I),dp)
[5037]	654	ENDIF
	655
	656	IF( ld_even(jp_J) )THEN
[6393]	657	dl_dy=1._dp/REAL(id_rho(jp_J)-1,dp)
[5037]	658	ELSE ! odd refinement
[6393]	659	dl_dy=1._dp/REAL(id_rho(jp_J),dp)
[5037]	660	ENDIF
	661
	662	ii=0
	663	DO jj=1,id_rho(jp_J)+1
	664
	665	IF( ld_even(jp_J) )THEN
[6393]	666	dl_y=REAL(jj-1,dp)dl_dy - dl_dy0.5_dp
[5037]	667	ELSE ! odd refinement
[13369]	668	dl_y=REAL(jj-1,dp)*dl_dy
[5037]	669	ENDIF
	670
	671	DO ji=1,id_rho(jp_I)+1
	672
	673	! iter
	674	ii=ii+1
	675
	676	IF( ld_even(jp_I) )THEN
[13369]	677	dl_x=REAL(ji-1,dp)dl_dx - dl_dx0.5_dp
[5037]	678	ELSE ! odd refinement
[13369]	679	dl_x=REAL(ji-1,dp)*dl_dx
[5037]	680	ENDIF
	681
	682	dd_weight(:,ii)=(/1._dp, dl_x, dl_y, dl_x*dl_y /)
	683
	684	ENDDO
	685	ENDDO
	686
	687	END SUBROUTINE interp_linear__get_weight2D
[12080]	688	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	689	SUBROUTINE interp_linear__get_weight1D(dd_weight, id_rho, ld_even)
[5037]	690	!-------------------------------------------------------------------
	691	!> @brief
[13369]	692	!> This subroutine compute interpoaltion weight for 1D array.
	693	!>
[5037]	694	!> @author J.Paul
[5617]	695	!> @date September, 2014 - Initial Version
[5037]	696	!>
	697	!> @param[in] dd_weight interpolation weight of 1D array
	698	!> @param[in] ld_even even refinment or not
[13369]	699	!> @param[in] id_rho refinement factor
[5037]	700	!-------------------------------------------------------------------
[12080]	701
[5037]	702	IMPLICIT NONE
[12080]	703
[5037]	704	! Argument
	705	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
	706	INTEGER(I4) , INTENT(IN ) :: id_rho
	707	LOGICAL , INTENT(IN ) :: ld_even
	708
	709	! local variable
	710	REAL(dp) :: dl_dx
[13369]	711	REAL(dp) :: dl_x
[5037]	712
	713	! loop indices
[13369]	714	INTEGER(i4) :: ji
[5037]	715	!----------------------------------------------------------------
[13369]	716
[5037]	717	IF( ld_even )THEN
[6393]	718	dl_dx=1._dp/REAL(id_rho-1,dp)
[5037]	719	ELSE ! odd refinement
[6393]	720	dl_dx=1._dp/REAL(id_rho,dp)
[5037]	721	ENDIF
	722
	723	DO ji=1,id_rho+1
	724	IF( ld_even )THEN
[13369]	725	dl_x=REAL(ji-1,dp)dl_dx - dl_dx0.5_dp
[5037]	726	ELSE ! odd refinement
[13369]	727	dl_x=REAL(ji-1,dp)*dl_dx
[5037]	728	ENDIF
	729
	730	dd_weight(:,ji)=(/1._dp, dl_x /)
	731	ENDDO
	732
	733	END SUBROUTINE interp_linear__get_weight1D
[12080]	734	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5037]	735	END MODULE interp_linear

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_linear.f90

Download in other formats: