prismtrs_bicubic_grad_2d.F90

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