- Timestamp:
- 2015-07-16T13:55:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90
r4213 r5602 8 8 !> @brief 9 9 !> This module manage interpolation on regular grid. 10 !> @note It is used to work on ORCA grid, as we work only with grid indices. 10 !> 11 !> @details Interpolation method to be used is specify inside variable 12 !> strcuture, as array of string character.<br/> 13 !> - td_var\%c_interp(1) string character is the interpolation name choose between: 14 !> - 'nearest' 15 !> - 'cubic ' 16 !> - 'linear ' 17 !> - td_var\%c_interp(2) string character is an operation to be used 18 !> on interpolated value.<br/> 19 !> operation have to be mulitplication '*' or division '/'.<br/> 20 !> coefficient have to be refinement factor following i-direction 'rhoi', 21 !> j-direction 'rhoj', or k-direction 'rhok'.<br/> 22 !> 23 !> Examples: '*rhoi', '/rhoj'. 24 !> 25 !> @note Those informations are read from namelist or variable configuration file (default).<br/> 26 !> Interplation method could be specify for each variable in namelist _namvar_, 27 !> defining string character _cn\_varinfo_.<br/> 28 !> Example: 29 !> - cn_varinfo='varname1:cubic/rhoi', 'varname2:linear' 30 !> 31 !> to create mixed grid (with coarse grid point needed to compute 32 !> interpolation):<br/> 33 !> @code 34 !> CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] ) 35 !> @endcode 36 !> - td_var is coarse grid variable (should be extrapolated) 37 !> - td_mix is mixed grid variable structure [output] 38 !> - id_rho is array of refinment factor [optional] 39 !> 40 !> to detected point to be interpolated:<br/> 41 !> @code 42 !> il_detect(:,:,:)=interp_detect( td_mix [,id_rho] ) 43 !> @endcode 44 !> - il_detect(:,:,:) is 3D array of detected point to be interpolated 45 !> - td_mix is mixed grid variable 46 !> - id_rho is array of refinement factor [optional] 47 !> 48 !> to interpolate variable value:<br/> 49 !> @code 50 !> CALL interp_fill_value( td_var [,id_rho] [,id_offset] ) 51 !> @endcode 52 !> - td_var is variable structure 53 !> - id_rho is array of refinement factor [optional] 54 !> - id_offset is array of offset between fine and coarse grid [optional] 55 !> 56 !> to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/> 57 !> @code 58 !> CALL interp_clean_mixed_grid( td_mix, td_var, id_rho ) 59 !> @endcode 60 !> - td_mix is mixed grid variable structure 61 !> - td_var is variable structure [output] 62 !> - id_rho is array of refinement factor [optional] 63 !> - id_offset is array of offset between fine and coarse grid [optional] 64 !> 65 !> @note It use to work on ORCA grid, as we work only with grid indices. 11 66 !> 12 67 !> @warning due to the use of second derivative when using cubic interpolation 13 !> you should add 2 extrabands68 !> you should add at least 2 extrabands. 14 69 !> 15 !> @details16 70 !> @author 17 71 !> J.Paul 18 72 ! REVISION HISTORY: 19 !> @date Nov, 2013 - Initial Version 20 ! 73 !> @date November, 2013 - Initial Version 74 !> @date September, 2014 75 !> - add header 76 !> - use interpolation method modules 77 !> 21 78 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 22 !> @todo23 !> - interp 3D24 !> - see issue when fill value is zero for check cubic_fill..25 79 !---------------------------------------------------------------------- 26 80 MODULE interp … … 29 83 USE global ! global variable 30 84 USE kind ! F90 kind parameter 31 USE logger 85 USE logger ! log file manager 32 86 USE fct ! basic useful function 33 !USE date ! date manager87 USE date ! date manager 34 88 USE att ! attribute manager 35 89 USE dim ! dimension manager 36 90 USE var ! variable manager 37 ! USE file ! file manager38 ! USE iom ! I/O manager39 ! USE dom ! domain manager40 91 USE grid ! grid manager 41 92 USE extrap ! extrapolation manager 42 ! USE interp ! interpolation manager 43 ! USE filter ! filter manager 44 ! USE mpp ! MPP manager 45 ! USE iom_mpp ! MPP I/O manager 93 USE interp_cubic ! cubic interpolation manager 94 USE interp_linear ! linear interpolation manager 95 USE interp_nearest ! nearest interpolation manager 46 96 47 97 IMPLICIT NONE 48 PRIVATE49 98 ! NOTE_avoid_public_variables_if_possible 50 99 … … 53 102 ! function and subroutine 54 103 PUBLIC :: interp_detect !< detected point to be interpolated 55 PUBLIC :: interp_fill_value !< interpolate value over detectected point104 PUBLIC :: interp_fill_value !< interpolate value 56 105 PUBLIC :: interp_create_mixed_grid !< create mixed grid 57 106 PUBLIC :: interp_clean_mixed_grid !< clean mixed grid 58 107 59 PRIVATE :: interp__detect !< detected point to be interpolated 60 PRIVATE :: interp__detect_wrapper !< detected point to be interpolated 61 PRIVATE :: interp__fill_value_wrapper !< interpolate value over detectected point 62 PRIVATE :: interp__fill_value !< interpolate value over detectected point 63 PRIVATE :: interp__clean_even_grid !< clean even mixed grid 64 PRIVATE :: interp__del_offset !< remove offset from interpolated grid 65 PRIVATE :: interp__check_method !< check if interpolation method available 66 ! PRIVATE :: interp__3D !< interpolate 3D grid 67 PRIVATE :: interp__2D !< interpolate 2D grid 68 PRIVATE :: interp__1D !< interpolate 1D grid 69 PRIVATE :: interp__2D_cubic_coef !< compute coefficient for bicubic interpolation 70 PRIVATE :: interp__2D_cubic_fill !< compute bicubic interpolation 71 PRIVATE :: interp__2D_linear_coef !< compute coefficient for bilinear interpolation 72 PRIVATE :: interp__2D_linear_fill !< compute bilinear interpolation 73 PRIVATE :: interp__2D_nearest_fill !< compute nearest interpolation 74 PRIVATE :: interp__1D_cubic_coef !< compute coefficient for cubic interpolation 75 PRIVATE :: interp__1D_cubic_fill !< compute cubic interpolation 76 PRIVATE :: interp__1D_linear_coef !< compute coefficient for linear interpolation 77 PRIVATE :: interp__1D_linear_fill !< compute linear interpolation 78 PRIVATE :: interp__1D_nearest_fill !< compute nearest interpolation 79 ! PRIVATE :: interp__longitude 108 PRIVATE :: interp__detect ! detected point to be interpolated 109 PRIVATE :: interp__detect_wrapper ! detected point to be interpolated 110 PRIVATE :: interp__fill_value_wrapper ! interpolate value over detectected point 111 PRIVATE :: interp__fill_value ! interpolate value over detectected point 112 PRIVATE :: interp__clean_even_grid ! clean even mixed grid 113 PRIVATE :: interp__check_method ! check if interpolation method available 80 114 81 115 TYPE TINTERP 82 !CHARACTER(LEN=lc) :: c_name = 'unknown' !< interpolation method name83 !CHARACTER(LEN=lc) :: c_factor = 'unknown' !< interpolation factor84 !CHARACTER(LEN=lc) :: c_divisor = 'unknown' !< interpolation divisor85 116 CHARACTER(LEN=lc) :: c_name = '' !< interpolation method name 86 117 CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor … … 89 120 90 121 INTERFACE interp_detect 91 MODULE PROCEDURE interp__detect_wrapper !< detected point to be interpolated122 MODULE PROCEDURE interp__detect_wrapper 92 123 END INTERFACE interp_detect 93 124 94 125 INTERFACE interp_fill_value 95 MODULE PROCEDURE interp__fill_value_wrapper !< detected point to be interpolated126 MODULE PROCEDURE interp__fill_value_wrapper 96 127 END INTERFACE interp_fill_value 97 128 … … 99 130 !------------------------------------------------------------------- 100 131 !> @brief 101 !> This function check if interpolation method available.132 !> This function check if interpolation method is available. 102 133 !> 103 134 !> @details 135 !> check if name of interpolation method is present in global list of string 136 !> character cp_interp_list (see global.f90). 104 137 !> 105 138 !> @author J.Paul 106 !> - Nov , 2013- Initial Version139 !> - November, 2013- Initial Version 107 140 ! 108 !> @param[in] cd_method :interpolation method141 !> @param[in] cd_method interpolation method 109 142 !> @return 110 !> @todo see extrap_detect 111 !------------------------------------------------------------------- 112 !> @code 143 !------------------------------------------------------------------- 113 144 FUNCTION interp__check_method( cd_method ) 114 145 IMPLICIT NONE … … 130 161 131 162 interp__check_method=.FALSE. 132 DO ji=1,i g_ninterp133 cl_interp=fct_lower(c g_interp_list(ji))163 DO ji=1,ip_ninterp 164 cl_interp=fct_lower(cp_interp_list(ji)) 134 165 IF( TRIM(cl_interp) == TRIM(cl_method) )THEN 135 166 interp__check_method=.TRUE. … … 139 170 140 171 END FUNCTION interp__check_method 141 !> @endcode142 172 !------------------------------------------------------------------- 143 173 !> @brief … … 149 179 !> 150 180 !> @author J.Paul 151 !> - Nov, 2013- Initial Version 152 ! 153 !> @param[in] td_mix : mixed grid variable (to interpolate) 154 !> @param[in] id_rho : table of refinement factor 155 !> @return table of detected point to be interpolated 156 !------------------------------------------------------------------- 157 !> @code 181 !> - November, 2013- Initial Version 182 !> 183 !> @param[in] td_mix mixed grid variable (to interpolate) 184 !> @param[in] id_rho array of refinement factor 185 !> @return 3D array of detected point to be interpolated 186 !------------------------------------------------------------------- 158 187 FUNCTION interp__detect_wrapper( td_mix, id_rho ) 159 188 IMPLICIT NONE … … 168 197 169 198 ! local variable 170 171 199 ! loop indices 172 200 !---------------------------------------------------------------- … … 208 236 209 237 END FUNCTION interp__detect_wrapper 210 !> @endcode211 238 !------------------------------------------------------------------- 212 239 !> @brief … … 217 244 !> 218 245 !> @author J.Paul 219 !> - Nov , 2013- Initial Version246 !> - November, 2013- Initial Version 220 247 ! 221 !> @param[in] td_mix : mixed grid variable (to interpolate) 222 !> @param[in] id_rho : table of refinement factor 223 !> @return table of detected point to be interpolated 224 !------------------------------------------------------------------- 225 !> @code 248 !> @param[in] td_mix mixed grid variable (to interpolate) 249 !> @param[in] id_rho array of refinement factor 250 !> @return 3D array of detected point to be interpolated 251 !------------------------------------------------------------------- 226 252 FUNCTION interp__detect( td_mix, id_rho ) 227 253 IMPLICIT NONE … … 253 279 ALLOCATE( il_rho(ip_maxdim) ) 254 280 il_rho(:)=1 255 IF( PRESENT(id_rho) ) il_rho( :)=id_rho(:)281 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 256 282 257 283 ! special case for even refinement on ARAKAWA-C grid … … 289 315 il_dim(:)=td_mix%t_dim(1:3)%i_len 290 316 317 ! init 291 318 interp__detect(:,:,:)=1 292 319 … … 297 324 298 325 ! do not compute point near fill value 299 FORALL( ji=1:il_dim(1):il_rho(jp_I), & 300 & jj=1:il_dim(2):il_rho(jp_J), & 301 & jk=1:il_dim(3):il_rho(jp_K), & 302 & td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill ) 303 304 ! i-direction 305 interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),& 306 & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 307 & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 308 ! j-direction 309 interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 310 & MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),& 311 & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 312 ! k-direction 313 interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 314 & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 315 & MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0 316 317 END FORALL 326 DO jk=1,il_dim(3),il_rho(jp_K) 327 DO jj=1,il_dim(2),il_rho(jp_J) 328 DO ji=1,il_dim(1),il_rho(jp_I) 329 330 IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN 331 332 ! i-direction 333 interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),& 334 & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 335 & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 336 ! j-direction 337 interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 338 & MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),& 339 & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 340 ! k-direction 341 interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 342 & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 343 & MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0 344 345 ENDIF 346 347 ENDDO 348 ENDDO 349 ENDDO 350 351 DEALLOCATE( il_rho ) 318 352 319 353 END FUNCTION interp__detect 320 !> @endcode321 354 !------------------------------------------------------------------- 322 355 !> @brief … … 330 363 !> 331 364 !> @author J.Paul 332 !> - Nov, 2013- Initial Version 333 ! 334 !> @param[in] td_var : coarse grid variable (should be extrapolated) 335 !> @param[out] td_mix : mixed grid variable 336 !> @param[in] id_rho : table of refinment factor 337 !> 338 !> @todo 339 !------------------------------------------------------------------- 340 !> @code 365 !> - November, 2013- Initial Version 366 !> 367 !> @param[in] td_var coarse grid variable (should be extrapolated) 368 !> @param[out] td_mix mixed grid variable 369 !> @param[in] id_rho array of refinment factor (default 1) 370 !------------------------------------------------------------------- 341 371 SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho ) 342 372 IMPLICIT NONE … … 359 389 ALLOCATE(il_rho(ip_maxdim)) 360 390 il_rho(:)=1 361 IF( PRESENT(id_rho) ) il_rho( :)=id_rho(:)391 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 362 392 363 393 ! special case for even refinement on ARAKAWA-C grid … … 381 411 382 412 ! copy variable 383 td_mix= td_var413 td_mix=var_copy(td_var) 384 414 385 415 ! compute new dimension length … … 408 438 & td_var%d_value(:,:,:,:) 409 439 440 DEALLOCATE(il_rho) 441 410 442 END SUBROUTINE interp_create_mixed_grid 411 !> @endcode412 443 !------------------------------------------------------------------- 413 444 !> @brief … … 418 449 !> 419 450 !> @author J.Paul 420 !> - Nov, 2013- Initial Version 421 ! 422 !> @param[in] td_mix : mixed grid variable 423 !> @param[in] id_rho : table of refinment factor 424 !> 425 !> @todo 426 !------------------------------------------------------------------- 427 !> @code 451 !> - November, 2013- Initial Version 452 !> 453 !> @param[inout] td_mix mixed grid variable 454 !> @param[in] id_rho array of refinment factor 455 !------------------------------------------------------------------- 428 456 SUBROUTINE interp__clean_even_grid( td_mix, id_rho ) 429 457 IMPLICIT NONE 430 458 ! Argument 431 TYPE(TVAR) , INTENT(INOUT) :: td_mix459 TYPE(TVAR) , INTENT(INOUT) :: td_mix 432 460 INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho 433 461 … … 468 496 END SELECT 469 497 470 ! copy variable471 tl_mix=td_mix472 473 498 ! remove some point only if refinement in some direction is even 474 499 IF( ANY(ll_even(:)) )THEN 500 501 ! copy variable 502 tl_mix=var_copy(td_mix) 475 503 476 504 ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, & … … 554 582 555 583 td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), & 556 & SHAPE=td_mix%t_dim(:)%i_len ) 584 & SHAPE=(/td_mix%t_dim(1)%i_len, & 585 & td_mix%t_dim(2)%i_len, & 586 & td_mix%t_dim(3)%i_len, & 587 & td_mix%t_dim(4)%i_len/) ) 557 588 558 589 DEALLOCATE( dl_vect ) … … 562 593 DEALLOCATE( ll_mask ) 563 594 595 ! clean 596 CALL var_clean(tl_mix) 597 564 598 ENDIF 565 599 566 CALL var_clean(tl_mix) 600 ! clean 601 DEALLOCATE(il_rho) 567 602 568 603 END SUBROUTINE interp__clean_even_grid 569 !> @endcode570 604 !------------------------------------------------------------------- 571 605 !> @brief 572 !> This subroutine save interpolated value over domain. 573 !> And so remove points added on mixed grid 574 !> to compute interpolation 606 !> This subroutine remove points added on mixed grid 607 !> to compute interpolation. And save interpolated value over domain. 575 608 !> 576 609 !> @details 577 610 !> 578 611 !> @author J.Paul 579 !> - Nov, 2013- Initial Version 580 ! 581 !> @param[in] td_var : table of mixed grid variable (to interpolate) 582 !> @param[in] id_rho : table of refinement factor 583 !> @return 584 !------------------------------------------------------------------- 585 !> @code 612 !> - November, 2013- Initial Version 613 !> @date September, 2014 614 !> - use offset to save useful domain 615 !> 616 !> @param[in] td_mix mixed grid variable structure 617 !> @param[out] td_var variable structure 618 !> @param[in] id_rho array of refinement factor (default 1) 619 !> @param[in] id_offset 2D array of offset between fine and coarse grid 620 !------------------------------------------------------------------- 586 621 SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, & 587 & id_rho, & 588 & id_offset ) 622 & id_rho, id_offset ) 589 623 IMPLICIT NONE 590 624 ! Argument 591 625 TYPE(TVAR) , INTENT(IN ) :: td_mix 592 626 TYPE(TVAR) , INTENT( OUT) :: td_var 593 INTEGER(I4), DIMENSION(:) , INTENT(IN ) , OPTIONAL:: id_rho594 INTEGER(I4), DIMENSION( :,:), INTENT(IN ), OPTIONAL:: id_offset627 INTEGER(I4), DIMENSION(:) , INTENT(IN ) :: id_rho 628 INTEGER(I4), DIMENSION(2,2), INTENT(IN ) :: id_offset 595 629 596 630 ! local variable … … 605 639 INTEGER(i4) :: il_jmax1 606 640 607 INTEGER(i4), DIMENSION(2,2) :: il_offset608 609 641 REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value 610 642 … … 613 645 ! loop indices 614 646 !---------------------------------------------------------------- 615 il_offset(:,:)=0616 IF( PRESENT(id_offset) )THEN617 IF( ANY( SHAPE(id_offset(:,:)) /= SHAPE(il_offset(:,:)) ) )THEN618 CALL logger_error("INTERP CLEAN MIXED GRID: invalid dimension of"//&619 & " offset table")620 ELSE621 il_offset(:,:)=id_offset(:,:)622 ENDIF623 ENDIF624 647 ! copy mixed variable in temporary structure 625 tl_mix= td_mix648 tl_mix=var_copy(td_mix) 626 649 627 650 ! remove unusefull points over mixed grid for even refinement … … 629 652 630 653 ! copy cleaned mixed variable 631 td_var=tl_mix 632 IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 633 634 ! delete table of value 635 CALL var_del_value(td_var) 636 637 ! compute domain indices 638 il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len 639 il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len 640 641 il_imin1=il_imin0+il_offset(1,1) 642 il_jmin1=il_jmin0+il_offset(2,1) 643 644 il_imax1=il_imax0-il_offset(1,2) 645 il_jmax1=il_jmax0-il_offset(2,2) 646 647 SELECT CASE(TRIM(td_var%c_point)) 648 CASE('U') 649 il_imin1=il_imin0+(il_offset(1,1)-1) 650 il_imax1=il_imax0-(il_offset(1,2)+1) 651 CASE('V') 652 il_jmin1=il_jmin0+(il_offset(2,1)-1) 653 il_jmax1=il_jmax0-(il_offset(2,2)+1) 654 CASE('F') 655 il_imin1=il_imin0+(il_offset(1,1)-1) 656 il_imax1=il_imax0-(il_offset(1,2)+1) 657 658 il_jmin1=il_jmin0+(il_offset(2,1)-1) 659 il_jmax1=il_jmax0-(il_offset(2,2)+1) 660 END SELECT 661 662 ! compute new dimension 663 td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 664 td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 654 td_var=var_copy(tl_mix) 655 656 ! delete array of value 657 CALL var_del_value(td_var) 658 659 ! compute domain indices in i-direction 660 il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len 665 661 666 ALLOCATE(dl_value(td_var%t_dim(1)%i_len, & 667 & td_var%t_dim(2)%i_len, & 668 & td_var%t_dim(3)%i_len, & 669 & td_var%t_dim(4)%i_len) ) 670 671 dl_value( 1:td_var%t_dim(1)%i_len, & 672 & :,:,:) = tl_mix%d_value( il_imin1:il_imax1, & 673 & il_jmin1:il_jmax1, & 674 & :, : ) 675 676 ! add variable value 677 CALL var_add_value(td_var,dl_value(:,:,:,:)) 678 679 ! save variable type 680 td_var%i_type=tl_mix%i_type 681 682 DEALLOCATE(dl_value) 662 IF( td_var%t_dim(1)%l_use )THEN 663 il_imin1=il_imin0+id_offset(Jp_I,1) 664 il_imax1=il_imax0-id_offset(Jp_I,2) 665 ELSE 666 667 il_imin1=il_imin0 668 il_imax1=il_imax0 683 669 684 670 ENDIF 685 671 672 ! compute domain indices in j-direction 673 il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len 674 675 IF( td_var%t_dim(2)%l_use )THEN 676 il_jmin1=il_jmin0+id_offset(Jp_J,1) 677 il_jmax1=il_jmax0-id_offset(Jp_J,2) 678 ELSE 679 680 il_jmin1=il_jmin0 681 il_jmax1=il_jmax0 682 683 ENDIF 684 685 ! compute new dimension 686 td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 687 td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 688 689 ALLOCATE(dl_value(td_var%t_dim(1)%i_len, & 690 & td_var%t_dim(2)%i_len, & 691 & td_var%t_dim(3)%i_len, & 692 & td_var%t_dim(4)%i_len) ) 693 694 dl_value( 1:td_var%t_dim(1)%i_len, & 695 & 1:td_var%t_dim(2)%i_len, & 696 & :,:) = tl_mix%d_value( il_imin1:il_imax1, & 697 & il_jmin1:il_jmax1, & 698 & :, : ) 699 700 ! add variable value 701 CALL var_add_value(td_var,dl_value(:,:,:,:)) 702 703 DEALLOCATE(dl_value) 704 705 ! clean 686 706 CALL var_clean(tl_mix) 687 707 688 708 END SUBROUTINE interp_clean_mixed_grid 689 !> @endcode690 709 !------------------------------------------------------------------- 691 710 !> @brief 692 !> This subroutine save interpolated value over domain. 693 !> And so remove points added on mixed grid 694 !> to compute interpolation 695 !> 696 !> @details 697 !> 698 !> @author J.Paul 699 !> - Nov, 2013- Initial Version 700 ! 701 !> @param[in] td_var : table of mixed grid variable (to interpolate) 702 !> @param[in] id_offset : table of offset 703 !> @return 704 !------------------------------------------------------------------- 705 !> @code 706 SUBROUTINE interp__del_offset( td_var, id_offset ) 707 IMPLICIT NONE 708 ! Argument 709 TYPE(TVAR) , INTENT(INOUT) :: td_var 710 INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset 711 712 ! local variable 713 INTEGER(i4) :: il_imin1 714 INTEGER(i4) :: il_jmin1 715 INTEGER(i4) :: il_imax1 716 INTEGER(i4) :: il_jmax1 717 718 REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value 719 720 ! loop indices 721 !---------------------------------------------------------------- 722 723 IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 724 725 il_imin1=1+id_offset(1,1) 726 il_jmin1=1+id_offset(2,1) 727 728 il_imax1=td_var%t_dim(1)%i_len-id_offset(2,1) 729 il_jmax1=td_var%t_dim(2)%i_len-id_offset(2,2) 730 731 ! compute new dimension 732 td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 733 td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 734 ALLOCATE( dl_value( td_var%t_dim(1)%i_len, & 735 & td_var%t_dim(2)%i_len, & 736 & td_var%t_dim(3)%i_len, & 737 & td_var%t_dim(4)%i_len) ) 738 739 dl_value(:,:,:,:)=td_var%d_value( il_imin1:il_imax1, & 740 & il_jmin1:il_jmax1, & 741 & :,: ) 742 743 ! delete table of value 744 CALL var_del_value(td_var) 745 746 ! add variable value 747 CALL var_add_value(td_var,dl_value(:,:,:,:)) 748 749 DEALLOCATE(dl_value) 750 751 ENDIF 752 753 END SUBROUTINE interp__del_offset 754 !> @endcode 755 !------------------------------------------------------------------- 756 !> @brief 757 !> This subroutine interpolate detected point. 711 !> This subroutine interpolate variable value. 758 712 !> 759 713 !> @details … … 762 716 !> 763 717 !> @author J.Paul 764 !> - Nov , 2013- Initial Version765 ! 766 !> @param[in ] td_var : variable to be interpolated767 !> @param[in] id_rho : tableof refinement factor768 ! -------------------------------------------------------------------769 ! > @code718 !> - November, 2013- Initial Version 719 !> 720 !> @param[inout] td_var variable structure 721 !> @param[in] id_rho array of refinement factor 722 !> @param[in] id_offset 2D array of offset between fine and coarse grid 723 !------------------------------------------------------------------- 770 724 SUBROUTINE interp__fill_value_wrapper( td_var, & 771 725 & id_rho, & … … 774 728 ! Argument 775 729 TYPE(TVAR) , INTENT(INOUT) :: td_var 776 INTEGER(I4), DIMENSION(:) ,INTENT(IN ), OPTIONAL :: id_rho730 INTEGER(I4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho 777 731 INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset 778 732 779 733 ! local variable 780 CHARACTER(LEN=lc) :: cl_method781 INTEGER(i4) , DIMENSION(:), ALLOCATABLE 782 INTEGER(i4) , DIMENSION(2,2) :: il_offset734 CHARACTER(LEN=lc) :: cl_method 735 INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho 736 INTEGER(i4) , DIMENSION(2,2) :: il_offset 783 737 784 738 ! loop indices 785 739 !---------------------------------------------------------------- 786 740 787 SELECT CASE(TRIM(td_var%c_interp(1)))788 CASE('cubic','linear','nearest')789 cl_method=TRIM(td_var%c_interp(1))790 CASE DEFAULT791 CALL logger_warn("INTERP FILL: interpolation method unknown."//&792 & " use linear interpolation")793 cl_method='linear'794 795 ! update variable structure value796 td_var%c_interp(1)='linear'797 END SELECT798 799 741 ALLOCATE( il_rho(ip_maxdim) ) 800 742 il_rho(:)=1 801 IF( PRESENT(id_rho) ) il_rho( :)=id_rho(:)743 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 802 744 IF( ANY(il_rho(:) < 0) )THEN 803 745 CALL logger_error("INTERP FILL VALUE: invalid "//& … … 808 750 IF( PRESENT(id_offset) )THEN 809 751 IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN 810 CALL logger_error("INTERP FILL VALUE: invalid tableof offset")752 CALL logger_error("INTERP FILL VALUE: invalid array of offset") 811 753 ELSE 812 754 il_offset(:,:)=id_offset(:,:) … … 817 759 & (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. & 818 760 & (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN 761 762 SELECT CASE(TRIM(td_var%c_interp(1))) 763 CASE('cubic','linear','nearest') 764 cl_method=TRIM(td_var%c_interp(1)) 765 CASE DEFAULT 766 CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//& 767 & " use linear interpolation") 768 cl_method='linear' 769 ! update variable structure value 770 td_var%c_interp(1)='linear' 771 END SELECT 819 772 820 773 CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//& … … 826 779 827 780 CALL interp__fill_value( td_var, cl_method, & 828 & il_rho(:), & 829 & il_offset(:,:) ) 781 & il_rho(:), il_offset(:,:) ) 830 782 831 783 SELECT CASE(TRIM(td_var%c_interp(2))) … … 864 816 END SELECT 865 817 818 ELSE 819 td_var%c_interp(:)='' 866 820 ENDIF 867 821 822 DEALLOCATE(il_rho) 823 868 824 END SUBROUTINE interp__fill_value_wrapper 869 !> @endcode870 825 !------------------------------------------------------------------- 871 826 !> @brief 872 827 !> This subroutine interpolate value over mixed grid. 873 828 !> 874 !> @details875 !>876 !>877 829 !> @author J.Paul 878 !> - Nov, 2013- Initial Version 879 ! 880 !> @param[inout] td_var : variable 881 !> @param[in] cd_method : interpolation method 882 !> @param[in] id_rho : table of refinment factor 883 !------------------------------------------------------------------- 884 !> @code 830 !> - November, 2013- Initial Version 831 !> @date September, 2014 832 !> - use interpolation method modules 833 !> 834 !> @param[inout] td_var variable structure 835 !> @param[in] cd_method interpolation method 836 !> @param[in] id_rho array of refinment factor 837 !> @param[in] id_offset 2D array of offset between fine and coarse grid 838 !------------------------------------------------------------------- 885 839 SUBROUTINE interp__fill_value( td_var, cd_method, & 886 & id_rho, & 887 & id_offset ) 840 & id_rho, id_offset ) 888 841 IMPLICIT NONE 889 842 ! Argument … … 891 844 CHARACTER(LEN=*) , INTENT(IN ) :: cd_method 892 845 INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho 893 INTEGER(I4) , DIMENSION( :,:), INTENT(IN ) :: id_offset846 INTEGER(I4) , DIMENSION(2,2), INTENT(IN ) :: id_offset 894 847 895 848 ! local variable … … 908 861 909 862 TYPE(TATT) :: tl_att 910 863 911 864 ! loop indices 912 INTEGER(i4) :: ji913 INTEGER(i4) :: jj914 INTEGER(i4) :: jk915 INTEGER(i4) :: jl916 865 !---------------------------------------------------------------- 917 866 … … 935 884 tl_att=att_init('interpolation',cl_interp) 936 885 CALL var_move_att(tl_mix, tl_att) 886 887 ! clean 888 CALL att_clean(tl_att) 937 889 938 890 ! special case for even refinement on ARAKAWA-C grid … … 972 924 973 925 !3- interpolate 974 DO jl=1,tl_mix%t_dim(4)%i_len 975 IF( il_rho(jp_K) /= 1 )THEN 976 ! CALL interp__3D(tl_mix%d_value(:,:,:,jl), tl_mix%d_fill, & 977 ! & il_detect(:,:,:), cd_method, & 978 ! & il_rhoi, il_rhoj, il_rhok, & 979 ! & ll_even(:), ll_discont ) 980 CALL logger_error("INTERP FILL: can not interpolate "//& 981 & "vertically for now ") 982 ENDIF 983 984 IF( ANY(il_detect(:,:,:)==1) )THEN 985 ! I-J plan 986 DO jk=1,tl_mix%t_dim(3)%i_len 987 CALL interp__2D(tl_mix%d_value(:,:,jk,jl), tl_mix%d_fill, & 988 & il_detect(:,:,jk), cd_method, & 989 & il_rho(jp_I), il_rho(jp_J), & 990 & ll_even(1:2), ll_discont) 991 ENDDO 992 IF( ALL(il_detect(:,:,:)==0) ) CYCLE 993 IF( il_rho(jp_K) /= 1 )THEN 994 ! I-K plan 995 DO jj=1,tl_mix%t_dim(2)%i_len 996 CALL interp__2D(tl_mix%d_value(:,jj,:,jl), tl_mix%d_fill, & 997 & il_detect(:,jj,:), cd_method, & 998 & il_rho(jp_J), il_rho(jp_K), & 999 & ll_even(1:3:2), ll_discont ) 1000 ENDDO 1001 IF( ALL(il_detect(:,:,:)==0) ) CYCLE 1002 ! J-K plan 1003 DO ji=1,tl_mix%t_dim(1)%i_len 1004 CALL interp__2D(tl_mix%d_value(ji,:,:,jl), tl_mix%d_fill, & 1005 & il_detect(ji,:,:), cd_method, & 1006 & il_rho(jp_J), il_rho(jp_K), & 1007 & ll_even(2:3), ll_discont ) 1008 ENDDO 1009 1010 ENDIF 1011 IF( ANY(il_detect(:,:,:)==1) )THEN 1012 ! I direction 1013 DO jk=1,tl_mix%t_dim(3)%i_len 1014 DO jj=1,tl_mix%t_dim(2)%i_len 1015 CALL interp__1D( tl_mix%d_value(:,jj,jk,jl), & 1016 & tl_mix%d_fill, & 1017 & il_detect(:,jj,jk), cd_method, & 1018 & il_rho(jp_I), & 1019 & ll_even(1), ll_discont ) 1020 ENDDO 1021 ENDDO 1022 IF( ALL(il_detect(:,:,:)==0) ) CYCLE 1023 ! J direction 1024 DO jk=1,tl_mix%t_dim(3)%i_len 1025 DO ji=1,tl_mix%t_dim(1)%i_len 1026 CALL interp__1D( tl_mix%d_value(ji,:,jk,jl), & 1027 & tl_mix%d_fill, & 1028 & il_detect(ji,:,jk), cd_method, & 1029 & il_rho(jp_J), & 1030 & ll_even(2), ll_discont ) 1031 ENDDO 1032 ENDDO 1033 IF( il_rho(jp_K) /= 1 )THEN 1034 IF( ALL(il_detect(:,:,:)==0) ) CYCLE 1035 ! K direction 1036 DO jj=1,tl_mix%t_dim(2)%i_len 1037 DO ji=1,tl_mix%t_dim(1)%i_len 1038 CALL interp__1D( tl_mix%d_value(ji,jj,:,jl), & 1039 & tl_mix%d_fill, & 1040 & il_detect(ji,jj,:), cd_method, & 1041 & il_rho(jp_K), & 1042 & ll_even(3), ll_discont ) 1043 ENDDO 1044 ENDDO 1045 ENDIF 1046 ENDIF 1047 ENDIF 1048 ENDDO 926 CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//& 927 & " discont "//TRIM(fct_str(ll_discont)) ) 928 SELECT CASE(TRIM(cd_method)) 929 CASE('cubic') 930 CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, & 931 & il_detect(:,:,:), & 932 & il_rho(:), ll_even(:), ll_discont ) 933 CASE('nearest') 934 CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), & 935 & il_detect(:,:,:), & 936 & il_rho(:) ) 937 CASE DEFAULT ! linear 938 CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, & 939 & il_detect(:,:,:), & 940 & il_rho(:), ll_even(:), ll_discont ) 941 END SELECT 1049 942 1050 943 IF( ANY(il_detect(:,:,:)==1) )THEN … … 1056 949 !4- save useful domain (remove offset) 1057 950 CALL interp_clean_mixed_grid( tl_mix, td_var, & 1058 & id_rho(:), & 1059 & id_offset(:,:) ) 951 & id_rho(:), id_offset(:,:) ) 1060 952 1061 953 ! clean variable structure 954 DEALLOCATE(il_rho) 1062 955 CALL var_clean(tl_mix) 1063 956 1064 957 END SUBROUTINE interp__fill_value 1065 !> @endcode1066 ! !-------------------------------------------------------------------1067 ! !> @brief1068 ! !> This subroutine1069 ! !>1070 ! !> @details1071 ! !>1072 ! !> @author J.Paul1073 ! !> - Nov, 2013- Initial Version1074 ! !1075 ! !> @param[in]1076 ! !> @param[out]1077 ! !-------------------------------------------------------------------1078 ! !> @code1079 ! SUBROUTINE interp__3D( dd_value, dd_fill, &1080 ! & id_detect, cd_method, &1081 ! & id_rhoi, id_rhoj, id_rhok, &1082 ! & ld_even, ld_discont )1083 ! IMPLICIT NONE1084 ! ! Argument1085 ! REAL(dp) , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value1086 ! REAL(dp) , INTENT(IN ) :: dd_fill1087 ! INTEGER(I4) , DIMENSION(:,:,:), INTENT(INOUT) :: id_detect1088 ! CHARACTER(LEN=*) , INTENT(IN ) :: cd_method1089 ! INTEGER(I4) , INTENT(IN ) :: id_rhoi1090 ! INTEGER(I4) , INTENT(IN ) :: id_rhoj1091 ! INTEGER(I4) , INTENT(IN ) :: id_rhok1092 ! LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even1093 ! LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont1094 !1095 ! ! local variable1096 ! INTEGER(i4), DIMENSION(3) :: il_shape1097 ! INTEGER(i4), DIMENSION(3) :: il_dim1098 !1099 ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coarse1100 ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx1101 ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy1102 ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_d2fdxy1103 !1104 ! LOGICAL :: ll_discont1105 !1106 ! ! loop indices1107 !! INTEGER(i4) :: ji1108 !! INTEGER(i4) :: jj1109 !! INTEGER(i4) :: ii1110 !! INTEGER(i4) :: ij1111 ! !----------------------------------------------------------------1112 ! ll_discont=.FALSE.1113 ! IF( PRESENT(ld_discont) ) ll_discont=ld_discont1114 !1115 ! il_shape(:)=SHAPE(dd_value)1116 ! il_dim(:)=(il_shape(:)-1)/21117 !1118 ! ALLOCATE( dl_coarse(il_dim(1),il_dim(2),il_dim(3)) )1119 ! ! value on coarse grid1120 ! dl_coarse(:,:,:)=dd_value( 1:il_shape(1):id_rhoi, &1121 ! & 1:il_shape(2):id_rhoj, &1122 ! & 1:il_shape(3):id_rhok )1123 ! SELECT CASE(TRIM(cd_method))1124 !1125 ! CASE('cubic')1126 !1127 ! ALLOCATE( dl_dfdx( il_dim(1),il_dim(2),il_dim(3)), &1128 ! & dl_dfdy( il_dim(1),il_dim(2),il_dim(3)), &1129 ! & dl_d2fdxy(il_dim(1),il_dim(2),il_dim(3)) )1130 !1131 !! ! compute derivative on coarse grid1132 !! dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'I')1133 !! dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'J')1134 !!1135 !! ! compute cross derivative on coarse grid1136 !! dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:,:), dd_fill, 'J')1137 !!1138 !! DO jj=1,il_shape(2)-1,id_rhoj1139 !! ij=((jj-1)/id_rhoj)+11140 !! DO ji=1,il_shape(1)-1,id_rhoi1141 !! ii=((ji-1)/id_rhoi)+11142 !!1143 !! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE1144 !! ! compute bicubic coefficient1145 !! dl_coef(:)=interp__2D_cubic_coef(dl_coarse(ii:ii+1,ij:ij+1),&1146 !! & dl_dfdx( ii:ii+1,ij:ij+1),&1147 !! & dl_dfdy( ii:ii+1,ij:ij+1),&1148 !! & dl_d2fdxy(ii:ii+1,ij:ij+1) )1149 !!1150 !! ! compute value on detetected point1151 !! CALL interp__2D_cubic_fill(dl_coef(:), &1152 !! & dd_value( ji:ji+id_rhoi, &1153 !! & jj:jj+id_rhoj ), &1154 !! & id_detect(ji:ji+id_rhoi, &1155 !! & jj:jj+id_rhoj ) )1156 !!1157 !! ENDDO1158 !! ENDDO1159 !1160 ! DEALLOCATE( dl_dfdx, &1161 ! & dl_dfdy, &1162 ! & dl_d2fdxy )1163 !1164 ! CASE('nearest')1165 !1166 ! CASE DEFAULT ! linear1167 !1168 !! DO jj=1,il_shape(2)-1,id_rhoj1169 !! ij=((jj-1)/id_rhoj)+11170 !! DO ji=1,il_shape(1)-1,id_rhoi1171 !! ii=((ji-1)/id_rhoi)+11172 !!1173 !! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE1174 !! ! compute bilinear coefficient1175 !! dl_coef(:)=interp__2D_linear_coef(dl_coarse(ii:ii+1,ij:ij+1))1176 !!1177 !! ! compute value on detetected point1178 !! CALL interp__2D_linear_fill(dl_coef(:), &1179 !! & dd_value( ji:ji+id_rhoi, &1180 !! & jj:jj+id_rhoj ), &1181 !! & id_detect(ji:ji+id_rhoi, &1182 !! & jj:jj+id_rhoj ) )1183 !!1184 !! ENDDO1185 !! ENDDO1186 !1187 ! END SELECT1188 !1189 ! DEALLOCATE( dl_coarse )1190 !1191 ! END SUBROUTINE interp__3D1192 ! !> @endcode1193 !-------------------------------------------------------------------1194 !> @brief1195 !> This subroutine1196 !>1197 !> @details1198 !>1199 !> @author J.Paul1200 !> - Nov, 2013- Initial Version1201 !1202 !> @param[in]1203 !> @param[out]1204 !-------------------------------------------------------------------1205 !> @code1206 SUBROUTINE interp__2D( dd_value, dd_fill, &1207 & id_detect, cd_method, &1208 & id_rhoi, id_rhoj, &1209 & ld_even, ld_discont )1210 1211 IMPLICIT NONE1212 ! Argument1213 REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value1214 REAL(dp) , INTENT(IN ) :: dd_fill1215 INTEGER(I4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect1216 CHARACTER(LEN=*) , INTENT(IN ) :: cd_method1217 INTEGER(I4) , INTENT(IN ) :: id_rhoi1218 INTEGER(I4) , INTENT(IN ) :: id_rhoj1219 LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even1220 LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont1221 1222 ! local variable1223 INTEGER(I4) :: il_xextra1224 INTEGER(I4) :: il_yextra1225 INTEGER(i4), DIMENSION(2) :: il_shape1226 INTEGER(i4), DIMENSION(2) :: il_dim1227 1228 REAL(dp) :: dl_min1229 REAL(dp) :: dl_max1230 REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_coef1231 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_coarse1232 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_tmp1233 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx1234 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy1235 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy1236 1237 LOGICAL :: ll_discont1238 ! loop indices1239 INTEGER(i4) :: ji1240 INTEGER(i4) :: jj1241 INTEGER(i4) :: ii1242 INTEGER(i4) :: ij1243 1244 !----------------------------------------------------------------1245 ll_discont=.FALSE.1246 IF( PRESENT(ld_discont) ) ll_discont=ld_discont1247 1248 CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&1249 & " discont "//TRIM(fct_str(ll_discont)) )1250 1251 il_shape(:)=SHAPE(dd_value)1252 1253 ! compute coarse grid dimension1254 il_xextra=id_rhoi-11255 il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi1256 1257 il_yextra=id_rhoj-11258 il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj1259 1260 ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )1261 1262 ! value on coarse grid1263 dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &1264 & 1:il_shape(2):id_rhoj )1265 1266 SELECT CASE(TRIM(cd_method))1267 1268 CASE('cubic')1269 1270 ALLOCATE( dl_dfdx( il_dim(1),il_dim(2)), &1271 & dl_dfdy( il_dim(1),il_dim(2)), &1272 & dl_d2fdxy(il_dim(1),il_dim(2)) )1273 1274 ! compute derivative on coarse grid1275 dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ll_discont)1276 dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ll_discont)1277 1278 ! compute cross derivative on coarse grid1279 dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ll_discont)1280 1281 ALLOCATE( dl_tmp(2,2) )1282 ALLOCATE( dl_coef(16) )1283 1284 DO jj=1,il_shape(2)-1,id_rhoj1285 ij=((jj-1)/id_rhoj)+11286 DO ji=1,il_shape(1)-1,id_rhoi1287 ii=((ji-1)/id_rhoi)+11288 1289 IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. &1290 & ANY( dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. &1291 & ANY( dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. &1292 & ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE1293 1294 dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)1295 IF( ll_discont )THEN1296 1297 dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )1298 dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )1299 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN1300 WHERE( dl_tmp(:,:) < 0_dp )1301 dl_tmp(:,:) = dl_tmp(:,:)+360._dp1302 END WHERE1303 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1304 WHERE( dl_tmp(:,:) > 180_dp )1305 dl_tmp(:,:) = dl_tmp(:,:)-180._dp1306 END WHERE1307 ENDIF1308 ENDIF1309 1310 ! compute bicubic coefficient1311 dl_coef(:)=interp__2D_cubic_coef(dl_tmp(:,:),&1312 & dl_dfdx( ii:ii+1,ij:ij+1),&1313 & dl_dfdy( ii:ii+1,ij:ij+1),&1314 & dl_d2fdxy( ii:ii+1,ij:ij+1),&1315 & dd_fill )1316 1317 ! compute value on detetected point1318 CALL interp__2D_cubic_fill(dd_value( ji:ji+id_rhoi, &1319 & jj:jj+id_rhoj ), &1320 & id_detect(ji:ji+id_rhoi, &1321 & jj:jj+id_rhoj ), &1322 & dl_coef(:), dd_fill, &1323 & ld_even(:), id_rhoi, id_rhoj )1324 1325 IF( ll_discont )THEN1326 WHERE( dd_value( ji:ji+id_rhoi, &1327 & jj:jj+id_rhoj ) >= 180._dp .AND. &1328 & dd_value( ji:ji+id_rhoi, &1329 & jj:jj+id_rhoj ) /= dd_fill )1330 dd_value( ji:ji+id_rhoi, &1331 & jj:jj+id_rhoj ) = &1332 & dd_value( ji:ji+id_rhoi, &1333 & jj:jj+id_rhoj ) - 360._dp1334 END WHERE1335 ENDIF1336 1337 ENDDO1338 ENDDO1339 1340 DEALLOCATE(dl_coef)1341 DEALLOCATE(dl_tmp )1342 1343 DEALLOCATE(dl_dfdx, &1344 & dl_dfdy, &1345 & dl_d2fdxy )1346 1347 CASE('nearest')1348 1349 DO jj=1,il_shape(2)-1,id_rhoj1350 ij=((jj-1)/id_rhoj)+11351 DO ji=1,il_shape(1)-1,id_rhoi1352 ii=((ji-1)/id_rhoi)+11353 1354 ! compute value on detetected point1355 CALL interp__2D_nearest_fill(dd_value( ji:ji+id_rhoi, &1356 & jj:jj+id_rhoj ), &1357 & id_detect(ji:ji+id_rhoi, &1358 & jj:jj+id_rhoj ) )1359 1360 ENDDO1361 ENDDO1362 1363 CASE DEFAULT ! linear1364 1365 ALLOCATE( dl_coef(4) )1366 ALLOCATE( dl_tmp(2,2) )1367 1368 DO jj=1,il_shape(2)-1,id_rhoj1369 ij=((jj-1)/id_rhoj)+11370 DO ji=1,il_shape(1)-1,id_rhoi1371 ii=((ji-1)/id_rhoi)+11372 1373 IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE1374 1375 dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)1376 IF( ll_discont )THEN1377 1378 dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )1379 dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )1380 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN1381 WHERE( dl_tmp(:,:) < 0_dp )1382 dl_tmp(:,:) = dl_tmp(:,:)+360._dp1383 END WHERE1384 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1385 WHERE( dl_tmp(:,:) > 180_dp )1386 dl_tmp(:,:) = dl_tmp(:,:)-180._dp1387 END WHERE1388 ENDIF1389 ENDIF1390 1391 ! compute bilinear coefficient1392 dl_coef(:)=interp__2D_linear_coef(dl_tmp(:,:), dd_fill)1393 1394 ! compute value on detetected point1395 CALL interp__2D_linear_fill(dd_value( ji:ji+id_rhoi, &1396 & jj:jj+id_rhoj ), &1397 & id_detect(ji:ji+id_rhoi, &1398 & jj:jj+id_rhoj ), &1399 & dl_coef(:), dd_fill, &1400 & ld_even(:), id_rhoi, id_rhoj )1401 1402 IF( ll_discont )THEN1403 WHERE( dd_value( ji:ji+id_rhoi, &1404 & jj:jj+id_rhoj ) >= 180._dp .AND. &1405 & dd_value( ji:ji+id_rhoi, &1406 & jj:jj+id_rhoj ) /= dd_fill )1407 dd_value( ji:ji+id_rhoi, &1408 & jj:jj+id_rhoj ) = &1409 & dd_value( ji:ji+id_rhoi, &1410 & jj:jj+id_rhoj ) - 360._dp1411 END WHERE1412 ENDIF1413 1414 ENDDO1415 ENDDO1416 1417 DEALLOCATE(dl_coef)1418 DEALLOCATE(dl_tmp )1419 1420 END SELECT1421 1422 DEALLOCATE( dl_coarse )1423 1424 END SUBROUTINE interp__2D1425 !> @endcode1426 !-------------------------------------------------------------------1427 !> @brief1428 !> This subroutine1429 !>1430 !> @details1431 !>1432 !> @author J.Paul1433 !> - Nov, 2013- Initial Version1434 !1435 !> @param[in]1436 !> @param[out]1437 !-------------------------------------------------------------------1438 !> @code1439 SUBROUTINE interp__1D( dd_value, dd_fill, &1440 & id_detect, cd_method, &1441 & id_rhoi, &1442 & ld_even, ld_discont )1443 IMPLICIT NONE1444 ! Argument1445 REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value1446 REAL(dp) , INTENT(IN ) :: dd_fill1447 INTEGER(I4) , DIMENSION(:), INTENT(INOUT) :: id_detect1448 CHARACTER(LEN=*) , INTENT(IN ) :: cd_method1449 INTEGER(I4) , INTENT(IN ) :: id_rhoi1450 LOGICAL , INTENT(IN ) :: ld_even1451 LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont1452 1453 ! local variable1454 INTEGER(i4), DIMENSION(1) :: il_shape1455 INTEGER(i4), DIMENSION(1) :: il_dim1456 INTEGER(I4) :: il_xextra1457 1458 REAL(dp) :: dl_min1459 REAL(dp) :: dl_max1460 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coarse1461 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_tmp1462 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_dfdx1463 1464 REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coef1465 1466 LOGICAL :: ll_discont1467 1468 ! loop indices1469 INTEGER(i4) :: ji1470 INTEGER(i4) :: ii1471 !----------------------------------------------------------------1472 ll_discont=.FALSE.1473 IF( PRESENT(ld_discont) ) ll_discont=ld_discont1474 1475 CALL logger_debug("INTERP 1D: interpolation method "//TRIM(cd_method)//&1476 & " discont "//TRIM(fct_str(ll_discont)) )1477 1478 il_shape(:)=SHAPE(dd_value)1479 1480 il_xextra=id_rhoi-11481 il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi1482 1483 ALLOCATE( dl_coarse(il_dim(1)) )1484 ! value on coarse grid1485 dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )1486 1487 SELECT CASE(TRIM(cd_method))1488 1489 CASE('cubic')1490 1491 ALLOCATE( dl_dfdx( il_dim(1)) )1492 1493 ! compute derivative on coarse grid1494 dl_dfdx(:)=extrap_deriv_1D(dl_coarse(:), dd_fill, ll_discont)1495 1496 ALLOCATE( dl_coef(4))1497 ALLOCATE( dl_tmp(2) )1498 1499 DO ji=1,il_shape(1)-1,id_rhoi1500 ii=((ji-1)/id_rhoi)+11501 1502 IF( ANY( dl_tmp(:)==dd_fill ) .OR. &1503 & ANY(dl_dfdx(:)==dd_fill ) ) CYCLE1504 1505 dl_tmp(:)=dl_coarse(ii:ii+1)1506 IF( ll_discont )THEN1507 1508 dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )1509 dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )1510 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN1511 WHERE( dl_tmp(:) < 0_dp )1512 dl_tmp(:) = dl_tmp(:)+360._dp1513 END WHERE1514 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1515 WHERE( dl_tmp(:) > 180_dp )1516 dl_tmp(:) = dl_tmp(:)-180._dp1517 END WHERE1518 ENDIF1519 ENDIF1520 1521 ! compute bicubic coefficient1522 dl_coef(:)=interp__1D_cubic_coef(dl_tmp(:),&1523 & dl_dfdx( ii:ii+1), &1524 & dd_fill )1525 1526 ! compute value on detetected point1527 CALL interp__1D_cubic_fill(dd_value( ji:ji+id_rhoi),&1528 & id_detect(ji:ji+id_rhoi),&1529 & dl_coef(:), dd_fill, &1530 & ld_even, id_rhoi )1531 1532 IF( ll_discont )THEN1533 WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &1534 & dd_value( ji:ji+id_rhoi ) /= dd_fill )1535 dd_value( ji:ji+id_rhoi ) = &1536 & dd_value( ji:ji+id_rhoi ) - 360._dp1537 END WHERE1538 ENDIF1539 1540 ENDDO1541 1542 DEALLOCATE(dl_coef)1543 DEALLOCATE(dl_tmp )1544 1545 CASE('nearest')1546 1547 DO ji=1,il_shape(1)-1,id_rhoi1548 ii=((ji-1)/id_rhoi)+11549 1550 ! compute value on detetected point1551 CALL interp__1D_nearest_fill(dd_value( ji:ji+id_rhoi), &1552 & id_detect(ji:ji+id_rhoi) )1553 1554 ENDDO1555 1556 CASE DEFAULT ! linear1557 1558 ALLOCATE(dl_coef(2))1559 ALLOCATE( dl_tmp(2) )1560 1561 DO ji=1,il_shape(1)-1,id_rhoi1562 ii=((ji-1)/id_rhoi)+11563 1564 IF( ANY(dl_coarse(ii:ii+1)==dd_fill ) ) CYCLE1565 1566 dl_tmp(:)=dl_coarse(ii:ii+1)1567 IF( ll_discont )THEN1568 1569 dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )1570 dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )1571 IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN1572 WHERE( dl_tmp(:) < 0_dp )1573 dl_tmp(:) = dl_tmp(:)+360._dp1574 END WHERE1575 ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN1576 WHERE( dl_tmp(:) > 180_dp )1577 dl_tmp(:) = dl_tmp(:)-180._dp1578 END WHERE1579 ENDIF1580 ENDIF1581 1582 ! compute bilinear coefficient1583 dl_coef(:)=interp__1D_linear_coef(dl_tmp(:), dd_fill)1584 1585 ! compute value on detetected point1586 CALL interp__1D_linear_fill( dd_value( ji:ji+id_rhoi),&1587 & id_detect(ji:ji+id_rhoi),&1588 & dl_coef(:), dd_fill, &1589 & ld_even, id_rhoi )1590 1591 IF( ll_discont )THEN1592 WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &1593 & dd_value( ji:ji+id_rhoi ) /= dd_fill )1594 dd_value( ji:ji+id_rhoi ) = &1595 & dd_value( ji:ji+id_rhoi ) - 360._dp1596 END WHERE1597 ENDIF1598 1599 ENDDO1600 1601 DEALLOCATE(dl_coef)1602 1603 END SELECT1604 1605 DEALLOCATE( dl_coarse )1606 1607 END SUBROUTINE interp__1D1608 !> @endcode1609 !-------------------------------------------------------------------1610 !> @brief1611 !> This subroutine1612 !>1613 !> @details1614 !>1615 !> @author J.Paul1616 !> - Nov, 2013- Initial Version1617 !1618 !> @param[in]1619 !> @param[out]1620 !-------------------------------------------------------------------1621 !> @code1622 FUNCTION interp__2D_cubic_coef( dd_value, &1623 & dd_dfdx, &1624 & dd_dfdy, &1625 & dd_d2fdxy,&1626 & dd_fill )1627 IMPLICIT NONE1628 ! Argument1629 REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value1630 REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx1631 REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy1632 REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy1633 REAL(dp) , INTENT(IN) :: dd_fill1634 1635 ! function1636 REAL(dp), DIMENSION(16) :: interp__2D_cubic_coef1637 1638 ! local variable1639 REAL(dp), DIMENSION(16,16), PARAMETER :: dp_matrix = RESHAPE( &1640 & (/ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&1641 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&1642 -3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&1643 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&1644 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&1645 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,&1646 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 ,&1647 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 ,&1648 -3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,&1649 0 , 0 , 0 , 0 ,-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 ,&1650 9 ,-9 ,-9 , 9 , 6 , 3 ,-6 ,-3 , 6 ,-6 , 3 ,-3 , 4 , 2 , 2 , 1 ,&1651 -6 , 6 , 6 ,-6 ,-3 ,-3 , 3 , 3 ,-4 , 4 ,-2 , 2 ,-2 ,-2 ,-1 ,-1 ,&1652 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,&1653 0 , 0 , 0 , 0 , 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 ,&1654 -6 , 6 , 6 ,-6 ,-4 ,-2 , 4 , 2 ,-3 , 3 ,-3 , 3 ,-2 ,-1 ,-2 ,-1 ,&1655 4 ,-4 ,-4 , 4 , 2 , 2 ,-2 ,-2 , 2 ,-2 , 2 ,-2 , 1 , 1 , 1 , 1 /), &1656 & (/ 16, 16 /) )1657 1658 REAL(dp), DIMENSION(16) :: dl_vect1659 1660 !----------------------------------------------------------------1661 ! init1662 interp__2D_cubic_coef(:)=dd_fill1663 1664 IF( ANY(SHAPE( dd_value(:,:))/= 2) .OR. &1665 & ANY(SHAPE( dd_dfdx(:,:))/= 2) .OR. &1666 & ANY(SHAPE( dd_dfdy(:,:))/= 2) .OR. &1667 & ANY(SHAPE(dd_d2fdxy(:,:))/= 2) )THEN1668 1669 CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&1670 & "input tables. shape should be (/2,2/)")1671 1672 ELSEIF( ANY( dd_value(:,:) == dd_fill) .OR. &1673 & ANY( dd_dfdx(:,:) == dd_fill) .OR. &1674 & ANY( dd_dfdy(:,:) == dd_fill) .OR. &1675 & ANY(dd_d2fdxy(:,:) == dd_fill) )THEN1676 1677 CALL logger_warn("INTERP CUBIC COEF: fill value detected. "//&1678 & "can not compute coefficient ")1679 1680 ELSE1681 1682 dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )1683 dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. )1684 dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. )1685 dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. )1686 1687 interp__2D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))1688 1689 ENDIF1690 END FUNCTION interp__2D_cubic_coef1691 !> @endcode1692 !-------------------------------------------------------------------1693 !> @brief1694 !> This subroutine1695 !>1696 !> @details1697 !>1698 !> @author J.Paul1699 !> - Nov, 2013- Initial Version1700 !1701 !> @param[in]1702 !> @param[out]1703 !-------------------------------------------------------------------1704 !> @code1705 SUBROUTINE interp__2D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &1706 & ld_even, id_rhoi, id_rhoj )1707 IMPLICIT NONE1708 ! Argument1709 REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value1710 INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect1711 REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef1712 REAL(dp) , INTENT(IN ) :: dd_fill1713 LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even1714 INTEGER(I4) , INTENT(IN ) :: id_rhoi1715 INTEGER(I4) , INTENT(IN ) :: id_rhoj1716 1717 ! local variable1718 INTEGER(i4), DIMENSION(2) :: il_shape1719 1720 REAL(dp) , DIMENSION(16) :: dl_vect1721 REAL(dp) :: dl_dx1722 REAL(dp) :: dl_dy1723 REAL(dp) :: dl_x1724 REAL(dp) :: dl_x21725 REAL(dp) :: dl_x31726 REAL(dp) :: dl_y1727 REAL(dp) :: dl_y21728 REAL(dp) :: dl_y31729 1730 ! loop indices1731 INTEGER(i4) :: ji1732 INTEGER(i4) :: jj1733 !----------------------------------------------------------------1734 1735 IF( SIZE(dd_coef(:)) /= 16 )THEN1736 CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&1737 & "coef table. shape should be (/16/)")1738 ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN1739 CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//&1740 & "can not compute interpolation.")1741 ELSE1742 1743 il_shape(:)=SHAPE(dd_value(:,:))1744 1745 dl_dx=1./REAL(id_rhoi)1746 dl_dy=1./REAL(id_rhoj)1747 1748 DO jj=1,il_shape(2)1749 1750 IF( ld_even(2) )THEN1751 dl_y=(jj-1)*dl_dy - dl_dy*0.51752 ELSE ! odd refinement1753 dl_y=(jj-1)*dl_dy1754 ENDIF1755 dl_y2=dl_y*dl_y1756 dl_y3=dl_y2*dl_y1757 1758 DO ji=1,il_shape(1)1759 1760 IF(id_detect(ji,jj)==1)THEN1761 1762 IF( ld_even(1) )THEN1763 dl_x=(ji-1)*dl_dx - dl_dx*0.51764 ELSE ! odd refinement1765 dl_x=(ji-1)*dl_dx1766 ENDIF1767 dl_x2=dl_x*dl_x1768 dl_x3=dl_x2*dl_x1769 1770 dl_vect(:)=(/1._dp, dl_x , dl_x2 , dl_x3 , &1771 & dl_y , dl_x*dl_y , dl_x2*dl_y , dl_x3*dl_y , &1772 & dl_y2, dl_x*dl_y2, dl_x2*dl_y2, dl_x3*dl_y2, &1773 & dl_y3, dl_x*dl_y3, dl_x2*dl_y3, dl_x3*dl_y3 /)1774 1775 dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))1776 id_detect(ji,jj)=01777 1778 ENDIF1779 1780 ENDDO1781 ENDDO1782 1783 ENDIF1784 1785 END SUBROUTINE interp__2D_cubic_fill1786 !> @endcode1787 !-------------------------------------------------------------------1788 !> @brief1789 !> This subroutine1790 !>1791 !> @details1792 !>1793 !> @author J.Paul1794 !> - Nov, 2013- Initial Version1795 !1796 !> @param[in]1797 !> @param[out]1798 !-------------------------------------------------------------------1799 !> @code1800 FUNCTION interp__2D_linear_coef( dd_value, dd_fill )1801 IMPLICIT NONE1802 ! Argument1803 REAL(dp), DIMENSION(:,:) , INTENT(IN) :: dd_value1804 REAL(dp) , INTENT(IN) :: dd_fill1805 1806 ! function1807 REAL(dp), DIMENSION(4) :: interp__2D_linear_coef1808 1809 ! local variable1810 1811 REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &1812 & (/ 1 , 0 , 0 , 0 ,&1813 -1 , 1 , 0 , 0 ,&1814 -1 , 0 , 1 , 0 ,&1815 1 ,-1 ,-1 , 1 /), &1816 & (/ 4, 4 /) )1817 1818 REAL(dp), DIMENSION(4) :: dl_vect1819 1820 !----------------------------------------------------------------1821 1822 IF( ANY(SHAPE(dd_value(:,:))/= 2) )THEN1823 CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&1824 & "input tables. shape should be (/2,2/)")1825 ELSEIF( ANY(dd_value(:,:)==dd_fill) )THEN1826 CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&1827 & "can not compute coefficient.")1828 ELSE1829 1830 dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )1831 1832 interp__2D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))1833 1834 ENDIF1835 END FUNCTION interp__2D_linear_coef1836 !> @endcode1837 !-------------------------------------------------------------------1838 !> @brief1839 !> This subroutine1840 !>1841 !> @details1842 !>1843 !> @author J.Paul1844 !> - Nov, 2013- Initial Version1845 !1846 !> @param[in]1847 !> @param[out]1848 !-------------------------------------------------------------------1849 !> @code1850 SUBROUTINE interp__2D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &1851 & ld_even, id_rhoi, id_rhoj )1852 IMPLICIT NONE1853 ! Argument1854 REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value1855 INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect1856 REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef1857 REAL(dp) , INTENT(IN ) :: dd_fill1858 LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even1859 INTEGER(I4) , INTENT(IN ) :: id_rhoi1860 INTEGER(I4) , INTENT(IN ) :: id_rhoj1861 1862 ! local variable1863 INTEGER(i4), DIMENSION(2) :: il_shape1864 1865 REAL(dp) , DIMENSION(4) :: dl_vect1866 REAL(dp) :: dl_dx1867 REAL(dp) :: dl_dy1868 REAL(dp) :: dl_x1869 REAL(dp) :: dl_y1870 1871 ! loop indices1872 INTEGER(i4) :: ji1873 INTEGER(i4) :: jj1874 !----------------------------------------------------------------1875 1876 IF( SIZE(dd_coef(:)) /= 4 )THEN1877 CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&1878 & "coef table. shape should be (/4/)")1879 ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN1880 CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//&1881 & "can not compute interpolation.")1882 ELSE1883 1884 il_shape(:)=SHAPE(dd_value(:,:))1885 1886 dl_dx=1./REAL(id_rhoi)1887 dl_dy=1./REAL(id_rhoj)1888 1889 DO jj=1,il_shape(2)1890 1891 IF( ld_even(2) )THEN1892 dl_y=(jj-1)*dl_dy - dl_dy*0.51893 ELSE ! odd refinement1894 dl_y=(jj-1)*dl_dy1895 ENDIF1896 1897 DO ji=1,il_shape(1)1898 1899 IF(id_detect(ji,jj)==1)THEN1900 1901 IF( ld_even(1) )THEN1902 dl_x=(ji-1)*dl_dx - dl_dx*0.51903 ELSE ! odd refinement1904 dl_x=(ji-1)*dl_dx1905 ENDIF1906 1907 dl_vect(:)=(/1._dp, dl_x, dl_y, dl_x*dl_y /)1908 1909 dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))1910 id_detect(ji,jj)=01911 1912 ENDIF1913 1914 ENDDO1915 ENDDO1916 1917 ENDIF1918 1919 END SUBROUTINE interp__2D_linear_fill1920 !> @endcode1921 !-------------------------------------------------------------------1922 !> @brief1923 !> This subroutine1924 !>1925 !> @details1926 !>1927 !> @author J.Paul1928 !> - Nov, 2013- Initial Version1929 !1930 !> @param[in]1931 !> @param[out]1932 !-------------------------------------------------------------------1933 !> @code1934 SUBROUTINE interp__2D_nearest_fill( dd_value, id_detect )1935 IMPLICIT NONE1936 ! Argument1937 REAL(dp) , DIMENSION(:,:) , INTENT(INOUT) :: dd_value1938 INTEGER(i4), DIMENSION(:,:) , INTENT(INOUT) :: id_detect1939 1940 ! local variable1941 INTEGER(i4), DIMENSION(2) :: il_shape1942 1943 INTEGER(i4) :: il_i11944 INTEGER(i4) :: il_i21945 INTEGER(i4) :: il_j11946 INTEGER(i4) :: il_j21947 1948 INTEGER(i4) :: il_half11949 INTEGER(i4) :: il_half21950 1951 ! loop indices1952 INTEGER(i4) :: ji1953 INTEGER(i4) :: jj1954 !----------------------------------------------------------------1955 1956 il_shape(:)=SHAPE(dd_value(:,:))1957 1958 il_i1=11959 il_i2=il_shape(1)1960 1961 il_j1=11962 il_j2=il_shape(2)1963 1964 il_half1=CEILING(il_shape(1)*0.5)1965 il_half2=CEILING(il_shape(2)*0.5)1966 1967 DO jj=1,il_half21968 1969 DO ji=1,il_half11970 1971 ! lower left point1972 IF(id_detect(ji,jj)==1)THEN1973 1974 dd_value( ji,jj)=dd_value(il_i1,il_j1)1975 id_detect(ji,jj)=01976 1977 ENDIF1978 1979 ! lower right point1980 IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN1981 1982 dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1)1983 id_detect(il_shape(1)-ji+1,jj)=01984 1985 ENDIF1986 1987 ! upper left point1988 IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN1989 1990 dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2)1991 id_detect(ji,il_shape(2)-jj+1)=01992 1993 ENDIF1994 1995 ! upper right point1996 IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN1997 1998 dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2)1999 id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=02000 2001 ENDIF2002 2003 ENDDO2004 2005 ENDDO2006 2007 END SUBROUTINE interp__2D_nearest_fill2008 !> @endcode2009 !-------------------------------------------------------------------2010 !> @brief2011 !> This subroutine2012 !>2013 !> @details2014 !>2015 !> @author J.Paul2016 !> - Nov, 2013- Initial Version2017 !2018 !> @param[in]2019 !> @param[out]2020 !-------------------------------------------------------------------2021 !> @code2022 FUNCTION interp__1D_cubic_coef( dd_value, &2023 & dd_dfdx, &2024 & dd_fill )2025 IMPLICIT NONE2026 ! Argument2027 REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value2028 REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_dfdx2029 REAL(dp) , INTENT(IN) :: dd_fill2030 2031 ! function2032 REAL(dp), DIMENSION(4) :: interp__1D_cubic_coef2033 2034 ! local variable2035 2036 REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &2037 & (/ 1 , 0 , 0 , 0 ,&2038 -1 , 1 , 0 , 0 ,&2039 -3 , 3 ,-2 ,-1 ,&2040 2 ,-2 , 1 , 1 /), &2041 & (/ 4, 4 /) )2042 2043 REAL(dp), DIMENSION(4) :: dl_vect2044 2045 !----------------------------------------------------------------2046 IF( SIZE(dd_value(:))/= 2 .OR. &2047 & SIZE(dd_dfdx(:) )/= 2 )THEN2048 2049 CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&2050 & "input tables. shape should be (/2,2/)")2051 2052 ELSEIF( ANY(dd_value(:)==dd_fill) .OR. &2053 & ANY( dd_dfdx(:)==dd_fill) )THEN2054 CALL logger_error("INTERP CUBIC COEF: fill value detected. "//&2055 & "can not compute coefficient.")2056 ELSE2057 2058 dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )2059 dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. )2060 2061 interp__1D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))2062 2063 ENDIF2064 2065 END FUNCTION interp__1D_cubic_coef2066 !> @endcode2067 !-------------------------------------------------------------------2068 !> @brief2069 !> This subroutine2070 !>2071 !> @details2072 !>2073 !> @author J.Paul2074 !> - Nov, 2013- Initial Version2075 !2076 !> @param[in]2077 !> @param[out]2078 !-------------------------------------------------------------------2079 !> @code2080 SUBROUTINE interp__1D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &2081 & ld_even, id_rhoi )2082 IMPLICIT NONE2083 ! Argument2084 REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value2085 INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect2086 REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef2087 REAL(dp) , INTENT(IN ) :: dd_fill2088 LOGICAL , INTENT(IN ) :: ld_even2089 INTEGER(I4) , INTENT(IN ) :: id_rhoi2090 2091 ! local variable2092 INTEGER(i4), DIMENSION(1) :: il_shape2093 2094 REAL(dp) , DIMENSION(4) :: dl_vect2095 REAL(dp) :: dl_dx2096 REAL(dp) :: dl_x2097 REAL(dp) :: dl_x22098 REAL(dp) :: dl_x32099 2100 ! loop indices2101 INTEGER(i4) :: ji2102 !----------------------------------------------------------------2103 2104 IF( SIZE(dd_coef(:)) /= 4 )THEN2105 CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&2106 & "coef table. shape should be (/4/)")2107 !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &2108 ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN2109 CALL logger_error("INTERP CUBIC FILL: fill value detected. "//&2110 & "can not compute interpolation")2111 ELSE2112 2113 il_shape(:)=SHAPE(dd_value(:))2114 2115 dl_dx=1./REAL(id_rhoi)2116 2117 DO ji=1,il_shape(1)2118 2119 IF(id_detect(ji)==1)THEN2120 2121 IF( ld_even )THEN2122 dl_x=(ji-1)*dl_dx - dl_dx*0.52123 ELSE ! odd refinement2124 dl_x=(ji-1)*dl_dx2125 ENDIF2126 dl_x2=dl_x*dl_x2127 dl_x3=dl_x2*dl_x2128 2129 dl_vect(:)=(/1._dp, dl_x, dl_x2, dl_x3 /)2130 2131 dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))2132 id_detect(ji)=02133 2134 ENDIF2135 2136 ENDDO2137 2138 ENDIF2139 2140 END SUBROUTINE interp__1D_cubic_fill2141 !> @endcode2142 !-------------------------------------------------------------------2143 !> @brief2144 !> This subroutine2145 !>2146 !> @details2147 !>2148 !> @author J.Paul2149 !> - Nov, 2013- Initial Version2150 !2151 !> @param[in]2152 !> @param[out]2153 !-------------------------------------------------------------------2154 !> @code2155 FUNCTION interp__1D_linear_coef( dd_value, dd_fill )2156 IMPLICIT NONE2157 ! Argument2158 REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value2159 REAL(dp) , INTENT(IN) :: dd_fill2160 2161 ! function2162 REAL(dp), DIMENSION(2) :: interp__1D_linear_coef2163 2164 ! local variable2165 2166 REAL(dp), DIMENSION(2,2), PARAMETER :: dp_matrix = RESHAPE( &2167 & (/ 1 , 0 ,&2168 -1 , 1 /), &2169 & (/ 2, 2 /) )2170 2171 REAL(dp), DIMENSION(2) :: dl_vect2172 2173 !----------------------------------------------------------------2174 2175 IF( ANY(SHAPE(dd_value(:))/= 2) )THEN2176 CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&2177 & "input tables. shape should be (/2/)")2178 ELSEIF( ANY(dd_value(:)==dd_fill) )THEN2179 CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&2180 & "can not compute coefficient.")2181 ELSE2182 2183 dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )2184 2185 interp__1D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))2186 2187 ENDIF2188 2189 END FUNCTION interp__1D_linear_coef2190 !> @endcode2191 !-------------------------------------------------------------------2192 !> @brief2193 !> This subroutine2194 !>2195 !> @details2196 !>2197 !> @author J.Paul2198 !> - Nov, 2013- Initial Version2199 !2200 !> @param[in]2201 !> @param[out]2202 !-------------------------------------------------------------------2203 !> @code2204 SUBROUTINE interp__1D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &2205 & ld_even, id_rhoi )2206 IMPLICIT NONE2207 ! Argument2208 REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value2209 INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect2210 REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef2211 REAL(dp) , INTENT(IN ) :: dd_fill2212 LOGICAL , INTENT(IN ) :: ld_even2213 INTEGER(I4) , INTENT(IN ) :: id_rhoi2214 2215 ! local variable2216 INTEGER(i4), DIMENSION(1) :: il_shape2217 2218 REAL(dp) , DIMENSION(2) :: dl_vect2219 REAL(dp) :: dl_dx2220 REAL(dp) :: dl_x2221 2222 ! loop indices2223 INTEGER(i4) :: ji2224 !----------------------------------------------------------------2225 2226 IF( SIZE(dd_coef(:)) /= 2 )THEN2227 CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&2228 & "coef table. shape should be (/2/)")2229 !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &2230 ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN2231 CALL logger_error("INTERP LINEAR FILL: fill value detected. "//&2232 & "can not compute interpolation")2233 ELSE2234 2235 il_shape(:)=SHAPE(dd_value)2236 2237 dl_dx=1./REAL(id_rhoi)2238 2239 DO ji=1,il_shape(1)2240 2241 IF(id_detect(ji)==1)THEN2242 2243 IF( ld_even )THEN2244 dl_x=(ji-1)*dl_dx - dl_dx*0.52245 ELSE ! odd refinement2246 dl_x=(ji-1)*dl_dx2247 ENDIF2248 2249 dl_vect(:)=(/1._dp, dl_x /)2250 2251 dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))2252 id_detect(ji)=02253 2254 ENDIF2255 2256 ENDDO2257 2258 ENDIF2259 2260 END SUBROUTINE interp__1D_linear_fill2261 !> @endcode2262 !-------------------------------------------------------------------2263 !> @brief2264 !> This subroutine2265 !>2266 !> @details2267 !>2268 !> @author J.Paul2269 !> - Nov, 2013- Initial Version2270 !2271 !> @param[in]2272 !> @param[out]2273 !-------------------------------------------------------------------2274 !> @code2275 SUBROUTINE interp__1D_nearest_fill( dd_value, id_detect )2276 IMPLICIT NONE2277 ! Argument2278 REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value2279 INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect2280 2281 ! local variable2282 INTEGER(i4), DIMENSION(1) :: il_shape2283 2284 INTEGER(i4) :: il_i12285 INTEGER(i4) :: il_i22286 2287 INTEGER(i4) :: il_half12288 2289 ! loop indices2290 INTEGER(i4) :: ji2291 !----------------------------------------------------------------2292 2293 il_shape(:)=SHAPE(dd_value)2294 2295 il_i1=12296 il_i2=il_shape(1)2297 2298 il_half1=CEILING(il_shape(1)*0.5)2299 2300 DO ji=1,il_half12301 2302 ! lower left point2303 IF(id_detect(ji)==1)THEN2304 2305 dd_value( ji)=dd_value(il_i1)2306 id_detect(ji)=02307 2308 ENDIF2309 2310 ! lower right point2311 IF(id_detect(il_shape(1)-ji+1)==1)THEN2312 2313 dd_value( il_shape(1)-ji+1)=dd_value(il_i2)2314 id_detect(il_shape(1)-ji+1)=02315 2316 ENDIF2317 2318 ENDDO2319 2320 END SUBROUTINE interp__1D_nearest_fill2321 !> @endcode2322 958 END MODULE interp
Note: See TracChangeset
for help on using the changeset viewer.