Changeset 5214 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/extrap.f90
- Timestamp:
- 2015-04-15T17:14:06+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/extrap.f90
r4213 r5214 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:dist_weight', 'varname2:min_error' 22 !> 23 !> to detect point to be extrapolated:<br/> 24 !> @code 25 !> il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext]) 26 !> @endcode 27 !> - il_detect(:,:,:) is 3D array of point to be extrapolated 28 !> - td_var is coarse grid variable to be extrapolated 29 !> - td_level is fine grid array of level (see vgrid_get_level) [optional] 30 !> - id_offset is array of offset between fine and coarse grid [optional] 31 !> - id_rho is array of refinment factor [optional] 32 !> - id_ext is array of number of points to be extrapolated [optional] 33 !> 34 !> to extrapolate variable:<br/> 35 !> @code 36 !> CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter]) 37 !> @endcode 38 !> - td_var is coarse grid variable to be extrapolated 39 !> - td_level is fine grid array of level (see vgrid_get_level) [optional] 40 !> - id_offset is array of offset between fine and coarse grid [optional] 41 !> - id_rho is array of refinment factor [optional] 42 !> - id_iext is number of points to be extrapolated in i-direction [optional] 43 !> - id_jext is number of points to be extrapolated in j-direction [optional] 44 !> - id_kext is number of points to be extrapolated in k-direction [optional] 45 !> - id_radius is radius of the halo used to compute extrapolation [optional] 46 !> - id_maxiter is maximum number of iteration [optional] 47 !> 48 !> to add extraband to the variable (to be extrapolated):<br/> 49 !> @code 50 !> CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] ) 51 !> @endcode 52 !> - td_var is variable structure 53 !> - id_isize : i-direction size of extra bands [optional] 54 !> - id_jsize : j-direction size of extra bands [optional] 55 !> 56 !> to delete extraband of a variable:<br/> 57 !> @code 58 !> CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] ) 59 !> @endcode 60 !> - td_var is variable structure 61 !> - id_isize : i-direction size of extra bands [optional] 62 !> - id_jsize : j-direction size of extra bands [optional] 63 !> 64 !> to compute first derivative of 1D array:<br/> 65 !> @code 66 !> dl_value(:)=extrap_deriv_1D( dd_value(:), dd_fill, [ld_discont] ) 67 !> @endcode 68 !> - dd_value is 1D array of variable 69 !> - dd_fill is FillValue of variable 70 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 71 !> 72 !> to compute first derivative of 2D array:<br/> 73 !> @code 74 !> dl_value(:,:)=extrap_deriv_2D( dd_value(:,:), dd_fill, cd_dim, [ld_discont] ) 75 !> @endcode 76 !> - dd_value is 2D array of variable 77 !> - dd_fill is FillValue of variable 78 !> - cd_dim is character to compute derivative on first (I) or second (J) dimension 79 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 80 !> 81 !> to compute first derivative of 3D array:<br/> 82 !> @code 83 !> dl_value(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim, [ld_discont] ) 84 !> @endcode 85 !> - dd_value is 3D array of variable 86 !> - dd_fill is FillValue of variable 87 !> - cd_dim is character to compute derivative on first (I), second (J), or third (K) dimension 88 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 89 !> 90 !> @warning _FillValue must not be zero (use var_chg_FillValue()) 12 91 !> 13 92 !> @author … … 15 94 ! REVISION HISTORY: 16 95 !> @date Nov, 2013 - Initial Version 96 !> @date September, 2014 97 !> - add header 17 98 !> 18 !> @note WARNING: FillValue must not be zero (use var_chg_FillValue) 99 !> @todo 100 !> - create module for each extrapolation method 19 101 !> 20 102 !> @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 103 !---------------------------------------------------------------------- 25 104 MODULE extrap 26 105 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 106 USE kind ! F90 kind parameter 107 USE phycst ! physical constant 108 USE global ! global variable 109 USE fct ! basic useful function 110 USE date ! date manager 111 USE logger ! log file manager 112 USE att ! attribute manager 113 USE dim ! dimension manager 114 USE var ! variable manager 35 115 36 116 IMPLICIT NONE 37 PRIVATE38 117 ! NOTE_avoid_public_variables_if_possible 39 118 40 119 ! type and variable 120 PRIVATE :: im_maxiter !< default maximum number of iteration 121 PRIVATE :: im_minext !< default minumum number of point to extrapolate 122 PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation 41 123 42 124 ! function and subroutine 43 125 PUBLIC :: extrap_detect !< detected point to be extrapolated 44 126 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 extrapolated52 PRIVATE :: extrap__ fill_value_wrapper !< extrapolate value over detected point53 PRIVATE :: extrap__fill_value !< extrapolate value over detected point54 PRIVATE :: extrap__ 3D55 PRIVATE :: extrap_ deriv_3D56 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 127 PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) 128 PUBLIC :: extrap_del_extrabands !< delete extraband of the variable 129 PUBLIC :: extrap_deriv_1D !< compute first derivative of 1D array 130 PUBLIC :: extrap_deriv_2D !< compute first derivative of 2D array 131 PUBLIC :: extrap_deriv_3D !< compute first derivative of 3D array 132 133 PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper 134 PRIVATE :: extrap__detect ! detected point to be extrapolated 135 PRIVATE :: extrap__fill_value_wrapper ! extrapolate value over detected point wrapper 136 PRIVATE :: extrap__fill_value ! extrapolate value over detected point 137 PRIVATE :: extrap__3D ! 138 PRIVATE :: extrap__3D_min_error_coef ! 139 PRIVATE :: extrap__3D_min_error_fill ! 140 PRIVATE :: extrap__3D_dist_weight_coef ! 141 PRIVATE :: extrap__3D_dist_weight_fill ! 60 142 61 143 INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration … … 64 146 65 147 INTERFACE extrap_detect 66 MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated148 MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated 67 149 END INTERFACE extrap_detect 68 150 … … 74 156 !------------------------------------------------------------------- 75 157 !> @brief 76 !> This function detected point to be extrapolated .158 !> This function detected point to be extrapolated, given variable structure. 77 159 !> 78 160 !> @details 161 !> optionaly, you could sepcify fine grid level, refinment factor (default 1), 162 !> offset between fine and coarse grid (default compute from refinment factor 163 !> as offset=(rho-1)/2), number of point to be extrapolated in each direction 164 !> (default im_minext).<br/> 79 165 !> 166 !> First coarsening fine grid level, if need be, then select point near 167 !> grid point already inform. 168 !> 169 !> @note point to be extrapolated are selected using FillValue, 170 !> so to avoid mistake FillValue should not be zero (use var_chg_FillValue) 171 !> 80 172 !> @author J.Paul 81 !> - Nov , 2013- Initial Version173 !> - November, 2013- Initial Version 82 174 ! 83 !> @param[in] td_var :variable to extrapolate84 !> @param[in] id_iext : number of points to be extrapolated in i-direction85 !> @param[in] id_ jext : number of points to be extrapolated in j-direction86 !> @param[in] id_ kext : number of points to be extrapolated in k-direction87 !> @ return88 ! -------------------------------------------------------------------89 ! > @code175 !> @param[in] td_var0 coarse grid variable to extrapolate 176 !> @param[in] td_level1 fine grid array of level 177 !> @param[in] id_offset array of offset between fine and coarse grid 178 !> @param[in] id_rho array of refinment factor 179 !> @param[in] id_ext array of number of points to be extrapolated 180 !> @return array of point to be extrapolated 181 !------------------------------------------------------------------- 90 182 FUNCTION extrap__detect( td_var0, td_level1, & 91 183 & id_offset, id_rho, id_ext ) 92 184 IMPLICIT NONE 93 185 ! Argument 94 TYPE(TVAR) , INTENT(IN OUT) :: td_var0186 TYPE(TVAR) , INTENT(IN ) :: td_var0 95 187 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1 96 188 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset … … 106 198 CHARACTER(LEN=lc) :: cl_level 107 199 108 INTEGER(i4) :: il_ varid200 INTEGER(i4) :: il_ind 109 201 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect 202 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp 110 203 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset 111 204 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1 … … 143 236 ALLOCATE( il_offset(ip_maxdim,2) ) 144 237 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 238 IF( PRESENT(id_offset) )THEN 148 239 il_offset(1:SIZE(id_offset(:,:),DIM=1),& 149 240 & 1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:) 241 ELSE 242 il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) 243 il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) 150 244 ENDIF 151 245 … … 160 254 161 255 ! select point already inform 162 WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1 256 DO jk0=1,td_var0%t_dim(3)%i_len 257 DO jj0=1,td_var0%t_dim(2)%i_len 258 DO ji0=1,td_var0%t_dim(1)%i_len 259 IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=1 260 ENDDO 261 ENDDO 262 ENDDO 163 263 164 264 IF( PRESENT(td_level1) )THEN 165 265 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'266 CASE DEFAULT !'T' 267 cl_level='tlevel' 268 CASE('U') 269 cl_level='ulevel' 270 CASE('V') 271 cl_level='vlevel' 272 CASE('F') 273 cl_level='flevel' 174 274 END SELECT 175 il_varid=var_get_id(td_level1(:),TRIM(cl_level)) 176 IF( il_varid == 0 )THEN 275 276 il_ind=var_get_index(td_level1(:),TRIM(cl_level)) 277 IF( il_ind == 0 )THEN 177 278 CALL logger_error("EXTRAP DETECT: can not compute point to be "//& 178 279 & "extrapolated for variable "//TRIM(td_var0%c_name)//& … … 180 281 & "level for variable point "//TRIM(TRIM(td_var0%c_point))) 181 282 ELSE 182 print *,'read ',trim(cl_level) 183 tl_var1=td_level1(il_varid) 283 tl_var1=var_copy(td_level1(il_ind)) 284 285 ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) ) 286 IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN 287 288 ! variable to be extrapolated use same resolution than level 289 il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4) 290 291 ELSE 292 ! variable to be extrapolated do not use same resolution than level 293 ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, & 294 & tl_var1%t_dim(2)%i_len) ) 295 ! match fine grid vertical level with coarse grid 296 il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K) 297 298 ALLOCATE( il_extra(ip_maxdim,2) ) 299 ! coarsening fine grid level 300 il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)*0.5_dp) 301 il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5_dp) 302 303 il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)*0.5_dp) 304 il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5_dp) 305 306 DO jj0=1,td_var0%t_dim(2)%i_len 307 308 jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1) 309 310 jj1m=MAX( jj1-il_extra(jp_J,1), 1 ) 311 jj1p=MIN( jj1+il_extra(jp_J,2), & 312 & tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) ) 313 314 DO ji0=1,td_var0%t_dim(1)%i_len 315 316 ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1) 317 318 ji1m=MAX( ji1-il_extra(jp_I,1), 1 ) 319 ji1p=MIN( ji1+il_extra(jp_I,2), & 320 & tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) ) 321 322 il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p)) 323 324 ENDDO 325 ENDDO 326 327 ! clean 328 DEALLOCATE( il_extra ) 329 DEALLOCATE( il_level1 ) 330 331 ENDIF 332 333 ! look for sea point 334 DO jk0=1,td_var0%t_dim(3)%i_len 335 WHERE( il_level1_G0(:,:) >= jk0) 336 il_detect(:,:,jk0)=1 337 END WHERE 338 ENDDO 339 340 ! clean 341 DEALLOCATE( il_level1_G0 ) 342 CALL var_clean(tl_var1) 343 184 344 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)) )THEN188 189 ! variable to be extrapolated use same resolution than level190 il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)191 192 ELSE193 ! variable to be extrapolated do not use same resolution than level194 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 grid197 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 level201 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_len208 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_len216 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 ENDDO225 ENDDO226 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 ENDDO258 259 ! clean260 DEALLOCATE( il_level1_G0 )261 262 345 ENDIF 263 346 … … 265 348 DEALLOCATE( il_offset ) 266 349 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 350 ALLOCATE( il_tmp(il_dim0(1),& 351 & il_dim0(2),& 352 & il_dim0(3)) ) 353 il_tmp(:,:,:)=il_detect(:,:,:) 354 ! select extra point depending on interpolation method 355 ! compute point near grid point already inform 356 DO jk0=1,il_dim0(3) 357 DO jj0=1,il_dim0(2) 358 DO ji0=1,il_dim0(1) 359 360 IF( il_tmp(ji0,jj0,jk0) == 1 )THEN 361 il_detect( & 362 & MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),& 363 & MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),& 364 & MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) & 365 & ) = 1 366 ENDIF 367 368 ENDDO 369 ENDDO 370 ENDDO 279 371 372 ! clean 373 DEALLOCATE( il_tmp ) 374 280 375 ! do not compute grid point already inform 281 WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0 376 DO jk0=1,td_var0%t_dim(3)%i_len 377 DO jj0=1,td_var0%t_dim(2)%i_len 378 DO ji0=1,td_var0%t_dim(1)%i_len 379 IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=0 380 ENDDO 381 ENDDO 382 ENDDO 282 383 283 384 ! save result … … 288 389 DEALLOCATE( il_ext ) 289 390 DEALLOCATE( il_detect ) 391 DEALLOCATE( il_rho ) 290 392 291 393 END FUNCTION extrap__detect 292 !> @endcode293 394 !------------------------------------------------------------------- 294 395 !> @brief 295 !> This function detected point to be extrapolated. 396 !> This function sort variable to be extrapolated, depending on number of 397 !> dimentsion, then detected point to be extrapolated. 296 398 !> 297 !> @details 399 !> @author J.Paul 400 !> - November, 2013- Initial Version 298 401 !> 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 402 !> @param[in] td_var coarse grid variable to extrapolate 403 !> @param[in] td_level fine grid array of level 404 !> @param[in] id_offset array of offset between fine and coarse grid 405 !> @param[in] id_rho array of refinment factor 406 !> @param[in] id_ext array of number of points to be extrapolated 407 !> @return 3D array of point to be extrapolated 408 !------------------------------------------------------------------- 309 409 FUNCTION extrap__detect_wrapper( td_var, td_level, & 310 410 & id_offset, id_rho, id_ext ) … … 312 412 IMPLICIT NONE 313 413 ! Argument 314 TYPE(TVAR) , INTENT(INOUT) :: td_var315 TYPE(TVAR) , DIMENSION(:) 414 TYPE(TVAR) , INTENT(IN ) :: td_var 415 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level 316 416 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 317 417 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho … … 335 435 ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN 336 436 337 ! detect point to be interpolated on I-J-K437 ! detect point to be extrapolated on I-J-K 338 438 CALL logger_debug(" EXTRAP DETECT: detect point "//& 339 439 & " for variable "//TRIM(td_var%c_name) ) … … 346 446 ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 347 447 348 ! detect point to be interpolated on I-J448 ! detect point to be extrapolated on I-J 349 449 CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//& 350 450 & " for variable "//TRIM(td_var%c_name) ) … … 357 457 ELSE IF( td_var%t_dim(3)%l_use )THEN 358 458 359 ! detect point to be interpolated on K459 ! detect point to be extrapolated on K 360 460 CALL logger_debug(" EXTRAP DETECT: detect vertical point "//& 361 461 & " for variable "//TRIM(td_var%c_name) ) … … 373 473 374 474 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 475 !------------------------------------------------------------------- 453 476 !> @brief 454 !> This subroutine 477 !> This subroutine select method to be used for extrapolation. 478 !> If need be, increase number of points to be extrapolated. 479 !> Finally launch extrap__fill_value. 455 480 !> 456 !> @details 481 !> @details 482 !> optionaly, you could specify :<br/> 483 !> - refinment factor (default 1) 484 !> - offset between fine and coarse grid (default compute from refinment factor 485 !> as offset=(rho-1)/2) 486 !> - number of point to be extrapolated in each direction (default im_minext) 487 !> - radius of the halo used to compute extrapolation 488 !> - maximum number of iteration 457 489 !> 458 490 !> @author J.Paul 459 491 !> - Nov, 2013- Initial Version 460 492 ! 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 493 !> @param[inout] td_var variable structure 494 !> @param[in] td_level fine grid array of level 495 !> @param[in] id_offset array of offset between fine and coarse grid 496 !> @param[in] id_rho array of refinment factor 497 !> @param[in] id_iext number of points to be extrapolated in i-direction 498 !> @param[in] id_jext number of points to be extrapolated in j-direction 499 !> @param[in] id_kext number of points to be extrapolated in k-direction 500 !> @param[in] id_radius radius of the halo used to compute extrapolation 501 !> @param[in] id_maxiter maximum number of iteration 502 !------------------------------------------------------------------- 469 503 SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & 470 504 & id_offset, & … … 496 530 !---------------------------------------------------------------- 497 531 IF( .NOT. ASSOCIATED(td_var%d_value) )THEN 498 CALL logger_error("EXTRAP FILL VALUE: no table ofvalue "//&532 CALL logger_error("EXTRAP FILL VALUE: no value "//& 499 533 & "associted to variable "//TRIM(td_var%c_name) ) 500 534 ELSE … … 542 576 IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & 543 577 & (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 578 & (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN 546 579 547 580 ! number of point use to compute box … … 577 610 578 611 END SUBROUTINE extrap__fill_value_wrapper 579 !> @endcode580 612 !------------------------------------------------------------------- 581 613 !> @brief 582 !> This subroutine 614 !> This subroutine compute point to be extrapolated, then extrapolate point. 583 615 !> 584 616 !> @details 617 !> optionaly, you could specify :<br/> 618 !> - refinment factor (default 1) 619 !> - offset between fine and coarse grid (default compute from refinment factor 620 !> as offset=(rho-1)/2) 585 621 !> 586 622 !> @author J.Paul 587 !> - Nov , 2013- Initial Version623 !> - November, 2013- Initial Version 588 624 ! 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 625 !> @param[inout] td_var variable structure 626 !> @param[in] cd_method extrapolation method 627 !> @param[in] id_iext number of points to be extrapolated in i-direction 628 !> @param[in] id_jext number of points to be extrapolated in j-direction 629 !> @param[in] id_kext number of points to be extrapolated in k-direction 630 !> @param[in] id_radius radius of the halo used to compute extrapolation 631 !> @param[in] id_maxiter maximum number of iteration 632 !> @param[in] td_level fine grid array of level 633 !> @param[in] id_offset array of offset between fine and coarse grid 634 !> @param[in] id_rho array of refinment factor 635 !------------------------------------------------------------------- 598 636 SUBROUTINE extrap__fill_value( td_var, cd_method, & 599 637 & id_iext, id_jext, id_kext, & … … 619 657 620 658 INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect 621 INTEGER(i4) :: il_radius622 INTEGER(i4) :: il_iter623 659 624 660 TYPE(TATT) :: tl_att 625 661 626 662 ! loop indices 627 INTEGER(i4) :: jl628 663 !---------------------------------------------------------------- 629 664 … … 637 672 & id_rho, & 638 673 & id_ext=(/id_iext, id_jext, id_kext/) ) 639 640 674 !2- add attribute to variable 641 675 cl_extrap=fct_concat(td_var%c_extrap(:)) … … 643 677 CALL var_move_att(td_var, tl_att) 644 678 645 CALL logger_warn(" EXTRAP FILL: "//& 679 CALL att_clean(tl_att) 680 681 CALL logger_info(" EXTRAP FILL: "//& 646 682 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 647 683 & " point(s) to extrapolate " ) 648 684 649 685 !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: "//& 659 & 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 " ) 675 ENDIF 686 CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, & 687 & il_detect(:,:,:), & 688 & cd_method, id_radius, id_maxiter ) 676 689 677 690 DEALLOCATE(il_detect) 678 691 679 692 END SUBROUTINE extrap__fill_value 680 !> @endcode681 693 !------------------------------------------------------------------- 682 694 !> @brief 683 !> This subroutine 695 !> This subroutine compute point to be extrapolated in 3D array. 684 696 !> 685 697 !> @details 698 !> in case of 'min_error' method:<br/> 699 !> - compute derivative in i-, j- and k- direction 700 !> - compute minimum error coefficient (distance to center of halo) 701 !> - compute extrapolatd values by calculated minimum error using taylor expansion 702 !> in case of 'dist_weight' method:<br/> 703 !> - compute distance weight coefficient (inverse of distance to center of halo) 704 !> - compute extrapolatd values using Inverse Distance Weighting 686 705 !> 687 706 !> @author J.Paul 688 707 !> - Nov, 2013- Initial Version 689 708 ! 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 709 !> @param[inout] dd_value 3D array of variable to be extrapolated 710 !> @param[in] dd_fill FillValue of variable 711 !> @param[inout] id_detect array of point to extrapolate 712 !> @param[in] cd_method extrapolation method 713 !> @param[in] id_radius radius of the halo used to compute extrapolation 714 !------------------------------------------------------------------- 697 715 SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 698 & cd_method, id_ ext)716 & cd_method, id_radius, id_maxiter ) 699 717 IMPLICIT NONE 700 718 ! Argument 701 REAL(dp) , DIMENSION(:,:,: ), INTENT(INOUT) :: dd_value719 REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 702 720 REAL(dp) , INTENT(IN ) :: dd_fill 703 721 INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 704 722 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 705 INTEGER(i4), INTENT(IN ) :: id_ext 723 INTEGER(i4), INTENT(IN ) :: id_radius 724 INTEGER(i4), INTENT(IN ) :: id_maxiter 706 725 707 726 ! local variable … … 712 731 INTEGER(i4) :: il_kmin 713 732 INTEGER(i4) :: il_kmax 714 715 INTEGER(i4), DIMENSION(3) :: il_shape 733 INTEGER(i4) :: il_iter 734 INTEGER(i4) :: il_radius 735 736 INTEGER(i4), DIMENSION(4) :: il_shape 716 737 INTEGER(i4), DIMENSION(3) :: il_dim 738 739 INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect 717 740 718 741 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx … … 725 748 INTEGER(i4) :: jj 726 749 INTEGER(i4) :: jk 750 INTEGER(i4) :: jl 727 751 !---------------------------------------------------------------- 728 752 729 753 il_shape(:)=SHAPE(dd_value) 730 754 755 ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) ) 756 731 757 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 758 CASE('min_error') 759 DO jl=1,il_shape(4) 760 761 ! intitialise table of poitn to be extrapolated 762 il_detect(:,:,:)=id_detect(:,:,:) 763 764 il_iter=1 765 DO WHILE( ANY(il_detect(:,:,:)==1) ) 766 ! change extend value to minimize number of iteration 767 il_radius=id_radius+(il_iter/id_maxiter) 768 769 ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) 770 ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) 771 ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) 772 773 ! compute derivative in i-direction 774 dl_dfdx(:,:,:)=dd_fill 775 IF( il_shape(1) > 1 )THEN 776 dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) 777 ENDIF 778 779 ! compute derivative in j-direction 780 dl_dfdy(:,:,:)=dd_fill 781 IF( il_shape(2) > 1 )THEN 782 dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) 783 ENDIF 784 785 ! compute derivative in k-direction 786 dl_dfdz(:,:,:)=dd_fill 787 IF( il_shape(3) > 1 )THEN 788 dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) 789 ENDIF 790 791 il_dim(1)=2*il_radius+1 792 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 793 il_dim(2)=2*il_radius+1 794 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 795 il_dim(3)=2*il_radius+1 796 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 797 798 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 799 800 dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & 801 & 1:il_dim(2), & 802 & 1:il_dim(3), & 803 & jl )) 804 805 DO jk=1,il_shape(3) 806 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 807 DO jj=1,il_shape(2) 808 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 809 DO ji=1,il_shape(1) 810 811 IF( il_detect(ji,jj,jk) == 1 )THEN 812 813 il_imin=MAX(ji-il_radius,1) 814 il_imax=MIN(ji+il_radius,il_shape(1)) 815 IF( il_dim(1) == 1 )THEN 816 il_imin=ji 817 il_imax=ji 818 ENDIF 819 820 il_jmin=MAX(jj-il_radius,1) 821 il_jmax=MIN(jj+il_radius,il_shape(2)) 822 IF( il_dim(2) == 1 )THEN 823 il_jmin=jj 824 il_jmax=jj 825 ENDIF 826 827 il_kmin=MAX(jk-il_radius,1) 828 il_kmax=MIN(jk+il_radius,il_shape(3)) 829 IF( il_dim(3) == 1 )THEN 830 il_kmin=jk 831 il_kmax=jk 832 ENDIF 833 834 dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 835 & dd_value( il_imin:il_imax, & 836 & il_jmin:il_jmax, & 837 & il_kmin:il_kmax,jl ), dd_fill, il_radius, & 838 & dl_dfdx( il_imin:il_imax, & 839 & il_jmin:il_jmax, & 840 & il_kmin:il_kmax ), & 841 & dl_dfdy( il_imin:il_imax, & 842 & il_jmin:il_jmax, & 843 & il_kmin:il_kmax ), & 844 & dl_dfdz( il_imin:il_imax, & 845 & il_jmin:il_jmax, & 846 & il_kmin:il_kmax ), & 847 & dl_coef(:,:,:) ) 848 849 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 850 il_detect(ji,jj,jk)= 0 851 ENDIF 852 783 853 ENDIF 784 854 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 855 ENDDO 820 856 ENDDO 821 857 ENDDO 858 859 DEALLOCATE( dl_dfdx ) 860 DEALLOCATE( dl_dfdy ) 861 DEALLOCATE( dl_dfdz ) 862 DEALLOCATE( dl_coef ) 863 864 il_iter=il_iter+1 822 865 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 866 ENDDO 867 868 CASE DEFAULT ! 'dist_weight' 869 DO jl=1,il_shape(4) 870 871 ! intitialise table of poitn to be extrapolated 872 il_detect(:,:,:)=id_detect(:,:,:) 873 874 il_iter=1 875 DO WHILE( ANY(il_detect(:,:,:)==1) ) 876 ! change extend value to minimize number of iteration 877 il_radius=id_radius+(il_iter/id_maxiter) 878 879 il_dim(1)=2*il_radius+1 880 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 881 il_dim(2)=2*il_radius+1 882 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 883 il_dim(3)=2*il_radius+1 884 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 885 886 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 887 888 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 889 & 1:il_dim(2), & 890 & 1:il_dim(3), & 891 & jl ) ) 892 893 DO jk=1,il_shape(3) 894 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 895 DO jj=1,il_shape(2) 896 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 897 DO ji=1,il_shape(1) 898 899 IF( il_detect(ji,jj,jk) == 1 )THEN 900 901 il_imin=MAX(ji-il_radius,1) 902 il_imax=MIN(ji+il_radius,il_shape(1)) 903 IF( il_dim(1) == 1 )THEN 904 il_imin=ji 905 il_imax=ji 906 ENDIF 907 908 il_jmin=MAX(jj-il_radius,1) 909 il_jmax=MIN(jj+il_radius,il_shape(2)) 910 IF( il_dim(2) == 1 )THEN 911 il_jmin=jj 912 il_jmax=jj 913 ENDIF 914 915 il_kmin=MAX(jk-il_radius,1) 916 il_kmax=MIN(jk+il_radius,il_shape(3)) 917 IF( il_dim(3) == 1 )THEN 918 il_kmin=jk 919 il_kmax=jk 920 ENDIF 921 922 dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 923 & dd_value( il_imin:il_imax, & 924 & il_jmin:il_jmax, & 925 & il_kmin:il_kmax, & 926 & jl), dd_fill, il_radius, & 927 & dl_coef(:,:,:) ) 928 929 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 930 il_detect(ji,jj,jk)= 0 931 ENDIF 932 857 933 ENDIF 858 934 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 935 ENDDO 885 936 ENDDO 886 937 ENDDO 938 939 DEALLOCATE( dl_coef ) 940 il_iter=il_iter+1 887 941 ENDDO 888 889 IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) 890 942 ENDDO 891 943 END SELECT 892 944 945 DEALLOCATE( il_detect ) 946 893 947 END SUBROUTINE extrap__3D 894 !> @endcode895 948 !------------------------------------------------------------------- 896 949 !> @brief 897 !> This function compute derivative of a function in i- and 898 !> j-direction 950 !> This function compute derivative of 1D array. 899 951 !> 900 952 !> @details 953 !> optionaly you could specify to take into account east west discontinuity 954 !> (-180° 180° or 0° 360° for longitude variable) 901 955 !> 902 956 !> @author J.Paul 903 !> - Nov , 2013- Initial Version957 !> - November, 2013- Initial Version 904 958 ! 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 959 !> @param[in] dd_value 1D array of variable to be extrapolated 960 !> @param[in] dd_fill FillValue of variable 961 !> @param[in] ld_discont logical to take into account east west discontinuity 962 !------------------------------------------------------------------- 910 963 PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) 911 964 … … 1003 1056 1004 1057 END FUNCTION extrap_deriv_1D 1005 !> @endcode1006 1058 !------------------------------------------------------------------- 1007 1059 !> @brief 1008 !> This function compute derivative of a function in i- and 1009 !> j-direction 1010 !> 1060 !> This function compute derivative of 2D array. 1061 !> you have to specify in which direction derivative have to be computed: 1062 !> first (I) or second (J) dimension. 1063 !> 1011 1064 !> @details 1065 !> optionaly you could specify to take into account east west discontinuity 1066 !> (-180° 180° or 0° 360° for longitude variable) 1012 1067 !> 1013 1068 !> @author J.Paul 1014 !> - Nov , 2013- Initial Version1069 !> - November, 2013- Initial Version 1015 1070 ! 1016 !> @param[in] dd_value : 2D tableof variable to be extrapolated1017 !> @param[in] dd_fill :FillValue of variable1018 !> @param[in] cd_dim :compute derivative on first (I) or second (J) dimension1019 ! -------------------------------------------------------------------1020 ! > @code1071 !> @param[in] dd_value 2D array of variable to be extrapolated 1072 !> @param[in] dd_fill FillValue of variable 1073 !> @param[in] cd_dim compute derivative on first (I) or second (J) dimension 1074 !> @param[in] ld_discont logical to take into account east west discontinuity 1075 !------------------------------------------------------------------- 1021 1076 FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) 1022 1077 … … 1123 1178 END WHERE 1124 1179 1125 ENDDO1180 ENDDO 1126 1181 1127 1182 CASE('J') … … 1187 1242 1188 1243 END FUNCTION extrap_deriv_2D 1189 !> @endcode1190 1244 !------------------------------------------------------------------- 1191 1245 !> @brief 1192 !> This function compute derivative of a function in i- and 1193 !> j-direction 1246 !> This function compute derivative of 3D array. 1247 !> you have to specify in which direction derivative have to be computed: 1248 !> first (I), second (J) or third (K) dimension. 1194 1249 !> 1195 1250 !> @details 1251 !> optionaly you could specify to take into account east west discontinuity 1252 !> (-180° 180° or 0° 360° for longitude variable) 1196 1253 !> 1197 1254 !> @author J.Paul 1198 !> - Nov , 2013- Initial Version1255 !> - November, 2013- Initial Version 1199 1256 ! 1200 !> @param[inout] dd_value : 3D tableof variable to be extrapolated1201 !> @param[in] dd_fill :FillValue of variable1202 !> @param[in] cd_dim :compute derivative on first (I) second (J) or third (K) dimension1203 ! -------------------------------------------------------------------1204 ! > @code1257 !> @param[inout] dd_value 3D array of variable to be extrapolated 1258 !> @param[in] dd_fill FillValue of variable 1259 !> @param[in] cd_dim compute derivative on first (I) second (J) or third (K) dimension 1260 !> @param[in] ld_discont logical to take into account east west discontinuity 1261 !------------------------------------------------------------------- 1205 1262 PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) 1206 1263 … … 1431 1488 1432 1489 END FUNCTION extrap_deriv_3D 1433 !> @endcode1434 1490 !------------------------------------------------------------------- 1435 1491 !> @brief 1436 !> This function compute extrapolatd values by calculated minimum error using 1437 !> taylor expansion 1492 !> This function compute coefficient for min_error extrapolation. 1438 1493 !> 1439 1494 !> @details 1495 !> coefficients are "grid distance" to the center of the box choosed to compute 1496 !> extrapolation. 1440 1497 !> 1441 1498 !> @author J.Paul 1442 !> - Nov , 2013- Initial Version1499 !> - November, 2013- Initial Version 1443 1500 ! 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 1501 !> @param[in] dd_value 3D array of variable to be extrapolated 1502 !> @return 3D array of coefficient for minimum error extrapolation 1503 !------------------------------------------------------------------- 1455 1504 PURE FUNCTION extrap__3D_min_error_coef( dd_value ) 1456 1505 … … 1514 1563 1515 1564 END FUNCTION extrap__3D_min_error_coef 1516 !> @endcode1517 1565 !------------------------------------------------------------------- 1518 1566 !> @brief 1519 !> This function compute extrapolatd value sby calculated minimum error using1567 !> This function compute extrapolatd value by calculated minimum error using 1520 1568 !> taylor expansion 1521 1569 !> 1522 !> @details 1570 !> @author J.Paul 1571 !> - November, 2013- Initial Version 1523 1572 !> 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, & 1573 !> @param[in] dd_value 3D array of variable to be extrapolated 1574 !> @param[in] dd_fill FillValue of variable 1575 !> @param[in] id_radius radius of the halo used to compute extrapolation 1576 !> @param[in] dd_dfdx derivative of function in i-direction 1577 !> @param[in] dd_dfdy derivative of function in j-direction 1578 !> @param[in] dd_dfdz derivative of function in k-direction 1579 !> @param[in] dd_coef array of coefficient for min_error extrapolation 1580 !> @return extrapolatd value 1581 !------------------------------------------------------------------- 1582 PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, & 1536 1583 & dd_dfdx, dd_dfdy, dd_dfdz, & 1537 1584 & dd_coef ) … … 1540 1587 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1541 1588 REAL(dp) , INTENT(IN) :: dd_fill 1542 INTEGER(i4), INTENT(IN) :: id_ ext1589 INTEGER(i4), INTENT(IN) :: id_radius 1543 1590 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx 1544 1591 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy … … 1564 1611 extrap__3D_min_error_fill=dd_fill 1565 1612 1566 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ ext*2))1613 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 1567 1614 1568 1615 IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN … … 1602 1649 1603 1650 END FUNCTION extrap__3D_min_error_fill 1604 !> @endcode1605 1651 !------------------------------------------------------------------- 1606 1652 !> @brief 1607 !> This function compute extrapolatd values by calculated minimum error using 1608 !> taylor expansion 1653 !> This function compute coefficient for inverse distance weighted method 1609 1654 !> 1610 1655 !> @details 1656 !> coefficients are inverse "grid distance" to the center of the box choosed to compute 1657 !> extrapolation. 1611 1658 !> 1612 1659 !> @author J.Paul 1613 !> - Nov , 2013- Initial Version1660 !> - November, 2013- Initial Version 1614 1661 ! 1615 !> @param[in] dd_value : 3D tableof variable to be extrapolated1616 ! -------------------------------------------------------------------1617 ! > @code1662 !> @param[in] dd_value 3D array of variable to be extrapolated 1663 !> @return 3D array of coefficient for inverse distance weighted extrapolation 1664 !------------------------------------------------------------------- 1618 1665 PURE FUNCTION extrap__3D_dist_weight_coef( dd_value ) 1619 1666 … … 1677 1724 1678 1725 END FUNCTION extrap__3D_dist_weight_coef 1679 !> @endcode1680 1726 !------------------------------------------------------------------- 1681 1727 !> @brief 1682 !> This function compute extrapolatd value s by calculated minimum error using1683 !> taylor expansion1728 !> This function compute extrapolatd value using inverse distance weighted 1729 !> method 1684 1730 !> 1685 1731 !> @details 1686 1732 !> 1687 1733 !> @author J.Paul 1688 !> - Nov , 2013- Initial Version1734 !> - November, 2013- Initial Version 1689 1735 ! 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 ) 1736 !> @param[in] dd_value 3D array of variable to be extrapolated 1737 !> @param[in] dd_fill FillValue of variable 1738 !> @param[in] id_radius radius of the halo used to compute extrapolation 1739 !> @param[in] dd_coef 3D array of coefficient for inverse distance weighted extrapolation 1740 !> @return extrapolatd value 1741 !------------------------------------------------------------------- 1742 FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, & 1743 & dd_coef ) 1697 1744 IMPLICIT NONE 1698 1745 ! Argument 1699 1746 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1700 1747 REAL(dp) , INTENT(IN) :: dd_fill 1701 INTEGER(i4), INTENT(IN) :: id_ ext1748 INTEGER(i4), INTENT(IN) :: id_radius 1702 1749 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef 1703 1750 … … 1722 1769 extrap__3D_dist_weight_fill=dd_fill 1723 1770 1724 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ ext*2))1771 il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 1725 1772 1726 1773 IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN … … 1733 1780 dl_coef(:,:,:)=0 1734 1781 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 1782 DO jk=1,il_shape(3) 1783 DO jj=1,il_shape(2) 1784 DO ji=1,il_shape(1) 1785 1786 IF( dd_value(ji,jj,jk) /= dd_fill )THEN 1787 ! compute factor 1788 dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) 1789 dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) 1790 ENDIF 1791 1792 ENDDO 1793 ENDDO 1794 ENDDO 1745 1795 1746 1796 ! return value 1747 1797 IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN 1748 extrap__3D_dist_weight_fill=SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 1798 extrap__3D_dist_weight_fill = & 1799 & SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 1749 1800 ENDIF 1750 1801 … … 1755 1806 1756 1807 END FUNCTION extrap__3D_dist_weight_fill 1757 !> @endcode1758 1808 !------------------------------------------------------------------- 1759 1809 !> @brief … … 1761 1811 !> extraband of N points at north,south,east and west boundaries. 1762 1812 !> 1813 !> @details 1814 !> optionaly you could specify size of extra bands in i- and j-direction 1815 !> 1763 1816 !> @author J.Paul 1764 !> - 2013-Initial version1817 !> - November, 2013-Initial version 1765 1818 ! 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)1819 !> @param[inout] td_var variable 1820 !> @param[in] id_isize i-direction size of extra bands (default=im_minext) 1821 !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) 1769 1822 !> @todo 1770 1823 !> - invalid special case for grid with north fold 1771 1824 !------------------------------------------------------------------- 1772 !> @code1773 1825 SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize ) 1774 1826 IMPLICIT NONE … … 1821 1873 dl_value(:,:,:,:)=td_var%d_value(:,:,:,:) 1822 1874 1875 1823 1876 td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize 1824 1877 td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize … … 1855 1908 1856 1909 END SUBROUTINE extrap_add_extrabands 1857 !> @endcode1858 1910 !------------------------------------------------------------------- 1859 1911 !> @brief … … 1861 1913 !> of N points at north,south,east and west boundaries. 1862 1914 !> 1915 !> @details 1916 !> optionaly you could specify size of extra bands in i- and j-direction 1917 !> 1863 1918 !> @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 1919 !> - November, 2013-Initial version 1920 !> 1921 !> @param[inout] td_var variable 1922 !> @param[in] id_isize i-direction size of extra bands (default=im_minext) 1923 !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) 1924 !------------------------------------------------------------------- 1871 1925 SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize ) 1872 1926 IMPLICIT NONE 1873 1927 ! Argument 1874 TYPE(TVAR) , INTENT(INOUT) 1928 TYPE(TVAR) , INTENT(INOUT) :: td_var 1875 1929 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize 1876 1930 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize … … 1881 1935 INTEGER(i4) :: il_isize 1882 1936 INTEGER(i4) :: il_jsize 1883 1937 1884 1938 INTEGER(i4) :: il_imin 1885 1939 INTEGER(i4) :: il_imax … … 1935 1989 1936 1990 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 1991 END MODULE extrap
Note: See TracChangeset
for help on using the changeset viewer.