- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/TOOLS/SIREN/src/extrap.f90
r4213 r6225 7 7 ! DESCRIPTION: 8 8 !> @brief 9 !> This module 9 !> This module manage extrapolation. 10 10 !> 11 11 !> @details 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' 17 !> 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: 21 !> - cn_varinfo='varname1:ext=dist_weight', 'varname2:ext=min_error' 22 !> 23 !> to detect point to be extrapolated:<br/> 24 !> @code 25 !> il_detect(:,:,:)=extrap_detect(td_var) 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 32 !> CALL extrap_fill_value( td_var, [id_radius]) 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()) 12 54 !> 13 55 !> @author 14 56 !> J.Paul 15 57 ! REVISION HISTORY: 16 !> @date Nov, 2013 - Initial Version 58 !> @date November, 2013 - Initial Version 59 !> @date September, 2014 60 !> - add header 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 17 67 !> 18 !> @note WARNING: FillValue must not be zero (use var_chg_FillValue) 68 !> @todo 69 !> - create module for each extrapolation method 70 !> - smooth extrapolated points 19 71 !> 20 72 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 21 !> @todo22 !> - something wrong when computing point to be extralopated23 !> - take care of ew value in variable structure24 73 !---------------------------------------------------------------------- 25 74 MODULE extrap 26 75 USE netcdf ! nf90 library 27 USE kind 28 USE phycst 29 USE global 30 USE fct 31 USE logger 32 USE dim 33 USE att 34 USE var 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 82 USE math ! mathematical function 83 USE att ! attribute manager 84 USE dim ! dimension manager 85 USE var ! variable manager 35 86 36 87 IMPLICIT NONE 37 PRIVATE38 88 ! NOTE_avoid_public_variables_if_possible 39 89 40 90 ! type and variable 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 41 93 42 94 ! function and subroutine 43 95 PUBLIC :: extrap_detect !< detected point to be extrapolated 44 96 PUBLIC :: extrap_fill_value !< extrapolate value over detected point 45 PUBLIC :: extrap_add_extrabands !< 46 PUBLIC :: extrap_del_extrabands !< 47 PUBLIC :: extrap_deriv_1D !< 48 PUBLIC :: extrap_deriv_2D !< 49 50 PRIVATE :: extrap__detect_wrapper !< detected point to be extrapolated 51 PRIVATE :: extrap__detect !< detected point to be extrapolated 52 PRIVATE :: extrap__fill_value_wrapper !< extrapolate value over detected point 53 PRIVATE :: extrap__fill_value !< extrapolate value over detected point 54 PRIVATE :: extrap__3D 55 PRIVATE :: extrap_deriv_3D 56 PRIVATE :: extrap__3D_min_error_coef 57 PRIVATE :: extrap__3D_min_error_fill 58 PRIVATE :: extrap__3D_dist_weight_coef 59 PRIVATE :: extrap__3D_dist_weight_fill 60 61 INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration 97 PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) 98 PUBLIC :: extrap_del_extrabands !< delete extraband of the variable 99 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 ! 109 62 110 INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate 63 111 INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation 64 112 65 113 INTERFACE extrap_detect 66 MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated114 MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated 67 115 END INTERFACE extrap_detect 68 116 … … 74 122 !------------------------------------------------------------------- 75 123 !> @brief 76 !> This function detected point to be extrapolated .124 !> This function detected point to be extrapolated, given variable structure. 77 125 !> 78 126 !> @details 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/> 79 131 !> 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 !> 80 138 !> @author J.Paul 81 !> - Nov, 2013- Initial Version 139 !> @date November, 2013 - Initial Version 140 !> @date June, 2015 141 !> - do not use level to select points to be extrapolated 82 142 ! 83 !> @param[in] td_var : variable to extrapolate 84 !> @param[in] id_iext : number of points to be extrapolated in i-direction 85 !> @param[in] id_jext : number of points to be extrapolated in j-direction 86 !> @param[in] id_kext : number of points to be extrapolated in k-direction 87 !> @return 88 !------------------------------------------------------------------- 89 !> @code 90 FUNCTION extrap__detect( td_var0, td_level1, & 91 & id_offset, id_rho, id_ext ) 143 !> @param[in] td_var0 coarse grid variable to extrapolate 144 !> @return array of point to be extrapolated 145 !------------------------------------------------------------------- 146 FUNCTION extrap__detect( td_var0 ) 92 147 IMPLICIT NONE 93 148 ! Argument 94 TYPE(TVAR) , INTENT(INOUT) :: td_var0 95 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1 96 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 97 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 98 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext 149 TYPE(TVAR) , INTENT(IN ) :: td_var0 99 150 100 151 ! function … … 104 155 105 156 ! local variable 106 CHARACTER(LEN=lc) :: cl_level107 108 INTEGER(i4) :: il_varid109 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect110 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset111 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1112 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1_G0113 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_extra114 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ext115 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_rho116 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_dim0117 118 TYPE(TVAR) :: tl_var1119 120 157 ! loop indices 121 158 INTEGER(i4) :: ji0 122 159 INTEGER(i4) :: jj0 123 160 INTEGER(i4) :: jk0 124 INTEGER(i4) :: ji1125 INTEGER(i4) :: jj1126 INTEGER(i4) :: ji1m127 INTEGER(i4) :: jj1m128 INTEGER(i4) :: ji1p129 INTEGER(i4) :: jj1p130 161 !---------------------------------------------------------------- 131 162 132 ! init 133 extrap__detect(:,:,:)=0 134 135 ALLOCATE( il_dim0(3) ) 136 il_dim0(:)=td_var0%t_dim(1:3)%i_len 137 138 ! optional argument 139 ALLOCATE( il_rho(ip_maxdim) ) 140 il_rho(:)=1 141 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 142 143 ALLOCATE( il_offset(ip_maxdim,2) ) 144 il_offset(:,:)=0 145 il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) 146 il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) 147 IF( PRESENT(id_offset) )THEN 148 il_offset(1:SIZE(id_offset(:,:),DIM=1),& 149 & 1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:) 150 ENDIF 151 152 ALLOCATE( il_ext(ip_maxdim) ) 153 il_ext(:)=im_minext 154 IF( PRESENT(id_ext) ) il_ext(1:SIZE(id_ext(:)))=id_ext(:) 155 156 ALLOCATE( il_detect(il_dim0(1),& 157 & il_dim0(2),& 158 & il_dim0(3)) ) 159 il_detect(:,:,:)=0 160 161 ! select point already inform 162 WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1 163 164 IF( PRESENT(td_level1) )THEN 165 SELECT CASE(TRIM(td_var0%c_point)) 166 CASE DEFAULT !'T' 167 cl_level='tlevel' 168 CASE('U') 169 cl_level='ulevel' 170 CASE('V') 171 cl_level='vlevel' 172 CASE('F') 173 cl_level='flevel' 174 END SELECT 175 il_varid=var_get_id(td_level1(:),TRIM(cl_level)) 176 IF( il_varid == 0 )THEN 177 CALL logger_error("EXTRAP DETECT: can not compute point to be "//& 178 & "extrapolated for variable "//TRIM(td_var0%c_name)//& 179 & ". can not find "//& 180 & "level for variable point "//TRIM(TRIM(td_var0%c_point))) 181 ELSE 182 print *,'read ',trim(cl_level) 183 tl_var1=td_level1(il_varid) 184 ENDIF 185 186 ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) ) 187 IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN 188 189 ! variable to be extrapolated use same resolution than level 190 il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4) 191 192 ELSE 193 ! variable to be extrapolated do not use same resolution than level 194 ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, & 195 & tl_var1%t_dim(2)%i_len) ) 196 ! match fine grid vertical level with coarse grid 197 il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K) 198 199 ALLOCATE( il_extra(ig_ndim,2) ) 200 ! coarsening fine grid level 201 il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)/2._dp) 202 il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)/2._dp) 203 204 il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)/2._dp) 205 il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)/2._dp) 206 207 DO jj0=1,td_var0%t_dim(2)%i_len 208 209 jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1) 210 211 jj1m=MAX( jj1-il_extra(jp_J,1), 1 ) 212 jj1p=MIN( jj1+il_extra(jp_J,2), & 213 & tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) ) 214 215 DO ji0=1,td_var0%t_dim(1)%i_len 216 217 ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1) 218 219 ji1m=MAX( ji1-il_extra(jp_I,1), 1 ) 220 ji1p=MIN( ji1+il_extra(jp_I,2), & 221 & tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) ) 222 223 il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p)) 224 ENDDO 163 ! force to extrapolated all points 164 extrap__detect(:,:,:)=1 165 166 ! do not compute grid point already inform 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 170 IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill )THEN 171 extrap__detect(ji0,jj0,jk0)=0 172 ENDIF 225 173 ENDDO 226 227 !il_level1_G0(:,:)=0228 !DO jj1=1,tl_var1%t_dim(2)%i_len229 230 ! jj0=INT(REAL((jj1+il_offset(jp_J,1)-1)-1,dp)/REAL(il_rho(jp_J),dp)) +1231 232 ! DO ji1=1,tl_var1%t_dim(1)%i_len233 234 ! ji0=INT(REAL((ji1+il_offset(jp_I,1)-1)-1,dp)/REAL(il_rho(jp_I),dp)) +1235 236 ! il_level1_G0(ji0,jj0)=MAX(il_level1_G0(ji0,jj0),il_level1(ji1,jj1))237 !238 ! ENDDO239 !ENDDO240 241 ! clean242 DEALLOCATE( il_extra )243 DEALLOCATE( il_level1 )244 245 ENDIF246 247 ! look for sea point248 !il_detect(:,:,1)=0249 DO jk0=1,td_var0%t_dim(3)%i_len250 WHERE( il_level1_G0(:,:) >= jk0)251 il_detect(:,:,jk0)=1252 END WHERE253 !il_detect(:,:,jk0)=il_level1_G0(:,:)254 !WHERE( td_var0%d_value(:,:,jk0,1) /= td_var0%d_fill )255 ! il_detect(:,:,1)=jk0-1256 !END WHERE257 174 ENDDO 258 259 ! clean 260 DEALLOCATE( il_level1_G0 ) 261 262 ENDIF 263 264 ! clean 265 DEALLOCATE( il_offset ) 266 267 !! select extra point depending on interpolation method 268 !! compute point near grid point already inform 269 !FORALL( ji0=1:il_dim0(1), & 270 !& jj0=1:il_dim0(2), & 271 !& jk0=1:il_dim0(3), & 272 !& il_detect(ji0,jj0,jk0) == 1 ) 273 274 ! il_detect(MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),& 275 ! & MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),& 276 ! & MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) )=1 277 278 !END FORALL 279 280 ! do not compute grid point already inform 281 WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0 282 283 ! save result 284 extrap__detect(:,:,:)=il_detect(:,:,:) 285 286 ! clean 287 DEALLOCATE( il_dim0 ) 288 DEALLOCATE( il_ext ) 289 DEALLOCATE( il_detect ) 175 ENDDO 290 176 291 177 END FUNCTION extrap__detect 292 !> @endcode293 178 !------------------------------------------------------------------- 294 179 !> @brief 295 !> This function detected point to be extrapolated. 180 !> This function sort variable to be extrapolated, depending on number of 181 !> dimentsion, then detected point to be extrapolated. 296 182 !> 297 !> @details 183 !> @author J.Paul 184 !> @date November, 2013 - Initial Version 185 !> @date June, 2015 186 !> - select all land points for extrapolation 298 187 !> 299 !> @author J.Paul 300 !> - Nov, 2013- Initial Version 301 ! 302 !> @param[in] td_var : variable to extrapolate 303 !> @param[in] id_iext : number of points to be extrapolated in i-direction 304 !> @param[in] id_jext : number of points to be extrapolated in j-direction 305 !> @param[in] id_kext : number of points to be extrapolated in k-direction 306 !> @return 307 !------------------------------------------------------------------- 308 !> @code 309 FUNCTION extrap__detect_wrapper( td_var, td_level, & 310 & id_offset, id_rho, id_ext ) 188 !> @param[in] td_var coarse grid variable to extrapolate 189 !> @return 3D array of point to be extrapolated 190 !------------------------------------------------------------------- 191 FUNCTION extrap__detect_wrapper( td_var ) 311 192 312 193 IMPLICIT NONE 313 194 ! Argument 314 TYPE(TVAR) , INTENT(INOUT) :: td_var 315 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level 316 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 317 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 318 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext 195 TYPE(TVAR) , INTENT(IN ) :: td_var 319 196 320 197 ! function … … 335 212 ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN 336 213 337 ! detect point to be interpolated on I-J-K214 ! detect point to be extrapolated on I-J-K 338 215 CALL logger_debug(" EXTRAP DETECT: detect point "//& 339 216 & " for variable "//TRIM(td_var%c_name) ) 340 217 341 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, & 342 & id_offset, & 343 & id_rho, & 344 & id_ext ) 218 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var ) 345 219 346 220 ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 347 221 348 ! detect point to be interpolated on I-J222 ! detect point to be extrapolated on I-J 349 223 CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//& 350 224 & " for variable "//TRIM(td_var%c_name) ) 351 225 352 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,& 353 & id_offset, & 354 & id_rho, & 355 & id_ext ) 226 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var ) 356 227 357 228 ELSE IF( td_var%t_dim(3)%l_use )THEN 358 229 359 ! detect point to be interpolated on K230 ! detect point to be extrapolated on K 360 231 CALL logger_debug(" EXTRAP DETECT: detect vertical point "//& 361 232 & " for variable "//TRIM(td_var%c_name) ) 362 233 363 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, & 364 & id_offset, & 365 & id_rho, & 366 & id_ext ) 234 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var ) 367 235 368 236 ENDIF … … 373 241 374 242 END FUNCTION extrap__detect_wrapper 375 !> @endcode376 ! !-------------------------------------------------------------------377 ! !> @brief378 ! !> This function detected point to be extrapolated.379 ! !>380 ! !> @details381 ! !>382 ! !> @author J.Paul383 ! !> - Nov, 2013- Initial Version384 ! !385 ! !> @param[in] td_var : variable to extrapolate386 ! !> @param[in] id_iext : number of points to be extrapolated in i-direction387 ! !> @param[in] id_jext : number of points to be extrapolated in j-direction388 ! !> @param[in] id_kext : number of points to be extrapolated in k-direction389 ! !> @return390 ! !391 ! !> @todo392 ! !-------------------------------------------------------------------393 ! !> @code394 ! FUNCTION extrap__detect(td_var, id_iext, id_jext, id_kext)395 ! IMPLICIT NONE396 ! ! Argument397 ! TYPE(TVAR) , INTENT(INOUT) :: td_var398 ! !INTEGER(i4), DIMENSION(:,:,:), INTENT(OUT ) :: id_detect399 ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext400 ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext401 ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext402 !403 ! ! function404 ! INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &405 ! & td_var%t_dim(2)%i_len, &406 ! & td_var%t_dim(3)%i_len) :: extrap__detect407 !408 ! ! local variable409 ! INTEGER(i4) :: il_iext410 ! INTEGER(i4) :: il_jext411 ! INTEGER(i4) :: il_kext412 !413 ! INTEGER(i4), DIMENSION(3) :: il_dim414 !415 ! ! loop indices416 ! INTEGER(i4) :: ji417 ! INTEGER(i4) :: jj418 ! INTEGER(i4) :: jk419 ! !----------------------------------------------------------------420 !421 ! ! optional argument422 ! il_iext=im_minext423 ! IF( PRESENT(id_iext) ) il_iext=id_iext424 ! il_jext=im_minext425 ! IF( PRESENT(id_jext) ) il_jext=id_jext426 ! il_kext=im_minext427 ! IF( PRESENT(id_kext) ) il_kext=id_kext428 !429 ! ! init430 ! extrap__detect(:,:,:)=0431 !432 ! il_dim(:)=td_var%t_dim(1:3)%i_len433 !434 ! ! compute point near grid point already inform435 ! FORALL( ji=1:il_dim(1), &436 ! & jj=1:il_dim(2), &437 ! & jk=1:il_dim(3), &438 ! & td_var%d_value(ji,jj,jk,1) /= td_var%d_fill )439 !440 ! extrap__detect(MAX(1,ji-il_iext):MIN(ji+il_iext,il_dim(1)),&441 ! & MAX(1,jj-il_jext):MIN(jj+il_jext,il_dim(2)),&442 ! & MAX(1,jk-il_kext):MIN(jk+il_kext,il_dim(3)) )=1443 !444 !445 ! END FORALL446 !447 ! ! do not compute grid point already inform448 ! WHERE( td_var%d_value(:,:,:,1) /= td_var%d_fill ) extrap__detect(:,:,:)=0449 !450 ! END FUNCTION extrap__detect451 ! !> @endcode452 243 !------------------------------------------------------------------- 453 244 !> @brief 454 !> This subroutine 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. 455 248 !> 456 !> @details 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 457 257 !> 458 258 !> @author J.Paul 459 !> - Nov, 2013- Initial Version 259 !> @date November, 2013 - Initial Version 260 !> @date June, 2015 261 !> - select all land points for extrapolation 460 262 ! 461 !> @param[inout] td_var : variable structure 462 !> @param[in] id_iext : number of points to be extrapolated in i-direction 463 !> @param[in] id_jext : number of points to be extrapolated in j-direction 464 !> @param[in] id_kext : number of points to be extrapolated in k-direction 465 !> @param[in] id_extend : radius of the box used to compute extrapolation 466 !> @param[in] id_maxiter : maximum nuber of iteration 467 !------------------------------------------------------------------- 468 !> @code 469 SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & 470 & id_offset, & 471 & id_rho, & 472 & id_iext, id_jext, id_kext, & 473 & id_radius, id_maxiter ) 263 !> @param[inout] td_var variable structure 264 !> @param[in] id_radius radius of the halo used to compute extrapolation 265 !------------------------------------------------------------------- 266 SUBROUTINE extrap__fill_value_wrapper( td_var, & 267 & id_radius ) 474 268 IMPLICIT NONE 475 269 ! Argument 476 270 TYPE(TVAR) , INTENT(INOUT) :: td_var 477 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level478 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset479 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho480 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext481 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext482 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext483 271 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius 484 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_maxiter485 272 486 273 ! local variable 487 INTEGER(i4) :: il_iext488 INTEGER(i4) :: il_jext489 INTEGER(i4) :: il_kext490 274 INTEGER(i4) :: il_radius 491 INTEGER(i4) :: il_maxiter492 275 493 276 CHARACTER(LEN=lc) :: cl_method … … 496 279 !---------------------------------------------------------------- 497 280 IF( .NOT. ASSOCIATED(td_var%d_value) )THEN 498 CALL logger_error("EXTRAP FILL VALUE: no table ofvalue "//&281 CALL logger_error("EXTRAP FILL VALUE: no value "//& 499 282 & "associted to variable "//TRIM(td_var%c_name) ) 500 283 ELSE … … 510 293 END SELECT 511 294 512 il_iext=im_minext 513 IF( PRESENT(id_iext) ) il_iext=id_iext 514 il_jext=im_minext 515 IF( PRESENT(id_jext) ) il_jext=id_jext 516 il_kext=0 517 IF( PRESENT(id_kext) ) il_kext=id_kext 518 519 IF( TRIM(td_var%c_interp(1)) == 'cubic')THEN 520 IF( il_iext > 0 .AND. il_iext < im_mincubic ) il_iext=im_mincubic 521 IF( il_jext > 0 .AND. il_jext < im_mincubic ) il_jext=im_mincubic 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 299 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 300 & " radius of the box used to compute extrapolation "//& 301 & "("//TRIM(fct_str(il_radius))//")") 522 302 ENDIF 523 303 524 IF( il_iext < 0 )THEN 525 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 526 & " number of points to be extrapolated in i-direction "//& 527 & "("//TRIM(fct_str(il_iext))//")") 528 ENDIF 529 530 IF( il_jext < 0 )THEN 531 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 532 & " number of points to be extrapolated in j-direction "//& 533 & "("//TRIM(fct_str(il_jext))//")") 534 ENDIF 535 536 IF( il_kext < 0 )THEN 537 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 538 & " number of points to be extrapolated in k-direction "//& 539 & "("//TRIM(fct_str(il_kext))//")") 540 ENDIF 541 542 IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & 543 & (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. & 544 & (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) .OR. & 545 & PRESENT(td_level) )THEN 546 547 ! number of point use to compute box 548 il_radius=1 549 IF( PRESENT(id_radius) ) il_radius=id_radius 550 IF( il_radius < 0 )THEN 551 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 552 & " radius of the box used to compute extrapolation "//& 553 & "("//TRIM(fct_str(il_radius))//")") 554 ENDIF 555 556 ! maximum number of iteration 557 il_maxiter=im_maxiter 558 IF( PRESENT(id_maxiter) ) il_maxiter=id_maxiter 559 IF( il_maxiter < 0 )THEN 560 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 561 & " maximum nuber of iteration "//& 562 & "("//TRIM(fct_str(il_maxiter))//")") 563 ENDIF 564 565 CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& 566 & " using "//TRIM(cl_method)//" method." ) 567 568 CALL extrap__fill_value( td_var, cl_method, & 569 & il_iext, il_jext, il_kext, & 570 & il_radius, il_maxiter, & 571 & td_level, & 572 & id_offset, id_rho ) 573 574 ENDIF 304 CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& 305 & " using "//TRIM(cl_method)//" method." ) 306 307 CALL extrap__fill_value( td_var, cl_method, & 308 & il_radius ) 575 309 576 310 ENDIF 577 311 578 312 END SUBROUTINE extrap__fill_value_wrapper 579 !> @endcode580 313 !------------------------------------------------------------------- 581 314 !> @brief 582 !> This subroutine 315 !> This subroutine compute point to be extrapolated, then extrapolate point. 583 316 !> 584 317 !> @details 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) 585 322 !> 586 323 !> @author J.Paul 587 !> - Nov, 2013- Initial Version 324 !> @date November, 2013 - Initial Version 325 !> @date June, 2015 326 !> - select all land points for extrapolation 588 327 ! 589 !> @param[inout] td_var : variable structure 590 !> @param[in] cd_method : extrapolation method 591 !> @param[in] id_iext : number of points to be extrapolated in i-direction 592 !> @param[in] id_jext : number of points to be extrapolated in j-direction 593 !> @param[in] id_kext : number of points to be extrapolated in k-direction 594 !> @param[in] id_radius : radius of the halo used to compute extrapolation 595 !> @param[in] id_maxiter : maximum nuber of iteration 596 !------------------------------------------------------------------- 597 !> @code 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 331 !------------------------------------------------------------------- 598 332 SUBROUTINE extrap__fill_value( td_var, cd_method, & 599 & id_iext, id_jext, id_kext, & 600 & id_radius, id_maxiter, & 601 & td_level, & 602 & id_offset, & 603 & id_rho ) 333 & id_radius ) 604 334 IMPLICIT NONE 605 335 ! Argument 606 336 TYPE(TVAR) , INTENT(INOUT) :: td_var 607 337 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 608 INTEGER(i4) , INTENT(IN ) :: id_iext609 INTEGER(i4) , INTENT(IN ) :: id_jext610 INTEGER(i4) , INTENT(IN ) :: id_kext611 338 INTEGER(i4) , INTENT(IN ) :: id_radius 612 INTEGER(i4) , INTENT(IN ) :: id_maxiter613 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level614 INTEGER(i4) , DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset615 INTEGER(i4) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho616 339 617 340 ! local variable … … 619 342 620 343 INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect 621 INTEGER(i4) :: il_radius622 INTEGER(i4) :: il_iter623 344 624 345 TYPE(TATT) :: tl_att 625 346 626 347 ! loop indices 627 INTEGER(i4) :: jl628 348 !---------------------------------------------------------------- 629 349 … … 633 353 & td_var%t_dim(3)%i_len) ) 634 354 635 il_detect(:,:,:) = extrap_detect( td_var, td_level, & 636 & id_offset, & 637 & id_rho, & 638 & id_ext=(/id_iext, id_jext, id_kext/) ) 355 il_detect(:,:,:) = extrap_detect( td_var ) 639 356 640 357 !2- add attribute to variable … … 643 360 CALL var_move_att(td_var, tl_att) 644 361 645 CALL logger_warn(" EXTRAP FILL: "//& 646 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 647 & " point(s) to extrapolate " ) 648 649 !3- extrapolate 650 DO jl=1,td_var%t_dim(4)%i_len 651 652 ! td_var%d_value(:,:,:,jl)=il_detect(:,:,:) 653 il_iter=1 654 DO WHILE( ANY(il_detect(:,:,:)==1) ) 655 ! change extend value to minimize number of iteration 656 il_radius=id_radius+(il_iter/id_maxiter) 657 658 CALL logger_debug(" EXTRAP FILL VALUE: "//& 362 CALL att_clean(tl_att) 363 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: "//& 659 370 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 660 & " points to extrapolate " ) 661 662 CALL extrap__3D(td_var%d_value(:,:,:,jl), td_var%d_fill, & 663 & il_detect(:,:,:), & 664 & cd_method, il_radius ) 665 666 il_iter=il_iter+1 667 ENDDO 668 669 ENDDO 670 671 IF( SUM(il_detect(:,:,:)) /= 0 )THEN 672 CALL logger_warn(" EXTRAP FILL: still "//& 673 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 674 & " point(s) to extrapolate " ) 371 & " point(s) to extrapolate " ) 372 373 CALL logger_info(" EXTRAP FILL: method "//& 374 & TRIM(cd_method) ) 375 376 !3- extrapolate 377 CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, & 378 & il_detect(:,:,:), & 379 & cd_method, id_radius ) 675 380 ENDIF 676 381 … … 678 383 679 384 END SUBROUTINE extrap__fill_value 680 !> @endcode681 385 !------------------------------------------------------------------- 682 386 !> @brief 683 !> This subroutine 387 !> This subroutine compute point to be extrapolated in 3D array. 684 388 !> 685 389 !> @details 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 686 397 !> 687 398 !> @author J.Paul 688 !> - Nov, 2013- Initial Version 399 !> @date November, 2013 - Initial Version 400 !> @date July, 2015 401 !> - compute coef indices to be used 402 !> - bug fix: force coef indice to 1, for dimension lenth equal to 1 689 403 ! 690 !> @param[inout] dd_value : 3D table of variable to be extrapolated 691 !> @param[in] dd_fill : FillValue of variable 692 !> @param[inout] id_detect : table of point to extrapolate 693 !> @param[in] id_ext : number of point use to compute box 694 !> @param[in] cd_method : extrapolation method 695 !------------------------------------------------------------------- 696 !> @code 404 !> @param[inout] dd_value 3D array of variable to be extrapolated 405 !> @param[in] dd_fill FillValue of variable 406 !> @param[inout] id_detect array of point to extrapolate 407 !> @param[in] cd_method extrapolation method 408 !> @param[in] id_radius radius of the halo used to compute extrapolation 409 !------------------------------------------------------------------- 697 410 SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 698 & cd_method, id_ ext)411 & cd_method, id_radius ) 699 412 IMPLICIT NONE 700 413 ! Argument 701 REAL(dp) , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value 702 REAL(dp) , INTENT(IN ) :: dd_fill 703 INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 704 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 705 INTEGER(i4), INTENT(IN ) :: id_ext 706 707 ! local variable 708 INTEGER(i4) :: il_imin 709 INTEGER(i4) :: il_imax 710 INTEGER(i4) :: il_jmin 711 INTEGER(i4) :: il_jmax 712 INTEGER(i4) :: il_kmin 713 INTEGER(i4) :: il_kmax 714 715 INTEGER(i4), DIMENSION(3) :: il_shape 716 INTEGER(i4), DIMENSION(3) :: il_dim 717 718 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx 719 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy 720 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz 721 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef 722 723 ! loop indices 724 INTEGER(i4) :: ji 725 INTEGER(i4) :: jj 726 INTEGER(i4) :: jk 727 !---------------------------------------------------------------- 728 729 il_shape(:)=SHAPE(dd_value) 730 731 SELECT CASE(TRIM(cd_method)) 732 CASE('min_error') 733 734 ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) 735 ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) 736 ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) 737 738 739 ! compute derivative in i-direction 740 dl_dfdx(:,:,:)=dd_fill 741 IF( il_shape(1) > 1 )THEN 742 dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'I' ) 743 ENDIF 744 745 ! compute derivative in i-direction 746 dl_dfdy(:,:,:)=dd_fill 747 IF( il_shape(2) > 1 )THEN 748 dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'J' ) 749 ENDIF 750 751 ! compute derivative in i-direction 752 dl_dfdz(:,:,:)=dd_fill 753 IF( il_shape(3) > 1 )THEN 754 dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'K' ) 755 ENDIF 756 757 il_dim(1)=2*id_ext+1 758 IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 759 il_dim(2)=2*id_ext+1 760 IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 761 il_dim(3)=2*id_ext+1 762 IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1 763 764 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 765 766 dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & 767 & 1:il_dim(2), & 768 & 1:il_dim(3))) 769 770 DO jk=1,il_shape(3) 771 IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE 772 DO jj=1,il_shape(2) 773 IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE 774 DO ji=1,il_shape(1) 775 776 IF( id_detect(ji,jj,jk) == 1 )THEN 777 778 il_imin=MAX(ji-id_ext,1) 779 il_imax=MIN(ji+id_ext,il_shape(1)) 780 IF( il_dim(1) == 1 )THEN 781 il_imin=ji 782 il_imax=ji 783 ENDIF 784 785 il_jmin=MAX(jj-id_ext,1) 786 il_jmax=MIN(jj+id_ext,il_shape(2)) 787 IF( il_dim(2) == 1 )THEN 788 il_jmin=jj 789 il_jmax=jj 790 ENDIF 791 792 il_kmin=MAX(jk-id_ext,1) 793 il_kmax=MIN(jk+id_ext,il_shape(3)) 794 IF( il_dim(3) == 1 )THEN 795 il_kmin=jk 796 il_kmax=jk 797 ENDIF 798 799 dd_value(ji,jj,jk)=extrap__3D_min_error_fill( & 800 & dd_value( il_imin:il_imax, & 801 & il_jmin:il_jmax, & 802 & il_kmin:il_kmax ), dd_fill, id_ext, & 803 & dl_dfdx( il_imin:il_imax, & 804 & il_jmin:il_jmax, & 805 & il_kmin:il_kmax ), & 806 & dl_dfdy( il_imin:il_imax, & 807 & il_jmin:il_jmax, & 808 & il_kmin:il_kmax ), & 809 & dl_dfdz( il_imin:il_imax, & 810 & il_jmin:il_jmax, & 811 & il_kmin:il_kmax ), & 812 & dl_coef(:,:,:) ) 813 814 IF( dd_value(ji,jj,jk) /= dd_fill )THEN 815 id_detect(ji,jj,jk)= 0 816 ENDIF 817 818 ENDIF 819 820 ENDDO 821 ENDDO 822 ENDDO 823 824 IF( ALLOCATED(dl_dfdx) ) DEALLOCATE( dl_dfdx ) 825 IF( ALLOCATED(dl_dfdy) ) DEALLOCATE( dl_dfdy ) 826 IF( ALLOCATED(dl_dfdz) ) DEALLOCATE( dl_dfdz ) 827 IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) 828 829 CASE DEFAULT ! 'dist_weight' 830 831 il_dim(1)=2*id_ext+1 832 IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 833 il_dim(2)=2*id_ext+1 834 IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 835 il_dim(3)=2*id_ext+1 836 IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1 837 838 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 839 840 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 841 & 1:il_dim(2), & 842 & 1:il_dim(3)) ) 843 844 DO jk=1,il_shape(3) 845 IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE 846 DO jj=1,il_shape(2) 847 IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE 848 DO ji=1,il_shape(1) 849 850 IF( id_detect(ji,jj,jk) == 1 )THEN 851 852 il_imin=MAX(ji-id_ext,1) 853 il_imax=MIN(ji+id_ext,il_shape(1)) 854 IF( il_dim(1) == 1 )THEN 855 il_imin=ji 856 il_imax=ji 857 ENDIF 858 859 il_jmin=MAX(jj-id_ext,1) 860 il_jmax=MIN(jj+id_ext,il_shape(2)) 861 IF( il_dim(2) == 1 )THEN 862 il_jmin=jj 863 il_jmax=jj 864 ENDIF 865 866 il_kmin=MAX(jk-id_ext,1) 867 il_kmax=MIN(jk+id_ext,il_shape(3)) 868 IF( il_dim(3) == 1 )THEN 869 il_kmin=jk 870 il_kmax=jk 871 ENDIF 872 873 dd_value(ji,jj,jk)=extrap__3D_dist_weight_fill( & 874 & dd_value( il_imin:il_imax, & 875 & il_jmin:il_jmax, & 876 & il_kmin:il_kmax ), dd_fill, id_ext, & 877 & dl_coef(:,:,:) ) 878 879 IF( dd_value(ji,jj,jk) /= dd_fill )THEN 880 id_detect(ji,jj,jk)= 0 881 ENDIF 882 883 ENDIF 884 885 ENDDO 886 ENDDO 887 ENDDO 888 889 IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) 890 891 END SELECT 892 893 END SUBROUTINE extrap__3D 894 !> @endcode 895 !------------------------------------------------------------------- 896 !> @brief 897 !> This function compute derivative of a function in i- and 898 !> j-direction 899 !> 900 !> @details 901 !> 902 !> @author J.Paul 903 !> - Nov, 2013- Initial Version 904 ! 905 !> @param[in] dd_value : 1D table of variable to be extrapolated 906 !> @param[in] dd_fill : FillValue of variable 907 !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension 908 !------------------------------------------------------------------- 909 !> @code 910 PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) 911 912 IMPLICIT NONE 913 ! Argument 914 REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value 915 REAL(dp) , INTENT(IN) :: dd_fill 916 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 917 918 ! function 919 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D 920 921 ! local variable 922 INTEGER(i4) :: il_imin 923 INTEGER(i4) :: il_imax 924 INTEGER(i4), DIMENSION(1) :: il_shape 925 926 REAL(dp) :: dl_min 927 REAL(dp) :: dl_max 928 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value 929 930 LOGICAL :: ll_discont 931 932 ! loop indices 933 INTEGER(i4) :: ji 934 935 INTEGER(i4) :: i1 936 INTEGER(i4) :: i2 937 !---------------------------------------------------------------- 938 ! init 939 extrap_deriv_1D(:)=dd_fill 940 941 ll_discont=.FALSE. 942 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 943 944 il_shape(:)=SHAPE(dd_value(:)) 945 946 ALLOCATE( dl_value(3)) 947 948 ! compute derivative in i-direction 949 DO ji=1,il_shape(1) 950 951 il_imin=MAX(ji-1,1) 952 il_imax=MIN(ji+1,il_shape(1)) 953 954 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 955 i1=1 ; i2=3 956 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 957 i1=1 ; i2=2 958 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 959 i1=2 ; i2=3 960 ENDIF 961 962 dl_value(i1:i2)=dd_value(il_imin:il_imax) 963 IF( il_imin == 1 )THEN 964 dl_value(:)=EOSHIFT( dl_value(:), & 965 & DIM=1, & 966 & SHIFT=-1, & 967 & BOUNDARY=dl_value(1) ) 968 ENDIF 969 IF( il_imax == il_shape(1) )THEN 970 dl_value(:)=EOSHIFT( dl_value(:), & 971 & DIM=1, & 972 & SHIFT=1, & 973 & BOUNDARY=dl_value(3)) 974 ENDIF 975 976 IF( ll_discont )THEN 977 dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill ) 978 dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill ) 979 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 980 WHERE( dl_value(:) < 0._dp ) 981 dl_value(:) = dl_value(:)+360._dp 982 END WHERE 983 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 984 WHERE( dl_value(:) > 180._dp ) 985 dl_value(:) = dl_value(:)-180._dp 986 END WHERE 987 ENDIF 988 ENDIF 989 990 IF( dl_value( 2) /= dd_fill .AND. & ! ji 991 & dl_value( 3) /= dd_fill .AND. & ! ji+1 992 & dl_value( 1) /= dd_fill )THEN ! ji-1 993 994 extrap_deriv_1D(ji)=& 995 & ( dl_value(3) - dl_value(1) ) / & 996 & REAL( il_imax-il_imin ,dp) 997 998 ENDIF 999 1000 ENDDO 1001 1002 DEALLOCATE( dl_value ) 1003 1004 END FUNCTION extrap_deriv_1D 1005 !> @endcode 1006 !------------------------------------------------------------------- 1007 !> @brief 1008 !> This function compute derivative of a function in i- and 1009 !> j-direction 1010 !> 1011 !> @details 1012 !> 1013 !> @author J.Paul 1014 !> - Nov, 2013- Initial Version 1015 ! 1016 !> @param[in] dd_value : 2D table of variable to be extrapolated 1017 !> @param[in] dd_fill : FillValue of variable 1018 !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension 1019 !------------------------------------------------------------------- 1020 !> @code 1021 FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) 1022 1023 IMPLICIT NONE 1024 ! Argument 1025 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value 1026 REAL(dp) , INTENT(IN) :: dd_fill 1027 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim 1028 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 1029 1030 ! function 1031 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 1032 & SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D 1033 1034 ! local variable 1035 INTEGER(i4) :: il_imin 1036 INTEGER(i4) :: il_imax 1037 INTEGER(i4) :: il_jmin 1038 INTEGER(i4) :: il_jmax 1039 INTEGER(i4), DIMENSION(2) :: il_shape 1040 1041 REAL(dp) :: dl_min 1042 REAL(dp) :: dl_max 1043 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value 1044 1045 LOGICAL :: ll_discont 1046 1047 ! loop indices 1048 INTEGER(i4) :: ji 1049 INTEGER(i4) :: jj 1050 1051 INTEGER(i4) :: i1 1052 INTEGER(i4) :: i2 1053 1054 INTEGER(i4) :: j1 1055 INTEGER(i4) :: j2 1056 !---------------------------------------------------------------- 1057 ! init 1058 extrap_deriv_2D(:,:)=dd_fill 1059 1060 ll_discont=.FALSE. 1061 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 1062 1063 il_shape(:)=SHAPE(dd_value(:,:)) 1064 1065 SELECT CASE(TRIM(fct_upper(cd_dim))) 1066 1067 CASE('I') 1068 1069 ALLOCATE( dl_value(3,il_shape(2)) ) 1070 ! compute derivative in i-direction 1071 DO ji=1,il_shape(1) 1072 1073 ! init 1074 dl_value(:,:)=dd_fill 1075 1076 il_imin=MAX(ji-1,1) 1077 il_imax=MIN(ji+1,il_shape(1)) 1078 1079 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 1080 i1=1 ; i2=3 1081 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 1082 i1=1 ; i2=2 1083 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 1084 i1=2 ; i2=3 1085 ENDIF 1086 1087 dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:) 1088 IF( il_imin == 1 )THEN 1089 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1090 & DIM=1, & 1091 & SHIFT=-1, & 1092 & BOUNDARY=dl_value(1,:) ) 1093 ENDIF 1094 IF( il_imax == il_shape(1) )THEN 1095 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1096 & DIM=1, & 1097 & SHIFT=1, & 1098 & BOUNDARY=dl_value(3,:)) 1099 ENDIF 1100 1101 IF( ll_discont )THEN 1102 dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1103 dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1104 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1105 WHERE( dl_value(:,:) < 0_dp ) 1106 dl_value(:,:) = dl_value(:,:)+360._dp 1107 END WHERE 1108 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1109 WHERE( dl_value(:,:) > 180 ) 1110 dl_value(:,:) = dl_value(:,:)-180._dp 1111 END WHERE 1112 ENDIF 1113 ENDIF 1114 1115 WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji 1116 & dl_value(3,:) /= dd_fill .AND. & ! ji+1 1117 & dl_value(1,:) /= dd_fill ) ! ji-1 1118 1119 extrap_deriv_2D(ji,:)=& 1120 & ( dl_value(3,:) - dl_value(1,:) ) / & 1121 & REAL( il_imax-il_imin,dp) 1122 1123 END WHERE 1124 1125 ENDDO 1126 1127 CASE('J') 1128 1129 ALLOCATE( dl_value(il_shape(1),3) ) 1130 ! compute derivative in j-direction 1131 DO jj=1,il_shape(2) 1132 1133 il_jmin=MAX(jj-1,1) 1134 il_jmax=MIN(jj+1,il_shape(2)) 1135 1136 IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 1137 j1=1 ; j2=3 1138 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 1139 j1=1 ; j2=2 1140 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 1141 j1=2 ; j2=3 1142 ENDIF 1143 1144 dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax) 1145 IF( il_jmin == 1 )THEN 1146 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1147 & DIM=2, & 1148 & SHIFT=-1, & 1149 & BOUNDARY=dl_value(:,1)) 1150 ENDIF 1151 IF( il_jmax == il_shape(2) )THEN 1152 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1153 & DIM=2, & 1154 & SHIFT=1, & 1155 & BOUNDARY=dl_value(:,3)) 1156 ENDIF 1157 1158 IF( ll_discont )THEN 1159 dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1160 dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1161 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1162 WHERE( dl_value(:,:) < 0_dp ) 1163 dl_value(:,:) = dl_value(:,:)+360._dp 1164 END WHERE 1165 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1166 WHERE( dl_value(:,:) > 180 ) 1167 dl_value(:,:) = dl_value(:,:)-180._dp 1168 END WHERE 1169 ENDIF 1170 ENDIF 1171 1172 WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj 1173 & dl_value(:, 3) /= dd_fill .AND. & ! jj+1 1174 & dl_value(:, 1) /= dd_fill ) ! jj-1 1175 1176 extrap_deriv_2D(:,jj)=& 1177 & ( dl_value(:,3) - dl_value(:,1) ) / & 1178 & REAL(il_jmax-il_jmin,dp) 1179 1180 END WHERE 1181 1182 ENDDO 1183 1184 END SELECT 1185 1186 DEALLOCATE( dl_value ) 1187 1188 END FUNCTION extrap_deriv_2D 1189 !> @endcode 1190 !------------------------------------------------------------------- 1191 !> @brief 1192 !> This function compute derivative of a function in i- and 1193 !> j-direction 1194 !> 1195 !> @details 1196 !> 1197 !> @author J.Paul 1198 !> - Nov, 2013- Initial Version 1199 ! 1200 !> @param[inout] dd_value : 3D table of variable to be extrapolated 1201 !> @param[in] dd_fill : FillValue of variable 1202 !> @param[in] cd_dim : compute derivative on first (I) second (J) or third (K) dimension 1203 !------------------------------------------------------------------- 1204 !> @code 1205 PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) 1206 1207 IMPLICIT NONE 1208 ! Argument 1209 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1210 REAL(dp) , INTENT(IN) :: dd_fill 1211 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim 1212 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 1213 1214 ! function 1215 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 1216 & SIZE(dd_value,DIM=2), & 1217 & SIZE(dd_value,DIM=3)) :: extrap_deriv_3D 414 REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 415 REAL(dp) , INTENT(IN ) :: dd_fill 416 INTEGER(i4), DIMENSION(:,:,:) , INTENT(INOUT) :: id_detect 417 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 418 INTEGER(i4), INTENT(IN ) :: id_radius 1218 419 1219 420 ! local variable … … 1224 425 INTEGER(i4) :: il_kmin 1225 426 INTEGER(i4) :: il_kmax 1226 INTEGER(i4), DIMENSION(3) :: il_shape 1227 1228 REAL(dp) :: dl_min 1229 REAL(dp) :: dl_max 1230 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value 1231 1232 LOGICAL :: ll_discont 427 INTEGER(i4) :: il_iter 428 INTEGER(i4) :: il_radius 429 INTEGER(i4) :: il_i1 430 INTEGER(i4) :: il_i2 431 INTEGER(i4) :: il_j1 432 INTEGER(i4) :: il_j2 433 INTEGER(i4) :: il_k1 434 INTEGER(i4) :: il_k2 435 436 INTEGER(i4), DIMENSION(4) :: il_shape 437 INTEGER(i4), DIMENSION(3) :: il_dim 438 439 INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect 440 441 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx 442 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy 443 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz 444 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef 445 446 LOGICAL :: ll_iter 1233 447 1234 448 ! loop indices … … 1236 450 INTEGER(i4) :: jj 1237 451 INTEGER(i4) :: jk 1238 1239 INTEGER(i4) :: i1 1240 INTEGER(i4) :: i2 1241 1242 INTEGER(i4) :: j1 1243 INTEGER(i4) :: j2 1244 1245 INTEGER(i4) :: k1 1246 INTEGER(i4) :: k2 452 INTEGER(i4) :: jl 1247 453 !---------------------------------------------------------------- 1248 ! init 1249 extrap_deriv_3D(:,:,:)=dd_fill 1250 1251 ll_discont=.FALSE. 1252 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 1253 1254 il_shape(:)=SHAPE(dd_value(:,:,:)) 1255 1256 1257 SELECT CASE(TRIM(fct_upper(cd_dim))) 1258 1259 CASE('I') 1260 1261 ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) ) 1262 ! compute derivative in i-direction 1263 DO ji=1,il_shape(1) 454 455 il_shape(:)=SHAPE(dd_value) 456 457 ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) ) 458 459 SELECT CASE(TRIM(cd_method)) 460 CASE('min_error') 461 DO jl=1,il_shape(4) 462 463 ! intitialise table of poitn to be extrapolated 464 il_detect(:,:,:)=id_detect(:,:,:) 465 466 il_iter=1 467 DO WHILE( ANY(il_detect(:,:,:)==1) ) 468 ! change extend value to minimize number of iteration 469 il_radius=id_radius+(il_iter-1) 470 ll_iter=.TRUE. 471 472 ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) 473 ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) 474 ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) 475 476 ! compute derivative in i-direction 477 dl_dfdx(:,:,:)=dd_fill 478 IF( il_shape(1) > 1 )THEN 479 dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 480 & dd_fill, 'I' ) 481 ENDIF 482 483 ! compute derivative in j-direction 484 dl_dfdy(:,:,:)=dd_fill 485 IF( il_shape(2) > 1 )THEN 486 dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 487 & dd_fill, 'J' ) 488 ENDIF 489 490 ! compute derivative in k-direction 491 dl_dfdz(:,:,:)=dd_fill 492 IF( il_shape(3) > 1 )THEN 493 dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 494 & dd_fill, 'K' ) 495 ENDIF 496 497 il_dim(1)=2*il_radius+1 498 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 499 il_dim(2)=2*il_radius+1 500 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 501 il_dim(3)=2*il_radius+1 502 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 503 504 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 505 506 dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & 507 & 1:il_dim(2), & 508 & 1:il_dim(3), & 509 & jl )) 510 511 DO jk=1,il_shape(3) 512 ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 513 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 514 DO jj=1,il_shape(2) 515 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 516 DO ji=1,il_shape(1) 517 518 IF( il_detect(ji,jj,jk) == 1 )THEN 519 520 il_imin=MAX(ji-il_radius,1) 521 il_imax=MIN(ji+il_radius,il_shape(1)) 522 ! coef indices to be used 523 il_i1 = il_radius-(ji-il_imin)+1 524 il_i2 = il_radius+(il_imax-ji)+1 525 IF( il_dim(1) == 1 )THEN 526 il_imin=ji 527 il_imax=ji 528 ! coef indices to be used 529 il_i1 = 1 530 il_i2 = 1 531 ENDIF 532 533 534 il_jmin=MAX(jj-il_radius,1) 535 il_jmax=MIN(jj+il_radius,il_shape(2)) 536 ! coef indices to be used 537 il_j1 = il_radius-(jj-il_jmin)+1 538 il_j2 = il_radius+(il_jmax-jj)+1 539 IF( il_dim(2) == 1 )THEN 540 il_jmin=jj 541 il_jmax=jj 542 ! coef indices to be used 543 il_j1 = 1 544 il_j2 = 1 545 ENDIF 546 547 il_kmin=MAX(jk-il_radius,1) 548 il_kmax=MIN(jk+il_radius,il_shape(3)) 549 ! coef indices to be used 550 il_k1 = il_radius-(jk-il_kmin)+1 551 il_k2 = il_radius+(il_kmax-jk)+1 552 IF( il_dim(3) == 1 )THEN 553 il_kmin=jk 554 il_kmax=jk 555 ! coef indices to be used 556 il_k1 = 1 557 il_k2 = 1 558 ENDIF 559 560 dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 561 & dd_value( il_imin:il_imax, & 562 & il_jmin:il_jmax, & 563 & il_kmin:il_kmax,jl ), dd_fill, il_radius, & 564 & dl_dfdx( il_imin:il_imax, & 565 & il_jmin:il_jmax, & 566 & il_kmin:il_kmax ), & 567 & dl_dfdy( il_imin:il_imax, & 568 & il_jmin:il_jmax, & 569 & il_kmin:il_kmax ), & 570 & dl_dfdz( il_imin:il_imax, & 571 & il_jmin:il_jmax, & 572 & il_kmin:il_kmax ), & 573 & dl_coef(il_i1:il_i2, & 574 & il_j1:il_j2, & 575 & il_k1:il_k2) ) 576 577 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 578 il_detect(ji,jj,jk)= 0 579 ll_iter=.FALSE. 580 ENDIF 581 582 ENDIF 583 584 ENDDO 585 ENDDO 586 ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 587 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 588 DO jj=il_shape(2),1,-1 589 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 590 DO ji=il_shape(1),1,-1 591 592 IF( il_detect(ji,jj,jk) == 1 )THEN 593 594 il_imin=MAX(ji-il_radius,1) 595 il_imax=MIN(ji+il_radius,il_shape(1)) 596 ! coef indices to be used 597 il_i1 = il_radius-(ji-il_imin)+1 598 il_i2 = il_radius+(il_imax-ji)+1 599 IF( il_dim(1) == 1 )THEN 600 il_imin=ji 601 il_imax=ji 602 ! coef indices to be used 603 il_i1 = 1 604 il_i2 = 1 605 ENDIF 606 607 608 il_jmin=MAX(jj-il_radius,1) 609 il_jmax=MIN(jj+il_radius,il_shape(2)) 610 ! coef indices to be used 611 il_j1 = il_radius-(jj-il_jmin)+1 612 il_j2 = il_radius+(il_jmax-jj)+1 613 IF( il_dim(2) == 1 )THEN 614 il_jmin=jj 615 il_jmax=jj 616 ! coef indices to be used 617 il_j1 = 1 618 il_j2 = 1 619 ENDIF 620 621 il_kmin=MAX(jk-il_radius,1) 622 il_kmax=MIN(jk+il_radius,il_shape(3)) 623 ! coef indices to be used 624 il_k1 = il_radius-(jk-il_kmin)+1 625 il_k2 = il_radius+(il_kmax-jk)+1 626 IF( il_dim(3) == 1 )THEN 627 il_kmin=jk 628 il_kmax=jk 629 ! coef indices to be used 630 il_k1 = 1 631 il_k2 = 1 632 ENDIF 633 634 dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 635 & dd_value( il_imin:il_imax, & 636 & il_jmin:il_jmax, & 637 & il_kmin:il_kmax,jl ), dd_fill, il_radius, & 638 & dl_dfdx( il_imin:il_imax, & 639 & il_jmin:il_jmax, & 640 & il_kmin:il_kmax ), & 641 & dl_dfdy( il_imin:il_imax, & 642 & il_jmin:il_jmax, & 643 & il_kmin:il_kmax ), & 644 & dl_dfdz( il_imin:il_imax, & 645 & il_jmin:il_jmax, & 646 & il_kmin:il_kmax ), & 647 & dl_coef(il_i1:il_i2, & 648 & il_j1:il_j2, & 649 & il_k1:il_k2) ) 650 651 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 652 il_detect(ji,jj,jk)= 0 653 ll_iter=.FALSE. 654 ENDIF 655 656 ENDIF 657 658 ENDDO 659 ENDDO 660 ENDDO 661 662 DEALLOCATE( dl_dfdx ) 663 DEALLOCATE( dl_dfdy ) 664 DEALLOCATE( dl_dfdz ) 665 DEALLOCATE( dl_coef ) 666 667 IF( ll_iter ) il_iter=il_iter+1 668 ENDDO 669 ENDDO 670 671 CASE DEFAULT ! 'dist_weight' 672 DO jl=1,il_shape(4) 673 674 ! intitialise table of poitn to be extrapolated 675 il_detect(:,:,:)=id_detect(:,:,:) 676 677 il_iter=1 678 DO WHILE( ANY(il_detect(:,:,:)==1) ) 679 ! change extend value to minimize number of iteration 680 il_radius=id_radius+(il_iter-1) 681 ll_iter=.TRUE. 682 683 il_dim(1)=2*il_radius+1 684 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 685 il_dim(2)=2*il_radius+1 686 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 687 il_dim(3)=2*il_radius+1 688 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 689 690 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 691 692 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1),& 693 & 1:il_dim(2),& 694 & 1:il_dim(3),& 695 & jl ) ) 696 697 DO jk=1,il_shape(3) 698 ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 699 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 700 DO jj=1,il_shape(2) 701 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 702 DO ji=1,il_shape(1) 703 704 IF( il_detect(ji,jj,jk) == 1 )THEN 705 706 il_imin=MAX(ji-il_radius,1) 707 il_imax=MIN(ji+il_radius,il_shape(1)) 708 ! coef indices to be used 709 il_i1 = il_radius-(ji-il_imin)+1 710 il_i2 = il_radius+(il_imax-ji)+1 711 IF( il_dim(1) == 1 )THEN 712 il_imin=ji 713 il_imax=ji 714 ! coef indices to be used 715 il_i1 = 1 716 il_i2 = 1 717 ENDIF 718 719 il_jmin=MAX(jj-il_radius,1) 720 il_jmax=MIN(jj+il_radius,il_shape(2)) 721 ! coef indices to be used 722 il_j1 = il_radius-(jj-il_jmin)+1 723 il_j2 = il_radius+(il_jmax-jj)+1 724 IF( il_dim(2) == 1 )THEN 725 il_jmin=jj 726 il_jmax=jj 727 ! coef indices to be used 728 il_j1 = 1 729 il_j2 = 1 730 ENDIF 731 732 il_kmin=MAX(jk-il_radius,1) 733 il_kmax=MIN(jk+il_radius,il_shape(3)) 734 ! coef indices to be used 735 il_k1 = il_radius-(jk-il_kmin)+1 736 il_k2 = il_radius+(il_kmax-jk)+1 737 IF( il_dim(3) == 1 )THEN 738 il_kmin=jk 739 il_kmax=jk 740 ! coef indices to be used 741 il_k1 = 1 742 il_k2 = 1 743 ENDIF 744 745 dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 746 & dd_value( il_imin:il_imax, & 747 & il_jmin:il_jmax, & 748 & il_kmin:il_kmax, & 749 & jl), dd_fill, il_radius, & 750 & dl_coef(il_i1:il_i2, & 751 & il_j1:il_j2, & 752 & il_k1:il_k2) ) 753 754 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 755 il_detect(ji,jj,jk)= 0 756 ll_iter=.FALSE. 757 ENDIF 758 759 ENDIF 760 761 ENDDO 762 ENDDO 763 ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 764 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 765 DO jj=il_shape(2),1,-1 766 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 767 DO ji=il_shape(1),1,-1 768 769 IF( il_detect(ji,jj,jk) == 1 )THEN 770 771 il_imin=MAX(ji-il_radius,1) 772 il_imax=MIN(ji+il_radius,il_shape(1)) 773 ! coef indices to be used 774 il_i1 = il_radius-(ji-il_imin)+1 775 il_i2 = il_radius+(il_imax-ji)+1 776 IF( il_dim(1) == 1 )THEN 777 il_imin=ji 778 il_imax=ji 779 ! coef indices to be used 780 il_i1 = 1 781 il_i2 = 1 782 ENDIF 783 784 il_jmin=MAX(jj-il_radius,1) 785 il_jmax=MIN(jj+il_radius,il_shape(2)) 786 ! coef indices to be used 787 il_j1 = il_radius-(jj-il_jmin)+1 788 il_j2 = il_radius+(il_jmax-jj)+1 789 IF( il_dim(2) == 1 )THEN 790 il_jmin=jj 791 il_jmax=jj 792 ! coef indices to be used 793 il_j1 = 1 794 il_j2 = 1 795 ENDIF 796 797 il_kmin=MAX(jk-il_radius,1) 798 il_kmax=MIN(jk+il_radius,il_shape(3)) 799 ! coef indices to be used 800 il_k1 = il_radius-(jk-il_kmin)+1 801 il_k2 = il_radius+(il_kmax-jk)+1 802 IF( il_dim(3) == 1 )THEN 803 il_kmin=jk 804 il_kmax=jk 805 ! coef indices to be used 806 il_k1 = 1 807 il_k2 = 1 808 ENDIF 809 810 dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 811 & dd_value( il_imin:il_imax, & 812 & il_jmin:il_jmax, & 813 & il_kmin:il_kmax, & 814 & jl), dd_fill, il_radius, & 815 & dl_coef(il_i1:il_i2, & 816 & il_j1:il_j2, & 817 & il_k1:il_k2) ) 818 819 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 820 il_detect(ji,jj,jk)= 0 821 ll_iter=.FALSE. 822 ENDIF 823 824 ENDIF 825 826 ENDDO 827 ENDDO 828 ENDDO 829 CALL logger_info(" EXTRAP 3D: "//& 830 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 831 & " point(s) to extrapolate " ) 1264 832 1265 il_imin=MAX(ji-1,1) 1266 il_imax=MIN(ji+1,il_shape(1)) 1267 1268 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 1269 i1=1 ; i2=3 1270 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 1271 i1=1 ; i2=2 1272 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 1273 i1=2 ; i2=3 1274 ENDIF 1275 1276 dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:) 1277 IF( il_imin == 1 )THEN 1278 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1279 & DIM=1, & 1280 & SHIFT=-1, & 1281 & BOUNDARY=dl_value(1,:,:) ) 1282 ENDIF 1283 IF( il_imax == il_shape(1) )THEN 1284 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1285 & DIM=1, & 1286 & SHIFT=1, & 1287 & BOUNDARY=dl_value(3,:,:)) 1288 ENDIF 1289 1290 IF( ll_discont )THEN 1291 dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1292 dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1293 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1294 WHERE( dl_value(:,:,:) < 0_dp ) 1295 dl_value(:,:,:) = dl_value(:,:,:)+360._dp 1296 END WHERE 1297 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1298 WHERE( dl_value(:,:,:) > 180 ) 1299 dl_value(:,:,:) = dl_value(:,:,:)-180._dp 1300 END WHERE 1301 ENDIF 1302 ENDIF 1303 1304 WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji 1305 & dl_value(3,:,:) /= dd_fill .AND. & !ji+1 1306 & dl_value(1,:,:) /= dd_fill ) !ji-1 1307 1308 extrap_deriv_3D(ji,:,:)= & 1309 & ( dl_value(3,:,:) - dl_value(1,:,:) ) / & 1310 & REAL( il_imax-il_imin ,dp) 1311 1312 END WHERE 1313 1314 ENDDO 1315 1316 CASE('J') 1317 1318 ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) ) 1319 ! compute derivative in j-direction 1320 DO jj=1,il_shape(2) 1321 1322 il_jmin=MAX(jj-1,1) 1323 il_jmax=MIN(jj+1,il_shape(2)) 1324 1325 IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 1326 j1=1 ; j2=3 1327 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 1328 j1=1 ; j2=2 1329 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 1330 j1=2 ; j2=3 1331 ENDIF 1332 1333 dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:) 1334 IF( il_jmin == 1 )THEN 1335 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1336 & DIM=2, & 1337 & SHIFT=-1, & 1338 & BOUNDARY=dl_value(:,1,:) ) 1339 ENDIF 1340 IF( il_jmax == il_shape(2) )THEN 1341 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1342 & DIM=2, & 1343 & SHIFT=1, & 1344 & BOUNDARY=dl_value(:,3,:)) 1345 ENDIF 1346 1347 IF( ll_discont )THEN 1348 dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1349 dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1350 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1351 WHERE( dl_value(:,:,:) < 0_dp ) 1352 dl_value(:,:,:) = dl_value(:,:,:)+360._dp 1353 END WHERE 1354 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1355 WHERE( dl_value(:,:,:) > 180 ) 1356 dl_value(:,:,:) = dl_value(:,:,:)-180._dp 1357 END WHERE 1358 ENDIF 1359 ENDIF 1360 1361 WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj 1362 & dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1 1363 & dl_value(:, 1,:) /= dd_fill ) ! jj-1 1364 1365 extrap_deriv_3D(:,jj,:)=& 1366 & ( dl_value(:,3,:) - dl_value(:,1,:) ) / & 1367 & REAL( il_jmax - il_jmin ,dp) 1368 1369 END WHERE 1370 1371 ENDDO 1372 1373 CASE('K') 1374 ! compute derivative in k-direction 1375 DO jk=1,il_shape(3) 1376 1377 il_kmin=MAX(jk-1,1) 1378 il_kmax=MIN(jk+1,il_shape(3)) 1379 1380 IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN 1381 k1=1 ; k2=3 1382 ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN 1383 k1=1 ; k2=2 1384 ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN 1385 k1=2 ; k2=3 1386 ENDIF 1387 1388 dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax) 1389 IF( il_kmin == 1 )THEN 1390 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1391 & DIM=3, & 1392 & SHIFT=-1, & 1393 & BOUNDARY=dl_value(:,:,1) ) 1394 ENDIF 1395 IF( il_kmax == il_shape(3) )THEN 1396 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1397 & DIM=3, & 1398 & SHIFT=1, & 1399 & BOUNDARY=dl_value(:,:,3)) 1400 ENDIF 1401 1402 IF( ll_discont )THEN 1403 dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1404 dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1405 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1406 WHERE( dl_value(:,:,:) < 0_dp ) 1407 dl_value(:,:,:) = dl_value(:,:,:)+360._dp 1408 END WHERE 1409 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1410 WHERE( dl_value(:,:,:) > 180 ) 1411 dl_value(:,:,:) = dl_value(:,:,:)-180._dp 1412 END WHERE 1413 ENDIF 1414 ENDIF 1415 1416 WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk 1417 & dl_value(:,:, 3) /= dd_fill .AND. & ! jk+1 1418 & dl_value(:,:, 1) /= dd_fill ) ! jk-1 1419 1420 extrap_deriv_3D(:,:,jk)=& 1421 & ( dl_value(:,:,3) - dl_value(:,:,1) ) / & 1422 & REAL( il_kmax-il_kmin,dp) 1423 1424 END WHERE 1425 1426 ENDDO 1427 833 DEALLOCATE( dl_coef ) 834 IF( ll_iter ) il_iter=il_iter+1 835 ENDDO 836 ENDDO 1428 837 END SELECT 1429 838 1430 DEALLOCATE( dl_value ) 1431 1432 END FUNCTION extrap_deriv_3D 1433 !> @endcode 839 DEALLOCATE( il_detect ) 840 841 END SUBROUTINE extrap__3D 1434 842 !------------------------------------------------------------------- 1435 843 !> @brief 1436 !> This function compute extrapolatd values by calculated minimum error using 1437 !> taylor expansion 844 !> This function compute coefficient for min_error extrapolation. 1438 845 !> 1439 846 !> @details 847 !> coefficients are "grid distance" to the center of the box 848 !> choosed to compute extrapolation. 1440 849 !> 1441 850 !> @author J.Paul 1442 !> - Nov, 2013- Initial Version 851 !> @date November, 2013 - Initial Version 852 !> @date July, 2015 853 !> - decrease weight of third dimension 1443 854 ! 1444 !> @param[in] dd_value : 3D table of variable to be extrapolated 1445 !> @param[in] dd_fill : FillValue of variable 1446 !> @param[in] dd_ideriv : derivative of function in i-direction 1447 !> @param[in] dd_jderiv : derivative of function in j-direction 1448 !> @param[in] dd_kderiv : derivative of function in k-direction 1449 !> @param[in] id_ji : i-direction indices table 1450 !> @param[in] id_jj : j-direction indices table 1451 !> @param[in] id_ii : i-direction indices of the point to extrapolate 1452 !> @param[in] id_ij : j-direction indices of the point to extrapolate 1453 !------------------------------------------------------------------- 1454 !> @code 855 !> @param[in] dd_value 3D array of variable to be extrapolated 856 !> @return 3D array of coefficient for minimum error extrapolation 857 !------------------------------------------------------------------- 1455 858 PURE FUNCTION extrap__3D_min_error_coef( dd_value ) 1456 859 … … 1495 898 1496 899 ! compute distance 900 ! "vertical weight" is lower than horizontal 1497 901 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 1498 902 & (jj-il_jmid)**2 + & 1499 & 903 & 3*(jk-il_kmid)**2 1500 904 1501 905 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 1514 918 1515 919 END FUNCTION extrap__3D_min_error_coef 1516 !> @endcode1517 920 !------------------------------------------------------------------- 1518 921 !> @brief 1519 !> This function compute extrapolatd value sby calculated minimum error using922 !> This function compute extrapolatd value by calculated minimum error using 1520 923 !> taylor expansion 1521 924 !> 1522 !> @details 925 !> @author J.Paul 926 !> @date November, 2013 - Initial Version 1523 927 !> 1524 !> @author J.Paul 1525 !> - Nov, 2013- Initial Version 1526 ! 1527 !> @param[in] dd_value : 3D table of variable to be extrapolated 1528 !> @param[in] dd_fill : FillValue of variable 1529 !> @param[in] dd_dfdx : derivative of function in i-direction 1530 !> @param[in] dd_dfdy : derivative of function in j-direction 1531 !> @param[in] dd_dfdz : derivative of function in k-direction 1532 !> @param[in] dd_coef : 1533 !------------------------------------------------------------------- 1534 !> @code 1535 PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_ext, & 928 !> @param[in] dd_value 3D array of variable to be extrapolated 929 !> @param[in] dd_fill FillValue of variable 930 !> @param[in] id_radius radius of the halo used to compute extrapolation 931 !> @param[in] dd_dfdx derivative of function in i-direction 932 !> @param[in] dd_dfdy derivative of function in j-direction 933 !> @param[in] dd_dfdz derivative of function in k-direction 934 !> @param[in] dd_coef array of coefficient for min_error extrapolation 935 !> @return extrapolatd value 936 !------------------------------------------------------------------- 937 PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, & 1536 938 & dd_dfdx, dd_dfdy, dd_dfdz, & 1537 939 & dd_coef ) … … 1540 942 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1541 943 REAL(dp) , INTENT(IN) :: dd_fill 1542 INTEGER(i4), INTENT(IN) :: id_ ext944 INTEGER(i4), INTENT(IN) :: id_radius 1543 945 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx 1544 946 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy … … 1564 966 extrap__3D_min_error_fill=dd_fill 1565 967 1566 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ ext*2))968 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 1567 969 1568 970 IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN … … 1602 1004 1603 1005 END FUNCTION extrap__3D_min_error_fill 1604 !> @endcode1605 1006 !------------------------------------------------------------------- 1606 1007 !> @brief 1607 !> This function compute extrapolatd values by calculated minimum error using 1608 !> taylor expansion 1008 !> This function compute coefficient for inverse distance weighted method 1609 1009 !> 1610 1010 !> @details 1011 !> coefficients are inverse "grid distance" to the center of the box choosed to compute 1012 !> extrapolation. 1611 1013 !> 1612 1014 !> @author J.Paul 1613 !> - Nov, 2013- Initial Version 1015 !> @date November, 2013 - Initial Version 1016 !> @date July, 2015 1017 !> - decrease weight of third dimension 1614 1018 ! 1615 !> @param[in] dd_value : 3D tableof variable to be extrapolated1616 ! -------------------------------------------------------------------1617 ! > @code1019 !> @param[in] dd_value 3D array of variable to be extrapolated 1020 !> @return 3D array of coefficient for inverse distance weighted extrapolation 1021 !------------------------------------------------------------------- 1618 1022 PURE FUNCTION extrap__3D_dist_weight_coef( dd_value ) 1619 1023 … … 1658 1062 1659 1063 ! compute distance 1064 ! "vertical weight" is lower than horizontal 1660 1065 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 1661 1066 & (jj-il_jmid)**2 + & 1662 & 1067 & 3*(jk-il_kmid)**2 1663 1068 1664 1069 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 1677 1082 1678 1083 END FUNCTION extrap__3D_dist_weight_coef 1679 !> @endcode1680 1084 !------------------------------------------------------------------- 1681 1085 !> @brief 1682 !> This function compute extrapolatd value s by calculated minimum error using1683 !> taylor expansion1086 !> This function compute extrapolatd value using inverse distance weighted 1087 !> method 1684 1088 !> 1685 1089 !> @details 1686 1090 !> 1687 1091 !> @author J.Paul 1688 !> - Nov, 2013- Initial Version1092 !> @date November, 2013 - Initial Version 1689 1093 ! 1690 !> @param[in] dd_value : 3D table of variable to be extrapolated 1691 !> @param[in] dd_fill : FillValue of variable 1692 !> @param[in] dd_coef : 1693 !------------------------------------------------------------------- 1694 !> @code 1695 FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_ext, & 1696 & dd_coef ) 1094 !> @param[in] dd_value 3D array of variable to be extrapolated 1095 !> @param[in] dd_fill FillValue of variable 1096 !> @param[in] id_radius radius of the halo used to compute extrapolation 1097 !> @param[in] dd_coef 3D array of coefficient for inverse distance weighted extrapolation 1098 !> @return extrapolatd value 1099 !------------------------------------------------------------------- 1100 FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, & 1101 & dd_coef ) 1697 1102 IMPLICIT NONE 1698 1103 ! Argument 1699 1104 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1700 1105 REAL(dp) , INTENT(IN) :: dd_fill 1701 INTEGER(i4), INTENT(IN) :: id_ ext1106 INTEGER(i4), INTENT(IN) :: id_radius 1702 1107 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef 1703 1108 … … 1716 1121 INTEGER(i4) :: jj 1717 1122 INTEGER(i4) :: jk 1718 1719 1123 !---------------------------------------------------------------- 1720 1124 … … 1722 1126 extrap__3D_dist_weight_fill=dd_fill 1723 1127 1724 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ ext*2))1128 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 1725 1129 1726 1130 IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN … … 1733 1137 dl_coef(:,:,:)=0 1734 1138 1735 FORALL( ji=1:il_shape(1), & 1736 & jj=1:il_shape(2), & 1737 & jk=1:il_shape(3), & 1738 & dd_value(ji,jj,jk) /= dd_fill ) 1739 1740 ! compute factor 1741 dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) 1742 dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) 1743 1744 END FORALL 1139 DO jk=1,il_shape(3) 1140 DO jj=1,il_shape(2) 1141 DO ji=1,il_shape(1) 1142 1143 IF( dd_value(ji,jj,jk) /= dd_fill )THEN 1144 ! compute factor 1145 dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) 1146 dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) 1147 ENDIF 1148 1149 ENDDO 1150 ENDDO 1151 ENDDO 1152 1745 1153 1746 1154 ! return value 1747 1155 IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN 1748 extrap__3D_dist_weight_fill=SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 1156 extrap__3D_dist_weight_fill = & 1157 & SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 1749 1158 ENDIF 1750 1159 … … 1755 1164 1756 1165 END FUNCTION extrap__3D_dist_weight_fill 1757 !> @endcode1758 1166 !------------------------------------------------------------------- 1759 1167 !> @brief … … 1761 1169 !> extraband of N points at north,south,east and west boundaries. 1762 1170 !> 1171 !> @details 1172 !> optionaly you could specify size of extra bands in i- and j-direction 1173 !> 1763 1174 !> @author J.Paul 1764 !> - 2013-Initial version1175 !> @date November, 2013 - Initial version 1765 1176 ! 1766 !> @param[inout] td_var :variable1767 !> @param[in] id_isize :i-direction size of extra bands (default=im_minext)1768 !> @param[in] id_jsize :j-direction size of extra bands (default=im_minext)1177 !> @param[inout] td_var variable 1178 !> @param[in] id_isize i-direction size of extra bands (default=im_minext) 1179 !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) 1769 1180 !> @todo 1770 1181 !> - invalid special case for grid with north fold 1771 1182 !------------------------------------------------------------------- 1772 !> @code1773 1183 SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize ) 1774 1184 IMPLICIT NONE … … 1821 1231 dl_value(:,:,:,:)=td_var%d_value(:,:,:,:) 1822 1232 1233 1823 1234 td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize 1824 1235 td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize … … 1855 1266 1856 1267 END SUBROUTINE extrap_add_extrabands 1857 !> @endcode1858 1268 !------------------------------------------------------------------- 1859 1269 !> @brief … … 1861 1271 !> of N points at north,south,east and west boundaries. 1862 1272 !> 1273 !> @details 1274 !> optionaly you could specify size of extra bands in i- and j-direction 1275 !> 1863 1276 !> @author J.Paul 1864 !> - 2013-Initial version 1865 ! 1866 !> @param[inout] td_var : variable 1867 !> @param[in] id_isize : i-direction size of extra bands (default=im_minext) 1868 !> @param[in] id_jsize : j-direction size of extra bands (default=im_minext) 1869 !------------------------------------------------------------------- 1870 !> @code 1277 !> @date November, 2013 - Initial version 1278 !> 1279 !> @param[inout] td_var variable 1280 !> @param[in] id_isize i-direction size of extra bands (default=im_minext) 1281 !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) 1282 !------------------------------------------------------------------- 1871 1283 SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize ) 1872 1284 IMPLICIT NONE 1873 1285 ! Argument 1874 TYPE(TVAR) , INTENT(INOUT) 1286 TYPE(TVAR) , INTENT(INOUT) :: td_var 1875 1287 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize 1876 1288 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize … … 1881 1293 INTEGER(i4) :: il_isize 1882 1294 INTEGER(i4) :: il_jsize 1883 1295 1884 1296 INTEGER(i4) :: il_imin 1885 1297 INTEGER(i4) :: il_imax … … 1935 1347 1936 1348 END SUBROUTINE extrap_del_extrabands 1937 !> @endcode1938 ! !-------------------------------------------------------------------1939 ! !> @brief1940 ! !> This function1941 ! !>1942 ! !> @details1943 ! !>1944 ! !> @author J.Paul1945 ! !> - Nov, 2013- Initial Version1946 ! !1947 ! !> @param[in]1948 ! !> @param[out]1949 ! !-------------------------------------------------------------------1950 ! !> @code1951 ! FUNCTION extrap_( )1952 ! IMPLICIT NONE1953 ! ! Argument1954 !1955 ! ! local variable1956 !1957 ! ! loop indices1958 ! !----------------------------------------------------------------1959 ! END FUNCTION extrap_1960 ! !> @endcode1961 ! !-------------------------------------------------------------------1962 ! !> @brief1963 ! !> This subroutine1964 ! !>1965 ! !> @details1966 ! !>1967 ! !> @author J.Paul1968 ! !> - Nov, 2013- Initial Version1969 ! !1970 ! !> @param[in]1971 ! !> @param[out]1972 ! !-------------------------------------------------------------------1973 ! !> @code1974 ! SUBROUTINE extrap_( )1975 ! IMPLICIT NONE1976 ! ! Argument1977 !1978 ! ! local variable1979 !1980 ! ! loop indices1981 ! !----------------------------------------------------------------1982 ! END SUBROUTINE extrap_1983 ! !> @endcode1984 1349 END MODULE extrap
Note: See TracChangeset
for help on using the changeset viewer.
