|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 MODULE remap_bicubic_reduced 00002 00003 !----------------------------------------------------------------------- 00004 ! BOP 00005 ! 00006 ! !MODULE: remap_bicubic_reduced 00007 ! 00008 ! !USES: 00009 USE mod_kinds_oasis ! defines common data types 00010 USE constants ! defines common constants 00011 USE grids ! module containing grid info 00012 USE remap_vars ! module containing remap info 00013 00014 ! !PUBLIC TYPES: 00015 USE mod_unit 00016 IMPLICIT NONE 00017 ! !PUBLIC MEMBER FUNCTIONS: 00018 ! 00019 ! !PUBLIC DATA MEMBERS: 00020 00021 ! !DESCRIPTION: 00022 ! This routine computes the weights for a bicubic interpolation 00023 ! with a reduced grid. Computes mappings from grid1 to grid2. 00024 ! 00025 ! !REVISION HISTORY: 00026 ! 2002.10.21 J.Ghattas created 00027 ! 00028 ! EOP 00029 !----------------------------------------------------------------------- 00030 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 00031 ! $Author: coquart $ 00032 !----------------------------------------------------------------------- 00033 00034 00035 00036 CONTAINS 00037 00038 !*********************************************************************** 00039 00040 00041 !----------------------------------------------------------------------- 00042 ! BOP 00043 ! 00044 ! !IROUTINE: remap_bicub_reduced 00045 ! 00046 ! !INTERFACE: 00047 00048 SUBROUTINE remap_bicub_reduced(ld_extrapdone) 00049 00050 ! !USES: 00051 00052 ! !RETURN VALUE: 00053 00054 ! !PARAMETERS: 00055 00056 LOGICAL, INTENT(in) :: 00057 ld_extrapdone ! logical, true if EXTRAP done on field 00058 LOGICAL :: 00059 ll_nnei ! true if extra search is done 00060 00061 INTEGER (KIND=ip_intwp_p), DIMENSION(4,4) :: 00062 ila_src_add ! address for source points non-masked 00063 00064 INTEGER (KIND=ip_intwp_p), DIMENSION(4) :: 00065 ila_nbr_found ! nrb of points found on each latitude 00066 00067 INTEGER (KIND=ip_intwp_p) :: 00068 ib_i, ! iter index 00069 ib_dst_add, ! destination address, target point 00070 il_count, ! nbr of latitudes with found points 00071 il_min, il_max, bin ! begin and end for distances calculation 00072 00073 REAL (KIND=ip_realwp_p), DIMENSION(4,4) :: 00074 rla_src_lons, ! longitudes for the points 'ila_src_add' 00075 rla_weight, ! bicubic weights for the points 'ila_src_add' 00076 rla_wght_lon ! temp. weights 00077 00078 REAL (KIND=ip_realwp_p), DIMENSION(4) :: 00079 rla_src_lats, ! latitudes for the points 'ila_src_add' 00080 rla_lats_temp, ! temp. latitudes 00081 rla_wght_lat, rla_wght_temp! temp. weights 00082 00083 REAL (KIND=ip_realwp_p) :: 00084 rl_plat, rl_plon ! latitude and longitude for destination address 00085 00086 REAL (KIND=ip_realwp_p) :: ! variables for distances calculation 00087 rl_coslat_dst, rl_sinlat_dst, 00088 rl_coslon_dst, rl_sinlon_dst, 00089 rl_distance, arg 00090 00091 REAL (KIND=ip_realwp_p), DIMENSION(2) :: 00092 rla_dist ! lat distances to point cible 00093 00094 INTEGER (KIND=ip_intwp_p), DIMENSION(4) :: 00095 ila_add_dist ! temporary addr. for distances 00096 00097 LOGICAL :: ll_linear ! flag 00098 00099 00100 ! !DESCRIPTION: 00101 ! This routine computes the weights for a bicubic interpolation 00102 ! with a reduced grid. Computes mappings from grid1 to grid2. 00103 ! 00104 ! !REVISION HISTORY: 00105 ! 2002.10.21 J.Ghattas created 00106 ! 00107 ! EOP 00108 !----------------------------------------------------------------------- 00109 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 00110 ! $Author: coquart $ 00111 !----------------------------------------------------------------------- 00112 ! 00113 IF (nlogprt .GE. 2) THEN 00114 WRITE (UNIT = nulou,FMT = *)' ' 00115 WRITE (UNIT = nulou,FMT = *) 'Entering routine remap_bicub_reduced' 00116 CALL FLUSH(nulou) 00117 ENDIF 00118 ! 00119 ll_nnei = .true. 00120 ! 00121 ! Loop over destination grid 00122 ! 00123 00124 DO ib_dst_add = 1, grid2_size ! for each target point 00125 ll_linear=.false. 00126 IF (.NOT. grid2_mask(ib_dst_add)) THEN 00127 CYCLE ! target point is masked 00128 END IF 00129 00130 rl_plat = grid2_center_lat(ib_dst_add) 00131 rl_plon = grid2_center_lon(ib_dst_add) 00132 00133 00134 ! 00135 ! Search for non-masked points among the closest 16 points 00136 ! on source grid (grid1) 00137 ! 00138 00139 CALL grid_search_16_points(ila_src_add, rla_src_lats, rla_src_lons,& 00140 ila_nbr_found, bin, rl_plat, & 00141 rl_plon, ld_extrapdone) 00142 00143 00144 ! 00145 ! If there is no point found, search the neaerst 00146 ! non-masked point 00147 ! 00148 00149 IF (SUM(ila_nbr_found)==0) THEN 00150 IF (ll_nnei .eqv. .TRUE. ) THEN 00151 IF (nlogprt .GE. 2) THEN 00152 WRITE(nulou,*) ' ' 00153 WRITE(nulou,*) & 00154 'All 16 surrounding source grid points are masked' 00155 WRITE(nulou,*) 'for target point ',ib_dst_add 00156 WRITE(nulou,*) 'with longitude and latitude', rl_plon, rl_plat 00157 WRITE(nulou,*) 'Using the nearest non-masked neighbour.' 00158 WRITE(nulou,*) ' ' 00159 CALL FLUSH(nulou) 00160 ENDIF 00161 00162 ! Search the nearest point in bin [il_min:il_max] 00163 IF (bin==0 .or. bin==1) THEN 00164 il_min=1 00165 il_max=bin_addr1_r(2,3) 00166 ELSE IF (bin==num_srch_red .or. bin==num_srch_red-1) THEN 00167 il_min=bin_addr1_r(1,num_srch_red-2) 00168 il_max=bin_addr1_r(2,num_srch_red) 00169 ELSE 00170 il_min=bin_addr1_r(1,bin-1)+1 00171 il_max=bin_addr1_r(2,bin+2) 00172 END IF 00173 00174 rl_coslat_dst = COS(rl_plat) 00175 rl_sinlat_dst = SIN(rl_plat) 00176 rl_coslon_dst = COS(rl_plon) 00177 rl_sinlon_dst = SIN(rl_plon) 00178 00179 rla_weight(1,1) = bignum 00180 ila_src_add(1,1) = 0 00181 !cdir novector 00182 DO ib_i=il_min, il_max 00183 IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN 00184 arg = rl_coslat_dst*COS(grid1_center_lat(ib_i))* & 00185 (rl_coslon_dst*COS(grid1_center_lon(ib_i)) + & 00186 rl_sinlon_dst*SIN(grid1_center_lon(ib_i)))+& 00187 rl_sinlat_dst*SIN(grid1_center_lat(ib_i)) 00188 IF (arg < -1.0d0) THEN 00189 arg = -1.0d0 00190 ELSE IF (arg > 1.0d0) THEN 00191 arg = 1.0d0 00192 END IF 00193 rl_distance = ACOS(arg) 00194 IF (rl_distance < rla_weight(1,1)) THEN 00195 rla_weight(1,1) = rl_distance 00196 ila_src_add(1,1) = ib_i 00197 END IF 00198 END IF 00199 END DO 00200 rla_weight(:,:) = 0 00201 rla_weight(1,1) = 1 00202 00203 CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight) 00204 IF (nlogprt .GE. 2) THEN 00205 WRITE(nulou,*) & 00206 'Nearest non masked neighbour is source point ', & 00207 ila_src_add(1,1) 00208 WRITE(nulou,*) 'with longitude and latitude', & 00209 grid1_center_lon(ila_src_add(1,1)), & 00210 grid1_center_lat(ila_src_add(1,1)) 00211 WRITE(nulou,*) ' ' 00212 ENDIF 00213 CYCLE 00214 END IF 00215 ENDIF 00216 00217 rla_weight(:,:) = 0 00218 ! if there is only one point found, save it 00219 IF (SUM(ila_nbr_found)==1) THEN 00220 DO ib_i=1,4 00221 IF (ila_nbr_found(ib_i)==1) THEN 00222 rla_weight(ib_i,1)=1 00223 EXIT 00224 END IF 00225 END DO 00226 CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight) 00227 CYCLE 00228 END IF 00229 00230 ! if there are only 2 points found, distance weighted average 00231 IF (SUM(ila_nbr_found)==2) THEN 00232 rl_coslat_dst = COS(rl_plat) 00233 rl_sinlat_dst = SIN(rl_plat) 00234 rl_coslon_dst = COS(rl_plon) 00235 rl_sinlon_dst = SIN(rl_plon) 00236 00237 rl_distance=0 ! count of total distance 00238 DO ib_i=1,4 00239 IF (ila_nbr_found(ib_i) > 0) THEN 00240 arg = rl_coslat_dst*COS(rla_src_lats(ib_i))* & 00241 (rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + & 00242 rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+& 00243 rl_sinlat_dst*SIN(rla_src_lats(ib_i)) 00244 IF (arg < -1.0d0) THEN 00245 arg = -1.0d0 00246 ELSE IF (arg > 1.0d0) THEN 00247 arg = 1.0d0 00248 END IF 00249 rla_weight(ib_i,1) = ACOS(arg) 00250 rl_distance = rl_distance+rla_weight(ib_i,1) 00251 IF (ila_nbr_found(ib_i)==2) THEN 00252 arg = rl_coslat_dst*COS(rla_src_lats(ib_i))* & 00253 (rl_coslon_dst*COS(rla_src_lons(ib_i,2)) + & 00254 rl_sinlon_dst*SIN(rla_src_lons(ib_i,2)))+& 00255 rl_sinlat_dst*SIN(rla_src_lats(ib_i)) 00256 IF (arg < -1.0d0) THEN 00257 arg = -1.0d0 00258 ELSE IF (arg > 1.0d0) THEN 00259 arg = 1.0d0 00260 END IF 00261 rla_weight(ib_i,2) = ACOS(arg) 00262 rl_distance = rl_distance+rla_weight(ib_i,2) 00263 END IF 00264 END IF 00265 END DO 00266 rla_weight=rla_weight/rl_distance 00267 00268 CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight) 00269 CYCLE 00270 END IF 00271 00272 ! Some case exceptional when just one point per line found 00273 00274 IF (ila_nbr_found(1)==1) THEN ! elimination of point 00275 ila_nbr_found(1)=0 00276 ila_src_add(1,1)=0 00277 END IF 00278 IF (ila_nbr_found(4)==1) THEN 00279 ila_nbr_found(4)=0 00280 ila_src_add(4,1)=0 00281 END IF 00282 00283 00284 00285 IF (ila_nbr_found(2)==1 .OR. ila_nbr_found(3)==1) THEN 00286 ila_add_dist(:)=4 00287 rla_dist(:)=bignum 00288 00289 ! search for the 2 nearest points or line of points 00290 DO ib_i=1,4 00291 IF (ila_nbr_found(ib_i) > 1) THEN 00292 rl_distance=ABS(rla_src_lats(ib_i)-rl_plat) 00293 ELSE IF (ila_nbr_found(ib_i)==1) THEN 00294 rl_coslat_dst = COS(rl_plat) 00295 rl_sinlat_dst = SIN(rl_plat) 00296 rl_coslon_dst = COS(rl_plon) 00297 rl_sinlon_dst = SIN(rl_plon) 00298 arg = rl_coslat_dst*COS(rla_src_lats(ib_i))* & 00299 (rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + & 00300 rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+& 00301 rl_sinlat_dst*SIN(rla_src_lats(ib_i)) 00302 IF (arg < -1.0d0) THEN 00303 arg = -1.0d0 00304 ELSE IF (arg > 1.0d0) THEN 00305 arg = 1.0d0 00306 END IF 00307 rl_distance= ACOS(arg) 00308 ELSE 00309 rl_distance=bignum 00310 END IF 00311 00312 IF (rl_distance < rla_dist(1)) THEN 00313 ila_add_dist(2)=ila_add_dist(1) 00314 ila_add_dist(1)=ib_i 00315 rla_dist(2)=rla_dist(1) 00316 rla_dist(1)=rl_distance 00317 ELSE IF (rl_distance < rla_dist(2)) THEN 00318 ila_add_dist(2)=ib_i 00319 rla_dist(2)=rl_distance 00320 END IF 00321 END DO 00322 00323 IF (ila_nbr_found(ila_add_dist(1))>1 .AND. & 00324 ila_nbr_found(ila_add_dist(2))>1) THEN 00325 ! linearie 00326 ll_linear=.true. 00327 ELSE 00328 ! do distance weighted averege 00329 rla_wght_lon(:,:)=0 00330 DO ib_i=1,2 00331 SELECT CASE (ila_nbr_found(ila_add_dist(ib_i))) 00332 CASE (4) 00333 CALL calcul_wght_irreg(rla_src_lons(ila_add_dist(ib_i),:),& 00334 rl_plon, rla_wght_lon(ila_add_dist(ib_i),:)) 00335 rla_wght_lon(ila_add_dist(ib_i),:)=& 00336 rla_wght_lon(ila_add_dist(ib_i),:)/& 00337 rla_dist(ib_i) 00338 CASE (3) 00339 CALL calcul_wght_3(rla_src_lons(ila_add_dist(ib_i),1:3),& 00340 rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:3)) 00341 rla_wght_lon(ila_add_dist(ib_i),1:3)=& 00342 rla_wght_lon(ila_add_dist(ib_i),1:3)/& 00343 rla_dist(ib_i) 00344 CASE (2) 00345 CALL calcul_wght_2(rla_src_lons(ila_add_dist(ib_i),1:2),& 00346 rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:2)) 00347 rla_wght_lon(ila_add_dist(ib_i),1:2)=& 00348 rla_wght_lon(ila_add_dist(ib_i),1:2)/& 00349 rla_dist(ib_i) 00350 CASE (1) 00351 rla_wght_lon(ila_add_dist(ib_i),1)=1/rla_dist(ib_i) 00352 END SELECT 00353 END DO 00354 rl_distance=0 00355 DO ib_i=1,4 00356 rl_distance=rl_distance + sum(rla_wght_lon(ib_i,:)) 00357 END DO 00358 rla_weight(:,:)=rla_wght_lon(:,:)/rl_distance 00359 00360 CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight) 00361 CYCLE 00362 END IF 00363 END IF 00364 00365 ! 00366 ! Calculation of weights for longitudes 00367 ! 00368 00369 rla_wght_lon(:,:)=0 00370 DO ib_i=1,4 00371 SELECT CASE (ila_nbr_found(ib_i)) 00372 CASE (4) 00373 CALL calcul_wght_irreg(rla_src_lons(ib_i,:), rl_plon, & 00374 rla_wght_lon(ib_i,:)) 00375 CASE (3) 00376 CALL calcul_wght_3(rla_src_lons(ib_i,1:3), rl_plon, & 00377 rla_wght_lon(ib_i,1:3)) 00378 CASE (2) 00379 CALL calcul_wght_2(rla_src_lons(ib_i,1:2), rl_plon, & 00380 rla_wght_lon(ib_i,1:2)) 00381 END SELECT 00382 END DO 00383 00384 00385 IF (ll_linear) THEN 00386 rla_wght_lat(:)=0 00387 CALL calcul_wght_2(rla_src_lats(ila_add_dist(:)), rl_plat, & 00388 rla_wght_temp(1:2)) 00389 rla_wght_lat(ila_add_dist(1))=rla_wght_temp(1) 00390 rla_wght_lat(ila_add_dist(2))=rla_wght_temp(2) 00391 DO ib_i=1,4 00392 rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:) 00393 END DO 00394 00395 CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight) 00396 CYCLE 00397 END IF 00398 00399 00400 ! 00401 ! Calculation of weights for latitudes 00402 ! 00403 00404 il_count=0 00405 DO ib_i=1,4 00406 IF (ila_nbr_found(ib_i)/=0) THEN 00407 il_count=il_count+1 00408 rla_lats_temp(il_count)=rla_src_lats(ib_i) 00409 END IF 00410 END DO 00411 00412 SELECT CASE (il_count) 00413 CASE (4) 00414 CALL calcul_wght_irreg(rla_lats_temp, rl_plat, rla_wght_temp(:)) 00415 CASE (3) 00416 CALL calcul_wght_3(rla_lats_temp(1:3), rl_plat, rla_wght_temp(1:3)) 00417 CASE (2) 00418 CALL calcul_wght_2(rla_lats_temp(1:2), rl_plat, rla_wght_temp(1:2)) 00419 CASE (1) 00420 rla_wght_temp(1)=1 00421 END SELECT 00422 00423 il_count=0 00424 DO ib_i=1,4 00425 IF (ila_nbr_found(ib_i)/=0) THEN 00426 il_count=il_count+1 00427 rla_wght_lat(ib_i)=rla_wght_temp(il_count) 00428 ELSE 00429 rla_wght_lat(ib_i)=0 00430 END IF 00431 END DO 00432 00433 ! 00434 ! Calculation of total weight, elementwise multiplication 00435 ! 00436 00437 DO ib_i=1,4 00438 rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:) 00439 END DO 00440 00441 CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight) 00442 00443 END DO 00444 ! 00445 IF (nlogprt .GE. 2) THEN 00446 WRITE (UNIT = nulou,FMT = *)' ' 00447 WRITE (UNIT = nulou,FMT = *) 'Leaving routine remap_bicub_reduced' 00448 CALL FLUSH(nulou) 00449 ENDIF 00450 ! 00451 END SUBROUTINE remap_bicub_reduced 00452 00453 00454 !----------------------------------------------------------------------- 00455 ! BOP 00456 ! 00457 ! !IROUTINE: grid_search_16_points 00458 ! 00459 ! !INTERFACE: 00460 ! 00461 SUBROUTINE grid_search_16_points(ida_src_add, rda_src_lats, rda_src_lons,& 00462 ida_nbr_found, bin, rd_plat, & 00463 rd_plon, ld_extrapdone) 00464 ! 00465 ! !USES: 00466 ! 00467 ! !RETURN VALUE: 00468 ! 00469 INTEGER (KIND=ip_intwp_p), DIMENSION(4,4), INTENT(out) :: 00470 ida_src_add ! searched addresses of the unmasked points enclosing 00471 ! target point 00472 00473 REAL (KIND=ip_realwp_p), DIMENSION(4,4), INTENT(out) :: 00474 rda_src_lons ! longitudes of the searched points 00475 00476 REAL (KIND=ip_realwp_p), DIMENSION(4), INTENT(out) :: 00477 rda_src_lats ! latitudes of the searched points 00478 00479 INTEGER (KIND=ip_intwp_p), DIMENSION(4), INTENT(out) :: 00480 ida_nbr_found ! indicates for each line how many points found 00481 00482 INTEGER (KIND=ip_intwp_p), INTENT(out) :: 00483 bin ! actual bin for target point 00484 ! 00485 ! !PARAMETERS: 00486 ! 00487 REAL (KIND=ip_realwp_p), INTENT(in) :: 00488 rd_plat, ! latitude of the target point 00489 rd_plon ! longitude of the target point 00490 00491 LOGICAL, INTENT(in) :: ld_extrapdone ! true if extrapolation done 00492 00493 INTEGER (KIND=ip_intwp_p) :: 00494 ib_k, ib_j, ib_i, ! iteration indices 00495 il_min, il_max, il_inter ! begin and end for actual bin 00496 00497 INTEGER (KIND=ip_intwp_p), DIMENSION(4,2) :: 00498 ila_corners ! temp addresses for bins 00499 00500 ! 00501 ! !DESCRIPTION: 00502 ! This routine finds the location of the target point in the source 00503 ! grid and returns those of the 16 nearest points that are unmasked. 00504 ! The source grid is a reduced grid. 00505 ! 00506 ! !REVISION HISTORY: 00507 ! 2002.10.21 J.Ghattas created 00508 ! 00509 ! EOP 00510 !----------------------------------------------------------------------- 00511 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 00512 ! $Author: coquart $ 00513 !----------------------------------------------------------------------- 00514 00515 00516 ! 00517 ! serch of actual bin 00518 ! 00519 00520 00521 IF (rd_plat > bin_lats_r(1,1)) THEN ! norther of the first bin 00522 bin=0 00523 ila_corners(1:2,1:2)= 0 00524 ila_corners(3,1)= bin_addr1_r(1,1)+1 00525 ila_corners(3,2)= bin_addr1_r(2,1) 00526 ila_corners(4,1)= bin_addr1_r(1,2) 00527 ila_corners(4,2)= bin_addr1_r(2,2) 00528 00529 ELSE IF (rd_plat > bin_lats_r(1,2)) THEN ! in the first bin 00530 bin=1 00531 ila_corners(1,1:2)= 0 00532 ila_corners(2,1)= bin_addr1_r(1,1)+1 00533 ila_corners(2,2)= bin_addr1_r(2,1) 00534 ila_corners(3,1)= bin_addr1_r(1,2) 00535 ila_corners(3,2)= bin_addr1_r(2,2) 00536 ila_corners(4,1)= bin_addr1_r(1,3) 00537 ila_corners(4,2)= bin_addr1_r(2,3) 00538 00539 ELSE IF (rd_plat < bin_lats_r(1,num_srch_red)) THEN 00540 ! South of the last complet bin 00541 bin=num_srch_red 00542 ila_corners(1,1) = bin_addr1_r(1,num_srch_red-1) 00543 ila_corners(1,2) = bin_addr1_r(2,num_srch_red-1) 00544 ila_corners(2,1) = bin_addr1_r(1,num_srch_red) 00545 ila_corners(2,2) = bin_addr1_r(2,num_srch_red) 00546 ila_corners(3:4,1:2) = 0 00547 00548 ELSE IF (rd_plat < bin_lats_r(1,num_srch_red-1)) THEN 00549 ! in the last bin which is complet 00550 ! the bin (num_srch_red-1) is the last bin which is complet 00551 bin=num_srch_red-1 00552 ila_corners(1,1) = bin_addr1_r(1,num_srch_red-2) 00553 ila_corners(1,2) = bin_addr1_r(2,num_srch_red-2) 00554 ila_corners(2,1) = bin_addr1_r(1,num_srch_red-1) 00555 ila_corners(2,2) = bin_addr1_r(2,num_srch_red-1) 00556 ila_corners(3,1) = bin_addr1_r(1,num_srch_red) 00557 ila_corners(3,2) = bin_addr1_r(2,num_srch_red) 00558 ila_corners(4,1:2) = 0 00559 ELSE 00560 il_min=2 00561 il_max=num_srch_red-1 00562 DO WHILE (il_min /= il_max-1) 00563 il_inter=(il_max-il_min)/2 + il_min 00564 IF (rd_plat <= bin_lats_r(1,il_min) .and. & 00565 rd_plat > bin_lats_r(1,il_inter)) THEN 00566 il_max=il_inter 00567 ELSE 00568 il_min=il_inter 00569 END IF 00570 END DO 00571 bin=il_min 00572 ila_corners(1,1) = bin_addr1_r(1,bin-1) 00573 ila_corners(1,2) = bin_addr1_r(2,bin-1) 00574 ila_corners(2,1) = bin_addr1_r(1,bin) 00575 ila_corners(2,2) = bin_addr1_r(2,bin) 00576 ila_corners(3,1) = bin_addr1_r(1,bin+1) 00577 ila_corners(3,2) = bin_addr1_r(2,bin+1) 00578 ila_corners(4,1) = bin_addr1_r(1,bin+2) 00579 ila_corners(4,2) = bin_addr1_r(2,bin+2) 00580 00581 IF (ila_corners(1,1)==0) THEN 00582 ila_corners(1,1)=1 00583 END IF 00584 END IF 00585 00586 DO ib_k=1,4 00587 IF (ila_corners(ib_k,1) .NE. 0) & 00588 rda_src_lats(ib_k)= grid1_center_lat(ila_corners(ib_k,1)) 00589 ENDDO 00590 00591 ! 00592 ! now perform a more detailed search for each line 00593 ! 00594 00595 ida_src_add(:,:)=0 00596 ida_nbr_found(:)=0 00597 rda_src_lons(:,:)=0 00598 00599 DO ib_k=1,4 ! for each line of found points 00600 IF (ila_corners(ib_k,1)==0) THEN 00601 CYCLE 00602 END IF 00603 00604 il_min=ila_corners(ib_k,1) 00605 il_max=ila_corners(ib_k,2) 00606 00607 IF (rd_plon < grid1_center_lon(il_min)) THEN 00608 DO ib_j=il_max-1, il_max 00609 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00610 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00611 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00612 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00613 grid1_center_lon(ib_j)-pi2 00614 END IF 00615 END DO 00616 DO ib_j=il_min, il_min+1 00617 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00618 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00619 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00620 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00621 grid1_center_lon(ib_j) 00622 END IF 00623 END DO 00624 00625 ELSE IF (rd_plon < grid1_center_lon(il_min+1)) THEN 00626 IF (grid1_mask(il_max) .or. ld_extrapdone) THEN 00627 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00628 ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_max 00629 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00630 grid1_center_lon(il_max)-pi2 00631 END IF 00632 DO ib_j=il_min, il_min+2 00633 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00634 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00635 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00636 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00637 grid1_center_lon(ib_j) 00638 END IF 00639 END DO 00640 00641 ELSE IF (rd_plon > grid1_center_lon(il_max)) THEN 00642 DO ib_j=il_max-1, il_max 00643 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00644 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00645 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00646 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00647 grid1_center_lon(ib_j) 00648 END IF 00649 END DO 00650 DO ib_j=il_min, il_min+1 00651 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00652 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00653 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00654 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00655 grid1_center_lon(ib_j)+pi2 00656 END IF 00657 END DO 00658 00659 ELSE IF (rd_plon > grid1_center_lon(il_max-1)) THEN 00660 DO ib_j=il_max-2, il_max 00661 IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN 00662 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00663 ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j 00664 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00665 grid1_center_lon(ib_j) 00666 END IF 00667 END DO 00668 IF (grid1_mask(il_min) .or. ld_extrapdone) THEN 00669 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00670 ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_min 00671 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00672 grid1_center_lon(il_min)+pi2 00673 END IF 00674 00675 ELSE 00676 00677 DO WHILE (il_min/=il_max-1) 00678 il_inter=(il_max-il_min)/2 + il_min 00679 IF (rd_plon >= grid1_center_lon(il_min) .and. & 00680 rd_plon < grid1_center_lon(il_inter)) THEN 00681 il_max=il_inter 00682 ELSE 00683 il_min=il_inter 00684 END IF 00685 END DO 00686 DO ib_i= il_min-1, il_min+2 00687 IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN 00688 ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1 00689 ida_src_add(ib_k,ida_nbr_found(ib_k))=ib_i 00690 rda_src_lons(ib_k,ida_nbr_found(ib_k))= & 00691 grid1_center_lon(ib_i) 00692 END IF 00693 END DO 00694 00695 END IF 00696 00697 END DO 00698 00699 00700 END SUBROUTINE grid_search_16_points 00701 00702 00703 00704 !----------------------------------------------------------------------- 00705 ! BOP 00706 ! 00707 ! !IROUTINE: calcul_wght_irreg 00708 ! 00709 ! !INTERFACE: 00710 ! 00711 SUBROUTINE calcul_wght_irreg(rda_x, rd_pt, rda_wght) 00712 ! 00713 ! !USES: 00714 ! 00715 ! !RETURN VALUE: 00716 ! 00717 REAL (KIND=ip_realwp_p), DIMENSION(4), INTENT(out) :: 00718 rda_wght ! array of weights for the points x 00719 ! 00720 ! !PARAMETERS: 00721 ! 00722 REAL (KIND=ip_realwp_p), DIMENSION(4), INTENT(in) :: 00723 rda_x ! array of positions on source grid, lat or lon 00724 00725 REAL (KIND=ip_realwp_p),INTENT(in) :: 00726 rd_pt ! position of target point to interpolate 00727 00728 REAL (KIND=ip_realwp_p) :: 00729 rl_t1, rl_t2, rl_t3, rl_t4, rl_t5, rl_t6, rl_t7, rl_t8, rl_t9, 00730 rl_u1, rl_u2, rl_u3, rl_u4, 00731 rl_k1, rl_k2, rl_k3, 00732 rl_d1, rl_d2, rl_d3, rl_d4, 00733 rl_c1, rl_c2, rl_c3, rl_c4, 00734 rl_b1, rl_b2, rl_b3, rl_b4, 00735 rl_a1, rl_a2, rl_a3, rl_a4, 00736 rl_y1, rl_y2, rl_y3, 00737 rl_a1_y, rl_a2_y, rl_a3_y, rl_a4_y, 00738 rl_b1_y, rl_b2_y, rl_b3_y, rl_b4_y, 00739 rl_c1_y, rl_c2_y, rl_c3_y, rl_c4_y 00740 ! 00741 ! !DESCRIPTION: 00742 ! Calculates a the weights of four points for a bicubic interpolation. 00743 ! The distances between the points can be irregulier. 00744 ! 00745 ! !REVISION HISTORY: 00746 ! 2002.10.21 J.Ghattas created 00747 ! 00748 ! EOP 00749 !----------------------------------------------------------------------- 00750 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 00751 ! $Author: coquart $ 00752 !----------------------------------------------------------------------- 00753 00754 00755 IF (rda_x(1)/=0.and. rda_x(2)/=0 .and. rda_x(3)/=0 .and. rda_x(4)/=0) THEN 00756 00757 rl_t1 = 1/rda_x(1) - 1/rda_x(2) 00758 rl_t2 = 1/rda_x(1)**2 - 1/rda_x(2)**2 00759 rl_t3 = 1/rda_x(1)**3 - 1/rda_x(2)**3 00760 rl_t4 = 1/rda_x(1) - 1/rda_x(3) 00761 rl_t5 = 1/rda_x(1)**2 - 1/rda_x(3)**2 00762 rl_t6 = 1/rda_x(1)**3 - 1/rda_x(3)**3 00763 rl_t7 = 1/rda_x(1) - 1/rda_x(4) 00764 rl_t8 = 1/rda_x(1)**2 - 1/rda_x(4)**2 00765 rl_t9 = 1/rda_x(1)**3 - 1/rda_x(4)**3 00766 00767 rl_u1 = rl_t2/rl_t1 - rl_t5/rl_t4 00768 rl_u2 = rl_t3/rl_t1 - rl_t6/rl_t4 00769 rl_u3 = rl_t2/rl_t1 - rl_t8/rl_t7 00770 rl_u4 = rl_t3/rl_t1 - rl_t9/rl_t7 00771 00772 rl_k1 = (1/(rl_t1*rl_u1)-1/(rl_t1*rl_u3)) / (rl_u2/rl_u1-rl_u4/rl_u3) 00773 rl_k2 = -1/(rl_t4*rl_u1) / (rl_u2/rl_u1-rl_u4/rl_u3) 00774 rl_k3 = 1/(rl_t7*rl_u3) / (rl_u2/rl_u1-rl_u4/rl_u3) 00775 00776 00777 rl_d1=(rl_k1+rl_k2+rl_k3)/rda_x(1)**3 00778 rl_d2 = -rl_k1/rda_x(2)**3 00779 rl_d3 = -rl_k2/rda_x(3)**3 00780 rl_d4 = -rl_k3/rda_x(4)**3 00781 00782 rl_c1 = 1/rl_u1*(1/(rl_t1*rda_x(1)**3)-1/(rl_t4*rda_x(1)**3)- & 00783 rl_u2*rl_d1) 00784 rl_c2 = 1/rl_u1*(1/(-rl_t1*rda_x(2)**3)-rl_u2*rl_d2) 00785 rl_c3 = 1/rl_u1*(1/(rl_t4*rda_x(3)**3)-rl_u2*rl_d3) 00786 rl_c4 = 1/rl_u1*(-rl_u2*rl_d4) 00787 00788 rl_b1 = 1/rl_t1/rda_x(1)**3-rl_t2/rl_t1*rl_c1-rl_t3/rl_t1*rl_d1 00789 rl_b2 = -1/rl_t1/rda_x(2)**3-rl_t2/rl_t1*rl_c2-rl_t3/rl_t1*rl_d2 00790 rl_b3 = -rl_t2/rl_t1*rl_c3-rl_t3/rl_t1*rl_d3 00791 rl_b4 = -rl_t2/rl_t1*rl_c4-rl_t3/rl_t1*rl_d4 00792 00793 rl_a1 = 1/rda_x(1)**3-1/rda_x(1)*rl_b1-1/rda_x(1)**2*rl_c1- & 00794 1/rda_x(1)**3*rl_d1 00795 rl_a2 = -1/rda_x(1)*rl_b2-1/rda_x(1)**2*rl_c2-1/rda_x(1)**3*rl_d2 00796 rl_a3 = -1/rda_x(1)*rl_b3-1/rda_x(1)**2*rl_c3-1/rda_x(1)**3*rl_d3 00797 rl_a4 = -1/rda_x(1)*rl_b4-1/rda_x(1)**2*rl_c4-1/rda_x(1)**3*rl_d4 00798 00799 ! Weights 00800 rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt + rl_d1 00801 rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt + rl_d2 00802 rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt + rl_d3 00803 rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt + rl_d4 00804 00805 ELSE ! there is one point = 0 00806 00807 rl_d1=0; rl_d2=0; rl_d3=0; rl_d4=0 00808 00809 ! Transformation for each case 00810 IF (rda_x(1)==0) THEN 00811 rl_y1=rda_x(2); rl_y2=rda_x(3); rl_y3=rda_x(4) 00812 rl_d1=1 00813 ELSE IF (rda_x(2)==0) THEN 00814 rl_y1=rda_x(1); rl_y2=rda_x(3); rl_y3=rda_x(4) 00815 rl_d2=1 00816 ELSE IF (rda_x(3)==0) THEN 00817 rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(4) 00818 rl_d3=1 00819 ELSE 00820 rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(3) 00821 rl_d4=1 00822 END IF 00823 00824 ! Solving the system 00825 rl_t1 = 1/rl_y1-1/rl_y2 00826 rl_t2 = 1/rl_y1**2-1/rl_y2**2 00827 rl_t3 = 1/rl_y1-1/rl_y3 00828 rl_t4 = 1/rl_y1**2-1/rl_y3**2 00829 00830 rl_c1_y =(1/rl_y1**3/rl_t1-1/rl_y1**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3) 00831 rl_c2_y = -1/rl_y2**3/rl_t1/(rl_t2/rl_t1-rl_t4/rl_t3) 00832 rl_c3_y = 1/rl_y3**3/rl_t3/(rl_t2/rl_t1-rl_t4/rl_t3) 00833 rl_c4_y=(-1/rl_y1**3/rl_t1+1/rl_y2**3/rl_t1+ & 00834 1/rl_y1**3/rl_t3-1/rl_y3**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3) 00835 00836 rl_b1_y = 1/rl_y1**3/rl_t1 - rl_c1_y*rl_t2/rl_t1 00837 rl_b2_y = -1/rl_y2**3/rl_t1 - rl_c2_y*rl_t2/rl_t1 00838 rl_b3_y = -rl_c3_y*rl_t2/rl_t1 00839 rl_b4_y = -1/rl_y1**3/rl_t1 + 1/rl_y2**3/rl_t1 - rl_c4_y*rl_t2/rl_t1 00840 00841 rl_a1_y = 1/rl_y1**3 - rl_b1_y/rl_y1 - rl_c1_y/rl_y1**2 00842 rl_a2_y = -rl_b2_y/rl_y1 - rl_c2_y/rl_y1**2 00843 rl_a3_y = -rl_b3_y/rl_y1 - rl_c3_y/rl_y1**2 00844 rl_a4_y = -1/rl_y1**3 - rl_b4_y/rl_y1 - rl_c4_y/rl_y1**2 00845 00846 ! Retransformation 00847 IF (rda_x(1)==0) THEN 00848 rl_a1=rl_a4_y; rl_a2=rl_a1_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y 00849 rl_b1=rl_b4_y; rl_b2=rl_b1_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y 00850 rl_c1=rl_c4_y; rl_c2=rl_c1_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y 00851 ELSE IF (rda_x(2)==0) THEN 00852 rl_a1=rl_a1_y; rl_a2=rl_a4_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y 00853 rl_b1=rl_b1_y; rl_b2=rl_b4_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y 00854 rl_c1=rl_c1_y; rl_c2=rl_c4_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y 00855 ELSE IF (rda_x(3)==0) THEN 00856 rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a4_y; rl_a4=rl_a3_y 00857 rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b4_y; rl_b4=rl_b3_y 00858 rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c4_y; rl_c4=rl_c3_y 00859 ELSE 00860 rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a3_y; rl_a4=rl_a4_y 00861 rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b3_y; rl_b4=rl_b4_y 00862 rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c3_y; rl_c4=rl_c4_y 00863 END IF 00864 00865 ! Weights 00866 rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt +rl_d1 00867 rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt +rl_d2 00868 rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt +rl_d3 00869 rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt +rl_d4 00870 00871 END IF 00872 00873 00874 END SUBROUTINE calcul_wght_irreg 00875 00876 !----------------------------------------------------------------------- 00877 ! BOP 00878 ! 00879 ! !IROUTINE: calcul_wght_3 00880 ! 00881 ! !INTERFACE: 00882 00883 SUBROUTINE calcul_wght_3(rda_x, rd_pt, rda_wght) 00884 00885 ! !USES: 00886 00887 ! !RETURN VALUE: 00888 00889 REAL (KIND=ip_realwp_p), DIMENSION(3), INTENT(out) :: 00890 rda_wght ! array of weights for the points x 00891 00892 ! !PARAMETERS: 00893 00894 REAL (KIND=ip_realwp_p), DIMENSION(3), INTENT(in) :: 00895 rda_x ! array of positions on source grid, lat or lon 00896 00897 REAL (KIND=ip_realwp_p), INTENT(in) :: 00898 rd_pt ! position of target point to interpolate 00899 00900 REAL (KIND=ip_realwp_p) :: 00901 rl_c1, rl_c2, rl_c3, 00902 rl_a1, rl_a2, rl_a3, 00903 rl_b1, rl_b2, rl_b3, 00904 rl_t1, rl_t2, rl_t3, rl_t4 00905 00906 ! !DESCRIPTION: 00907 ! Calculates a the weights of 3 points for a parabolic interpolation. 00908 ! 00909 ! !REVISION HISTORY: 00910 ! 2002.10.21 J.Ghattas created 00911 ! 00912 ! EOP 00913 !----------------------------------------------------------------------- 00914 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 00915 ! $Author: coquart $ 00916 !----------------------------------------------------------------------- 00917 00918 00919 IF (rda_x(1)/=0 .and. rda_x(2)/=0 .and. rda_x(3)/=0) THEN 00920 rl_t1 = 1/rda_x(1)-1/rda_x(2) 00921 rl_t2 = 1/rda_x(1)**2-1/rda_x(2)**2 00922 rl_t3 = 1/rda_x(1)-1/rda_x(3) 00923 rl_t4 = 1/rda_x(1)**2-1/rda_x(3)**2 00924 00925 rl_c1 = (1/rda_x(1)**2/rl_t1-1/rda_x(1)**2/rl_t3) / & 00926 (rl_t2/rl_t1-rl_t4/rl_t3) 00927 rl_c2 = -1/rda_x(2)**2/rl_t1 / (rl_t2/rl_t1-rl_t4/rl_t3) 00928 rl_c3 = 1/rda_x(3)**2/rl_t3 / (rl_t2/rl_t1-rl_t4/rl_t3) 00929 00930 rl_b1 = 1/rda_x(1)**2/rl_t1 - rl_c1*rl_t2/rl_t1 00931 rl_b2 = -1/rda_x(2)**2/rl_t1 - rl_c2*rl_t2/rl_t1 00932 rl_b3 = - rl_c3*rl_t2/rl_t1 00933 00934 rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1) - rl_c1/rda_x(1)**2 00935 rl_a2 = - rl_b2/rda_x(1) - rl_c2/rda_x(1)**2 00936 rl_a3 = - rl_b3/rda_x(1) - rl_c3/rda_x(1)**2 00937 00938 00939 ELSE IF (rda_x(1)==0) THEN 00940 rl_c1 = 1; rl_c2 = 0; rl_c3 = 0 00941 rl_b1 = (-1/rda_x(2)**2+1/rda_x(3)**2) / (1/rda_x(2)-1/rda_x(3)) 00942 rl_b2 = 1/rda_x(2)**2 / (1/rda_x(2)-1/rda_x(3)) 00943 rl_b3 = -1/rda_x(3)**2 / (1/rda_x(2)-1/rda_x(3)) 00944 00945 rl_a1 = -1/rda_x(2)**2 - rl_b1/rda_x(2) 00946 rl_a2 = 1/rda_x(2)**2 - rl_b2/rda_x(2) 00947 rl_a3 = - rl_b3/rda_x(2) 00948 00949 ELSE IF (rda_x(2)==0) THEN 00950 00951 rl_c1 = 0; rl_c2 = 1; rl_c3 = 0 00952 rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(3)) 00953 rl_b2 = (-1/rda_x(1)**2+1/rda_x(3)**2) / (1/rda_x(1)-1/rda_x(3)) 00954 rl_b3 = -1/rda_x(3)**2 / (1/rda_x(1)-1/rda_x(3)) 00955 00956 rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1) 00957 rl_a2 = -1/rda_x(1)**2 - rl_b2/rda_x(1) 00958 rl_a3 = - rl_b3/rda_x(1) 00959 00960 ELSE !rda_x(3)==0 00961 rl_c1 = 0; rl_c2 = 0; rl_c3 = 1 00962 rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(2)) 00963 rl_b2 = -1/rda_x(2)**2 / (1/rda_x(1)-1/rda_x(2)) 00964 rl_b3 = (-1/rda_x(1)**2+1/rda_x(2)**2) / (1/rda_x(1)-1/rda_x(2)) 00965 00966 rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1) 00967 rl_a2 = - rl_b2/rda_x(1) 00968 rl_a3 = -1/rda_x(1)**2 - rl_b3/rda_x(1) 00969 00970 00971 END IF 00972 00973 ! Weights 00974 rda_wght(1) = rl_a1*rd_pt**2 + rl_b1*rd_pt + rl_c1 00975 rda_wght(2) = rl_a2*rd_pt**2 + rl_b2*rd_pt + rl_c2 00976 rda_wght(3) = rl_a3*rd_pt**2 + rl_b3*rd_pt + rl_c3 00977 00978 00979 END SUBROUTINE calcul_wght_3 00980 00981 00982 !----------------------------------------------------------------------- 00983 ! BOP 00984 ! 00985 ! !IROUTINE: calcul_wght_2 00986 ! 00987 ! !INTERFACE: 00988 00989 SUBROUTINE calcul_wght_2(rda_x, rd_pt, rda_wght) 00990 00991 ! !USES: 00992 00993 ! !RETURN VALUE: 00994 00995 REAL (KIND=ip_realwp_p), DIMENSION(2), INTENT(out) :: 00996 rda_wght ! array of weights for the points x 00997 00998 ! !PARAMETERS: 00999 01000 REAL (KIND=ip_realwp_p), DIMENSION(2), INTENT(in) :: 01001 rda_x ! array of positions on source grid, lat or lon 01002 01003 REAL (KIND=ip_realwp_p), INTENT(in) :: 01004 rd_pt ! position of target point to interpolate 01005 01006 REAL (KIND=ip_realwp_p) :: rl_b1, rl_b2, rl_a1, rl_a2 01007 01008 ! !DESCRIPTION: 01009 ! Calculates a the weights of 2 points for a linair interpolation. 01010 ! 01011 ! !REVISION HISTORY: 01012 ! 2002.10.21 J.Ghattas created 01013 ! 01014 ! EOP 01015 !----------------------------------------------------------------------- 01016 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 01017 ! $Author: coquart $ 01018 !----------------------------------------------------------------------- 01019 01020 01021 IF (rda_x(1)/=0 .and. rda_x(2)/=0) THEN 01022 rl_b1 = 1/(1-rda_x(1)/rda_x(2)) 01023 rl_b2 = -1/(rda_x(2)/rda_x(1)-1) 01024 rl_a1 = 1/rda_x(1) - rl_b1/rda_x(1) 01025 rl_a2 = - rl_b2/rda_x(1) 01026 01027 ELSE IF (rda_x(1)==0) THEN 01028 rl_b1=1; rl_b2=0 01029 rl_a1=-1/rda_x(2) 01030 rl_a2=1/rda_x(2) 01031 ELSE 01032 rl_b1=0; rl_b2=1 01033 rl_a1=1/rda_x(1) 01034 rl_a2=-1/rda_x(1) 01035 END IF 01036 01037 rda_wght(1) = rl_a1*rd_pt + rl_b1 01038 rda_wght(2) = rl_a2*rd_pt + rl_b2 01039 01040 END SUBROUTINE calcul_wght_2 01041 01042 01043 !----------------------------------------------------------------------- 01044 ! BOP 01045 ! 01046 ! !IROUTINE: store_link_bicub 01047 ! 01048 ! !INTERFACE: 01049 01050 SUBROUTINE store_link_bicub(id_dst_add, ida_src_add, rda_wght) 01051 01052 ! !USES: 01053 01054 ! !RETURN VALUE: 01055 01056 ! !PARAMETERS: 01057 01058 INTEGER (KIND=ip_intwp_p), INTENT(in) :: 01059 id_dst_add ! address on destination grid 01060 01061 INTEGER (KIND=ip_intwp_p), DIMENSION(4,4), INTENT(in) :: 01062 ida_src_add ! addresses for links on source grid 01063 01064 REAL (KIND=ip_realwp_p), DIMENSION(4,4), INTENT(in) :: 01065 rda_wght ! array of remapping weights for these links 01066 01067 INTEGER (KIND=ip_intwp_p) :: ib_i, 01068 il_num_links_old ! placeholder for old link number 01069 01070 INTEGER (KIND=ip_intwp_p), DIMENSION(16) :: 01071 ila_src_add ! reshaped addresses 01072 01073 REAL (KIND=ip_realwp_p), DIMENSION(16) :: 01074 rla_wght ! reshaped weights 01075 01076 ! !DESCRIPTION: 01077 ! This routine stores the addresses and weights for 16 links associated 01078 ! with one destination point in the appropriate address. 01079 ! 01080 ! !REVISION HISTORY: 01081 ! 2002.10.21 J.Ghattas created 01082 ! 01083 ! EOP 01084 !----------------------------------------------------------------------- 01085 ! $Id: remap_bicubic_reduced.F90 3500 2012-03-22 14:58:30Z coquart $ 01086 ! $Author: coquart $ 01087 !----------------------------------------------------------------------- 01088 01089 01090 ! 01091 ! Increment number of links and check if remap arrays need 01092 ! to be increased to accomodate the new link. then store the link. 01093 ! 01094 01095 il_num_links_old = num_links_map1 01096 num_links_map1 = il_num_links_old + 16 01097 01098 IF (num_links_map1 > max_links_map1) THEN 01099 CALL resize_remap_vars(1,MAX(resize_increment,16)) 01100 END IF 01101 01102 ila_src_add=RESHAPE(ida_src_add,(/16/)) 01103 rla_wght=RESHAPE(rda_wght,(/16/)) 01104 01105 DO ib_i=1,16 01106 grid1_add_map1(il_num_links_old+ib_i) = ila_src_add(ib_i) 01107 grid2_add_map1(il_num_links_old+ib_i) = id_dst_add 01108 wts_map1(1,il_num_links_old+ib_i) = rla_wght(ib_i) 01109 END DO 01110 01111 END SUBROUTINE store_link_bicub 01112 01113 01114 END MODULE remap_bicubic_reduced 01115 01116 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01117
1.7.4 
