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 3241 2011-06-23 09:55:11Z 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 3241 2011-06-23 09:55:11Z 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 #ifdef DEBUGX
00126   DOUBLE PRECISION                    :: dbl_rad2deg
00127 #endif
00128 
00129 ! ---------------------------------------------------------------------
00130 !
00131 #ifdef VERBOSE
00132   PRINT *, '| | | | | | | Enter PRISMTrs_Bilinear_weight_2d'
00133   PRINT*, 'Size of target epio :',id_tgt_size
00134   call psmile_flushstd
00135 #endif
00136   id_err = 0
00137   dda_weights = 0.0
00138   common_grid_range(1) = 0.
00139   common_grid_range(2) =  dble_pi2
00140   shiftSize = common_grid_range(2) - common_grid_range(1)
00141 !
00142 ! 0. Treat the lat and lon arrays
00143 !
00144 ! 0.1. convert longitudes to 0,2pi interval
00145 !
00146   DO ib = 1, id_src_size
00147     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00148       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00149     ENDDO
00150     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00151       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00152     ENDDO
00153   ENDDO
00154   DO ib = 1, id_tgt_size
00155     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00156       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00157     ENDDO
00158     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00159       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00160     ENDDO
00161   ENDDO
00162 !
00163 ! 0.2. make sure input latitude range is within the machine values
00164 !      for +/- pi/2
00165 !
00166   DO ib = 1, id_src_size
00167     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00168         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00169         CALL psmile_abort
00170     ENDIF
00171   ENDDO
00172   DO ib = 1, id_tgt_size
00173     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00174         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00175         CALL psmile_abort
00176     ENDIF
00177   ENDDO
00178 
00179   nearest_neighbor = .false.
00180 
00181 #ifndef SX
00182 !
00183 ! 1. In the loop
00184 !
00185   DO ib = 1, id_tgt_size
00186 
00187     ! 1.1. Set the lat and lon of the TARGET point
00188     plat = dda_tgt_lat(ib)
00189     plon = dda_tgt_lon(ib)
00190 
00191     IF ( ida_neighbor_index(ib, 1) > 0 ) THEN
00192 
00193        iter = 0
00194 
00195 ! 1.2. For the TARGET point ib, check that the 4 neighours were found
00196 
00197        IF ( ida_neighbor_index(ib, 4) <= 0 ) nearest_neighbor = .true.
00198 
00199        IF ( .not. nearest_neighbor ) THEN
00200 
00201 ! 1.2.1. For the TARGET point ib, set the lat and lon of the four neighours
00202 
00203           DO ib_bis = 1, 4
00204              src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00205              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00206           END DO
00207 
00208           DO ib_bis = 1, 4
00209             vec_lon(ib_bis) = src_lons (ib_bis) - plon
00210             IF (vec_lon(ib_bis) > dble_pi) THEN
00211                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00212              ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00213                  src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00214              ENDIF
00215           ENDDO
00216 
00217 #ifdef DEBUGX
00218             dbl_rad2deg = 360.0/dble_pi2
00219             IF ( (ib == 4388) .OR. (ib == 4389) ) THEN
00220                 PRINT* 
00221                 PRINT*, '+++++++++++++++++++++++++++++++'
00222                 PRINT*, 'EPIOT number :',ib
00223                 PRINT*, '+++++++++++++++++++++++++++++++'
00224                 PRINT*
00225                 PRINT*, 'lat and lon of TARGET point :',dbl_rad2deg*dda_tgt_lat(ib),dbl_rad2deg*dda_tgt_lon(ib)
00226                 DO ib_bis=1,4
00227                   PRINT*, 'EPIOS number, lat and lon of neighbours:',ida_neighbor_index(ib, ib_bis),&
00228                      dbl_rad2deg*dda_src_lat (ida_neighbor_index(ib, ib_bis)), &
00229                      dbl_rad2deg*dda_src_lon (ida_neighbor_index(ib, ib_bis))
00230                 ENDDO
00231                 PRINT*
00232             ENDIF
00233 #endif
00234 
00235 ! 1.2.2. Compute the weights for a bilinear method
00236 
00237           dth1 = src_lats(2) - src_lats(1)
00238           dth2 = src_lats(4) - src_lats(1)
00239           dth3 = src_lats(3) - src_lats(2) - dth2
00240 
00241           dph1 = src_lons(2) - src_lons(1)
00242           dph2 = src_lons(4) - src_lons(1)
00243           dph3 = src_lons(3) - src_lons(2)
00244 
00245           IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00246           IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00247           IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
00248           IF (dph1 < -three*dble_pih) dph1 = dph1 + dble_pi2
00249           IF (dph2 < -three*dble_pih) dph2 = dph2 + dble_pi2
00250           IF (dph3 < -three*dble_pih) dph3 = dph3 + dble_pi2
00251 
00252           dph3 = dph3 - dph2
00253 
00254           iguess = half
00255           jguess = half
00256 
00257           ! Direct methods exist to calculate the weights. We need to
00258           ! investigate whether this is more performant. rr. 
00259           DO iter= 1, PSMILe_trans_Max_iter
00260 
00261              dthp = plat - src_lats(1)   &
00262                          - dth1*iguess   &
00263                          - dth2*jguess   &
00264                          - dth3*iguess*jguess
00265 
00266              dphp = plon - src_lons(1)
00267 
00268              IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00269              IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00270 
00271              dphp = dphp - dph1*iguess &
00272                          - dph2*jguess &
00273                          - dph3*iguess*jguess
00274 
00275              mat1 = dth1 + dth3*jguess
00276              mat2 = dth2 + dth3*iguess
00277              mat3 = dph1 + dph3*jguess
00278              mat4 = dph2 + dph3*iguess
00279 
00280              determinant = mat1*mat4 - mat2*mat3
00281 
00282              IF ( determinant == 0.0 ) EXIT
00283 
00284              deli = (dthp*mat4 - mat2*dphp)/determinant
00285              delj = (mat1*dphp - dthp*mat3)/determinant
00286 
00287              IF ( ABS(deli) < converge .AND. &
00288                   ABS(delj) < converge) EXIT
00289 
00290              iguess = iguess + deli
00291              jguess = jguess + delj
00292 
00293           END DO
00294 
00295           IF (iter <= PSMILe_trans_Max_iter) THEN
00296 
00297              !***
00298              !*** successfully found i,j - compute weights
00299              !***
00300 
00301              dda_weights(ib,1) = ABS((one-iguess)*(one-jguess))
00302              dda_weights(ib,2) = ABS(iguess*(one-jguess))
00303              dda_weights(ib,3) = ABS(iguess*jguess)
00304              dda_weights(ib,4) = ABS((one-iguess)*jguess)
00305 
00306              IF ( ABS(SUM(dda_weights(ib,1:4))-1.0D0) > tolerance ) THEN
00307                ! Something went wrong with the weights,
00308                ! reset and try again with nearest neighbor.
00309 #ifdef VERBOSE
00310                PRINT *, '| | | | | | | | We switch to nearest neighbor.',ABS(SUM(dda_weights(ib,1:4)))
00311 #endif
00312                dda_weights(ib,:) = 0.0
00313                nearest_neighbor = .true.
00314              ENDIF
00315 
00316           ELSE
00317 
00318 #ifdef VERBOSE
00319              PRINT *, '| | | | | | | | Point coords  : ',plat,plon
00320              PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00321              PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00322              PRINT *, '| | | | | | | | Current i,j   : ',iguess, jguess
00323              PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00324              PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00325 #endif
00326              nearest_neighbor = .true.
00327 
00328           ENDIF
00329 
00330        ENDIF ! switch to nearest neighbor
00331 
00332 ! We switch to nearest neighbor in case the iteration did not converge or produced
00333 ! wrong weights or in case the number of neighbors found was not sufficient to
00334 ! perform bilinear interpolation.
00335 
00336        IF ( nearest_neighbor ) THEN
00337 
00338           dist2v = 0.0d0
00339 
00340           DO ib_bis = 1, 4
00341              IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00342           ENDDO
00343 
00344           nb = ib_bis - 1
00345 
00346           DO ib_bis = 1, nb
00347 
00348               srch_add = ida_neighbor_index(ib, ib_bis)
00349 
00350               arg = SIN(plat)*SIN(dda_src_lat(srch_add)) +  &
00351                     COS(plat)*COS(dda_src_lat(srch_add)) *  &
00352                    (COS(plon)*COS(dda_src_lon(srch_add)) +  &
00353                     SIN(plon)*SIN(dda_src_lon(srch_add)))
00354 
00355               IF (arg < -1.0d0) THEN
00356                  arg = -1.0d0
00357               ELSE IF (arg > 1.0d0) THEN
00358                  arg = 1.0d0
00359               END IF
00360 
00361               dist2v(ib_bis) = ACOS(arg)
00362               dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00363 
00364           END DO
00365 
00366           DO ib_bis = 1, nb
00367              IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00368           END DO
00369 
00370           IF (ib_bis <= nb) THEN
00371 !             ... Point to be interpolated is located on Point i
00372              dda_weights(ib,ib_bis) = 1.0d0
00373 
00374           ELSE IF ( nb == 1 ) THEN
00375              dda_weights(ib,nb) = 1.0d0
00376 
00377           ELSE IF ( nb > 1 ) THEN
00378              dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00379              rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00380              dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00381           ENDIF
00382 
00383        END IF
00384 
00385        nearest_neighbor = .false.
00386 
00387      END IF
00388 
00389 #ifdef DEBUGX
00390      IF ( (ib == 4388) .OR. (ib == 4389) ) THEN
00391          PRINT*
00392          PRINT*, 'EPIOT number, EPIOS numbers of neighbours:',ib,ida_neighbor_index(ib,:)
00393          PRINT*, 'dda_weights(ib,:) :',dda_weights(ib,:)
00394          PRINT*
00395      ENDIF
00396 #endif
00397 
00398    END DO
00399 
00400 #else
00401 
00402 ! Vector length is close to maximum. Probably
00403 ! there is not much more to optimize.
00404 !
00405 ! 1.1. Set the lat and lon of the TARGET point
00406 !
00407   plat = dda_tgt_lat
00408   plon = dda_tgt_lon
00409 !
00410 ! 1. In the loop
00411 !
00412   DO ib = 1, id_tgt_size
00413      if ( ida_neighbor_index(ib, 4) <= 0  .AND. &
00414           ida_neighbor_index(ib, 1)  > 0  .AND. &
00415           ida_tgt_mask(ib) == 1 ) nearest_neighbor(ib) = .TRUE.
00416   END DO
00417 
00418   DO ib_part = 1, id_tgt_size, 5000
00419   DO ib = ib_part, min(id_tgt_size,ib_part+5000-1)
00420 
00421 ! 1.2. For the TARGET point ib, check that the 4 neighours were found
00422 
00423      IF ( .not. nearest_neighbor(ib) .AND. &
00424           ida_tgt_mask(ib) == 1      .AND. &
00425           ida_neighbor_index(ib, 1) > 0 ) THEN
00426 
00427 ! 1.2.1. For the TARGET point ib, set the lat and lon of the four neighours
00428 
00429           DO ib_bis = 1, 4
00430              src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00431              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00432           END DO
00433 
00434           DO ib_bis = 1, 4
00435             vec_lon(ib_bis) = src_lons (ib_bis) - plon(ib)
00436             IF (vec_lon(ib_bis) > dble_pi) THEN
00437                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00438             ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00439                 src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00440             ENDIF
00441           ENDDO
00442 
00443 ! 1.2.2. Compute the weights for a bilinear method
00444 
00445           dth1 = src_lats(2) - src_lats(1)
00446           dth2 = src_lats(4) - src_lats(1)
00447           dth3 = src_lats(3) - src_lats(2) - dth2
00448 
00449           dph1 = src_lons(2) - src_lons(1)
00450           dph2 = src_lons(4) - src_lons(1)
00451           dph3 = src_lons(3) - src_lons(2)
00452 
00453           IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00454           IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00455           IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
00456           IF (dph1 < -three*dble_pih) dph1 = dph1 + dble_pi2
00457           IF (dph2 < -three*dble_pih) dph2 = dph2 + dble_pi2
00458           IF (dph3 < -three*dble_pih) dph3 = dph3 + dble_pi2
00459 
00460           dph3 = dph3 - dph2
00461 
00462           iguess = half
00463           jguess = half
00464 
00465           ! Direct methods exist to calculate the weights. We need to
00466           ! investigate whether this is more performant. rr. 
00467           DO iter= 1, PSMILe_trans_Max_iter
00468 
00469              dthp = plat(ib) - src_lats(1)   &
00470                              - dth1*iguess   &
00471                              - dth2*jguess   &
00472                              - dth3*iguess*jguess
00473 
00474              dphp = plon(ib) - src_lons(1)
00475 
00476              IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00477              IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00478 
00479              dphp = dphp - dph1*iguess &
00480                          - dph2*jguess &
00481                          - dph3*iguess*jguess
00482 
00483              mat1 = dth1 + dth3*jguess
00484              mat2 = dth2 + dth3*iguess
00485              mat3 = dph1 + dph3*jguess
00486              mat4 = dph2 + dph3*iguess
00487 
00488              determinant = mat1*mat4 - mat2*mat3
00489 
00490              IF ( determinant == 0.0 ) EXIT
00491 
00492              deli = (dthp*mat4 - mat2*dphp)/determinant
00493              delj = (mat1*dphp - dthp*mat3)/determinant
00494 
00495              IF ( ABS(deli) < converge .AND. &
00496                   ABS(delj) < converge) EXIT
00497 
00498              iguess = iguess + deli
00499              jguess = jguess + delj
00500 
00501           END DO
00502 
00503           IF (iter <= PSMILe_trans_Max_iter) THEN
00504 
00505              !***
00506              !*** successfully found i,j - compute weights
00507              !***
00508 
00509              dda_weights(ib,1) = ABS((one-iguess)*(one-jguess))
00510              dda_weights(ib,2) = ABS(iguess*(one-jguess))
00511              dda_weights(ib,3) = ABS(iguess*jguess)
00512              dda_weights(ib,4) = ABS((one-iguess)*jguess)
00513 
00514              IF ( ABS(SUM(dda_weights(ib,1:4))-1.0D0) > tolerance ) THEN
00515                ! Something went wrong with the weights,
00516                ! reset and try again with nearest neighbor.
00517 #ifdef VERBOSE
00518                PRINT *, '| | | | | | | | We switch to nearest neighbor.',ABS(SUM(dda_weights(ib,1:4)))
00519 #endif
00520                dda_weights(ib,:) = 0.0
00521                nearest_neighbor(ib) = .true.
00522              ENDIF
00523 
00524           ELSE
00525 
00526 #ifdef VERBOSE
00527              PRINT *, '| | | | | | | | Point coords  : ',plat(ib),plon(ib)
00528              PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00529              PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00530              PRINT *, '| | | | | | | | Current i,j   : ',iguess, jguess
00531              PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00532              PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00533 #endif
00534              nearest_neighbor(ib) = .true.
00535 
00536           ENDIF
00537 
00538        ENDIF ! switch to nearest neighbor
00539 
00540   END DO
00541   END DO
00542 
00543 ! We switch to nearest neighbor in case the iteration did not converge or produced
00544 ! wrong weights or in case the number of neighbors found was not sufficient to
00545 ! perform bilinear interpolation.
00546 
00547   DO ib = 1, id_tgt_size
00548 
00549        IF ( nearest_neighbor(ib) ) THEN
00550 
00551           dist2v = 0.0d0
00552 
00553           DO ib_bis = 1, 4
00554              IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00555           ENDDO
00556 
00557           nb = ib_bis - 1
00558 
00559           DO ib_bis = 1, nb
00560 
00561               srch_add = ida_neighbor_index(ib, ib_bis)
00562 
00563               arg = SIN(plat(ib))*SIN(dda_src_lat(srch_add)) +  &
00564                     COS(plat(ib))*COS(dda_src_lat(srch_add)) *  &
00565                    (COS(plon(ib))*COS(dda_src_lon(srch_add)) +  &
00566                     SIN(plon(ib))*SIN(dda_src_lon(srch_add)))
00567 
00568               IF (arg < -1.0d0) THEN
00569                  arg = -1.0d0
00570               ELSE IF (arg > 1.0d0) THEN
00571                  arg = 1.0d0
00572               END IF
00573 
00574               dist2v(ib_bis) = ACOS(arg)
00575               dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00576 
00577           END DO
00578 
00579           DO ib_bis = 1, nb
00580              IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00581           END DO
00582 
00583           IF (ib_bis <= nb) THEN
00584 !             ... Point to be interpolated is located on Point i
00585              dda_weights(ib,ib_bis) = 1.0d0
00586 
00587           ELSE IF ( nb == 1 ) THEN
00588              dda_weights(ib,nb) = 1.0d0
00589 
00590           ELSE IF ( nb > 0 ) THEN
00591              dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00592              rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00593              dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00594           ENDIF
00595 
00596        END IF
00597 
00598 !       nearest_neighbor = .false.
00599 
00600 
00601   END DO
00602 
00603 #endif
00604 !
00605 #ifdef VERBOSE
00606   PRINT *, '| | | | | | | Quit PRISMTrs_Bilinear_weight_2d'
00607   call psmile_flushstd
00608 #endif
00609 
00610 END SUBROUTINE PRISMTrs_Bilinear_weight_2d
00611 
00612 
00613 
00614 
00615