prismtrs_bilinear_weight_2d.F90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------
00002 ! Copyright 2006-2010, CERFACS, Toulouse, France.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PRISMTrs_Bilinear_weight_2d
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_bilinear_weight_2d(id_src_size,         &
00011                                        dda_src_lat,         &
00012                                        dda_src_lon,         &
00013                                        id_tgt_size,         &
00014                                        dda_tgt_lat,         &
00015                                        dda_tgt_lon,         &
00016                                        ida_tgt_mask,        &
00017                                        ida_neighbor_index,  & 
00018                                        dda_weights,         &
00019                                        id_err)
00020 
00021 !
00022 ! !USES:
00023 !
00024   USE PRISMDrv, dummy_INTERFACE => PRISMTrs_Bilinear_weight_2d
00025 
00026   IMPLICIT NONE
00027 
00028 !
00029 ! !PARAMETERS:
00030 !
00031 ! size of the source and TARGET grid to interpolate
00032   INTEGER, INTENT (IN)                         :: id_src_size
00033   INTEGER, INTENT (IN)                         :: id_tgt_size
00034  
00035 ! latitudes and longitudes of the source grid
00036   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lat
00037   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lon
00038 
00039 ! latitudes and longitudes of the TARGET grid
00040   DOUBLE PRECISION, DIMENSION(id_tgt_size)     :: dda_tgt_lat
00041   DOUBLE PRECISION, DIMENSION(id_tgt_size)     :: dda_tgt_lon
00042 
00043   INTEGER, DIMENSION(id_tgt_size), INTENT (IN) :: ida_tgt_mask
00044 
00045 ! indices of the foud neighbors on the source grid for each TARGET point
00046   INTEGER, DIMENSION(id_tgt_size,4), INTENT (IN) :: ida_neighbor_index
00047 
00048 !
00049 ! ! RETURN VALUE
00050 !
00051 ! bilinear weights for four corners
00052   DOUBLE PRECISION, DIMENSION(id_tgt_size,4), INTENT (Out) :: dda_weights
00053 
00054   INTEGER, INTENT (Out)                        :: id_err   ! error value
00055 
00056 ! !DESCRIPTION
00057 ! SUBROUTINE "PRISMTrs_Bilinear_weight_2d" computes the weights with 
00058 ! the bilinear interpolation method, derived from the SCRIP 1.4
00059 ! library available from Los Alamos National Laboratory.
00060 ! The inputs are the neighboring adresses computed in the PSMILe.
00061 !
00062 ! !REVISED HISTORY
00063 !   Date      Programmer   Description
00064 ! ----------  ----------   -----------
00065 ! 06/12/2002  D. Declat    Creation
00066 ! 22/02/2005  R. Redler    Switch to nearest neighbor interpolation corrected.
00067 !
00068 ! EOP
00069 !----------------------------------------------------------------------
00070 ! $Id: prismtrs_bilinear_weight_2d.F90 2963 2011-02-17 14:52:53Z coquart $
00071 ! $Author: coquart $
00072 !----------------------------------------------------------------------
00073 !
00074 ! Local declarations
00075 !
00076   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00077      '$Id: prismtrs_bilinear_weight_2d.F90 2963 2011-02-17 14:52:53Z coquart $'
00078 
00079 ! loop indices
00080   INTEGER                             :: ib, ib_bis, iter
00081   DOUBLE PRECISION, PARAMETER         :: converge  = 1.0D-2
00082   DOUBLE PRECISION, PARAMETER         :: tolerance = 1.0D-12
00083 
00084 ! vectors of lon and lat of the source points used to compute the weights
00085   DOUBLE PRECISION, DIMENSION(4)      :: src_lats
00086   DOUBLE PRECISION, DIMENSION(4)      :: src_lons
00087   DOUBLE PRECISION, DIMENSION(4)      :: vec_lon
00088   DOUBLE PRECISION                    :: common_grid_range(2)
00089   DOUBLE PRECISION                    :: shiftSize
00090 
00091 ! lat/lon coords of destination point
00092 #ifndef SX
00093   DOUBLE PRECISION                    :: plat, plon
00094 #else
00095   DOUBLE PRECISION, DIMENSION(id_tgt_size) :: plat, plon
00096 #endif
00097 ! current guess for bilinear coordinate
00098   DOUBLE PRECISION                    :: iguess, jguess
00099 ! corrections to i,j
00100   DOUBLE PRECISION                    :: deli, delj
00101 ! some latitude  differences
00102   DOUBLE PRECISION                    :: dth1, dth2, dth3
00103 ! some longitude differences
00104   DOUBLE PRECISION                    :: dph1, dph2, dph3
00105  ! difference between point and sw corner
00106   DOUBLE PRECISION                    :: dthp, dphp
00107 ! matrix elements
00108   DOUBLE PRECISION                    :: mat1, mat2, mat3, mat4
00109 ! matrix determinant
00110   DOUBLE PRECISION                    :: determinant
00111 
00112 ! for computing dist-weighted avg
00113   DOUBLE PRECISION                    :: dist2v(4)
00114   DOUBLE PRECISION                    :: rdistance, arg
00115 !
00116   INTEGER                             :: nb
00117   INTEGER                             :: srch_add
00118 #ifndef SX
00119   LOGICAL                             :: nearest_neighbor
00120 #else
00121   INTEGER                             :: ib_part
00122   LOGICAL, DIMENSION(id_tgt_size)     :: nearest_neighbor
00123 #endif
00124 
00125 ! ---------------------------------------------------------------------
00126 !
00127 #ifdef VERBOSE
00128   PRINT *, '| | | | | | | Enter PRISMTrs_Bilinear_weight_2d'
00129   call psmile_flushstd
00130 #endif
00131   id_err = 0
00132   dda_weights = 0.0
00133   common_grid_range(1) = -dble_pi
00134   common_grid_range(2) =  dble_pi
00135   shiftSize = common_grid_range(2) - common_grid_range(1)
00136 !
00137 ! 0. Treat the lat and lon arrays
00138 !
00139 ! 0.1. convert longitudes to 0,2pi interval
00140 !
00141   DO ib = 1, id_src_size
00142     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00143       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00144     ENDDO
00145     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00146       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00147     ENDDO
00148   ENDDO
00149   DO ib = 1, id_tgt_size
00150     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00151       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00152     ENDDO
00153     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00154       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00155     ENDDO
00156   ENDDO
00157 !
00158 ! 0.2. make sure input latitude range is within the machine values
00159 !      for +/- pi/2
00160 !
00161   DO ib = 1, id_src_size
00162     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00163         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00164         CALL psmile_abort
00165     ENDIF
00166   ENDDO
00167   DO ib = 1, id_tgt_size
00168     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00169         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00170         CALL psmile_abort
00171     ENDIF
00172   ENDDO
00173 
00174   nearest_neighbor = .false.
00175 
00176 #ifndef SX
00177 !
00178 ! 1. In the loop
00179 !
00180   DO ib = 1, id_tgt_size
00181 
00182     ! 1.1. Set the lat and lon of the TARGET point
00183     plat = dda_tgt_lat(ib)
00184     plon = dda_tgt_lon(ib)
00185 
00186     IF ( ida_neighbor_index(ib, 1) > 0 ) THEN
00187 
00188        iter = 0
00189 
00190 ! 1.2. For the TARGET point ib, check that the 4 neighours were found
00191 
00192        IF ( ida_neighbor_index(ib, 4) <= 0 ) nearest_neighbor = .true.
00193 
00194        IF ( .not. nearest_neighbor ) THEN
00195 
00196 ! 1.2.1. For the TARGET point ib, set the lat and lon of the four neighours
00197 
00198           DO ib_bis = 1, 4
00199              src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00200              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00201           END DO
00202 
00203           DO ib_bis = 1, 4
00204             vec_lon(ib_bis) = src_lons (ib_bis) - plon
00205             IF (vec_lon(ib_bis) > dble_pi) THEN
00206                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00207              ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00208                  src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00209              ENDIF
00210           ENDDO
00211 
00212 ! 1.2.2. Compute the weights for a bilinear method
00213 
00214           dth1 = src_lats(2) - src_lats(1)
00215           dth2 = src_lats(4) - src_lats(1)
00216           dth3 = src_lats(3) - src_lats(2) - dth2
00217 
00218           dph1 = src_lons(2) - src_lons(1)
00219           dph2 = src_lons(4) - src_lons(1)
00220           dph3 = src_lons(3) - src_lons(2)
00221 
00222           IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00223           IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00224           IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
00225           IF (dph1 < -three*dble_pih) dph1 = dph1 + dble_pi2
00226           IF (dph2 < -three*dble_pih) dph2 = dph2 + dble_pi2
00227           IF (dph3 < -three*dble_pih) dph3 = dph3 + dble_pi2
00228 
00229           dph3 = dph3 - dph2
00230 
00231           iguess = half
00232           jguess = half
00233 
00234           ! Direct methods exist to calculate the weights. We need to
00235           ! investigate whether this is more performant. rr. 
00236           DO iter= 1, PSMILe_trans_Max_iter
00237 
00238              dthp = plat - src_lats(1)   &
00239                          - dth1*iguess   &
00240                          - dth2*jguess   &
00241                          - dth3*iguess*jguess
00242 
00243              dphp = plon - src_lons(1)
00244 
00245              IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00246              IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00247 
00248              dphp = dphp - dph1*iguess &
00249                          - dph2*jguess &
00250                          - dph3*iguess*jguess
00251 
00252              mat1 = dth1 + dth3*jguess
00253              mat2 = dth2 + dth3*iguess
00254              mat3 = dph1 + dph3*jguess
00255              mat4 = dph2 + dph3*iguess
00256 
00257              determinant = mat1*mat4 - mat2*mat3
00258 
00259              IF ( determinant == 0.0 ) EXIT
00260 
00261              deli = (dthp*mat4 - mat2*dphp)/determinant
00262              delj = (mat1*dphp - dthp*mat3)/determinant
00263 
00264              IF ( ABS(deli) < converge .AND. &
00265                   ABS(delj) < converge) EXIT
00266 
00267              iguess = iguess + deli
00268              jguess = jguess + delj
00269 
00270           END DO
00271 
00272           IF (iter <= PSMILe_trans_Max_iter) THEN
00273 
00274              !***
00275              !*** successfully found i,j - compute weights
00276              !***
00277 
00278              dda_weights(ib,1) = ABS((one-iguess)*(one-jguess))
00279              dda_weights(ib,2) = ABS(iguess*(one-jguess))
00280              dda_weights(ib,3) = ABS(iguess*jguess)
00281              dda_weights(ib,4) = ABS((one-iguess)*jguess)
00282 
00283              IF ( ABS(SUM(dda_weights(ib,1:4))-1.0D0) > tolerance ) THEN
00284                ! Something went wrong with the weights,
00285                ! reset and try again with nearest neighbor.
00286 #ifdef VERBOSE
00287                PRINT *, '| | | | | | | | We switch to nearest neighbor.',ABS(SUM(dda_weights(ib,1:4)))
00288 #endif
00289                dda_weights(ib,:) = 0.0
00290                nearest_neighbor = .true.
00291              ENDIF
00292 
00293           ELSE
00294 
00295 #ifdef VERBOSE
00296              PRINT *, '| | | | | | | | Point coords  : ',plat,plon
00297              PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00298              PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00299              PRINT *, '| | | | | | | | Current i,j   : ',iguess, jguess
00300              PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00301              PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00302 #endif
00303              nearest_neighbor = .true.
00304 
00305           ENDIF
00306 
00307        ENDIF ! switch to nearest neighbor
00308 
00309 ! We switch to nearest neighbor in case the iteration did not converge or produced
00310 ! wrong weights or in case the number of neighbors found was not sufficient to
00311 ! perform bilinear interpolation.
00312 
00313        IF ( nearest_neighbor ) THEN
00314 
00315           dist2v = 0.0d0
00316 
00317           DO ib_bis = 1, 4
00318              IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00319           ENDDO
00320 
00321           nb = ib_bis - 1
00322 
00323           DO ib_bis = 1, nb
00324 
00325               srch_add = ida_neighbor_index(ib, ib_bis)
00326 
00327               arg = SIN(plat)*SIN(dda_src_lat(srch_add)) +  &
00328                     COS(plat)*COS(dda_src_lat(srch_add)) *  &
00329                    (COS(plon)*COS(dda_src_lon(srch_add)) +  &
00330                     SIN(plon)*SIN(dda_src_lon(srch_add)))
00331 
00332               IF (arg < -1.0d0) THEN
00333                  arg = -1.0d0
00334               ELSE IF (arg > 1.0d0) THEN
00335                  arg = 1.0d0
00336               END IF
00337 
00338               dist2v(ib_bis) = ACOS(arg)
00339               dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00340 
00341           END DO
00342 
00343           DO ib_bis = 1, nb
00344              IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00345           END DO
00346 
00347           IF (ib_bis <= nb) THEN
00348 !             ... Point to be interpolated is located on Point i
00349              dda_weights(ib,ib_bis) = 1.0d0
00350 
00351           ELSE IF ( nb == 1 ) THEN
00352              dda_weights(ib,nb) = 1.0d0
00353 
00354           ELSE IF ( nb > 1 ) THEN
00355              dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00356              rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00357              dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00358           ENDIF
00359 
00360        END IF
00361 
00362        nearest_neighbor = .false.
00363 
00364      END IF
00365 
00366   END DO
00367 
00368 #else
00369 
00370 ! Vector length is close to maximum. Probably
00371 ! there is not much more to optimize.
00372 !
00373 ! 1.1. Set the lat and lon of the TARGET point
00374 !
00375   plat = dda_tgt_lat
00376   plon = dda_tgt_lon
00377 !
00378 ! 1. In the loop
00379 !
00380   DO ib = 1, id_tgt_size
00381      if ( ida_neighbor_index(ib, 4) <= 0  .AND. &
00382           ida_neighbor_index(ib, 1)  > 0  .AND. &
00383           ida_tgt_mask(ib) == 1 ) nearest_neighbor(ib) = .TRUE.
00384   END DO
00385 
00386   DO ib_part = 1, id_tgt_size, 5000
00387   DO ib = ib_part, min(id_tgt_size,ib_part+5000-1)
00388 
00389 ! 1.2. For the TARGET point ib, check that the 4 neighours were found
00390 
00391      IF ( .not. nearest_neighbor(ib) .AND. &
00392           ida_tgt_mask(ib) == 1      .AND. &
00393           ida_neighbor_index(ib, 1) > 0 ) THEN
00394 
00395 ! 1.2.1. For the TARGET point ib, set the lat and lon of the four neighours
00396 
00397           DO ib_bis = 1, 4
00398              src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00399              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00400           END DO
00401 
00402           DO ib_bis = 1, 4
00403             vec_lon(ib_bis) = src_lons (ib_bis) - plon(ib)
00404             IF (vec_lon(ib_bis) > dble_pi) THEN
00405                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00406             ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00407                 src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00408             ENDIF
00409           ENDDO
00410 
00411 ! 1.2.2. Compute the weights for a bilinear method
00412 
00413           dth1 = src_lats(2) - src_lats(1)
00414           dth2 = src_lats(4) - src_lats(1)
00415           dth3 = src_lats(3) - src_lats(2) - dth2
00416 
00417           dph1 = src_lons(2) - src_lons(1)
00418           dph2 = src_lons(4) - src_lons(1)
00419           dph3 = src_lons(3) - src_lons(2)
00420 
00421           IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00422           IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00423           IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
00424           IF (dph1 < -three*dble_pih) dph1 = dph1 + dble_pi2
00425           IF (dph2 < -three*dble_pih) dph2 = dph2 + dble_pi2
00426           IF (dph3 < -three*dble_pih) dph3 = dph3 + dble_pi2
00427 
00428           dph3 = dph3 - dph2
00429 
00430           iguess = half
00431           jguess = half
00432 
00433           ! Direct methods exist to calculate the weights. We need to
00434           ! investigate whether this is more performant. rr. 
00435           DO iter= 1, PSMILe_trans_Max_iter
00436 
00437              dthp = plat(ib) - src_lats(1)   &
00438                              - dth1*iguess   &
00439                              - dth2*jguess   &
00440                              - dth3*iguess*jguess
00441 
00442              dphp = plon(ib) - src_lons(1)
00443 
00444              IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00445              IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00446 
00447              dphp = dphp - dph1*iguess &
00448                          - dph2*jguess &
00449                          - dph3*iguess*jguess
00450 
00451              mat1 = dth1 + dth3*jguess
00452              mat2 = dth2 + dth3*iguess
00453              mat3 = dph1 + dph3*jguess
00454              mat4 = dph2 + dph3*iguess
00455 
00456              determinant = mat1*mat4 - mat2*mat3
00457 
00458              IF ( determinant == 0.0 ) EXIT
00459 
00460              deli = (dthp*mat4 - mat2*dphp)/determinant
00461              delj = (mat1*dphp - dthp*mat3)/determinant
00462 
00463              IF ( ABS(deli) < converge .AND. &
00464                   ABS(delj) < converge) EXIT
00465 
00466              iguess = iguess + deli
00467              jguess = jguess + delj
00468 
00469           END DO
00470 
00471           IF (iter <= PSMILe_trans_Max_iter) THEN
00472 
00473              !***
00474              !*** successfully found i,j - compute weights
00475              !***
00476 
00477              dda_weights(ib,1) = ABS((one-iguess)*(one-jguess))
00478              dda_weights(ib,2) = ABS(iguess*(one-jguess))
00479              dda_weights(ib,3) = ABS(iguess*jguess)
00480              dda_weights(ib,4) = ABS((one-iguess)*jguess)
00481 
00482              IF ( ABS(SUM(dda_weights(ib,1:4))-1.0D0) > tolerance ) THEN
00483                ! Something went wrong with the weights,
00484                ! reset and try again with nearest neighbor.
00485 #ifdef VERBOSE
00486                PRINT *, '| | | | | | | | We switch to nearest neighbor.',ABS(SUM(dda_weights(ib,1:4)))
00487 #endif
00488                dda_weights(ib,:) = 0.0
00489                nearest_neighbor(ib) = .true.
00490              ENDIF
00491 
00492           ELSE
00493 
00494 #ifdef VERBOSE
00495              PRINT *, '| | | | | | | | Point coords  : ',plat(ib),plon(ib)
00496              PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00497              PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00498              PRINT *, '| | | | | | | | Current i,j   : ',iguess, jguess
00499              PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00500              PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00501 #endif
00502              nearest_neighbor(ib) = .true.
00503 
00504           ENDIF
00505 
00506        ENDIF ! switch to nearest neighbor
00507 
00508   END DO
00509   END DO
00510 
00511 ! We switch to nearest neighbor in case the iteration did not converge or produced
00512 ! wrong weights or in case the number of neighbors found was not sufficient to
00513 ! perform bilinear interpolation.
00514 
00515   DO ib = 1, id_tgt_size
00516 
00517        IF ( nearest_neighbor(ib) ) THEN
00518 
00519           dist2v = 0.0d0
00520 
00521           DO ib_bis = 1, 4
00522              IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00523           ENDDO
00524 
00525           nb = ib_bis - 1
00526 
00527           DO ib_bis = 1, nb
00528 
00529               srch_add = ida_neighbor_index(ib, ib_bis)
00530 
00531               arg = SIN(plat(ib))*SIN(dda_src_lat(srch_add)) +  &
00532                     COS(plat(ib))*COS(dda_src_lat(srch_add)) *  &
00533                    (COS(plon(ib))*COS(dda_src_lon(srch_add)) +  &
00534                     SIN(plon(ib))*SIN(dda_src_lon(srch_add)))
00535 
00536               IF (arg < -1.0d0) THEN
00537                  arg = -1.0d0
00538               ELSE IF (arg > 1.0d0) THEN
00539                  arg = 1.0d0
00540               END IF
00541 
00542               dist2v(ib_bis) = ACOS(arg)
00543               dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00544 
00545           END DO
00546 
00547           DO ib_bis = 1, nb
00548              IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00549           END DO
00550 
00551           IF (ib_bis <= nb) THEN
00552 !             ... Point to be interpolated is located on Point i
00553              dda_weights(ib,ib_bis) = 1.0d0
00554 
00555           ELSE IF ( nb == 1 ) THEN
00556              dda_weights(ib,nb) = 1.0d0
00557 
00558           ELSE IF ( nb > 0 ) THEN
00559              dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00560              rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00561              dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00562           ENDIF
00563 
00564        END IF
00565 
00566 !       nearest_neighbor = .false.
00567 
00568 
00569   END DO
00570 
00571 #endif
00572 !
00573 #ifdef VERBOSE
00574   PRINT *, '| | | | | | | Quit PRISMTrs_Bilinear_weight_2d'
00575   call psmile_flushstd
00576 #endif
00577 
00578 END SUBROUTINE PRISMTrs_Bilinear_weight_2d
00579 
00580 
00581 
00582 
00583