- Timestamp:
- 2018-10-29T15:20:26+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/TOOLS/SIREN/src/extrap.f90
r10248 r10251 19 19 !> defining string character _cn\_varinfo_. By default _dist_weight_.<br/> 20 20 !> Example: 21 !> - cn_varinfo='varname1: ext=dist_weight', 'varname2:ext=min_error'21 !> - cn_varinfo='varname1:dist_weight', 'varname2:min_error' 22 22 !> 23 23 !> to detect point to be extrapolated:<br/> 24 24 !> @code 25 !> il_detect(:,:,:)=extrap_detect(td_var )25 !> il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext]) 26 26 !> @endcode 27 27 !> - il_detect(:,:,:) is 3D array of point to be extrapolated 28 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] 29 33 !> 30 34 !> to extrapolate variable:<br/> 31 35 !> @code 32 !> CALL extrap_fill_value( td_var, [ id_radius])36 !> CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter]) 33 37 !> @endcode 34 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] 35 45 !> - id_radius is radius of the halo used to compute extrapolation [optional] 46 !> - id_maxiter is maximum number of iteration [optional] 36 47 !> 37 48 !> to add extraband to the variable (to be extrapolated):<br/> … … 51 62 !> - id_jsize : j-direction size of extra bands [optional] 52 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 !> 53 90 !> @warning _FillValue must not be zero (use var_chg_FillValue()) 54 91 !> … … 56 93 !> J.Paul 57 94 ! REVISION HISTORY: 58 !> @date Nov ember, 2013 - Initial Version95 !> @date Nov, 2013 - Initial Version 59 96 !> @date September, 2014 60 97 !> - add header 61 !> @date June, 201562 !> - extrapolate all land points (_FillValue)63 !> - move deriv function to math module64 !> @date July, 201565 !> - compute extrapolation from north west to south east,66 !> and from south east to north west67 98 !> 68 99 !> @todo 69 100 !> - create module for each extrapolation method 70 !> - smooth extrapolated points71 101 !> 72 102 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 80 110 USE date ! date manager 81 111 USE logger ! log file manager 82 USE math ! mathematical function83 112 USE att ! attribute manager 84 113 USE dim ! dimension manager … … 89 118 90 119 ! type and variable 120 PRIVATE :: im_maxiter !< default maximum number of iteration 91 121 PRIVATE :: im_minext !< default minumum number of point to extrapolate 92 122 PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation … … 97 127 PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) 98 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 99 132 100 133 PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper … … 108 141 PRIVATE :: extrap__3D_dist_weight_fill ! 109 142 143 INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration 110 144 INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate 111 145 INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation … … 137 171 !> 138 172 !> @author J.Paul 139 !> @date November, 2013 - Initial Version 140 !> @date June, 2015 141 !> - do not use level to select points to be extrapolated 173 !> - November, 2013- Initial Version 142 174 ! 143 175 !> @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 144 180 !> @return array of point to be extrapolated 145 181 !------------------------------------------------------------------- 146 FUNCTION extrap__detect( td_var0 ) 182 FUNCTION extrap__detect( td_var0, td_level1, & 183 & id_offset, id_rho, id_ext ) 147 184 IMPLICIT NONE 148 185 ! Argument 149 186 TYPE(TVAR) , INTENT(IN ) :: td_var0 187 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1 188 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 189 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 190 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext 150 191 151 192 ! function … … 155 196 156 197 ! local variable 198 CHARACTER(LEN=lc) :: cl_level 199 200 INTEGER(i4) :: il_ind 201 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect 202 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp 203 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset 204 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1 205 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1_G0 206 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_extra 207 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ext 208 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_rho 209 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_dim0 210 211 TYPE(TVAR) :: tl_var1 212 157 213 ! loop indices 158 214 INTEGER(i4) :: ji0 159 215 INTEGER(i4) :: jj0 160 216 INTEGER(i4) :: jk0 217 INTEGER(i4) :: ji1 218 INTEGER(i4) :: jj1 219 INTEGER(i4) :: ji1m 220 INTEGER(i4) :: jj1m 221 INTEGER(i4) :: ji1p 222 INTEGER(i4) :: jj1p 161 223 !---------------------------------------------------------------- 162 224 163 ! force to extrapolated all points 164 extrap__detect(:,:,:)=1 225 ! init 226 extrap__detect(:,:,:)=0 227 228 ALLOCATE( il_dim0(3) ) 229 il_dim0(:)=td_var0%t_dim(1:3)%i_len 230 231 ! optional argument 232 ALLOCATE( il_rho(ip_maxdim) ) 233 il_rho(:)=1 234 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 235 236 ALLOCATE( il_offset(ip_maxdim,2) ) 237 il_offset(:,:)=0 238 IF( PRESENT(id_offset) )THEN 239 il_offset(1:SIZE(id_offset(:,:),DIM=1),& 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) 244 ENDIF 245 246 ALLOCATE( il_ext(ip_maxdim) ) 247 il_ext(:)=im_minext 248 IF( PRESENT(id_ext) ) il_ext(1:SIZE(id_ext(:)))=id_ext(:) 249 250 ALLOCATE( il_detect(il_dim0(1),& 251 & il_dim0(2),& 252 & il_dim0(3)) ) 253 il_detect(:,:,:)=0 254 255 ! select point already inform 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 263 264 IF( PRESENT(td_level1) )THEN 265 SELECT CASE(TRIM(td_var0%c_point)) 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' 274 END SELECT 275 276 il_ind=var_get_index(td_level1(:),TRIM(cl_level)) 277 IF( il_ind == 0 )THEN 278 CALL logger_error("EXTRAP DETECT: can not compute point to be "//& 279 & "extrapolated for variable "//TRIM(td_var0%c_name)//& 280 & ". can not find "//& 281 & "level for variable point "//TRIM(TRIM(td_var0%c_point))) 282 ELSE 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 344 ENDIF 345 ENDIF 346 347 ! clean 348 DEALLOCATE( il_offset ) 349 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 371 372 ! clean 373 DEALLOCATE( il_tmp ) 165 374 166 375 ! do not compute grid point already inform … … 168 377 DO jj0=1,td_var0%t_dim(2)%i_len 169 378 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 379 IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=0 173 380 ENDDO 174 381 ENDDO 175 382 ENDDO 383 384 ! save result 385 extrap__detect(:,:,:)=il_detect(:,:,:) 386 387 ! clean 388 DEALLOCATE( il_dim0 ) 389 DEALLOCATE( il_ext ) 390 DEALLOCATE( il_detect ) 391 DEALLOCATE( il_rho ) 176 392 177 393 END FUNCTION extrap__detect … … 182 398 !> 183 399 !> @author J.Paul 184 !> @date November, 2013 - Initial Version 185 !> @date June, 2015 186 !> - select all land points for extrapolation 400 !> - November, 2013- Initial Version 187 401 !> 188 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 189 407 !> @return 3D array of point to be extrapolated 190 408 !------------------------------------------------------------------- 191 FUNCTION extrap__detect_wrapper( td_var ) 409 FUNCTION extrap__detect_wrapper( td_var, td_level, & 410 & id_offset, id_rho, id_ext ) 192 411 193 412 IMPLICIT NONE 194 413 ! Argument 195 414 TYPE(TVAR) , INTENT(IN ) :: td_var 415 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level 416 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 417 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 418 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext 196 419 197 420 ! function … … 216 439 & " for variable "//TRIM(td_var%c_name) ) 217 440 218 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var ) 441 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, & 442 & id_offset, & 443 & id_rho, & 444 & id_ext ) 219 445 220 446 ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN … … 224 450 & " for variable "//TRIM(td_var%c_name) ) 225 451 226 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var ) 452 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,& 453 & id_offset, & 454 & id_rho, & 455 & id_ext ) 227 456 228 457 ELSE IF( td_var%t_dim(3)%l_use )THEN … … 232 461 & " for variable "//TRIM(td_var%c_name) ) 233 462 234 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var ) 463 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, & 464 & id_offset, & 465 & id_rho, & 466 & id_ext ) 235 467 236 468 ENDIF … … 257 489 !> 258 490 !> @author J.Paul 259 !> @date November, 2013 - Initial Version 260 !> @date June, 2015 261 !> - select all land points for extrapolation 491 !> - Nov, 2013- Initial Version 262 492 ! 263 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 264 500 !> @param[in] id_radius radius of the halo used to compute extrapolation 265 !------------------------------------------------------------------- 266 SUBROUTINE extrap__fill_value_wrapper( td_var, & 267 & id_radius ) 501 !> @param[in] id_maxiter maximum number of iteration 502 !------------------------------------------------------------------- 503 SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & 504 & id_offset, & 505 & id_rho, & 506 & id_iext, id_jext, id_kext, & 507 & id_radius, id_maxiter ) 268 508 IMPLICIT NONE 269 509 ! Argument 270 510 TYPE(TVAR) , INTENT(INOUT) :: td_var 511 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level 512 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 513 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 514 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext 515 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext 516 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext 271 517 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius 518 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_maxiter 272 519 273 520 ! local variable 521 INTEGER(i4) :: il_iext 522 INTEGER(i4) :: il_jext 523 INTEGER(i4) :: il_kext 274 524 INTEGER(i4) :: il_radius 525 INTEGER(i4) :: il_maxiter 275 526 276 527 CHARACTER(LEN=lc) :: cl_method … … 293 544 END SELECT 294 545 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 546 il_iext=im_minext 547 IF( PRESENT(id_iext) ) il_iext=id_iext 548 il_jext=im_minext 549 IF( PRESENT(id_jext) ) il_jext=id_jext 550 il_kext=0 551 IF( PRESENT(id_kext) ) il_kext=id_kext 552 553 IF( TRIM(td_var%c_interp(1)) == 'cubic')THEN 554 IF( il_iext > 0 .AND. il_iext < im_mincubic ) il_iext=im_mincubic 555 IF( il_jext > 0 .AND. il_jext < im_mincubic ) il_jext=im_mincubic 556 ENDIF 557 558 IF( il_iext < 0 )THEN 299 559 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 300 & " radius of the box used to compute extrapolation "//&301 & "("//TRIM(fct_str(il_ radius))//")")560 & " number of points to be extrapolated in i-direction "//& 561 & "("//TRIM(fct_str(il_iext))//")") 302 562 ENDIF 303 563 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 ) 564 IF( il_jext < 0 )THEN 565 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 566 & " number of points to be extrapolated in j-direction "//& 567 & "("//TRIM(fct_str(il_jext))//")") 568 ENDIF 569 570 IF( il_kext < 0 )THEN 571 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 572 & " number of points to be extrapolated in k-direction "//& 573 & "("//TRIM(fct_str(il_kext))//")") 574 ENDIF 575 576 IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & 577 & (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. & 578 & (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN 579 580 ! number of point use to compute box 581 il_radius=1 582 IF( PRESENT(id_radius) ) il_radius=id_radius 583 IF( il_radius < 0 )THEN 584 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 585 & " radius of the box used to compute extrapolation "//& 586 & "("//TRIM(fct_str(il_radius))//")") 587 ENDIF 588 589 ! maximum number of iteration 590 il_maxiter=im_maxiter 591 IF( PRESENT(id_maxiter) ) il_maxiter=id_maxiter 592 IF( il_maxiter < 0 )THEN 593 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 594 & " maximum nuber of iteration "//& 595 & "("//TRIM(fct_str(il_maxiter))//")") 596 ENDIF 597 598 CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& 599 & " using "//TRIM(cl_method)//" method." ) 600 601 CALL extrap__fill_value( td_var, cl_method, & 602 & il_iext, il_jext, il_kext, & 603 & il_radius, il_maxiter, & 604 & td_level, & 605 & id_offset, id_rho ) 606 607 ENDIF 309 608 310 609 ENDIF … … 322 621 !> 323 622 !> @author J.Paul 324 !> @date November, 2013 - Initial Version 325 !> @date June, 2015 326 !> - select all land points for extrapolation 623 !> - November, 2013- Initial Version 327 624 ! 328 625 !> @param[inout] td_var variable structure 329 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 330 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 331 635 !------------------------------------------------------------------- 332 636 SUBROUTINE extrap__fill_value( td_var, cd_method, & 333 & id_radius ) 637 & id_iext, id_jext, id_kext, & 638 & id_radius, id_maxiter, & 639 & td_level, & 640 & id_offset, & 641 & id_rho ) 334 642 IMPLICIT NONE 335 643 ! Argument 336 644 TYPE(TVAR) , INTENT(INOUT) :: td_var 337 645 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 646 INTEGER(i4) , INTENT(IN ) :: id_iext 647 INTEGER(i4) , INTENT(IN ) :: id_jext 648 INTEGER(i4) , INTENT(IN ) :: id_kext 338 649 INTEGER(i4) , INTENT(IN ) :: id_radius 650 INTEGER(i4) , INTENT(IN ) :: id_maxiter 651 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level 652 INTEGER(i4) , DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 653 INTEGER(i4) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 339 654 340 655 ! local variable … … 353 668 & td_var%t_dim(3)%i_len) ) 354 669 355 il_detect(:,:,:) = extrap_detect( td_var ) 356 670 il_detect(:,:,:) = extrap_detect( td_var, td_level, & 671 & id_offset, & 672 & id_rho, & 673 & id_ext=(/id_iext, id_jext, id_kext/) ) 357 674 !2- add attribute to variable 358 675 cl_extrap=fct_concat(td_var%c_extrap(:)) … … 362 679 CALL att_clean(tl_att) 363 680 364 IF( ALL(il_detect(:,:,:)==1) )THEN 365 CALL logger_warn(" EXTRAP FILL: "//& 366 & " can not extrapolate "//TRIM(td_var%c_name)//& 367 & ". no value inform." ) 368 ELSE 369 CALL logger_info(" EXTRAP FILL: "//& 370 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 371 & " point(s) to extrapolate " ) 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 ) 380 ENDIF 681 CALL logger_info(" EXTRAP FILL: "//& 682 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 683 & " point(s) to extrapolate " ) 684 685 !3- extrapolate 686 CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, & 687 & il_detect(:,:,:), & 688 & cd_method, id_radius, id_maxiter ) 381 689 382 690 DEALLOCATE(il_detect) … … 397 705 !> 398 706 !> @author J.Paul 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 707 !> - Nov, 2013- Initial Version 403 708 ! 404 709 !> @param[inout] dd_value 3D array of variable to be extrapolated … … 409 714 !------------------------------------------------------------------- 410 715 SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 411 & cd_method, id_radius )716 & cd_method, id_radius, id_maxiter ) 412 717 IMPLICIT NONE 413 718 ! Argument 414 719 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 720 REAL(dp) , INTENT(IN ) :: dd_fill 721 INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 722 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 723 INTEGER(i4), INTENT(IN ) :: id_radius 724 INTEGER(i4), INTENT(IN ) :: id_maxiter 419 725 420 726 ! local variable 421 INTEGER(i4) :: il_imin 422 INTEGER(i4) :: il_imax 423 INTEGER(i4) :: il_jmin 424 INTEGER(i4) :: il_jmax 425 INTEGER(i4) :: il_kmin 426 INTEGER(i4) :: il_kmax 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 727 INTEGER(i4) :: il_imin 728 INTEGER(i4) :: il_imax 729 INTEGER(i4) :: il_jmin 730 INTEGER(i4) :: il_jmax 731 INTEGER(i4) :: il_kmin 732 INTEGER(i4) :: il_kmax 733 INTEGER(i4) :: il_iter 734 INTEGER(i4) :: il_radius 735 736 INTEGER(i4), DIMENSION(4) :: il_shape 737 INTEGER(i4), DIMENSION(3) :: il_dim 438 738 439 739 INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect … … 443 743 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz 444 744 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef 445 446 LOGICAL :: ll_iter447 745 448 746 ! loop indices … … 467 765 DO WHILE( ANY(il_detect(:,:,:)==1) ) 468 766 ! change extend value to minimize number of iteration 469 il_radius=id_radius+(il_iter-1) 470 ll_iter=.TRUE. 767 il_radius=id_radius+(il_iter/id_maxiter) 471 768 472 769 ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) … … 477 774 dl_dfdx(:,:,:)=dd_fill 478 775 IF( il_shape(1) > 1 )THEN 479 dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 480 & dd_fill, 'I' ) 776 dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) 481 777 ENDIF 482 778 … … 484 780 dl_dfdy(:,:,:)=dd_fill 485 781 IF( il_shape(2) > 1 )THEN 486 dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 487 & dd_fill, 'J' ) 782 dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) 488 783 ENDIF 489 784 … … 491 786 dl_dfdz(:,:,:)=dd_fill 492 787 IF( il_shape(3) > 1 )THEN 493 dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 494 & dd_fill, 'K' ) 788 dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) 495 789 ENDIF 496 790 … … 510 804 511 805 DO jk=1,il_shape(3) 512 ! from North West(1,1) to South East(il_shape(1),il_shape(2))513 806 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 514 807 DO jj=1,il_shape(2) … … 520 813 il_imin=MAX(ji-il_radius,1) 521 814 il_imax=MIN(ji+il_radius,il_shape(1)) 522 ! coef indices to be used523 il_i1 = il_radius-(ji-il_imin)+1524 il_i2 = il_radius+(il_imax-ji)+1525 815 IF( il_dim(1) == 1 )THEN 526 816 il_imin=ji 527 817 il_imax=ji 528 ! coef indices to be used529 il_i1 = 1530 il_i2 = 1531 818 ENDIF 532 533 819 534 820 il_jmin=MAX(jj-il_radius,1) 535 821 il_jmax=MIN(jj+il_radius,il_shape(2)) 536 ! coef indices to be used537 il_j1 = il_radius-(jj-il_jmin)+1538 il_j2 = il_radius+(il_jmax-jj)+1539 822 IF( il_dim(2) == 1 )THEN 540 823 il_jmin=jj 541 824 il_jmax=jj 542 ! coef indices to be used543 il_j1 = 1544 il_j2 = 1545 825 ENDIF 546 826 547 827 il_kmin=MAX(jk-il_radius,1) 548 828 il_kmax=MIN(jk+il_radius,il_shape(3)) 549 ! coef indices to be used550 il_k1 = il_radius-(jk-il_kmin)+1551 il_k2 = il_radius+(il_kmax-jk)+1552 829 IF( il_dim(3) == 1 )THEN 553 830 il_kmin=jk 554 831 il_kmax=jk 555 ! coef indices to be used556 il_k1 = 1557 il_k2 = 1558 832 ENDIF 559 833 … … 571 845 & il_jmin:il_jmax, & 572 846 & il_kmin:il_kmax ), & 573 & dl_coef(il_i1:il_i2, & 574 & il_j1:il_j2, & 575 & il_k1:il_k2) ) 847 & dl_coef(:,:,:) ) 576 848 577 849 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 578 850 il_detect(ji,jj,jk)= 0 579 ll_iter=.FALSE.580 ENDIF581 582 ENDIF583 584 ENDDO585 ENDDO586 ! from South East(il_shape(1),il_shape(2)) to North West(1,1)587 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE588 DO jj=il_shape(2),1,-1589 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE590 DO ji=il_shape(1),1,-1591 592 IF( il_detect(ji,jj,jk) == 1 )THEN593 594 il_imin=MAX(ji-il_radius,1)595 il_imax=MIN(ji+il_radius,il_shape(1))596 ! coef indices to be used597 il_i1 = il_radius-(ji-il_imin)+1598 il_i2 = il_radius+(il_imax-ji)+1599 IF( il_dim(1) == 1 )THEN600 il_imin=ji601 il_imax=ji602 ! coef indices to be used603 il_i1 = 1604 il_i2 = 1605 ENDIF606 607 608 il_jmin=MAX(jj-il_radius,1)609 il_jmax=MIN(jj+il_radius,il_shape(2))610 ! coef indices to be used611 il_j1 = il_radius-(jj-il_jmin)+1612 il_j2 = il_radius+(il_jmax-jj)+1613 IF( il_dim(2) == 1 )THEN614 il_jmin=jj615 il_jmax=jj616 ! coef indices to be used617 il_j1 = 1618 il_j2 = 1619 ENDIF620 621 il_kmin=MAX(jk-il_radius,1)622 il_kmax=MIN(jk+il_radius,il_shape(3))623 ! coef indices to be used624 il_k1 = il_radius-(jk-il_kmin)+1625 il_k2 = il_radius+(il_kmax-jk)+1626 IF( il_dim(3) == 1 )THEN627 il_kmin=jk628 il_kmax=jk629 ! coef indices to be used630 il_k1 = 1631 il_k2 = 1632 ENDIF633 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 )THEN652 il_detect(ji,jj,jk)= 0653 ll_iter=.FALSE.654 851 ENDIF 655 852 … … 665 862 DEALLOCATE( dl_coef ) 666 863 667 IF( ll_iter )il_iter=il_iter+1864 il_iter=il_iter+1 668 865 ENDDO 669 866 ENDDO … … 678 875 DO WHILE( ANY(il_detect(:,:,:)==1) ) 679 876 ! change extend value to minimize number of iteration 680 il_radius=id_radius+(il_iter-1) 681 ll_iter=.TRUE. 877 il_radius=id_radius+(il_iter/id_maxiter) 682 878 683 879 il_dim(1)=2*il_radius+1 … … 690 886 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 691 887 692 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), &693 & 1:il_dim(2), &694 & 1:il_dim(3), &888 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 889 & 1:il_dim(2), & 890 & 1:il_dim(3), & 695 891 & jl ) ) 696 892 697 893 DO jk=1,il_shape(3) 698 ! from North West(1,1) to South East(il_shape(1),il_shape(2))699 894 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 700 895 DO jj=1,il_shape(2) … … 706 901 il_imin=MAX(ji-il_radius,1) 707 902 il_imax=MIN(ji+il_radius,il_shape(1)) 708 ! coef indices to be used709 il_i1 = il_radius-(ji-il_imin)+1710 il_i2 = il_radius+(il_imax-ji)+1711 903 IF( il_dim(1) == 1 )THEN 712 904 il_imin=ji 713 905 il_imax=ji 714 ! coef indices to be used715 il_i1 = 1716 il_i2 = 1717 906 ENDIF 718 907 719 908 il_jmin=MAX(jj-il_radius,1) 720 909 il_jmax=MIN(jj+il_radius,il_shape(2)) 721 ! coef indices to be used722 il_j1 = il_radius-(jj-il_jmin)+1723 il_j2 = il_radius+(il_jmax-jj)+1724 910 IF( il_dim(2) == 1 )THEN 725 911 il_jmin=jj 726 912 il_jmax=jj 727 ! coef indices to be used728 il_j1 = 1729 il_j2 = 1730 913 ENDIF 731 914 732 915 il_kmin=MAX(jk-il_radius,1) 733 916 il_kmax=MIN(jk+il_radius,il_shape(3)) 734 ! coef indices to be used735 il_k1 = il_radius-(jk-il_kmin)+1736 il_k2 = il_radius+(il_kmax-jk)+1737 917 IF( il_dim(3) == 1 )THEN 738 918 il_kmin=jk 739 919 il_kmax=jk 740 ! coef indices to be used741 il_k1 = 1742 il_k2 = 1743 920 ENDIF 744 921 … … 748 925 & il_kmin:il_kmax, & 749 926 & jl), dd_fill, il_radius, & 750 & dl_coef(il_i1:il_i2, & 751 & il_j1:il_j2, & 752 & il_k1:il_k2) ) 927 & dl_coef(:,:,:) ) 753 928 754 929 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 755 930 il_detect(ji,jj,jk)= 0 756 ll_iter=.FALSE.757 ENDIF758 759 ENDIF760 761 ENDDO762 ENDDO763 ! from South East(il_shape(1),il_shape(2)) to North West(1,1)764 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE765 DO jj=il_shape(2),1,-1766 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE767 DO ji=il_shape(1),1,-1768 769 IF( il_detect(ji,jj,jk) == 1 )THEN770 771 il_imin=MAX(ji-il_radius,1)772 il_imax=MIN(ji+il_radius,il_shape(1))773 ! coef indices to be used774 il_i1 = il_radius-(ji-il_imin)+1775 il_i2 = il_radius+(il_imax-ji)+1776 IF( il_dim(1) == 1 )THEN777 il_imin=ji778 il_imax=ji779 ! coef indices to be used780 il_i1 = 1781 il_i2 = 1782 ENDIF783 784 il_jmin=MAX(jj-il_radius,1)785 il_jmax=MIN(jj+il_radius,il_shape(2))786 ! coef indices to be used787 il_j1 = il_radius-(jj-il_jmin)+1788 il_j2 = il_radius+(il_jmax-jj)+1789 IF( il_dim(2) == 1 )THEN790 il_jmin=jj791 il_jmax=jj792 ! coef indices to be used793 il_j1 = 1794 il_j2 = 1795 ENDIF796 797 il_kmin=MAX(jk-il_radius,1)798 il_kmax=MIN(jk+il_radius,il_shape(3))799 ! coef indices to be used800 il_k1 = il_radius-(jk-il_kmin)+1801 il_k2 = il_radius+(il_kmax-jk)+1802 IF( il_dim(3) == 1 )THEN803 il_kmin=jk804 il_kmax=jk805 ! coef indices to be used806 il_k1 = 1807 il_k2 = 1808 ENDIF809 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 )THEN820 il_detect(ji,jj,jk)= 0821 ll_iter=.FALSE.822 931 ENDIF 823 932 … … 827 936 ENDDO 828 937 ENDDO 829 CALL logger_info(" EXTRAP 3D: "//& 830 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 831 & " point(s) to extrapolate " ) 832 938 833 939 DEALLOCATE( dl_coef ) 834 IF( ll_iter )il_iter=il_iter+1940 il_iter=il_iter+1 835 941 ENDDO 836 942 ENDDO … … 840 946 841 947 END SUBROUTINE extrap__3D 948 !------------------------------------------------------------------- 949 !> @brief 950 !> This function compute derivative of 1D array. 951 !> 952 !> @details 953 !> optionaly you could specify to take into account east west discontinuity 954 !> (-180° 180° or 0° 360° for longitude variable) 955 !> 956 !> @author J.Paul 957 !> - November, 2013- Initial Version 958 ! 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 !------------------------------------------------------------------- 963 PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) 964 965 IMPLICIT NONE 966 ! Argument 967 REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value 968 REAL(dp) , INTENT(IN) :: dd_fill 969 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 970 971 ! function 972 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D 973 974 ! local variable 975 INTEGER(i4) :: il_imin 976 INTEGER(i4) :: il_imax 977 INTEGER(i4), DIMENSION(1) :: il_shape 978 979 REAL(dp) :: dl_min 980 REAL(dp) :: dl_max 981 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value 982 983 LOGICAL :: ll_discont 984 985 ! loop indices 986 INTEGER(i4) :: ji 987 988 INTEGER(i4) :: i1 989 INTEGER(i4) :: i2 990 !---------------------------------------------------------------- 991 ! init 992 extrap_deriv_1D(:)=dd_fill 993 994 ll_discont=.FALSE. 995 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 996 997 il_shape(:)=SHAPE(dd_value(:)) 998 999 ALLOCATE( dl_value(3)) 1000 1001 ! compute derivative in i-direction 1002 DO ji=1,il_shape(1) 1003 1004 il_imin=MAX(ji-1,1) 1005 il_imax=MIN(ji+1,il_shape(1)) 1006 1007 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 1008 i1=1 ; i2=3 1009 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 1010 i1=1 ; i2=2 1011 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 1012 i1=2 ; i2=3 1013 ENDIF 1014 1015 dl_value(i1:i2)=dd_value(il_imin:il_imax) 1016 IF( il_imin == 1 )THEN 1017 dl_value(:)=EOSHIFT( dl_value(:), & 1018 & DIM=1, & 1019 & SHIFT=-1, & 1020 & BOUNDARY=dl_value(1) ) 1021 ENDIF 1022 IF( il_imax == il_shape(1) )THEN 1023 dl_value(:)=EOSHIFT( dl_value(:), & 1024 & DIM=1, & 1025 & SHIFT=1, & 1026 & BOUNDARY=dl_value(3)) 1027 ENDIF 1028 1029 IF( ll_discont )THEN 1030 dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill ) 1031 dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill ) 1032 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1033 WHERE( dl_value(:) < 0._dp ) 1034 dl_value(:) = dl_value(:)+360._dp 1035 END WHERE 1036 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1037 WHERE( dl_value(:) > 180._dp ) 1038 dl_value(:) = dl_value(:)-180._dp 1039 END WHERE 1040 ENDIF 1041 ENDIF 1042 1043 IF( dl_value( 2) /= dd_fill .AND. & ! ji 1044 & dl_value( 3) /= dd_fill .AND. & ! ji+1 1045 & dl_value( 1) /= dd_fill )THEN ! ji-1 1046 1047 extrap_deriv_1D(ji)=& 1048 & ( dl_value(3) - dl_value(1) ) / & 1049 & REAL( il_imax-il_imin ,dp) 1050 1051 ENDIF 1052 1053 ENDDO 1054 1055 DEALLOCATE( dl_value ) 1056 1057 END FUNCTION extrap_deriv_1D 1058 !------------------------------------------------------------------- 1059 !> @brief 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 !> 1064 !> @details 1065 !> optionaly you could specify to take into account east west discontinuity 1066 !> (-180° 180° or 0° 360° for longitude variable) 1067 !> 1068 !> @author J.Paul 1069 !> - November, 2013- Initial Version 1070 ! 1071 !> @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 !------------------------------------------------------------------- 1076 FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) 1077 1078 IMPLICIT NONE 1079 ! Argument 1080 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value 1081 REAL(dp) , INTENT(IN) :: dd_fill 1082 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim 1083 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 1084 1085 ! function 1086 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 1087 & SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D 1088 1089 ! local variable 1090 INTEGER(i4) :: il_imin 1091 INTEGER(i4) :: il_imax 1092 INTEGER(i4) :: il_jmin 1093 INTEGER(i4) :: il_jmax 1094 INTEGER(i4), DIMENSION(2) :: il_shape 1095 1096 REAL(dp) :: dl_min 1097 REAL(dp) :: dl_max 1098 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value 1099 1100 LOGICAL :: ll_discont 1101 1102 ! loop indices 1103 INTEGER(i4) :: ji 1104 INTEGER(i4) :: jj 1105 1106 INTEGER(i4) :: i1 1107 INTEGER(i4) :: i2 1108 1109 INTEGER(i4) :: j1 1110 INTEGER(i4) :: j2 1111 !---------------------------------------------------------------- 1112 ! init 1113 extrap_deriv_2D(:,:)=dd_fill 1114 1115 ll_discont=.FALSE. 1116 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 1117 1118 il_shape(:)=SHAPE(dd_value(:,:)) 1119 1120 SELECT CASE(TRIM(fct_upper(cd_dim))) 1121 1122 CASE('I') 1123 1124 ALLOCATE( dl_value(3,il_shape(2)) ) 1125 ! compute derivative in i-direction 1126 DO ji=1,il_shape(1) 1127 1128 ! init 1129 dl_value(:,:)=dd_fill 1130 1131 il_imin=MAX(ji-1,1) 1132 il_imax=MIN(ji+1,il_shape(1)) 1133 1134 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 1135 i1=1 ; i2=3 1136 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 1137 i1=1 ; i2=2 1138 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 1139 i1=2 ; i2=3 1140 ENDIF 1141 1142 dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:) 1143 IF( il_imin == 1 )THEN 1144 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1145 & DIM=1, & 1146 & SHIFT=-1, & 1147 & BOUNDARY=dl_value(1,:) ) 1148 ENDIF 1149 IF( il_imax == il_shape(1) )THEN 1150 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1151 & DIM=1, & 1152 & SHIFT=1, & 1153 & BOUNDARY=dl_value(3,:)) 1154 ENDIF 1155 1156 IF( ll_discont )THEN 1157 dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1158 dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1159 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1160 WHERE( dl_value(:,:) < 0_dp ) 1161 dl_value(:,:) = dl_value(:,:)+360._dp 1162 END WHERE 1163 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1164 WHERE( dl_value(:,:) > 180 ) 1165 dl_value(:,:) = dl_value(:,:)-180._dp 1166 END WHERE 1167 ENDIF 1168 ENDIF 1169 1170 WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji 1171 & dl_value(3,:) /= dd_fill .AND. & ! ji+1 1172 & dl_value(1,:) /= dd_fill ) ! ji-1 1173 1174 extrap_deriv_2D(ji,:)=& 1175 & ( dl_value(3,:) - dl_value(1,:) ) / & 1176 & REAL( il_imax-il_imin,dp) 1177 1178 END WHERE 1179 1180 ENDDO 1181 1182 CASE('J') 1183 1184 ALLOCATE( dl_value(il_shape(1),3) ) 1185 ! compute derivative in j-direction 1186 DO jj=1,il_shape(2) 1187 1188 il_jmin=MAX(jj-1,1) 1189 il_jmax=MIN(jj+1,il_shape(2)) 1190 1191 IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 1192 j1=1 ; j2=3 1193 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 1194 j1=1 ; j2=2 1195 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 1196 j1=2 ; j2=3 1197 ENDIF 1198 1199 dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax) 1200 IF( il_jmin == 1 )THEN 1201 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1202 & DIM=2, & 1203 & SHIFT=-1, & 1204 & BOUNDARY=dl_value(:,1)) 1205 ENDIF 1206 IF( il_jmax == il_shape(2) )THEN 1207 dl_value(:,:)=EOSHIFT( dl_value(:,:), & 1208 & DIM=2, & 1209 & SHIFT=1, & 1210 & BOUNDARY=dl_value(:,3)) 1211 ENDIF 1212 1213 IF( ll_discont )THEN 1214 dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1215 dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 1216 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1217 WHERE( dl_value(:,:) < 0_dp ) 1218 dl_value(:,:) = dl_value(:,:)+360._dp 1219 END WHERE 1220 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1221 WHERE( dl_value(:,:) > 180 ) 1222 dl_value(:,:) = dl_value(:,:)-180._dp 1223 END WHERE 1224 ENDIF 1225 ENDIF 1226 1227 WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj 1228 & dl_value(:, 3) /= dd_fill .AND. & ! jj+1 1229 & dl_value(:, 1) /= dd_fill ) ! jj-1 1230 1231 extrap_deriv_2D(:,jj)=& 1232 & ( dl_value(:,3) - dl_value(:,1) ) / & 1233 & REAL(il_jmax-il_jmin,dp) 1234 1235 END WHERE 1236 1237 ENDDO 1238 1239 END SELECT 1240 1241 DEALLOCATE( dl_value ) 1242 1243 END FUNCTION extrap_deriv_2D 1244 !------------------------------------------------------------------- 1245 !> @brief 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. 1249 !> 1250 !> @details 1251 !> optionaly you could specify to take into account east west discontinuity 1252 !> (-180° 180° or 0° 360° for longitude variable) 1253 !> 1254 !> @author J.Paul 1255 !> - November, 2013- Initial Version 1256 ! 1257 !> @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 !------------------------------------------------------------------- 1262 PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) 1263 1264 IMPLICIT NONE 1265 ! Argument 1266 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value 1267 REAL(dp) , INTENT(IN) :: dd_fill 1268 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim 1269 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont 1270 1271 ! function 1272 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 1273 & SIZE(dd_value,DIM=2), & 1274 & SIZE(dd_value,DIM=3)) :: extrap_deriv_3D 1275 1276 ! local variable 1277 INTEGER(i4) :: il_imin 1278 INTEGER(i4) :: il_imax 1279 INTEGER(i4) :: il_jmin 1280 INTEGER(i4) :: il_jmax 1281 INTEGER(i4) :: il_kmin 1282 INTEGER(i4) :: il_kmax 1283 INTEGER(i4), DIMENSION(3) :: il_shape 1284 1285 REAL(dp) :: dl_min 1286 REAL(dp) :: dl_max 1287 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value 1288 1289 LOGICAL :: ll_discont 1290 1291 ! loop indices 1292 INTEGER(i4) :: ji 1293 INTEGER(i4) :: jj 1294 INTEGER(i4) :: jk 1295 1296 INTEGER(i4) :: i1 1297 INTEGER(i4) :: i2 1298 1299 INTEGER(i4) :: j1 1300 INTEGER(i4) :: j2 1301 1302 INTEGER(i4) :: k1 1303 INTEGER(i4) :: k2 1304 !---------------------------------------------------------------- 1305 ! init 1306 extrap_deriv_3D(:,:,:)=dd_fill 1307 1308 ll_discont=.FALSE. 1309 IF( PRESENT(ld_discont) ) ll_discont=ld_discont 1310 1311 il_shape(:)=SHAPE(dd_value(:,:,:)) 1312 1313 1314 SELECT CASE(TRIM(fct_upper(cd_dim))) 1315 1316 CASE('I') 1317 1318 ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) ) 1319 ! compute derivative in i-direction 1320 DO ji=1,il_shape(1) 1321 1322 il_imin=MAX(ji-1,1) 1323 il_imax=MIN(ji+1,il_shape(1)) 1324 1325 IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 1326 i1=1 ; i2=3 1327 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 1328 i1=1 ; i2=2 1329 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 1330 i1=2 ; i2=3 1331 ENDIF 1332 1333 dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:) 1334 IF( il_imin == 1 )THEN 1335 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1336 & DIM=1, & 1337 & SHIFT=-1, & 1338 & BOUNDARY=dl_value(1,:,:) ) 1339 ENDIF 1340 IF( il_imax == il_shape(1) )THEN 1341 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1342 & DIM=1, & 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. & ! ji 1362 & dl_value(3,:,:) /= dd_fill .AND. & !ji+1 1363 & dl_value(1,:,:) /= dd_fill ) !ji-1 1364 1365 extrap_deriv_3D(ji,:,:)= & 1366 & ( dl_value(3,:,:) - dl_value(1,:,:) ) / & 1367 & REAL( il_imax-il_imin ,dp) 1368 1369 END WHERE 1370 1371 ENDDO 1372 1373 CASE('J') 1374 1375 ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) ) 1376 ! compute derivative in j-direction 1377 DO jj=1,il_shape(2) 1378 1379 il_jmin=MAX(jj-1,1) 1380 il_jmax=MIN(jj+1,il_shape(2)) 1381 1382 IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 1383 j1=1 ; j2=3 1384 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 1385 j1=1 ; j2=2 1386 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 1387 j1=2 ; j2=3 1388 ENDIF 1389 1390 dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:) 1391 IF( il_jmin == 1 )THEN 1392 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1393 & DIM=2, & 1394 & SHIFT=-1, & 1395 & BOUNDARY=dl_value(:,1,:) ) 1396 ENDIF 1397 IF( il_jmax == il_shape(2) )THEN 1398 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1399 & DIM=2, & 1400 & SHIFT=1, & 1401 & BOUNDARY=dl_value(:,3,:)) 1402 ENDIF 1403 1404 IF( ll_discont )THEN 1405 dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1406 dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1407 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1408 WHERE( dl_value(:,:,:) < 0_dp ) 1409 dl_value(:,:,:) = dl_value(:,:,:)+360._dp 1410 END WHERE 1411 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1412 WHERE( dl_value(:,:,:) > 180 ) 1413 dl_value(:,:,:) = dl_value(:,:,:)-180._dp 1414 END WHERE 1415 ENDIF 1416 ENDIF 1417 1418 WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj 1419 & dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1 1420 & dl_value(:, 1,:) /= dd_fill ) ! jj-1 1421 1422 extrap_deriv_3D(:,jj,:)=& 1423 & ( dl_value(:,3,:) - dl_value(:,1,:) ) / & 1424 & REAL( il_jmax - il_jmin ,dp) 1425 1426 END WHERE 1427 1428 ENDDO 1429 1430 CASE('K') 1431 ! compute derivative in k-direction 1432 DO jk=1,il_shape(3) 1433 1434 il_kmin=MAX(jk-1,1) 1435 il_kmax=MIN(jk+1,il_shape(3)) 1436 1437 IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN 1438 k1=1 ; k2=3 1439 ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN 1440 k1=1 ; k2=2 1441 ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN 1442 k1=2 ; k2=3 1443 ENDIF 1444 1445 dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax) 1446 IF( il_kmin == 1 )THEN 1447 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1448 & DIM=3, & 1449 & SHIFT=-1, & 1450 & BOUNDARY=dl_value(:,:,1) ) 1451 ENDIF 1452 IF( il_kmax == il_shape(3) )THEN 1453 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 1454 & DIM=3, & 1455 & SHIFT=1, & 1456 & BOUNDARY=dl_value(:,:,3)) 1457 ENDIF 1458 1459 IF( ll_discont )THEN 1460 dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1461 dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 1462 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 1463 WHERE( dl_value(:,:,:) < 0_dp ) 1464 dl_value(:,:,:) = dl_value(:,:,:)+360._dp 1465 END WHERE 1466 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 1467 WHERE( dl_value(:,:,:) > 180 ) 1468 dl_value(:,:,:) = dl_value(:,:,:)-180._dp 1469 END WHERE 1470 ENDIF 1471 ENDIF 1472 1473 WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk 1474 & dl_value(:,:, 3) /= dd_fill .AND. & ! jk+1 1475 & dl_value(:,:, 1) /= dd_fill ) ! jk-1 1476 1477 extrap_deriv_3D(:,:,jk)=& 1478 & ( dl_value(:,:,3) - dl_value(:,:,1) ) / & 1479 & REAL( il_kmax-il_kmin,dp) 1480 1481 END WHERE 1482 1483 ENDDO 1484 1485 END SELECT 1486 1487 DEALLOCATE( dl_value ) 1488 1489 END FUNCTION extrap_deriv_3D 842 1490 !------------------------------------------------------------------- 843 1491 !> @brief … … 845 1493 !> 846 1494 !> @details 847 !> coefficients are "grid distance" to the center of the box 848 !> choosed to computeextrapolation.1495 !> coefficients are "grid distance" to the center of the box choosed to compute 1496 !> extrapolation. 849 1497 !> 850 1498 !> @author J.Paul 851 !> @date November, 2013 - Initial Version 852 !> @date July, 2015 853 !> - decrease weight of third dimension 1499 !> - November, 2013- Initial Version 854 1500 ! 855 1501 !> @param[in] dd_value 3D array of variable to be extrapolated … … 898 1544 899 1545 ! compute distance 900 ! "vertical weight" is lower than horizontal901 1546 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 902 1547 & (jj-il_jmid)**2 + & 903 & 3*(jk-il_kmid)**21548 & (jk-il_kmid)**2 904 1549 905 1550 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 924 1569 !> 925 1570 !> @author J.Paul 926 !> @date November, 2013- Initial Version1571 !> - November, 2013- Initial Version 927 1572 !> 928 1573 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1013 1658 !> 1014 1659 !> @author J.Paul 1015 !> @date November, 2013 - Initial Version 1016 !> @date July, 2015 1017 !> - decrease weight of third dimension 1660 !> - November, 2013- Initial Version 1018 1661 ! 1019 1662 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1062 1705 1063 1706 ! compute distance 1064 ! "vertical weight" is lower than horizontal1065 1707 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 1066 1708 & (jj-il_jmid)**2 + & 1067 & 3*(jk-il_kmid)**21709 & (jk-il_kmid)**2 1068 1710 1069 1711 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 1090 1732 !> 1091 1733 !> @author J.Paul 1092 !> @date November, 2013- Initial Version1734 !> - November, 2013- Initial Version 1093 1735 ! 1094 1736 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1121 1763 INTEGER(i4) :: jj 1122 1764 INTEGER(i4) :: jk 1765 1123 1766 !---------------------------------------------------------------- 1124 1767 … … 1150 1793 ENDDO 1151 1794 ENDDO 1152 1153 1795 1154 1796 ! return value … … 1173 1815 !> 1174 1816 !> @author J.Paul 1175 !> @date November, 2013 -Initial version1817 !> - November, 2013-Initial version 1176 1818 ! 1177 1819 !> @param[inout] td_var variable … … 1275 1917 !> 1276 1918 !> @author J.Paul 1277 !> @date November, 2013 -Initial version1919 !> - November, 2013-Initial version 1278 1920 !> 1279 1921 !> @param[inout] td_var variable
Note: See TracChangeset
for help on using the changeset viewer.
