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.

extrap.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/extrap.f90 @ 5608

Last change on this file since 5608 was 5608, checked in by jpaul, 9 years ago
commit changes/bugfix/... for SIREN; see ticket #1580
File size: 50.7 KB

Rev	Line
[4213]	1	!----------------------------------------------------------------------
	2	! NEMO system team, System and Interface for oceanic RElocable Nesting
	3	!----------------------------------------------------------------------
	4	!
	5	! MODULE: extrap
	6	!
	7	! DESCRIPTION:
	8	!> @brief
[5037]	9	!> This module manage extrapolation.
[4213]	10	!>
	11	!> @details
[5037]	12	!> Extrapolation method to be used is specify inside variable
	13	!> strcuture, as array of string character.<br/>
	14	!> - td_var\%c_extrap(1) string character is the interpolation name choose between:
	15	!> - 'dist_weight'
	16	!> - 'min_error'
[4213]	17	!>
[5037]	18	!> @note Extrapolation method could be specify for each variable in namelist _namvar_,
	19	!> defining string character _cn\_varinfo_. By default _dist_weight_.<br/>
	20	!> Example:
[5608]	21	!> - cn_varinfo='varname1:ext=dist_weight', 'varname2:ext=min_error'
[5037]	22	!>
	23	!> to detect point to be extrapolated:<br/>
	24	!> @code
[5608]	25	!> il_detect(:,:,:)=extrap_detect(td_var)
[5037]	26	!> @endcode
	27	!> - il_detect(:,:,:) is 3D array of point to be extrapolated
	28	!> - td_var is coarse grid variable to be extrapolated
	29	!>
	30	!> to extrapolate variable:<br/>
	31	!> @code
[5608]	32	!> CALL extrap_fill_value( td_var, [id_radius])
[5037]	33	!> @endcode
	34	!> - td_var is coarse grid variable to be extrapolated
	35	!> - id_radius is radius of the halo used to compute extrapolation [optional]
	36	!>
	37	!> to add extraband to the variable (to be extrapolated):<br/>
	38	!> @code
	39	!> CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] )
	40	!> @endcode
	41	!> - td_var is variable structure
	42	!> - id_isize : i-direction size of extra bands [optional]
	43	!> - id_jsize : j-direction size of extra bands [optional]
	44	!>
	45	!> to delete extraband of a variable:<br/>
	46	!> @code
	47	!> CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] )
	48	!> @endcode
	49	!> - td_var is variable structure
	50	!> - id_isize : i-direction size of extra bands [optional]
	51	!> - id_jsize : j-direction size of extra bands [optional]
	52	!>
	53	!> @warning _FillValue must not be zero (use var_chg_FillValue())
	54	!>
[4213]	55	!> @author
	56	!> J.Paul
	57	! REVISION HISTORY:
[5608]	58	!> @date November, 2013 - Initial Version
[5037]	59	!> @date September, 2014
	60	!> - add header
[5608]	61	!> @date June, 2015
	62	!> - extrapolate all land points (_FillValue)
	63	!> - move deriv function to math module
	64	!> @date July, 2015
	65	!> - compute extrapolation from north west to south east,
	66	!> and from south east to north west
[4213]	67	!>
[5037]	68	!> @todo
	69	!> - create module for each extrapolation method
[5608]	70	!> - smooth extrapolated points
[4213]	71	!>
	72	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	73	!----------------------------------------------------------------------
	74	MODULE extrap
	75	USE netcdf ! nf90 library
[5037]	76	USE kind ! F90 kind parameter
	77	USE phycst ! physical constant
	78	USE global ! global variable
	79	USE fct ! basic useful function
	80	USE date ! date manager
	81	USE logger ! log file manager
[5608]	82	USE math ! mathematical function
[5037]	83	USE att ! attribute manager
	84	USE dim ! dimension manager
	85	USE var ! variable manager
[4213]	86
	87	IMPLICIT NONE
	88	! NOTE_avoid_public_variables_if_possible
	89
	90	! type and variable
[5037]	91	PRIVATE :: im_minext !< default minumum number of point to extrapolate
	92	PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation
[4213]	93
	94	! function and subroutine
	95	PUBLIC :: extrap_detect !< detected point to be extrapolated
	96	PUBLIC :: extrap_fill_value !< extrapolate value over detected point
[5037]	97	PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated)
	98	PUBLIC :: extrap_del_extrabands !< delete extraband of the variable
[4213]	99
[5037]	100	PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper
	101	PRIVATE :: extrap__detect ! detected point to be extrapolated
	102	PRIVATE :: extrap__fill_value_wrapper ! extrapolate value over detected point wrapper
	103	PRIVATE :: extrap__fill_value ! extrapolate value over detected point
	104	PRIVATE :: extrap__3D !
	105	PRIVATE :: extrap__3D_min_error_coef !
	106	PRIVATE :: extrap__3D_min_error_fill !
	107	PRIVATE :: extrap__3D_dist_weight_coef !
	108	PRIVATE :: extrap__3D_dist_weight_fill !
[4213]	109
	110	INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate
	111	INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation
	112
	113	INTERFACE extrap_detect
[5037]	114	MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated
[4213]	115	END INTERFACE extrap_detect
	116
	117	INTERFACE extrap_fill_value
	118	MODULE PROCEDURE extrap__fill_value_wrapper !< detected point to be interpolated
	119	END INTERFACE extrap_fill_value
	120
	121	CONTAINS
	122	!-------------------------------------------------------------------
	123	!> @brief
[5037]	124	!> This function detected point to be extrapolated, given variable structure.
[4213]	125	!>
	126	!> @details
[5037]	127	!> optionaly, you could sepcify fine grid level, refinment factor (default 1),
	128	!> offset between fine and coarse grid (default compute from refinment factor
	129	!> as offset=(rho-1)/2), number of point to be extrapolated in each direction
	130	!> (default im_minext).<br/>
[4213]	131	!>
[5037]	132	!> First coarsening fine grid level, if need be, then select point near
	133	!> grid point already inform.
	134	!>
	135	!> @note point to be extrapolated are selected using FillValue,
	136	!> so to avoid mistake FillValue should not be zero (use var_chg_FillValue)
	137	!>
[4213]	138	!> @author J.Paul
[5608]	139	!> @date November, 2013 - Initial Version
	140	!> @date June, 2015
	141	!> - do not use level to select points to be extrapolated
[4213]	142	!
[5037]	143	!> @param[in] td_var0 coarse grid variable to extrapolate
	144	!> @return array of point to be extrapolated
[4213]	145	!-------------------------------------------------------------------
[5608]	146	FUNCTION extrap__detect( td_var0 )
[4213]	147	IMPLICIT NONE
	148	! Argument
[5037]	149	TYPE(TVAR) , INTENT(IN ) :: td_var0
[4213]	150
	151	! function
	152	INTEGER(i4), DIMENSION(td_var0%t_dim(1)%i_len,&
	153	& td_var0%t_dim(2)%i_len,&
	154	& td_var0%t_dim(3)%i_len ) :: extrap__detect
	155
	156	! local variable
	157	! loop indices
	158	INTEGER(i4) :: ji0
	159	INTEGER(i4) :: jj0
	160	INTEGER(i4) :: jk0
	161	!----------------------------------------------------------------
	162
[5608]	163	! force to extrapolated all points
	164	extrap__detect(:,:,:)=1
[4213]	165
[5608]	166	! do not compute grid point already inform
[5037]	167	DO jk0=1,td_var0%t_dim(3)%i_len
	168	DO jj0=1,td_var0%t_dim(2)%i_len
	169	DO ji0=1,td_var0%t_dim(1)%i_len
[5608]	170	IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill )THEN
	171	extrap__detect(ji0,jj0,jk0)=0
[5037]	172	ENDIF
	173	ENDDO
	174	ENDDO
	175	ENDDO
	176
[4213]	177	END FUNCTION extrap__detect
	178	!-------------------------------------------------------------------
	179	!> @brief
[5037]	180	!> This function sort variable to be extrapolated, depending on number of
	181	!> dimentsion, then detected point to be extrapolated.
[4213]	182	!>
[5037]	183	!> @author J.Paul
[5608]	184	!> @date November, 2013 - Initial Version
	185	!> @date June, 2015
	186	!> - select all land points for extrapolation
[4213]	187	!>
[5037]	188	!> @param[in] td_var coarse grid variable to extrapolate
	189	!> @return 3D array of point to be extrapolated
[4213]	190	!-------------------------------------------------------------------
[5608]	191	FUNCTION extrap__detect_wrapper( td_var )
[4213]	192
	193	IMPLICIT NONE
	194	! Argument
[5037]	195	TYPE(TVAR) , INTENT(IN ) :: td_var
[4213]	196
	197	! function
	198	INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len,&
	199	& td_var%t_dim(2)%i_len,&
	200	& td_var%t_dim(3)%i_len ) :: extrap__detect_wrapper
	201
	202	! local variable
	203	! loop indices
	204	!----------------------------------------------------------------
	205	! init
	206	extrap__detect_wrapper(:,:,:)=0
	207
	208	IF( .NOT. ANY(td_var%t_dim(1:3)%l_use) )THEN
	209	! no dimension I-J-K used
	210	CALL logger_debug(" EXTRAP DETECT: nothing done for variable"//&
	211	& TRIM(td_var%c_name) )
	212	ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN
	213
[5037]	214	! detect point to be extrapolated on I-J-K
[4213]	215	CALL logger_debug(" EXTRAP DETECT: detect point "//&
	216	& " for variable "//TRIM(td_var%c_name) )
	217
[5608]	218	extrap__detect_wrapper(:,:,:)=extrap__detect( td_var )
[4213]	219
	220	ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
	221
[5037]	222	! detect point to be extrapolated on I-J
[4213]	223	CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//&
	224	& " for variable "//TRIM(td_var%c_name) )
	225
[5608]	226	extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var )
[4213]	227
	228	ELSE IF( td_var%t_dim(3)%l_use )THEN
	229
[5037]	230	! detect point to be extrapolated on K
[4213]	231	CALL logger_debug(" EXTRAP DETECT: detect vertical point "//&
	232	& " for variable "//TRIM(td_var%c_name) )
	233
[5608]	234	extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var )
[4213]	235
	236	ENDIF
	237
	238	CALL logger_debug(" EXTRAP DETECT: "//&
	239	& TRIM(fct_str(SUM(extrap__detect_wrapper(:,:,:))))//&
	240	& " points to be extrapolated" )
	241
	242	END FUNCTION extrap__detect_wrapper
	243	!-------------------------------------------------------------------
	244	!> @brief
[5037]	245	!> This subroutine select method to be used for extrapolation.
	246	!> If need be, increase number of points to be extrapolated.
	247	!> Finally launch extrap__fill_value.
[4213]	248	!>
[5037]	249	!> @details
	250	!> optionaly, you could specify :<br/>
	251	!> - refinment factor (default 1)
	252	!> - offset between fine and coarse grid (default compute from refinment factor
	253	!> as offset=(rho-1)/2)
	254	!> - number of point to be extrapolated in each direction (default im_minext)
	255	!> - radius of the halo used to compute extrapolation
	256	!> - maximum number of iteration
[4213]	257	!>
	258	!> @author J.Paul
[5608]	259	!> @date November, 2013 - Initial Version
	260	!> @date June, 2015
	261	!> - select all land points for extrapolation
[4213]	262	!
[5037]	263	!> @param[inout] td_var variable structure
	264	!> @param[in] id_radius radius of the halo used to compute extrapolation
[4213]	265	!-------------------------------------------------------------------
[5608]	266	SUBROUTINE extrap__fill_value_wrapper( td_var, &
	267	& id_radius )
[4213]	268	IMPLICIT NONE
	269	! Argument
	270	TYPE(TVAR) , INTENT(INOUT) :: td_var
	271	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius
	272
	273	! local variable
	274	INTEGER(i4) :: il_radius
	275
	276	CHARACTER(LEN=lc) :: cl_method
	277
	278	! loop indices
	279	!----------------------------------------------------------------
	280	IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
[5037]	281	CALL logger_error("EXTRAP FILL VALUE: no value "//&
[4213]	282	& "associted to variable "//TRIM(td_var%c_name) )
	283	ELSE
	284
	285	SELECT CASE(TRIM(td_var%c_extrap(1)))
	286	CASE('min_error')
	287	cl_method='min_error'
	288	CASE DEFAULT
	289	cl_method='dist_weight'
	290
	291	!update variable structure
	292	td_var%c_extrap(1)='dist_weight'
	293	END SELECT
	294
[5608]	295	! number of point use to compute box
	296	il_radius=1
	297	IF( PRESENT(id_radius) ) il_radius=id_radius
	298	IF( il_radius < 0 )THEN
[4213]	299	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
[5608]	300	& " radius of the box used to compute extrapolation "//&
	301	& "("//TRIM(fct_str(il_radius))//")")
[4213]	302	ENDIF
	303
[5608]	304	CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//&
	305	& " using "//TRIM(cl_method)//" method." )
[4213]	306
[5608]	307	CALL extrap__fill_value( td_var, cl_method, &
	308	& il_radius )
[4213]	309
	310	ENDIF
	311
	312	END SUBROUTINE extrap__fill_value_wrapper
	313	!-------------------------------------------------------------------
	314	!> @brief
[5037]	315	!> This subroutine compute point to be extrapolated, then extrapolate point.
[4213]	316	!>
	317	!> @details
[5037]	318	!> optionaly, you could specify :<br/>
	319	!> - refinment factor (default 1)
	320	!> - offset between fine and coarse grid (default compute from refinment factor
	321	!> as offset=(rho-1)/2)
[4213]	322	!>
	323	!> @author J.Paul
[5608]	324	!> @date November, 2013 - Initial Version
	325	!> @date June, 2015
	326	!> - select all land points for extrapolation
[4213]	327	!
[5037]	328	!> @param[inout] td_var variable structure
	329	!> @param[in] cd_method extrapolation method
	330	!> @param[in] id_radius radius of the halo used to compute extrapolation
[4213]	331	!-------------------------------------------------------------------
	332	SUBROUTINE extrap__fill_value( td_var, cd_method, &
[5608]	333	& id_radius )
[4213]	334	IMPLICIT NONE
	335	! Argument
	336	TYPE(TVAR) , INTENT(INOUT) :: td_var
	337	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
	338	INTEGER(i4) , INTENT(IN ) :: id_radius
	339
	340	! local variable
	341	CHARACTER(LEN=lc) :: cl_extrap
	342
	343	INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
	344
	345	TYPE(TATT) :: tl_att
	346
	347	! loop indices
	348	!----------------------------------------------------------------
	349
	350	!1- detect point to be extrapolated
	351	ALLOCATE( il_detect( td_var%t_dim(1)%i_len, &
	352	& td_var%t_dim(2)%i_len, &
	353	& td_var%t_dim(3)%i_len) )
	354
[5608]	355	il_detect(:,:,:) = extrap_detect( td_var )
	356
[4213]	357	!2- add attribute to variable
	358	cl_extrap=fct_concat(td_var%c_extrap(:))
	359	tl_att=att_init('extrapolation',cl_extrap)
	360	CALL var_move_att(td_var, tl_att)
	361
[5037]	362	CALL att_clean(tl_att)
	363
[5608]	364	IF( ALL(il_detect(:,:,:)==1) )THEN
	365	CALL logger_warn(" EXTRAP FILL: "//&
	366	& " can not extrapolate "//TRIM(td_var%c_name)//&
	367	& ". no value inform." )
	368	ELSE
	369	CALL logger_info(" EXTRAP FILL: "//&
	370	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
	371	& " point(s) to extrapolate " )
[4213]	372
[5608]	373	CALL logger_info(" EXTRAP FILL: method "//&
	374	& TRIM(cd_method) )
[4213]	375
[5608]	376	!3- extrapolate
	377	CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, &
	378	& il_detect(:,:,:), &
	379	& cd_method, id_radius )
	380	ENDIF
	381
[4213]	382	DEALLOCATE(il_detect)
	383
	384	END SUBROUTINE extrap__fill_value
	385	!-------------------------------------------------------------------
	386	!> @brief
[5037]	387	!> This subroutine compute point to be extrapolated in 3D array.
[4213]	388	!>
	389	!> @details
[5037]	390	!> in case of 'min_error' method:<br/>
	391	!> - compute derivative in i-, j- and k- direction
	392	!> - compute minimum error coefficient (distance to center of halo)
	393	!> - compute extrapolatd values by calculated minimum error using taylor expansion
	394	!> in case of 'dist_weight' method:<br/>
	395	!> - compute distance weight coefficient (inverse of distance to center of halo)
	396	!> - compute extrapolatd values using Inverse Distance Weighting
[4213]	397	!>
	398	!> @author J.Paul
[5608]	399	!> @date November, 2013 - Initial Version
	400	!> @date July, 2015
	401	!> - compute coef indices to be used
[4213]	402	!
[5037]	403	!> @param[inout] dd_value 3D array of variable to be extrapolated
	404	!> @param[in] dd_fill FillValue of variable
	405	!> @param[inout] id_detect array of point to extrapolate
	406	!> @param[in] cd_method extrapolation method
	407	!> @param[in] id_radius radius of the halo used to compute extrapolation
[4213]	408	!-------------------------------------------------------------------
	409	SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,&
[5608]	410	& cd_method, id_radius )
[4213]	411	IMPLICIT NONE
	412	! Argument
[5037]	413	REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value
[5608]	414	REAL(dp) , INTENT(IN ) :: dd_fill
	415	INTEGER(i4), DIMENSION(:,:,:) , INTENT(INOUT) :: id_detect
	416	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
	417	INTEGER(i4), INTENT(IN ) :: id_radius
[4213]	418
	419	! local variable
[5608]	420	INTEGER(i4) :: il_imin
	421	INTEGER(i4) :: il_imax
	422	INTEGER(i4) :: il_jmin
	423	INTEGER(i4) :: il_jmax
	424	INTEGER(i4) :: il_kmin
	425	INTEGER(i4) :: il_kmax
	426	INTEGER(i4) :: il_iter
	427	INTEGER(i4) :: il_radius
	428	INTEGER(i4) :: il_i1
	429	INTEGER(i4) :: il_i2
	430	INTEGER(i4) :: il_j1
	431	INTEGER(i4) :: il_j2
	432	INTEGER(i4) :: il_k1
	433	INTEGER(i4) :: il_k2
[4213]	434
[5608]	435	INTEGER(i4), DIMENSION(4) :: il_shape
	436	INTEGER(i4), DIMENSION(3) :: il_dim
[4213]	437
[5037]	438	INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect
	439
[4213]	440	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
	441	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy
	442	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz
	443	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
	444
[5608]	445	LOGICAL :: ll_iter
	446
[4213]	447	! loop indices
	448	INTEGER(i4) :: ji
	449	INTEGER(i4) :: jj
	450	INTEGER(i4) :: jk
[5037]	451	INTEGER(i4) :: jl
[4213]	452	!----------------------------------------------------------------
	453
	454	il_shape(:)=SHAPE(dd_value)
	455
[5037]	456	ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) )
	457
[4213]	458	SELECT CASE(TRIM(cd_method))
[5037]	459	CASE('min_error')
	460	DO jl=1,il_shape(4)
[4213]	461
[5037]	462	! intitialise table of poitn to be extrapolated
	463	il_detect(:,:,:)=id_detect(:,:,:)
[4213]	464
[5037]	465	il_iter=1
	466	DO WHILE( ANY(il_detect(:,:,:)==1) )
	467	! change extend value to minimize number of iteration
[5608]	468	il_radius=id_radius+(il_iter-1)
	469	ll_iter=.TRUE.
[4213]	470
[5037]	471	ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )
	472	ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) )
	473	ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) )
[4213]	474
[5037]	475	! compute derivative in i-direction
	476	dl_dfdx(:,:,:)=dd_fill
	477	IF( il_shape(1) > 1 )THEN
[5608]	478	dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
	479	& dd_fill, 'I' )
[5037]	480	ENDIF
[4213]	481
[5037]	482	! compute derivative in j-direction
	483	dl_dfdy(:,:,:)=dd_fill
	484	IF( il_shape(2) > 1 )THEN
[5608]	485	dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
	486	& dd_fill, 'J' )
[5037]	487	ENDIF
	488
	489	! compute derivative in k-direction
	490	dl_dfdz(:,:,:)=dd_fill
	491	IF( il_shape(3) > 1 )THEN
[5608]	492	dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
	493	& dd_fill, 'K' )
[5037]	494	ENDIF
	495
	496	il_dim(1)=2*il_radius+1
	497	IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
	498	il_dim(2)=2*il_radius+1
	499	IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
	500	il_dim(3)=2*il_radius+1
	501	IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
[4213]	502
[5037]	503	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
[4213]	504
[5037]	505	dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), &
	506	& 1:il_dim(2), &
	507	& 1:il_dim(3), &
	508	& jl ))
[4213]	509
[5037]	510	DO jk=1,il_shape(3)
[5608]	511	! from North West(1,1) to South East(il_shape(1),il_shape(2))
[5037]	512	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
	513	DO jj=1,il_shape(2)
	514	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
	515	DO ji=1,il_shape(1)
[4213]	516
[5037]	517	IF( il_detect(ji,jj,jk) == 1 )THEN
	518
	519	il_imin=MAX(ji-il_radius,1)
	520	il_imax=MIN(ji+il_radius,il_shape(1))
[5608]	521	! coef indices to be used
	522	il_i1 = il_radius-(ji-il_imin)+1
	523	il_i2 = il_radius+(il_imax-ji)+1
[5037]	524	IF( il_dim(1) == 1 )THEN
	525	il_imin=ji
	526	il_imax=ji
[5608]	527	! coef indices to be used
	528	il_i1 = 1
	529	il_i2 = 2
[5037]	530	ENDIF
[4213]	531
[5608]	532
[5037]	533	il_jmin=MAX(jj-il_radius,1)
	534	il_jmax=MIN(jj+il_radius,il_shape(2))
[5608]	535	! coef indices to be used
	536	il_j1 = il_radius-(jj-il_jmin)+1
	537	il_j2 = il_radius+(il_jmax-jj)+1
[5037]	538	IF( il_dim(2) == 1 )THEN
	539	il_jmin=jj
	540	il_jmax=jj
[5608]	541	! coef indices to be used
	542	il_j1 = 1
	543	il_j2 = 2
[5037]	544	ENDIF
[4213]	545
[5037]	546	il_kmin=MAX(jk-il_radius,1)
	547	il_kmax=MIN(jk+il_radius,il_shape(3))
[5608]	548	! coef indices to be used
	549	il_k1 = il_radius-(jk-il_kmin)+1
	550	il_k2 = il_radius+(il_kmax-jk)+1
[5037]	551	IF( il_dim(3) == 1 )THEN
	552	il_kmin=jk
	553	il_kmax=jk
[5608]	554	! coef indices to be used
	555	il_k1 = 1
	556	il_k2 = 2
[5037]	557	ENDIF
[4213]	558
[5037]	559	dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( &
	560	& dd_value( il_imin:il_imax, &
	561	& il_jmin:il_jmax, &
	562	& il_kmin:il_kmax,jl ), dd_fill, il_radius, &
	563	& dl_dfdx( il_imin:il_imax, &
	564	& il_jmin:il_jmax, &
	565	& il_kmin:il_kmax ), &
	566	& dl_dfdy( il_imin:il_imax, &
	567	& il_jmin:il_jmax, &
	568	& il_kmin:il_kmax ), &
	569	& dl_dfdz( il_imin:il_imax, &
	570	& il_jmin:il_jmax, &
	571	& il_kmin:il_kmax ), &
[5608]	572	& dl_coef(il_i1:il_i2, &
	573	& il_j1:il_j2, &
	574	& il_k1:il_k2) )
[5037]	575
	576	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
	577	il_detect(ji,jj,jk)= 0
[5608]	578	ll_iter=.FALSE.
[5037]	579	ENDIF
	580
[4213]	581	ENDIF
	582
[5037]	583	ENDDO
[4213]	584	ENDDO
[5608]	585	! from South East(il_shape(1),il_shape(2)) to North West(1,1)
	586	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
	587	DO jj=il_shape(2),1,-1
	588	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
	589	DO ji=il_shape(1),1,-1
	590
	591	IF( il_detect(ji,jj,jk) == 1 )THEN
	592
	593	il_imin=MAX(ji-il_radius,1)
	594	il_imax=MIN(ji+il_radius,il_shape(1))
	595	! coef indices to be used
	596	il_i1 = il_radius-(ji-il_imin)+1
	597	il_i2 = il_radius+(il_imax-ji)+1
	598	IF( il_dim(1) == 1 )THEN
	599	il_imin=ji
	600	il_imax=ji
	601	! coef indices to be used
	602	il_i1 = 1
	603	il_i2 = 2
	604	ENDIF
	605
	606
	607	il_jmin=MAX(jj-il_radius,1)
	608	il_jmax=MIN(jj+il_radius,il_shape(2))
	609	! coef indices to be used
	610	il_j1 = il_radius-(jj-il_jmin)+1
	611	il_j2 = il_radius+(il_jmax-jj)+1
	612	IF( il_dim(2) == 1 )THEN
	613	il_jmin=jj
	614	il_jmax=jj
	615	! coef indices to be used
	616	il_j1 = 1
	617	il_j2 = 2
	618	ENDIF
	619
	620	il_kmin=MAX(jk-il_radius,1)
	621	il_kmax=MIN(jk+il_radius,il_shape(3))
	622	! coef indices to be used
	623	il_k1 = il_radius-(jk-il_kmin)+1
	624	il_k2 = il_radius+(il_kmax-jk)+1
	625	IF( il_dim(3) == 1 )THEN
	626	il_kmin=jk
	627	il_kmax=jk
	628	! coef indices to be used
	629	il_k1 = 1
	630	il_k2 = 2
	631	ENDIF
	632
	633	dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( &
	634	& dd_value( il_imin:il_imax, &
	635	& il_jmin:il_jmax, &
	636	& il_kmin:il_kmax,jl ), dd_fill, il_radius, &
	637	& dl_dfdx( il_imin:il_imax, &
	638	& il_jmin:il_jmax, &
	639	& il_kmin:il_kmax ), &
	640	& dl_dfdy( il_imin:il_imax, &
	641	& il_jmin:il_jmax, &
	642	& il_kmin:il_kmax ), &
	643	& dl_dfdz( il_imin:il_imax, &
	644	& il_jmin:il_jmax, &
	645	& il_kmin:il_kmax ), &
	646	& dl_coef(il_i1:il_i2, &
	647	& il_j1:il_j2, &
	648	& il_k1:il_k2) )
	649
	650	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
	651	il_detect(ji,jj,jk)= 0
	652	ll_iter=.FALSE.
	653	ENDIF
	654
	655	ENDIF
	656
	657	ENDDO
	658	ENDDO
[4213]	659	ENDDO
[5037]	660
	661	DEALLOCATE( dl_dfdx )
	662	DEALLOCATE( dl_dfdy )
	663	DEALLOCATE( dl_dfdz )
	664	DEALLOCATE( dl_coef )
	665
[5608]	666	IF( ll_iter ) il_iter=il_iter+1
[4213]	667	ENDDO
[5037]	668	ENDDO
[4213]	669
[5037]	670	CASE DEFAULT ! 'dist_weight'
	671	DO jl=1,il_shape(4)
[4213]	672
[5037]	673	! intitialise table of poitn to be extrapolated
	674	il_detect(:,:,:)=id_detect(:,:,:)
[4213]	675
[5037]	676	il_iter=1
	677	DO WHILE( ANY(il_detect(:,:,:)==1) )
	678	! change extend value to minimize number of iteration
[5608]	679	il_radius=id_radius+(il_iter-1)
	680	ll_iter=.TRUE.
[4213]	681
[5037]	682	il_dim(1)=2*il_radius+1
	683	IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
	684	il_dim(2)=2*il_radius+1
	685	IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
	686	il_dim(3)=2*il_radius+1
	687	IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
	688
	689	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
[4213]	690
[5608]	691	dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1),&
	692	& 1:il_dim(2),&
	693	& 1:il_dim(3),&
[5037]	694	& jl ) )
[5608]	695
[5037]	696	DO jk=1,il_shape(3)
[5608]	697	! from North West(1,1) to South East(il_shape(1),il_shape(2))
[5037]	698	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
	699	DO jj=1,il_shape(2)
	700	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
	701	DO ji=1,il_shape(1)
[4213]	702
[5037]	703	IF( il_detect(ji,jj,jk) == 1 )THEN
	704
	705	il_imin=MAX(ji-il_radius,1)
	706	il_imax=MIN(ji+il_radius,il_shape(1))
[5608]	707	! coef indices to be used
	708	il_i1 = il_radius-(ji-il_imin)+1
	709	il_i2 = il_radius+(il_imax-ji)+1
[5037]	710	IF( il_dim(1) == 1 )THEN
	711	il_imin=ji
	712	il_imax=ji
[5608]	713	! coef indices to be used
	714	il_i1 = 1
	715	il_i2 = 2
[5037]	716	ENDIF
[4213]	717
[5037]	718	il_jmin=MAX(jj-il_radius,1)
	719	il_jmax=MIN(jj+il_radius,il_shape(2))
[5608]	720	! coef indices to be used
	721	il_j1 = il_radius-(jj-il_jmin)+1
	722	il_j2 = il_radius+(il_jmax-jj)+1
[5037]	723	IF( il_dim(2) == 1 )THEN
	724	il_jmin=jj
	725	il_jmax=jj
[5608]	726	! coef indices to be used
	727	il_j1 = 1
	728	il_j2 = 2
[5037]	729	ENDIF
[4213]	730
[5037]	731	il_kmin=MAX(jk-il_radius,1)
	732	il_kmax=MIN(jk+il_radius,il_shape(3))
[5608]	733	! coef indices to be used
	734	il_k1 = il_radius-(jk-il_kmin)+1
	735	il_k2 = il_radius+(il_kmax-jk)+1
[5037]	736	IF( il_dim(3) == 1 )THEN
	737	il_kmin=jk
	738	il_kmax=jk
[5608]	739	! coef indices to be used
	740	il_k1 = 1
	741	il_k2 = 2
[5037]	742	ENDIF
[4213]	743
[5037]	744	dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( &
	745	& dd_value( il_imin:il_imax, &
	746	& il_jmin:il_jmax, &
	747	& il_kmin:il_kmax, &
	748	& jl), dd_fill, il_radius, &
[5608]	749	& dl_coef(il_i1:il_i2, &
	750	& il_j1:il_j2, &
	751	& il_k1:il_k2) )
[5037]	752
	753	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
	754	il_detect(ji,jj,jk)= 0
[5608]	755	ll_iter=.FALSE.
[5037]	756	ENDIF
	757
[4213]	758	ENDIF
	759
[5037]	760	ENDDO
[4213]	761	ENDDO
[5608]	762	! from South East(il_shape(1),il_shape(2)) to North West(1,1)
	763	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
	764	DO jj=il_shape(2),1,-1
	765	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
	766	DO ji=il_shape(1),1,-1
[5037]	767
[5608]	768	IF( il_detect(ji,jj,jk) == 1 )THEN
	769
	770	il_imin=MAX(ji-il_radius,1)
	771	il_imax=MIN(ji+il_radius,il_shape(1))
	772	! coef indices to be used
	773	il_i1 = il_radius-(ji-il_imin)+1
	774	il_i2 = il_radius+(il_imax-ji)+1
	775	IF( il_dim(1) == 1 )THEN
	776	il_imin=ji
	777	il_imax=ji
	778	! coef indices to be used
	779	il_i1 = 1
	780	il_i2 = 2
	781	ENDIF
[4213]	782
[5608]	783	il_jmin=MAX(jj-il_radius,1)
	784	il_jmax=MIN(jj+il_radius,il_shape(2))
	785	! coef indices to be used
	786	il_j1 = il_radius-(jj-il_jmin)+1
	787	il_j2 = il_radius+(il_jmax-jj)+1
	788	IF( il_dim(2) == 1 )THEN
	789	il_jmin=jj
	790	il_jmax=jj
	791	! coef indices to be used
	792	il_j1 = 1
	793	il_j2 = 2
	794	ENDIF
[5037]	795
[5608]	796	il_kmin=MAX(jk-il_radius,1)
	797	il_kmax=MIN(jk+il_radius,il_shape(3))
	798	! coef indices to be used
	799	il_k1 = il_radius-(jk-il_kmin)+1
	800	il_k2 = il_radius+(il_kmax-jk)+1
	801	IF( il_dim(3) == 1 )THEN
	802	il_kmin=jk
	803	il_kmax=jk
	804	! coef indices to be used
	805	il_k1 = 1
	806	il_k2 = 2
	807	ENDIF
[4213]	808
[5608]	809	dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( &
	810	& dd_value( il_imin:il_imax, &
	811	& il_jmin:il_jmax, &
	812	& il_kmin:il_kmax, &
	813	& jl), dd_fill, il_radius, &
	814	& dl_coef(il_i1:il_i2, &
	815	& il_j1:il_j2, &
	816	& il_k1:il_k2) )
[4213]	817
[5608]	818	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
	819	il_detect(ji,jj,jk)= 0
	820	ll_iter=.FALSE.
	821	ENDIF
[4213]	822
[5608]	823	ENDIF
[4213]	824
[5608]	825	ENDDO
	826	ENDDO
	827	ENDDO
	828	CALL logger_info(" EXTRAP 3D: "//&
	829	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
	830	& " point(s) to extrapolate " )
[4213]	831
[5608]	832	DEALLOCATE( dl_coef )
	833	IF( ll_iter ) il_iter=il_iter+1
	834	ENDDO
	835	ENDDO
[4213]	836	END SELECT
	837
[5608]	838	DEALLOCATE( il_detect )
[4213]	839
[5608]	840	END SUBROUTINE extrap__3D
[4213]	841	!-------------------------------------------------------------------
	842	!> @brief
[5037]	843	!> This function compute coefficient for min_error extrapolation.
[4213]	844	!>
	845	!> @details
[5608]	846	!> coefficients are "grid distance" to the center of the box
	847	!> choosed to compute extrapolation.
[4213]	848	!>
	849	!> @author J.Paul
[5608]	850	!> @date November, 2013 - Initial Version
	851	!> @date July, 2015
	852	!> - decrease weight of third dimension
[4213]	853	!
[5037]	854	!> @param[in] dd_value 3D array of variable to be extrapolated
	855	!> @return 3D array of coefficient for minimum error extrapolation
[4213]	856	!-------------------------------------------------------------------
	857	PURE FUNCTION extrap__3D_min_error_coef( dd_value )
	858
	859	IMPLICIT NONE
	860	! Argument
	861	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
	862
	863	! function
	864	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
	865	& SIZE(dd_value(:,:,:),DIM=2), &
	866	& SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_min_error_coef
	867
	868	! local variable
	869	INTEGER(i4), DIMENSION(3) :: il_shape
	870
	871	INTEGER(i4) :: il_imid
	872	INTEGER(i4) :: il_jmid
	873	INTEGER(i4) :: il_kmid
	874
	875	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
	876
	877	! loop indices
	878	INTEGER(i4) :: ji
	879	INTEGER(i4) :: jj
	880	INTEGER(i4) :: jk
	881	!----------------------------------------------------------------
	882
	883	! init
	884	extrap__3D_min_error_coef(:,:,:)=0
	885
	886	il_shape(:)=SHAPE(dd_value(:,:,:))
	887
	888	il_imid=INT(REAL(il_shape(1),sp)*0.5+1)
	889	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1)
	890	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1)
	891
	892	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
	893
	894	DO jk=1,il_shape(3)
	895	DO jj=1,il_shape(2)
	896	DO ji=1,il_shape(1)
	897
	898	! compute distance
[5608]	899	! "vertical weight" is lower than horizontal
[4213]	900	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
	901	& (jj-il_jmid)**2 + &
[5608]	902	& 3(jk-il_kmid)*2
[4213]	903
	904	IF( dl_dist(ji,jj,jk) /= 0 )THEN
	905	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
	906	ENDIF
	907
	908	ENDDO
	909	ENDDO
	910	ENDDO
	911
	912	WHERE( dl_dist(:,:,:) /= 0 )
	913	extrap__3D_min_error_coef(:,:,:)=dl_dist(:,:,:)
	914	END WHERE
	915
	916	DEALLOCATE( dl_dist )
	917
	918	END FUNCTION extrap__3D_min_error_coef
	919	!-------------------------------------------------------------------
	920	!> @brief
[5037]	921	!> This function compute extrapolatd value by calculated minimum error using
[4213]	922	!> taylor expansion
	923	!>
[5037]	924	!> @author J.Paul
	925	!> - November, 2013- Initial Version
[4213]	926	!>
[5037]	927	!> @param[in] dd_value 3D array of variable to be extrapolated
	928	!> @param[in] dd_fill FillValue of variable
	929	!> @param[in] id_radius radius of the halo used to compute extrapolation
	930	!> @param[in] dd_dfdx derivative of function in i-direction
	931	!> @param[in] dd_dfdy derivative of function in j-direction
	932	!> @param[in] dd_dfdz derivative of function in k-direction
	933	!> @param[in] dd_coef array of coefficient for min_error extrapolation
	934	!> @return extrapolatd value
[4213]	935	!-------------------------------------------------------------------
[5037]	936	PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, &
[4213]	937	& dd_dfdx, dd_dfdy, dd_dfdz, &
	938	& dd_coef )
	939	IMPLICIT NONE
	940	! Argument
	941	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
	942	REAL(dp) , INTENT(IN) :: dd_fill
[5037]	943	INTEGER(i4), INTENT(IN) :: id_radius
[4213]	944	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx
	945	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy
	946	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdz
	947	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
	948
	949	! function
	950	REAL(dp) :: extrap__3d_min_error_fill
	951
	952	! local variable
	953	INTEGER(i4), DIMENSION(3) :: il_shape
	954	INTEGER(i4), DIMENSION(3) :: il_ind
	955
	956	INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_mask
	957	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_error
	958
	959	INTEGER(i4) :: il_min
	960	! loop indices
	961
	962	!----------------------------------------------------------------
	963
	964	! init
	965	extrap__3D_min_error_fill=dd_fill
	966
[5037]	967	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
[4213]	968
	969	IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN
	970
	971	il_shape(:)=SHAPE(dd_value(:,:,:))
	972	ALLOCATE( il_mask( il_shape(1),il_shape(2),il_shape(3)) )
	973	ALLOCATE( dl_error(il_shape(1),il_shape(2),il_shape(3)) )
	974
	975	! compute error
	976	dl_error(:,:,:)=0.
	977	il_mask(:,:,:)=0
	978	WHERE( dd_dfdx(:,:,:) /= dd_fill )
	979	dl_error(:,:,:)=dd_coef(:,:,:)*dd_dfdx(:,:,:)
	980	il_mask(:,:,:)=1
	981	END WHERE
	982	WHERE( dd_dfdy(:,:,:) /= dd_fill )
	983	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdy(:,:,:))
	984	il_mask(:,:,:)=1
	985	END WHERE
	986	WHERE( dd_dfdz(:,:,:) /= dd_fill )
	987	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdz(:,:,:))
	988	il_mask(:,:,:)=1
	989	END WHERE
	990
	991	! get minimum error indices
	992	il_ind(:)=MINLOC(dl_error(:,:,:),il_mask(:,:,:)==1)
	993
	994	! return value
	995	IF( ALL(il_ind(:)/=0) )THEN
	996	extrap__3D_min_error_fill=dd_value(il_ind(1),il_ind(2),il_ind(3))
	997	ENDIF
	998
	999	DEALLOCATE( il_mask )
	1000	DEALLOCATE( dl_error )
	1001
	1002	ENDIF
	1003
	1004	END FUNCTION extrap__3D_min_error_fill
	1005	!-------------------------------------------------------------------
	1006	!> @brief
[5037]	1007	!> This function compute coefficient for inverse distance weighted method
[4213]	1008	!>
	1009	!> @details
[5037]	1010	!> coefficients are inverse "grid distance" to the center of the box choosed to compute
	1011	!> extrapolation.
[4213]	1012	!>
	1013	!> @author J.Paul
[5608]	1014	!> @date November, 2013 - Initial Version
	1015	!> @date July, 2015
	1016	!> - decrease weight of third dimension
[4213]	1017	!
[5037]	1018	!> @param[in] dd_value 3D array of variable to be extrapolated
	1019	!> @return 3D array of coefficient for inverse distance weighted extrapolation
[4213]	1020	!-------------------------------------------------------------------
	1021	PURE FUNCTION extrap__3D_dist_weight_coef( dd_value )
	1022
	1023	IMPLICIT NONE
	1024	! Argument
	1025	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
	1026
	1027	! function
	1028	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
	1029	& SIZE(dd_value(:,:,:),DIM=2), &
	1030	& SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_dist_weight_coef
	1031
	1032	! local variable
	1033	INTEGER(i4), DIMENSION(3) :: il_shape
	1034
	1035	INTEGER(i4) :: il_imid
	1036	INTEGER(i4) :: il_jmid
	1037	INTEGER(i4) :: il_kmid
	1038
	1039	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
	1040
	1041	! loop indices
	1042	INTEGER(i4) :: ji
	1043	INTEGER(i4) :: jj
	1044	INTEGER(i4) :: jk
	1045	!----------------------------------------------------------------
	1046
	1047	! init
	1048	extrap__3D_dist_weight_coef(:,:,:)=0
	1049
	1050	il_shape(:)=SHAPE(dd_value(:,:,:))
	1051
	1052	il_imid=INT(REAL(il_shape(1),sp)*0.5+1,i4)
	1053	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1,i4)
	1054	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1,i4)
	1055
	1056	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
	1057
	1058	DO jk=1,il_shape(3)
	1059	DO jj=1,il_shape(2)
	1060	DO ji=1,il_shape(1)
	1061
	1062	! compute distance
[5608]	1063	! "vertical weight" is lower than horizontal
[4213]	1064	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
	1065	& (jj-il_jmid)**2 + &
[5608]	1066	& 3(jk-il_kmid)*2
[4213]	1067
	1068	IF( dl_dist(ji,jj,jk) /= 0 )THEN
	1069	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
	1070	ENDIF
	1071
	1072	ENDDO
	1073	ENDDO
	1074	ENDDO
	1075
	1076	WHERE( dl_dist(:,:,:) /= 0 )
	1077	extrap__3D_dist_weight_coef(:,:,:)=1./dl_dist(:,:,:)
	1078	END WHERE
	1079
	1080	DEALLOCATE( dl_dist )
	1081
	1082	END FUNCTION extrap__3D_dist_weight_coef
	1083	!-------------------------------------------------------------------
	1084	!> @brief
[5037]	1085	!> This function compute extrapolatd value using inverse distance weighted
	1086	!> method
[4213]	1087	!>
	1088	!> @details
	1089	!>
	1090	!> @author J.Paul
[5608]	1091	!> - November, 2013 - Initial Version
[4213]	1092	!
[5037]	1093	!> @param[in] dd_value 3D array of variable to be extrapolated
	1094	!> @param[in] dd_fill FillValue of variable
	1095	!> @param[in] id_radius radius of the halo used to compute extrapolation
	1096	!> @param[in] dd_coef 3D array of coefficient for inverse distance weighted extrapolation
	1097	!> @return extrapolatd value
[4213]	1098	!-------------------------------------------------------------------
[5037]	1099	FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, &
	1100	& dd_coef )
[4213]	1101	IMPLICIT NONE
	1102	! Argument
	1103	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
	1104	REAL(dp) , INTENT(IN) :: dd_fill
[5037]	1105	INTEGER(i4), INTENT(IN) :: id_radius
[4213]	1106	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
	1107
	1108	! function
	1109	REAL(dp) :: extrap__3D_dist_weight_fill
	1110
	1111	! local variable
	1112	INTEGER(i4), DIMENSION(3) :: il_shape
	1113
	1114	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value
	1115	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
	1116
	1117	INTEGER(i4) :: il_min
	1118	! loop indices
	1119	INTEGER(i4) :: ji
	1120	INTEGER(i4) :: jj
	1121	INTEGER(i4) :: jk
	1122	!----------------------------------------------------------------
	1123
	1124	! init
	1125	extrap__3D_dist_weight_fill=dd_fill
	1126
[5037]	1127	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
[4213]	1128
	1129	IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN
	1130
	1131	il_shape(:)=SHAPE(dd_value(:,:,:))
	1132	ALLOCATE( dl_value( il_shape(1),il_shape(2),il_shape(3)) )
	1133	ALLOCATE( dl_coef( il_shape(1),il_shape(2),il_shape(3)) )
	1134
	1135	dl_value(:,:,:)=0
	1136	dl_coef(:,:,:)=0
	1137
[5037]	1138	DO jk=1,il_shape(3)
	1139	DO jj=1,il_shape(2)
	1140	DO ji=1,il_shape(1)
[4213]	1141
[5037]	1142	IF( dd_value(ji,jj,jk) /= dd_fill )THEN
	1143	! compute factor
	1144	dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk)
	1145	dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk)
	1146	ENDIF
[4213]	1147
[5037]	1148	ENDDO
	1149	ENDDO
	1150	ENDDO
[4213]	1151
[5608]	1152
[4213]	1153	! return value
	1154	IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN
[5037]	1155	extrap__3D_dist_weight_fill = &
	1156	& SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) )
[4213]	1157	ENDIF
	1158
	1159	DEALLOCATE( dl_value )
	1160	DEALLOCATE( dl_coef )
	1161
	1162	ENDIF
	1163
	1164	END FUNCTION extrap__3D_dist_weight_fill
	1165	!-------------------------------------------------------------------
	1166	!> @brief
	1167	!> This subroutine add to the variable (to be extrapolated) an
	1168	!> extraband of N points at north,south,east and west boundaries.
	1169	!>
[5037]	1170	!> @details
	1171	!> optionaly you could specify size of extra bands in i- and j-direction
	1172	!>
[4213]	1173	!> @author J.Paul
[5037]	1174	!> - November, 2013-Initial version
[4213]	1175	!
[5037]	1176	!> @param[inout] td_var variable
	1177	!> @param[in] id_isize i-direction size of extra bands (default=im_minext)
	1178	!> @param[in] id_jsize j-direction size of extra bands (default=im_minext)
[4213]	1179	!> @todo
	1180	!> - invalid special case for grid with north fold
	1181	!-------------------------------------------------------------------
	1182	SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize )
	1183	IMPLICIT NONE
	1184	! Argument
	1185	TYPE(TVAR) , INTENT(INOUT) :: td_var
	1186	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
	1187	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
	1188
	1189	! local variable
	1190	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
	1191
	1192	INTEGER(i4) :: il_isize
	1193	INTEGER(i4) :: il_jsize
	1194	INTEGER(i4) :: il_tmp
	1195
	1196	! loop indices
	1197	INTEGER(i4) :: ji
	1198	INTEGER(i4) :: ij
	1199	!----------------------------------------------------------------
	1200	il_isize=im_minext
	1201	IF(PRESENT(id_isize)) il_isize=id_isize
	1202	IF( il_isize < im_minext .AND. &
	1203	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
	1204	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
	1205	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
	1206	& " cubic interpolation ")
	1207	ENDIF
	1208
	1209	il_jsize=im_minext
	1210	IF(PRESENT(id_jsize)) il_jsize=id_jsize
	1211	IF( il_jsize < im_minext .AND. &
	1212	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
	1213	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
	1214	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
	1215	& " cubic interpolation ")
	1216	ENDIF
	1217
	1218	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
	1219	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
	1220
	1221	CALL logger_trace( "EXTRAP ADD EXTRABANDS: dimension change "//&
	1222	& "in variable "//TRIM(td_var%c_name) )
	1223
	1224	! add extrabands in variable
	1225	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
	1226	& td_var%t_dim(2)%i_len, &
	1227	& td_var%t_dim(3)%i_len, &
	1228	& td_var%t_dim(4)%i_len ))
	1229
	1230	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
	1231
[5037]	1232
[4213]	1233	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize
	1234	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize
	1235
	1236	DEALLOCATE(td_var%d_value)
	1237	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
	1238	& td_var%t_dim(2)%i_len, &
	1239	& td_var%t_dim(3)%i_len, &
	1240	& td_var%t_dim(4)%i_len ) )
	1241
	1242	! intialise
	1243	td_var%d_value(:,:,:,:)=td_var%d_fill
	1244
	1245	! fill center
	1246	td_var%d_value( 1+il_isize:td_var%t_dim(1)%i_len-il_isize, &
	1247	& 1+il_jsize:td_var%t_dim(2)%i_len-il_jsize, &
	1248	& :,:) = dl_value(:,:,:,:)
	1249
	1250	! special case for overlap
	1251	IF( td_var%i_ew >= 0 .AND. il_isize /= 0 )THEN
	1252	DO ji=1,il_isize
	1253	! from east to west
	1254	il_tmp=td_var%t_dim(1)%i_len-td_var%i_ew+ji-2*il_isize
	1255	td_var%d_value(ji,:,:,:) = td_var%d_value(il_tmp,:,:,:)
	1256
	1257	! from west to east
	1258	ij=td_var%t_dim(1)%i_len-ji+1
	1259	il_tmp=td_var%i_ew-ji+2*il_isize+1
	1260	td_var%d_value(ij,:,:,:) = td_var%d_value(il_tmp,:,:,:)
	1261	ENDDO
	1262	ENDIF
	1263
	1264	DEALLOCATE( dl_value )
	1265
	1266	END SUBROUTINE extrap_add_extrabands
	1267	!-------------------------------------------------------------------
	1268	!> @brief
	1269	!> This subroutine remove of the variable an extraband
	1270	!> of N points at north,south,east and west boundaries.
	1271	!>
[5037]	1272	!> @details
	1273	!> optionaly you could specify size of extra bands in i- and j-direction
	1274	!>
[4213]	1275	!> @author J.Paul
[5608]	1276	!> - November, 2013 - Initial version
[5037]	1277	!>
	1278	!> @param[inout] td_var variable
	1279	!> @param[in] id_isize i-direction size of extra bands (default=im_minext)
	1280	!> @param[in] id_jsize j-direction size of extra bands (default=im_minext)
[4213]	1281	!-------------------------------------------------------------------
	1282	SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize )
	1283	IMPLICIT NONE
	1284	! Argument
[5037]	1285	TYPE(TVAR) , INTENT(INOUT) :: td_var
[4213]	1286	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
	1287	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
	1288
	1289	! local variable
	1290	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
	1291
	1292	INTEGER(i4) :: il_isize
	1293	INTEGER(i4) :: il_jsize
[5037]	1294
[4213]	1295	INTEGER(i4) :: il_imin
	1296	INTEGER(i4) :: il_imax
	1297	INTEGER(i4) :: il_jmin
	1298	INTEGER(i4) :: il_jmax
	1299
	1300	! loop indices
	1301	!----------------------------------------------------------------
	1302	il_isize=im_minext
	1303	IF(PRESENT(id_isize)) il_isize=id_isize
	1304
	1305	il_jsize=im_minext
	1306	IF(PRESENT(id_jsize)) il_jsize=id_jsize
	1307
	1308	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
	1309	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
	1310
	1311	CALL logger_trace( "EXTRAP DEL EXTRABANDS: dimension change "//&
	1312	& "in variable "//TRIM(td_var%c_name) )
	1313
	1314	! add extrabands in variable
	1315	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
	1316	& td_var%t_dim(2)%i_len, &
	1317	& td_var%t_dim(3)%i_len, &
	1318	& td_var%t_dim(4)%i_len ))
	1319
	1320	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
	1321
	1322	! fill center
	1323	il_imin=1+il_isize
	1324	il_imax=td_var%t_dim(1)%i_len-il_isize
	1325
	1326	il_jmin=1+il_jsize
	1327	il_jmax=td_var%t_dim(2)%i_len-il_jsize
	1328
	1329	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len - 2*il_isize
	1330	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len - 2*il_jsize
	1331
	1332	DEALLOCATE(td_var%d_value)
	1333	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
	1334	& td_var%t_dim(2)%i_len, &
	1335	& td_var%t_dim(3)%i_len, &
	1336	& td_var%t_dim(4)%i_len ) )
	1337
	1338	! intialise
	1339	td_var%d_value(:,:,:,:)=td_var%d_fill
	1340
	1341	td_var%d_value(:,:,:,:)=dl_value(il_imin:il_imax,&
	1342	& il_jmin:il_jmax,&
	1343	& :,:)
	1344
	1345	DEALLOCATE( dl_value )
	1346
	1347	END SUBROUTINE extrap_del_extrabands
	1348	END MODULE extrap

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

Download in other formats: