Changeset 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src/extrap.f90
- Timestamp:
- 2015-12-03T09:10:32+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src/extrap.f90
r5260 r5989 19 19 !> defining string character _cn\_varinfo_. By default _dist_weight_.<br/> 20 20 !> Example: 21 !> - cn_varinfo='varname1: dist_weight', 'varname2:min_error'21 !> - cn_varinfo='varname1:ext=dist_weight', 'varname2:ext=min_error' 22 22 !> 23 23 !> to detect point to be extrapolated:<br/> 24 24 !> @code 25 !> il_detect(:,:,:)=extrap_detect(td_var , [td_level], [id_offset,] [id_rho,] [id_ext])25 !> il_detect(:,:,:)=extrap_detect(td_var) 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]33 29 !> 34 30 !> to extrapolate variable:<br/> 35 31 !> @code 36 !> CALL extrap_fill_value( td_var, [ td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter])32 !> CALL extrap_fill_value( td_var, [id_radius]) 37 33 !> @endcode 38 34 !> - 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 35 !> - id_radius is radius of the halo used to compute extrapolation [optional] 46 !> - id_maxiter is maximum number of iteration [optional]47 36 !> 48 37 !> to add extraband to the variable (to be extrapolated):<br/> … … 62 51 !> - id_jsize : j-direction size of extra bands [optional] 63 52 !> 64 !> to compute first derivative of 1D array:<br/>65 !> @code66 !> dl_value(:)=extrap_deriv_1D( dd_value(:), dd_fill, [ld_discont] )67 !> @endcode68 !> - dd_value is 1D array of variable69 !> - dd_fill is FillValue of variable70 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]71 !>72 !> to compute first derivative of 2D array:<br/>73 !> @code74 !> dl_value(:,:)=extrap_deriv_2D( dd_value(:,:), dd_fill, cd_dim, [ld_discont] )75 !> @endcode76 !> - dd_value is 2D array of variable77 !> - dd_fill is FillValue of variable78 !> - cd_dim is character to compute derivative on first (I) or second (J) dimension79 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]80 !>81 !> to compute first derivative of 3D array:<br/>82 !> @code83 !> dl_value(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim, [ld_discont] )84 !> @endcode85 !> - dd_value is 3D array of variable86 !> - dd_fill is FillValue of variable87 !> - cd_dim is character to compute derivative on first (I), second (J), or third (K) dimension88 !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]89 !>90 53 !> @warning _FillValue must not be zero (use var_chg_FillValue()) 91 54 !> … … 93 56 !> J.Paul 94 57 ! REVISION HISTORY: 95 !> @date Nov , 2013 - Initial Version58 !> @date November, 2013 - Initial Version 96 59 !> @date September, 2014 97 60 !> - add header 61 !> @date June, 2015 62 !> - extrapolate all land points (_FillValue) 63 !> - move deriv function to math module 64 !> @date July, 2015 65 !> - compute extrapolation from north west to south east, 66 !> and from south east to north west 98 67 !> 99 68 !> @todo 100 69 !> - create module for each extrapolation method 70 !> - smooth extrapolated points 101 71 !> 102 72 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 110 80 USE date ! date manager 111 81 USE logger ! log file manager 82 USE math ! mathematical function 112 83 USE att ! attribute manager 113 84 USE dim ! dimension manager … … 118 89 119 90 ! type and variable 120 PRIVATE :: im_maxiter !< default maximum number of iteration121 91 PRIVATE :: im_minext !< default minumum number of point to extrapolate 122 92 PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation … … 127 97 PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) 128 98 PUBLIC :: extrap_del_extrabands !< delete extraband of the variable 129 PUBLIC :: extrap_deriv_1D !< compute first derivative of 1D array130 PUBLIC :: extrap_deriv_2D !< compute first derivative of 2D array131 PUBLIC :: extrap_deriv_3D !< compute first derivative of 3D array132 99 133 100 PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper … … 141 108 PRIVATE :: extrap__3D_dist_weight_fill ! 142 109 143 INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration144 110 INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate 145 111 INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation … … 171 137 !> 172 138 !> @author J.Paul 173 !> - November, 2013- Initial Version 139 !> @date November, 2013 - Initial Version 140 !> @date June, 2015 141 !> - do not use level to select points to be extrapolated 174 142 ! 175 143 !> @param[in] td_var0 coarse grid variable to extrapolate 176 !> @param[in] td_level1 fine grid array of level177 !> @param[in] id_offset array of offset between fine and coarse grid178 !> @param[in] id_rho array of refinment factor179 !> @param[in] id_ext array of number of points to be extrapolated180 144 !> @return array of point to be extrapolated 181 145 !------------------------------------------------------------------- 182 FUNCTION extrap__detect( td_var0, td_level1, & 183 & id_offset, id_rho, id_ext ) 146 FUNCTION extrap__detect( td_var0 ) 184 147 IMPLICIT NONE 185 148 ! Argument 186 149 TYPE(TVAR) , INTENT(IN ) :: td_var0 187 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1188 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset189 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho190 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext191 150 192 151 ! function … … 196 155 197 156 ! local variable 198 CHARACTER(LEN=lc) :: cl_level199 200 INTEGER(i4) :: il_ind201 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect202 INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp203 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset204 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1205 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1_G0206 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_extra207 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ext208 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_rho209 INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_dim0210 211 TYPE(TVAR) :: tl_var1212 213 157 ! loop indices 214 158 INTEGER(i4) :: ji0 215 159 INTEGER(i4) :: jj0 216 160 INTEGER(i4) :: jk0 217 INTEGER(i4) :: ji1218 INTEGER(i4) :: jj1219 INTEGER(i4) :: ji1m220 INTEGER(i4) :: jj1m221 INTEGER(i4) :: ji1p222 INTEGER(i4) :: jj1p223 161 !---------------------------------------------------------------- 224 162 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 ) 163 ! force to extrapolated all points 164 extrap__detect(:,:,:)=1 374 165 375 166 ! do not compute grid point already inform … … 377 168 DO jj0=1,td_var0%t_dim(2)%i_len 378 169 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 170 IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill )THEN 171 extrap__detect(ji0,jj0,jk0)=0 172 ENDIF 380 173 ENDDO 381 174 ENDDO 382 175 ENDDO 383 384 ! save result385 extrap__detect(:,:,:)=il_detect(:,:,:)386 387 ! clean388 DEALLOCATE( il_dim0 )389 DEALLOCATE( il_ext )390 DEALLOCATE( il_detect )391 DEALLOCATE( il_rho )392 176 393 177 END FUNCTION extrap__detect … … 398 182 !> 399 183 !> @author J.Paul 400 !> - November, 2013- Initial Version 184 !> @date November, 2013 - Initial Version 185 !> @date June, 2015 186 !> - select all land points for extrapolation 401 187 !> 402 188 !> @param[in] td_var coarse grid variable to extrapolate 403 !> @param[in] td_level fine grid array of level404 !> @param[in] id_offset array of offset between fine and coarse grid405 !> @param[in] id_rho array of refinment factor406 !> @param[in] id_ext array of number of points to be extrapolated407 189 !> @return 3D array of point to be extrapolated 408 190 !------------------------------------------------------------------- 409 FUNCTION extrap__detect_wrapper( td_var, td_level, & 410 & id_offset, id_rho, id_ext ) 191 FUNCTION extrap__detect_wrapper( td_var ) 411 192 412 193 IMPLICIT NONE 413 194 ! Argument 414 195 TYPE(TVAR) , INTENT(IN ) :: td_var 415 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level416 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset417 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho418 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext419 196 420 197 ! function … … 439 216 & " for variable "//TRIM(td_var%c_name) ) 440 217 441 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, & 442 & id_offset, & 443 & id_rho, & 444 & id_ext ) 218 extrap__detect_wrapper(:,:,:)=extrap__detect( td_var ) 445 219 446 220 ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN … … 450 224 & " for variable "//TRIM(td_var%c_name) ) 451 225 452 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,& 453 & id_offset, & 454 & id_rho, & 455 & id_ext ) 226 extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var ) 456 227 457 228 ELSE IF( td_var%t_dim(3)%l_use )THEN … … 461 232 & " for variable "//TRIM(td_var%c_name) ) 462 233 463 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, & 464 & id_offset, & 465 & id_rho, & 466 & id_ext ) 234 extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var ) 467 235 468 236 ENDIF … … 489 257 !> 490 258 !> @author J.Paul 491 !> - Nov, 2013- Initial Version 259 !> @date November, 2013 - Initial Version 260 !> @date June, 2015 261 !> - select all land points for extrapolation 492 262 ! 493 263 !> @param[inout] td_var variable structure 494 !> @param[in] td_level fine grid array of level495 !> @param[in] id_offset array of offset between fine and coarse grid496 !> @param[in] id_rho array of refinment factor497 !> @param[in] id_iext number of points to be extrapolated in i-direction498 !> @param[in] id_jext number of points to be extrapolated in j-direction499 !> @param[in] id_kext number of points to be extrapolated in k-direction500 264 !> @param[in] id_radius radius of the halo used to compute extrapolation 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 ) 265 !------------------------------------------------------------------- 266 SUBROUTINE extrap__fill_value_wrapper( td_var, & 267 & id_radius ) 508 268 IMPLICIT NONE 509 269 ! Argument 510 270 TYPE(TVAR) , INTENT(INOUT) :: td_var 511 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level512 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset513 INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho514 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext515 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext516 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext517 271 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius 518 INTEGER(i4), INTENT(IN ), OPTIONAL :: id_maxiter519 272 520 273 ! local variable 521 INTEGER(i4) :: il_iext522 INTEGER(i4) :: il_jext523 INTEGER(i4) :: il_kext524 274 INTEGER(i4) :: il_radius 525 INTEGER(i4) :: il_maxiter526 275 527 276 CHARACTER(LEN=lc) :: cl_method … … 544 293 END SELECT 545 294 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 295 ! number of point use to compute box 296 il_radius=1 297 IF( PRESENT(id_radius) ) il_radius=id_radius 298 IF( il_radius < 0 )THEN 299 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 300 & " radius of the box used to compute extrapolation "//& 301 & "("//TRIM(fct_str(il_radius))//")") 556 302 ENDIF 557 303 558 IF( il_iext < 0 )THEN 559 CALL logger_error("EXTRAP FILL VALUE: invalid "//& 560 & " number of points to be extrapolated in i-direction "//& 561 & "("//TRIM(fct_str(il_iext))//")") 562 ENDIF 563 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 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 ) 608 309 609 310 ENDIF … … 621 322 !> 622 323 !> @author J.Paul 623 !> - November, 2013- Initial Version 324 !> @date November, 2013 - Initial Version 325 !> @date June, 2015 326 !> - select all land points for extrapolation 624 327 ! 625 328 !> @param[inout] td_var variable structure 626 329 !> @param[in] cd_method extrapolation method 627 !> @param[in] id_iext number of points to be extrapolated in i-direction628 !> @param[in] id_jext number of points to be extrapolated in j-direction629 !> @param[in] id_kext number of points to be extrapolated in k-direction630 330 !> @param[in] id_radius radius of the halo used to compute extrapolation 631 !> @param[in] id_maxiter maximum number of iteration632 !> @param[in] td_level fine grid array of level633 !> @param[in] id_offset array of offset between fine and coarse grid634 !> @param[in] id_rho array of refinment factor635 331 !------------------------------------------------------------------- 636 332 SUBROUTINE extrap__fill_value( td_var, cd_method, & 637 & id_iext, id_jext, id_kext, & 638 & id_radius, id_maxiter, & 639 & td_level, & 640 & id_offset, & 641 & id_rho ) 333 & id_radius ) 642 334 IMPLICIT NONE 643 335 ! Argument 644 336 TYPE(TVAR) , INTENT(INOUT) :: td_var 645 337 CHARACTER(LEN=*), INTENT(IN ) :: cd_method 646 INTEGER(i4) , INTENT(IN ) :: id_iext647 INTEGER(i4) , INTENT(IN ) :: id_jext648 INTEGER(i4) , INTENT(IN ) :: id_kext649 338 INTEGER(i4) , INTENT(IN ) :: id_radius 650 INTEGER(i4) , INTENT(IN ) :: id_maxiter651 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level652 INTEGER(i4) , DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset653 INTEGER(i4) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho654 339 655 340 ! local variable … … 668 353 & td_var%t_dim(3)%i_len) ) 669 354 670 il_detect(:,:,:) = extrap_detect( td_var, td_level, & 671 & id_offset, & 672 & id_rho, & 673 & id_ext=(/id_iext, id_jext, id_kext/) ) 355 il_detect(:,:,:) = extrap_detect( td_var ) 356 674 357 !2- add attribute to variable 675 358 cl_extrap=fct_concat(td_var%c_extrap(:)) … … 679 362 CALL att_clean(tl_att) 680 363 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 ) 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 689 381 690 382 DEALLOCATE(il_detect) … … 705 397 !> 706 398 !> @author J.Paul 707 !> - Nov, 2013- Initial Version 399 !> @date November, 2013 - Initial Version 400 !> @date July, 2015 401 !> - compute coef indices to be used 402 !> - bug fix: force coef indice to 1, for dimension lenth equal to 1 708 403 ! 709 404 !> @param[inout] dd_value 3D array of variable to be extrapolated … … 714 409 !------------------------------------------------------------------- 715 410 SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 716 & cd_method, id_radius , id_maxiter)411 & cd_method, id_radius ) 717 412 IMPLICIT NONE 718 413 ! Argument 719 414 REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 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 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 725 419 726 420 ! local variable 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 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 738 438 739 439 INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect … … 743 443 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz 744 444 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef 445 446 LOGICAL :: ll_iter 745 447 746 448 ! loop indices … … 765 467 DO WHILE( ANY(il_detect(:,:,:)==1) ) 766 468 ! change extend value to minimize number of iteration 767 il_radius=id_radius+(il_iter/id_maxiter) 469 il_radius=id_radius+(il_iter-1) 470 ll_iter=.TRUE. 768 471 769 472 ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) … … 774 477 dl_dfdx(:,:,:)=dd_fill 775 478 IF( il_shape(1) > 1 )THEN 776 dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) 479 dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 480 & dd_fill, 'I' ) 777 481 ENDIF 778 482 … … 780 484 dl_dfdy(:,:,:)=dd_fill 781 485 IF( il_shape(2) > 1 )THEN 782 dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) 486 dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 487 & dd_fill, 'J' ) 783 488 ENDIF 784 489 … … 786 491 dl_dfdz(:,:,:)=dd_fill 787 492 IF( il_shape(3) > 1 )THEN 788 dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) 493 dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 494 & dd_fill, 'K' ) 789 495 ENDIF 790 496 … … 804 510 805 511 DO jk=1,il_shape(3) 512 ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 806 513 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 807 514 DO jj=1,il_shape(2) … … 813 520 il_imin=MAX(ji-il_radius,1) 814 521 il_imax=MIN(ji+il_radius,il_shape(1)) 522 ! coef indices to be used 523 il_i1 = il_radius-(ji-il_imin)+1 524 il_i2 = il_radius+(il_imax-ji)+1 815 525 IF( il_dim(1) == 1 )THEN 816 526 il_imin=ji 817 527 il_imax=ji 818 ENDIF 528 ! coef indices to be used 529 il_i1 = 1 530 il_i2 = 1 531 ENDIF 532 819 533 820 534 il_jmin=MAX(jj-il_radius,1) 821 535 il_jmax=MIN(jj+il_radius,il_shape(2)) 536 ! coef indices to be used 537 il_j1 = il_radius-(jj-il_jmin)+1 538 il_j2 = il_radius+(il_jmax-jj)+1 822 539 IF( il_dim(2) == 1 )THEN 823 540 il_jmin=jj 824 541 il_jmax=jj 542 ! coef indices to be used 543 il_j1 = 1 544 il_j2 = 1 825 545 ENDIF 826 546 827 547 il_kmin=MAX(jk-il_radius,1) 828 548 il_kmax=MIN(jk+il_radius,il_shape(3)) 549 ! coef indices to be used 550 il_k1 = il_radius-(jk-il_kmin)+1 551 il_k2 = il_radius+(il_kmax-jk)+1 829 552 IF( il_dim(3) == 1 )THEN 830 553 il_kmin=jk 831 554 il_kmax=jk 555 ! coef indices to be used 556 il_k1 = 1 557 il_k2 = 1 832 558 ENDIF 833 559 … … 845 571 & il_jmin:il_jmax, & 846 572 & il_kmin:il_kmax ), & 847 & dl_coef(:,:,:) ) 573 & dl_coef(il_i1:il_i2, & 574 & il_j1:il_j2, & 575 & il_k1:il_k2) ) 848 576 849 577 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 850 578 il_detect(ji,jj,jk)= 0 579 ll_iter=.FALSE. 580 ENDIF 581 582 ENDIF 583 584 ENDDO 585 ENDDO 586 ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 587 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 588 DO jj=il_shape(2),1,-1 589 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 590 DO ji=il_shape(1),1,-1 591 592 IF( il_detect(ji,jj,jk) == 1 )THEN 593 594 il_imin=MAX(ji-il_radius,1) 595 il_imax=MIN(ji+il_radius,il_shape(1)) 596 ! coef indices to be used 597 il_i1 = il_radius-(ji-il_imin)+1 598 il_i2 = il_radius+(il_imax-ji)+1 599 IF( il_dim(1) == 1 )THEN 600 il_imin=ji 601 il_imax=ji 602 ! coef indices to be used 603 il_i1 = 1 604 il_i2 = 1 605 ENDIF 606 607 608 il_jmin=MAX(jj-il_radius,1) 609 il_jmax=MIN(jj+il_radius,il_shape(2)) 610 ! coef indices to be used 611 il_j1 = il_radius-(jj-il_jmin)+1 612 il_j2 = il_radius+(il_jmax-jj)+1 613 IF( il_dim(2) == 1 )THEN 614 il_jmin=jj 615 il_jmax=jj 616 ! coef indices to be used 617 il_j1 = 1 618 il_j2 = 1 619 ENDIF 620 621 il_kmin=MAX(jk-il_radius,1) 622 il_kmax=MIN(jk+il_radius,il_shape(3)) 623 ! coef indices to be used 624 il_k1 = il_radius-(jk-il_kmin)+1 625 il_k2 = il_radius+(il_kmax-jk)+1 626 IF( il_dim(3) == 1 )THEN 627 il_kmin=jk 628 il_kmax=jk 629 ! coef indices to be used 630 il_k1 = 1 631 il_k2 = 1 632 ENDIF 633 634 dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 635 & dd_value( il_imin:il_imax, & 636 & il_jmin:il_jmax, & 637 & il_kmin:il_kmax,jl ), dd_fill, il_radius, & 638 & dl_dfdx( il_imin:il_imax, & 639 & il_jmin:il_jmax, & 640 & il_kmin:il_kmax ), & 641 & dl_dfdy( il_imin:il_imax, & 642 & il_jmin:il_jmax, & 643 & il_kmin:il_kmax ), & 644 & dl_dfdz( il_imin:il_imax, & 645 & il_jmin:il_jmax, & 646 & il_kmin:il_kmax ), & 647 & dl_coef(il_i1:il_i2, & 648 & il_j1:il_j2, & 649 & il_k1:il_k2) ) 650 651 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 652 il_detect(ji,jj,jk)= 0 653 ll_iter=.FALSE. 851 654 ENDIF 852 655 … … 862 665 DEALLOCATE( dl_coef ) 863 666 864 il_iter=il_iter+1667 IF( ll_iter ) il_iter=il_iter+1 865 668 ENDDO 866 669 ENDDO … … 875 678 DO WHILE( ANY(il_detect(:,:,:)==1) ) 876 679 ! change extend value to minimize number of iteration 877 il_radius=id_radius+(il_iter/id_maxiter) 680 il_radius=id_radius+(il_iter-1) 681 ll_iter=.TRUE. 878 682 879 683 il_dim(1)=2*il_radius+1 … … 886 690 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 887 691 888 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), 889 & 1:il_dim(2), 890 & 1:il_dim(3), 692 dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1),& 693 & 1:il_dim(2),& 694 & 1:il_dim(3),& 891 695 & jl ) ) 892 696 893 697 DO jk=1,il_shape(3) 698 ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 894 699 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 895 700 DO jj=1,il_shape(2) … … 901 706 il_imin=MAX(ji-il_radius,1) 902 707 il_imax=MIN(ji+il_radius,il_shape(1)) 708 ! coef indices to be used 709 il_i1 = il_radius-(ji-il_imin)+1 710 il_i2 = il_radius+(il_imax-ji)+1 903 711 IF( il_dim(1) == 1 )THEN 904 712 il_imin=ji 905 713 il_imax=ji 714 ! coef indices to be used 715 il_i1 = 1 716 il_i2 = 1 906 717 ENDIF 907 718 908 719 il_jmin=MAX(jj-il_radius,1) 909 720 il_jmax=MIN(jj+il_radius,il_shape(2)) 721 ! coef indices to be used 722 il_j1 = il_radius-(jj-il_jmin)+1 723 il_j2 = il_radius+(il_jmax-jj)+1 910 724 IF( il_dim(2) == 1 )THEN 911 725 il_jmin=jj 912 726 il_jmax=jj 727 ! coef indices to be used 728 il_j1 = 1 729 il_j2 = 1 913 730 ENDIF 914 731 915 732 il_kmin=MAX(jk-il_radius,1) 916 733 il_kmax=MIN(jk+il_radius,il_shape(3)) 734 ! coef indices to be used 735 il_k1 = il_radius-(jk-il_kmin)+1 736 il_k2 = il_radius+(il_kmax-jk)+1 917 737 IF( il_dim(3) == 1 )THEN 918 738 il_kmin=jk 919 739 il_kmax=jk 740 ! coef indices to be used 741 il_k1 = 1 742 il_k2 = 1 920 743 ENDIF 921 744 … … 925 748 & il_kmin:il_kmax, & 926 749 & jl), dd_fill, il_radius, & 927 & dl_coef(:,:,:) ) 750 & dl_coef(il_i1:il_i2, & 751 & il_j1:il_j2, & 752 & il_k1:il_k2) ) 928 753 929 754 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 930 755 il_detect(ji,jj,jk)= 0 756 ll_iter=.FALSE. 757 ENDIF 758 759 ENDIF 760 761 ENDDO 762 ENDDO 763 ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 764 IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 765 DO jj=il_shape(2),1,-1 766 IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 767 DO ji=il_shape(1),1,-1 768 769 IF( il_detect(ji,jj,jk) == 1 )THEN 770 771 il_imin=MAX(ji-il_radius,1) 772 il_imax=MIN(ji+il_radius,il_shape(1)) 773 ! coef indices to be used 774 il_i1 = il_radius-(ji-il_imin)+1 775 il_i2 = il_radius+(il_imax-ji)+1 776 IF( il_dim(1) == 1 )THEN 777 il_imin=ji 778 il_imax=ji 779 ! coef indices to be used 780 il_i1 = 1 781 il_i2 = 1 782 ENDIF 783 784 il_jmin=MAX(jj-il_radius,1) 785 il_jmax=MIN(jj+il_radius,il_shape(2)) 786 ! coef indices to be used 787 il_j1 = il_radius-(jj-il_jmin)+1 788 il_j2 = il_radius+(il_jmax-jj)+1 789 IF( il_dim(2) == 1 )THEN 790 il_jmin=jj 791 il_jmax=jj 792 ! coef indices to be used 793 il_j1 = 1 794 il_j2 = 1 795 ENDIF 796 797 il_kmin=MAX(jk-il_radius,1) 798 il_kmax=MIN(jk+il_radius,il_shape(3)) 799 ! coef indices to be used 800 il_k1 = il_radius-(jk-il_kmin)+1 801 il_k2 = il_radius+(il_kmax-jk)+1 802 IF( il_dim(3) == 1 )THEN 803 il_kmin=jk 804 il_kmax=jk 805 ! coef indices to be used 806 il_k1 = 1 807 il_k2 = 1 808 ENDIF 809 810 dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 811 & dd_value( il_imin:il_imax, & 812 & il_jmin:il_jmax, & 813 & il_kmin:il_kmax, & 814 & jl), dd_fill, il_radius, & 815 & dl_coef(il_i1:il_i2, & 816 & il_j1:il_j2, & 817 & il_k1:il_k2) ) 818 819 IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 820 il_detect(ji,jj,jk)= 0 821 ll_iter=.FALSE. 931 822 ENDIF 932 823 … … 936 827 ENDDO 937 828 ENDDO 938 829 CALL logger_info(" EXTRAP 3D: "//& 830 & TRIM(fct_str(SUM(il_detect(:,:,:))))//& 831 & " point(s) to extrapolate " ) 832 939 833 DEALLOCATE( dl_coef ) 940 il_iter=il_iter+1834 IF( ll_iter ) il_iter=il_iter+1 941 835 ENDDO 942 836 ENDDO … … 946 840 947 841 END SUBROUTINE extrap__3D 948 !-------------------------------------------------------------------949 !> @brief950 !> This function compute derivative of 1D array.951 !>952 !> @details953 !> optionaly you could specify to take into account east west discontinuity954 !> (-180° 180° or 0° 360° for longitude variable)955 !>956 !> @author J.Paul957 !> - November, 2013- Initial Version958 !959 !> @param[in] dd_value 1D array of variable to be extrapolated960 !> @param[in] dd_fill FillValue of variable961 !> @param[in] ld_discont logical to take into account east west discontinuity962 !-------------------------------------------------------------------963 PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont )964 965 IMPLICIT NONE966 ! Argument967 REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value968 REAL(dp) , INTENT(IN) :: dd_fill969 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont970 971 ! function972 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D973 974 ! local variable975 INTEGER(i4) :: il_imin976 INTEGER(i4) :: il_imax977 INTEGER(i4), DIMENSION(1) :: il_shape978 979 REAL(dp) :: dl_min980 REAL(dp) :: dl_max981 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value982 983 LOGICAL :: ll_discont984 985 ! loop indices986 INTEGER(i4) :: ji987 988 INTEGER(i4) :: i1989 INTEGER(i4) :: i2990 !----------------------------------------------------------------991 ! init992 extrap_deriv_1D(:)=dd_fill993 994 ll_discont=.FALSE.995 IF( PRESENT(ld_discont) ) ll_discont=ld_discont996 997 il_shape(:)=SHAPE(dd_value(:))998 999 ALLOCATE( dl_value(3))1000 1001 ! compute derivative in i-direction1002 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 )THEN1008 i1=1 ; i2=31009 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN1010 i1=1 ; i2=21011 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN1012 i1=2 ; i2=31013 ENDIF1014 1015 dl_value(i1:i2)=dd_value(il_imin:il_imax)1016 IF( il_imin == 1 )THEN1017 dl_value(:)=EOSHIFT( dl_value(:), &1018 & DIM=1, &1019 & SHIFT=-1, &1020 & BOUNDARY=dl_value(1) )1021 ENDIF1022 IF( il_imax == il_shape(1) )THEN1023 dl_value(:)=EOSHIFT( dl_value(:), &1024 & DIM=1, &1025 & SHIFT=1, &1026 & BOUNDARY=dl_value(3))1027 ENDIF1028 1029 IF( ll_discont )THEN1030 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 )THEN1033 WHERE( dl_value(:) < 0._dp )1034 dl_value(:) = dl_value(:)+360._dp1035 END WHERE1036 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1037 WHERE( dl_value(:) > 180._dp )1038 dl_value(:) = dl_value(:)-180._dp1039 END WHERE1040 ENDIF1041 ENDIF1042 1043 IF( dl_value( 2) /= dd_fill .AND. & ! ji1044 & dl_value( 3) /= dd_fill .AND. & ! ji+11045 & dl_value( 1) /= dd_fill )THEN ! ji-11046 1047 extrap_deriv_1D(ji)=&1048 & ( dl_value(3) - dl_value(1) ) / &1049 & REAL( il_imax-il_imin ,dp)1050 1051 ENDIF1052 1053 ENDDO1054 1055 DEALLOCATE( dl_value )1056 1057 END FUNCTION extrap_deriv_1D1058 !-------------------------------------------------------------------1059 !> @brief1060 !> 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 !> @details1065 !> optionaly you could specify to take into account east west discontinuity1066 !> (-180° 180° or 0° 360° for longitude variable)1067 !>1068 !> @author J.Paul1069 !> - November, 2013- Initial Version1070 !1071 !> @param[in] dd_value 2D array of variable to be extrapolated1072 !> @param[in] dd_fill FillValue of variable1073 !> @param[in] cd_dim compute derivative on first (I) or second (J) dimension1074 !> @param[in] ld_discont logical to take into account east west discontinuity1075 !-------------------------------------------------------------------1076 FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont )1077 1078 IMPLICIT NONE1079 ! Argument1080 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value1081 REAL(dp) , INTENT(IN) :: dd_fill1082 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim1083 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont1084 1085 ! function1086 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &1087 & SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D1088 1089 ! local variable1090 INTEGER(i4) :: il_imin1091 INTEGER(i4) :: il_imax1092 INTEGER(i4) :: il_jmin1093 INTEGER(i4) :: il_jmax1094 INTEGER(i4), DIMENSION(2) :: il_shape1095 1096 REAL(dp) :: dl_min1097 REAL(dp) :: dl_max1098 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value1099 1100 LOGICAL :: ll_discont1101 1102 ! loop indices1103 INTEGER(i4) :: ji1104 INTEGER(i4) :: jj1105 1106 INTEGER(i4) :: i11107 INTEGER(i4) :: i21108 1109 INTEGER(i4) :: j11110 INTEGER(i4) :: j21111 !----------------------------------------------------------------1112 ! init1113 extrap_deriv_2D(:,:)=dd_fill1114 1115 ll_discont=.FALSE.1116 IF( PRESENT(ld_discont) ) ll_discont=ld_discont1117 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-direction1126 DO ji=1,il_shape(1)1127 1128 ! init1129 dl_value(:,:)=dd_fill1130 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 )THEN1135 i1=1 ; i2=31136 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN1137 i1=1 ; i2=21138 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN1139 i1=2 ; i2=31140 ENDIF1141 1142 dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:)1143 IF( il_imin == 1 )THEN1144 dl_value(:,:)=EOSHIFT( dl_value(:,:), &1145 & DIM=1, &1146 & SHIFT=-1, &1147 & BOUNDARY=dl_value(1,:) )1148 ENDIF1149 IF( il_imax == il_shape(1) )THEN1150 dl_value(:,:)=EOSHIFT( dl_value(:,:), &1151 & DIM=1, &1152 & SHIFT=1, &1153 & BOUNDARY=dl_value(3,:))1154 ENDIF1155 1156 IF( ll_discont )THEN1157 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 )THEN1160 WHERE( dl_value(:,:) < 0_dp )1161 dl_value(:,:) = dl_value(:,:)+360._dp1162 END WHERE1163 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1164 WHERE( dl_value(:,:) > 180 )1165 dl_value(:,:) = dl_value(:,:)-180._dp1166 END WHERE1167 ENDIF1168 ENDIF1169 1170 WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji1171 & dl_value(3,:) /= dd_fill .AND. & ! ji+11172 & dl_value(1,:) /= dd_fill ) ! ji-11173 1174 extrap_deriv_2D(ji,:)=&1175 & ( dl_value(3,:) - dl_value(1,:) ) / &1176 & REAL( il_imax-il_imin,dp)1177 1178 END WHERE1179 1180 ENDDO1181 1182 CASE('J')1183 1184 ALLOCATE( dl_value(il_shape(1),3) )1185 ! compute derivative in j-direction1186 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 )THEN1192 j1=1 ; j2=31193 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN1194 j1=1 ; j2=21195 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN1196 j1=2 ; j2=31197 ENDIF1198 1199 dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax)1200 IF( il_jmin == 1 )THEN1201 dl_value(:,:)=EOSHIFT( dl_value(:,:), &1202 & DIM=2, &1203 & SHIFT=-1, &1204 & BOUNDARY=dl_value(:,1))1205 ENDIF1206 IF( il_jmax == il_shape(2) )THEN1207 dl_value(:,:)=EOSHIFT( dl_value(:,:), &1208 & DIM=2, &1209 & SHIFT=1, &1210 & BOUNDARY=dl_value(:,3))1211 ENDIF1212 1213 IF( ll_discont )THEN1214 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 )THEN1217 WHERE( dl_value(:,:) < 0_dp )1218 dl_value(:,:) = dl_value(:,:)+360._dp1219 END WHERE1220 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1221 WHERE( dl_value(:,:) > 180 )1222 dl_value(:,:) = dl_value(:,:)-180._dp1223 END WHERE1224 ENDIF1225 ENDIF1226 1227 WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj1228 & dl_value(:, 3) /= dd_fill .AND. & ! jj+11229 & dl_value(:, 1) /= dd_fill ) ! jj-11230 1231 extrap_deriv_2D(:,jj)=&1232 & ( dl_value(:,3) - dl_value(:,1) ) / &1233 & REAL(il_jmax-il_jmin,dp)1234 1235 END WHERE1236 1237 ENDDO1238 1239 END SELECT1240 1241 DEALLOCATE( dl_value )1242 1243 END FUNCTION extrap_deriv_2D1244 !-------------------------------------------------------------------1245 !> @brief1246 !> 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 !> @details1251 !> optionaly you could specify to take into account east west discontinuity1252 !> (-180° 180° or 0° 360° for longitude variable)1253 !>1254 !> @author J.Paul1255 !> - November, 2013- Initial Version1256 !1257 !> @param[inout] dd_value 3D array of variable to be extrapolated1258 !> @param[in] dd_fill FillValue of variable1259 !> @param[in] cd_dim compute derivative on first (I) second (J) or third (K) dimension1260 !> @param[in] ld_discont logical to take into account east west discontinuity1261 !-------------------------------------------------------------------1262 PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont )1263 1264 IMPLICIT NONE1265 ! Argument1266 REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value1267 REAL(dp) , INTENT(IN) :: dd_fill1268 CHARACTER(LEN=*) , INTENT(IN) :: cd_dim1269 LOGICAL , INTENT(IN), OPTIONAL :: ld_discont1270 1271 ! function1272 REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &1273 & SIZE(dd_value,DIM=2), &1274 & SIZE(dd_value,DIM=3)) :: extrap_deriv_3D1275 1276 ! local variable1277 INTEGER(i4) :: il_imin1278 INTEGER(i4) :: il_imax1279 INTEGER(i4) :: il_jmin1280 INTEGER(i4) :: il_jmax1281 INTEGER(i4) :: il_kmin1282 INTEGER(i4) :: il_kmax1283 INTEGER(i4), DIMENSION(3) :: il_shape1284 1285 REAL(dp) :: dl_min1286 REAL(dp) :: dl_max1287 REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value1288 1289 LOGICAL :: ll_discont1290 1291 ! loop indices1292 INTEGER(i4) :: ji1293 INTEGER(i4) :: jj1294 INTEGER(i4) :: jk1295 1296 INTEGER(i4) :: i11297 INTEGER(i4) :: i21298 1299 INTEGER(i4) :: j11300 INTEGER(i4) :: j21301 1302 INTEGER(i4) :: k11303 INTEGER(i4) :: k21304 !----------------------------------------------------------------1305 ! init1306 extrap_deriv_3D(:,:,:)=dd_fill1307 1308 ll_discont=.FALSE.1309 IF( PRESENT(ld_discont) ) ll_discont=ld_discont1310 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-direction1320 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 )THEN1326 i1=1 ; i2=31327 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN1328 i1=1 ; i2=21329 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN1330 i1=2 ; i2=31331 ENDIF1332 1333 dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:)1334 IF( il_imin == 1 )THEN1335 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1336 & DIM=1, &1337 & SHIFT=-1, &1338 & BOUNDARY=dl_value(1,:,:) )1339 ENDIF1340 IF( il_imax == il_shape(1) )THEN1341 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1342 & DIM=1, &1343 & SHIFT=1, &1344 & BOUNDARY=dl_value(3,:,:))1345 ENDIF1346 1347 IF( ll_discont )THEN1348 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 )THEN1351 WHERE( dl_value(:,:,:) < 0_dp )1352 dl_value(:,:,:) = dl_value(:,:,:)+360._dp1353 END WHERE1354 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1355 WHERE( dl_value(:,:,:) > 180 )1356 dl_value(:,:,:) = dl_value(:,:,:)-180._dp1357 END WHERE1358 ENDIF1359 ENDIF1360 1361 WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji1362 & dl_value(3,:,:) /= dd_fill .AND. & !ji+11363 & dl_value(1,:,:) /= dd_fill ) !ji-11364 1365 extrap_deriv_3D(ji,:,:)= &1366 & ( dl_value(3,:,:) - dl_value(1,:,:) ) / &1367 & REAL( il_imax-il_imin ,dp)1368 1369 END WHERE1370 1371 ENDDO1372 1373 CASE('J')1374 1375 ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) )1376 ! compute derivative in j-direction1377 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 )THEN1383 j1=1 ; j2=31384 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN1385 j1=1 ; j2=21386 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN1387 j1=2 ; j2=31388 ENDIF1389 1390 dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:)1391 IF( il_jmin == 1 )THEN1392 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1393 & DIM=2, &1394 & SHIFT=-1, &1395 & BOUNDARY=dl_value(:,1,:) )1396 ENDIF1397 IF( il_jmax == il_shape(2) )THEN1398 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1399 & DIM=2, &1400 & SHIFT=1, &1401 & BOUNDARY=dl_value(:,3,:))1402 ENDIF1403 1404 IF( ll_discont )THEN1405 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 )THEN1408 WHERE( dl_value(:,:,:) < 0_dp )1409 dl_value(:,:,:) = dl_value(:,:,:)+360._dp1410 END WHERE1411 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1412 WHERE( dl_value(:,:,:) > 180 )1413 dl_value(:,:,:) = dl_value(:,:,:)-180._dp1414 END WHERE1415 ENDIF1416 ENDIF1417 1418 WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj1419 & dl_value(:, 3,:) /= dd_fill .AND. & ! jj+11420 & dl_value(:, 1,:) /= dd_fill ) ! jj-11421 1422 extrap_deriv_3D(:,jj,:)=&1423 & ( dl_value(:,3,:) - dl_value(:,1,:) ) / &1424 & REAL( il_jmax - il_jmin ,dp)1425 1426 END WHERE1427 1428 ENDDO1429 1430 CASE('K')1431 ! compute derivative in k-direction1432 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 )THEN1438 k1=1 ; k2=31439 ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN1440 k1=1 ; k2=21441 ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN1442 k1=2 ; k2=31443 ENDIF1444 1445 dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax)1446 IF( il_kmin == 1 )THEN1447 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1448 & DIM=3, &1449 & SHIFT=-1, &1450 & BOUNDARY=dl_value(:,:,1) )1451 ENDIF1452 IF( il_kmax == il_shape(3) )THEN1453 dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &1454 & DIM=3, &1455 & SHIFT=1, &1456 & BOUNDARY=dl_value(:,:,3))1457 ENDIF1458 1459 IF( ll_discont )THEN1460 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 )THEN1463 WHERE( dl_value(:,:,:) < 0_dp )1464 dl_value(:,:,:) = dl_value(:,:,:)+360._dp1465 END WHERE1466 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1467 WHERE( dl_value(:,:,:) > 180 )1468 dl_value(:,:,:) = dl_value(:,:,:)-180._dp1469 END WHERE1470 ENDIF1471 ENDIF1472 1473 WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk1474 & dl_value(:,:, 3) /= dd_fill .AND. & ! jk+11475 & dl_value(:,:, 1) /= dd_fill ) ! jk-11476 1477 extrap_deriv_3D(:,:,jk)=&1478 & ( dl_value(:,:,3) - dl_value(:,:,1) ) / &1479 & REAL( il_kmax-il_kmin,dp)1480 1481 END WHERE1482 1483 ENDDO1484 1485 END SELECT1486 1487 DEALLOCATE( dl_value )1488 1489 END FUNCTION extrap_deriv_3D1490 842 !------------------------------------------------------------------- 1491 843 !> @brief … … 1493 845 !> 1494 846 !> @details 1495 !> coefficients are "grid distance" to the center of the box choosed to compute1496 !> extrapolation.847 !> coefficients are "grid distance" to the center of the box 848 !> choosed to compute extrapolation. 1497 849 !> 1498 850 !> @author J.Paul 1499 !> - November, 2013- Initial Version 851 !> @date November, 2013 - Initial Version 852 !> @date July, 2015 853 !> - decrease weight of third dimension 1500 854 ! 1501 855 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1544 898 1545 899 ! compute distance 900 ! "vertical weight" is lower than horizontal 1546 901 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 1547 902 & (jj-il_jmid)**2 + & 1548 & 903 & 3*(jk-il_kmid)**2 1549 904 1550 905 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 1569 924 !> 1570 925 !> @author J.Paul 1571 !> - November, 2013- Initial Version926 !> @date November, 2013 - Initial Version 1572 927 !> 1573 928 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1658 1013 !> 1659 1014 !> @author J.Paul 1660 !> - November, 2013- Initial Version 1015 !> @date November, 2013 - Initial Version 1016 !> @date July, 2015 1017 !> - decrease weight of third dimension 1661 1018 ! 1662 1019 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1705 1062 1706 1063 ! compute distance 1064 ! "vertical weight" is lower than horizontal 1707 1065 dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 1708 1066 & (jj-il_jmid)**2 + & 1709 & 1067 & 3*(jk-il_kmid)**2 1710 1068 1711 1069 IF( dl_dist(ji,jj,jk) /= 0 )THEN … … 1732 1090 !> 1733 1091 !> @author J.Paul 1734 !> - November, 2013- Initial Version1092 !> @date November, 2013 - Initial Version 1735 1093 ! 1736 1094 !> @param[in] dd_value 3D array of variable to be extrapolated … … 1763 1121 INTEGER(i4) :: jj 1764 1122 INTEGER(i4) :: jk 1765 1766 1123 !---------------------------------------------------------------- 1767 1124 … … 1793 1150 ENDDO 1794 1151 ENDDO 1152 1795 1153 1796 1154 ! return value … … 1815 1173 !> 1816 1174 !> @author J.Paul 1817 !> - November, 2013-Initial version1175 !> @date November, 2013 - Initial version 1818 1176 ! 1819 1177 !> @param[inout] td_var variable … … 1917 1275 !> 1918 1276 !> @author J.Paul 1919 !> - November, 2013-Initial version1277 !> @date November, 2013 - Initial version 1920 1278 !> 1921 1279 !> @param[inout] td_var variable
Note: See TracChangeset
for help on using the changeset viewer.