prismtrs_trilinear_weight.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 !
00008 ! !ROUTINE: PRISMTrs_Trilinear_weight
00009 !
00010 ! !INTERFACE
00011 subroutine prismtrs_trilinear_weight(id_src_size,         &
00012                                      dda_src_lat,         &
00013                                      dda_src_lon,         &
00014                                      dda_src_z,           &
00015                                      id_tgt_size,         &
00016                                      dda_tgt_lat,         &
00017                                      dda_tgt_lon,         &
00018                                      dda_tgt_z,           &
00019                                      ida_tgt_mask,        &
00020                                      ida_neighbor_index,  & 
00021                                      dda_weights,         &
00022                                      id_err)
00023 !
00024 ! !USES:
00025 !
00026   USE PRISMDrv, dummy_interface => PRISMTrs_Trilinear_weight
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, longitudes and z 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   DOUBLE PRECISION, DIMENSION(id_src_size)             :: dda_src_z
00041 
00042 ! latitudes, longitudes and z of the target grid
00043   DOUBLE PRECISION, DIMENSION(id_tgt_size)             :: dda_tgt_lat
00044   DOUBLE PRECISION, DIMENSION(id_tgt_size)             :: dda_tgt_lon
00045   DOUBLE PRECISION, DIMENSION(id_tgt_size)             :: dda_tgt_z
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,8), INTENT(IN) :: ida_neighbor_index
00051 
00052 !
00053 ! ! RETURN VALUE
00054 !
00055 ! trilinear weights for 8 corners
00056   DOUBLE PRECISION, DIMENSION(id_tgt_size,8), INTENT(Out) :: dda_weights
00057 
00058   INTEGER, INTENT (Out)                      :: id_err
00059 
00060 ! !DESCRIPTION
00061 ! Subroutine "PRISMTrs_Trilinear_weight" computes the weights for the 
00062 ! trilinear interpolation method. 
00063 ! The inputs are the results of the search neighboring done in the PSMILe.
00064 ! The weights are stored in files which reference is given to the 
00065 ! transformer as output of the routine.
00066 !
00067 ! !REVISED HISTORY
00068 !   Date      Programmer   Description
00069 ! ----------  ----------   -----------
00070 ! 25/07/2003  C. Cylly     Creation
00071 ! 23/072005   R. Redler    Corrected the switch to nearest neighbor.
00072 !
00073 ! EOP
00074 !----------------------------------------------------------------------
00075 ! $Id: prismtrs_trilinear_weight.F90 2963 2011-02-17 14:52:53Z coquart $
00076 ! $Author: coquart $
00077 !----------------------------------------------------------------------
00078 !
00079 ! Local declarations
00080 !
00081   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00082      '$Id: prismtrs_trilinear_weight.F90 2963 2011-02-17 14:52:53Z coquart $'
00083 
00084 ! Earth radius
00085   DOUBLE PRECISION, PARAMETER :: a = 6400.0
00086 
00087 ! loop indices
00088   INTEGER                     :: ib, ib_bis, iter
00089   DOUBLE PRECISION, PARAMETER :: converge = 1.0D-2
00090   DOUBLE PRECISION, PARAMETER :: tolerance = 1.0D-8
00091 
00092 ! vectors of lon and lat of the source points used to compute the weights
00093   DOUBLE PRECISION                  :: common_grid_range(2)
00094   DOUBLE PRECISION                  :: shiftSize
00095 #ifndef SX
00096   DOUBLE PRECISION, DIMENSION(8) :: src_lats
00097   DOUBLE PRECISION, DIMENSION(8) :: src_lons,vec_lon
00098   DOUBLE PRECISION, DIMENSION(8) :: src_zs
00099 ! lat/lon/z coords of destination point
00100   DOUBLE PRECISION               :: plat, plon, pz
00101 ! some latitude  differences
00102   DOUBLE PRECISION :: dth1, dth2, dth3, dth4, dth5, dth6, dth7
00103 ! some longitude differences
00104   DOUBLE PRECISION :: dph1, dph2, dph3, dph4, dph5, dph6, dph7
00105 ! some z  differences
00106   DOUBLE PRECISION :: dz1, dz2, dz3, dz4, dz5, dz6, dz7
00107 #else
00108   INTEGER                                  :: ib_part
00109 ! lon coords of the 8 source points
00110   DOUBLE PRECISION, DIMENSION(8)           :: src_lons,vec_lon
00111 ! lat/lon/z coords of destination point
00112   DOUBLE PRECISION, DIMENSION(id_tgt_size) :: plat, plon, pz
00113 ! some latitude  differences
00114   DOUBLE PRECISION, DIMENSION(id_tgt_size) :: dth1, dth2, dth3, dth4, dth5, dth6, dth7
00115 ! some longitude differences
00116   DOUBLE PRECISION, DIMENSION(id_tgt_size) :: dph1, dph2, dph3, dph4, dph5, dph6, dph7
00117 ! some z  differences
00118   DOUBLE PRECISION, DIMENSION(id_tgt_size) :: dz1, dz2, dz3, dz4, dz5, dz6, dz7
00119 #endif
00120 ! current guess for trilinear coordinate
00121   DOUBLE PRECISION :: iguess, jguess, kguess
00122 ! corrections to i,j,k
00123   DOUBLE PRECISION :: deli, delj, delk
00124   DOUBLE PRECISION :: dphtmp1, dphtmp2, dphtmp3
00125  ! difference between point and sw corner of lower box-face
00126   DOUBLE PRECISION :: dthp, dphp, dzp
00127 ! matrix elements
00128   DOUBLE PRECISION :: mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9 
00129 ! matrix determinant
00130   DOUBLE PRECISION :: determinant , deti, detj, detk     
00131 
00132 ! for computing dist-weighted avg
00133   DOUBLE PRECISION                    :: dist2v(8)
00134   DOUBLE PRECISION                    :: rdistance, distance, arg
00135 !
00136   INTEGER                             :: nb
00137   INTEGER                             :: srch_add
00138 #ifndef SX
00139   LOGICAL                             :: nearest_neighbor
00140 #else
00141   LOGICAL, DIMENSION(id_tgt_size)     :: nearest_neighbor
00142 #endif
00143 ! ---------------------------------------------------------------------
00144 !
00145 #ifdef VERBOSE
00146   PRINT *, '| | | | | | | Enter PRISMTrs_Trilinear_weight'
00147   call psmile_flushstd
00148 #endif
00149 
00150   id_err = 0
00151   dda_weights = 0.0
00152   common_grid_range(1) = -dble_pi
00153   common_grid_range(2) =  dble_pi
00154   shiftSize = common_grid_range(2) - common_grid_range(1)
00155 !
00156 ! 0. Treat the lat and lon arrays
00157 !
00158 ! 0.1. convert longitudes to 0,2pi interval
00159 !
00160   DO ib = 1, id_src_size
00161     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00162       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00163     ENDDO
00164     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00165       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00166     ENDDO
00167   ENDDO
00168   DO ib = 1, id_tgt_size
00169     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00170       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00171     ENDDO
00172     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00173       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00174     ENDDO
00175   ENDDO
00176 !
00177 ! 0.2. make sure input latitude range is within the machine values
00178 !      for +/- pi/2
00179 !  
00180   DO ib = 1, id_src_size
00181     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00182         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00183         CALL psmile_abort
00184     ENDIF
00185   ENDDO
00186   DO ib = 1, id_tgt_size
00187     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00188         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00189         CALL psmile_abort
00190     ENDIF
00191   ENDDO
00192 
00193   nearest_neighbor = .false.
00194 
00195 #ifndef SX
00196 !
00197 ! 1. In the loop
00198 !
00199   DO ib = 1, id_tgt_size
00200 
00201 #ifdef DEBUG
00202     PRINT *, 'Target index and neighbours = ', ib, ida_neighbor_index(ib, :)
00203 #endif
00204 
00205     IF ( ida_tgt_mask(ib) == 1 .AND. ida_neighbor_index(ib, 1) > 0 ) THEN
00206 
00207 ! 1.1. Set the lat and lon of the TARGET point
00208        plat = dda_tgt_lat(ib)
00209        plon = dda_tgt_lon(ib)
00210        pz   = dda_tgt_z(ib)
00211 
00212 ! 1.2. For the target point ib, check that the eight neighours were found
00213 
00214        IF ( ida_neighbor_index(ib, 8) <= 0 ) nearest_neighbor = .true.
00215 
00216        IF ( .not. nearest_neighbor ) THEN
00217 #ifdef DEBUG
00218            PRINT *, 'not nearest_neighbor, ib, lat, lon, z', ib, plat, plon, pz
00219 #endif
00220     
00221 ! 1.2.1. For the TARGET point ib, set the lat and lon of the eight neighours
00222 
00223           DO ib_bis = 1, 8
00224              src_zs (ib_bis)   = dda_src_z (ida_neighbor_index(ib, ib_bis))
00225              src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00226              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00227           END DO
00228 
00229           DO ib_bis = 1, 8
00230             vec_lon(ib_bis) = src_lons (ib_bis) - plon
00231             IF (vec_lon(ib_bis) > dble_pi) THEN
00232                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00233             ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00234                 src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00235             ENDIF
00236           ENDDO
00237 
00238 ! 1.2.2. Compute the weights for a trilinear method
00239 
00240           dth1 = src_lats(2) - src_lats(1)
00241           dth2 = src_lats(4) - src_lats(1)
00242           dth3 = src_lats(5) - src_lats(1)
00243 
00244           dth4 = src_lats(3) - src_lats(4) - dth1
00245           dth5 = src_lats(6) - src_lats(5) - dth1
00246           dth6 = src_lats(8) - src_lats(5) - dth2
00247 
00248           dth7 = src_lats(2) - src_lats(1) + &
00249                  src_lats(4) - src_lats(3) + &
00250                  src_lats(5) - src_lats(6) + &
00251                  src_lats(7) - src_lats(8)
00252 
00253           dph1 = src_lons(2) - src_lons(1)
00254           dph2 = src_lons(4) - src_lons(1)
00255           dph3 = src_lons(5) - src_lons(1)
00256 
00257           IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00258           IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00259           IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
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 
00264           dphtmp1 = src_lons(3) - src_lons(4)
00265           dphtmp2 = src_lons(6) - src_lons(5)
00266           dphtmp3 = src_lons(8) - src_lons(5)
00267 
00268           IF (dphtmp1 >  three*dble_pih) dphtmp1 = dphtmp1 - dble_pi2
00269           IF (dphtmp2 >  three*dble_pih) dphtmp2 = dphtmp2 - dble_pi2
00270           IF (dphtmp3 >  three*dble_pih) dphtmp3 = dphtmp3 - dble_pi2
00271           IF (dphtmp1 < -three*dble_pih) dphtmp1 = dphtmp1 + dble_pi2
00272           IF (dphtmp2 < -three*dble_pih) dphtmp2 = dphtmp2 + dble_pi2
00273           IF (dphtmp3 < -three*dble_pih) dphtmp3 = dphtmp3 + dble_pi2
00274 
00275           dph4 = dphtmp1 - dph1
00276           dph5 = dphtmp2 - dph1
00277           dph6 = dphtmp3 - dph2
00278 
00279           dphtmp3 = src_lons(7) - src_lons(8)
00280 
00281           IF (dphtmp3 >  three*dble_pih) dphtmp3 = dphtmp3 - dble_pi2
00282           IF (dphtmp3 < -three*dble_pih) dphtmp3 = dphtmp3 + dble_pi2
00283 
00284           dph7 = dph1 - dphtmp1 - dphtmp2 + dphtmp3
00285 
00286           dz1 = src_zs(2) - src_zs(1)
00287           dz2 = src_zs(4) - src_zs(1)
00288           dz3 = src_zs(5) - src_zs(1)
00289 
00290           dz4 = src_zs(3) - src_zs(4) - dz1
00291           dz5 = src_zs(6) - src_zs(5) - dz1
00292           dz6 = src_zs(8) - src_zs(5) - dz2
00293 
00294           dz7 = src_zs(2) - src_zs(1) + &
00295                 src_zs(4) - src_zs(3) + &
00296                 src_zs(5) - src_zs(6) + &
00297                 src_zs(7) - src_zs(8)
00298 
00299           iguess = half
00300           jguess = half
00301           kguess = half
00302        
00303           ! Direct methods exist to calculate the weights. We need to
00304           ! investigate whether this is more performant. rr. 
00305           DO iter = 1, PSMILe_trans_Max_iter
00306 
00307              dthp = plat - src_lats(1)        &
00308                          - dth1*iguess        &
00309                          - dth2*jguess        &
00310                          - dth3*kguess        &
00311                          - dth4*iguess*jguess &
00312                          - dth5*iguess*kguess &
00313                          - dth6*jguess*kguess &
00314                          - dth7*iguess*jguess*kguess
00315 
00316              dphp = plon - src_lons(1)
00317 
00318              IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00319              IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00320 
00321              dphp = dphp - dph1*iguess        &
00322                          - dph2*jguess        &
00323                          - dph3*kguess        &
00324                          - dph4*iguess*jguess &
00325                          - dph5*iguess*kguess &
00326                          - dph6*jguess*kguess &
00327                          - dph7*iguess*jguess*kguess
00328 
00329 
00330              dzp = pz - src_zs(1)         &
00331                       - dz1*iguess        &
00332                       - dz2*jguess        &
00333                       - dz3*kguess        &
00334                       - dz4*iguess*jguess &
00335                       - dz5*iguess*kguess &
00336                       - dz6*jguess*kguess &
00337                       - dz7*iguess*jguess*kguess
00338 
00339 
00340              mat1 = dth1 + dth4*jguess + dth5*kguess + dth7*jguess*kguess
00341              mat2 = dth2 + dth4*iguess + dth6*kguess + dth7*iguess*kguess
00342              mat3 = dth3 + dth5*iguess + dth6*jguess + dth7*iguess*jguess
00343 
00344              mat4 = dph1 + dph4*jguess + dph5*kguess + dph7*jguess*kguess
00345              mat5 = dph2 + dph4*iguess + dph6*kguess + dph7*iguess*kguess
00346              mat6 = dph3 + dph5*iguess + dph6*jguess + dph7*iguess*jguess
00347 
00348              mat7 = dz1 + dz4*jguess + dz5*kguess + dz7*jguess*kguess
00349              mat8 = dz2 + dz4*iguess + dz6*kguess + dz7*iguess*kguess
00350              mat9 = dz3 + dz5*iguess + dz6*jguess + dz7*iguess*jguess
00351 
00352 
00353              determinant = mat1*mat5*mat9 &
00354                          + mat2*mat6*mat7 &
00355                          + mat3*mat4*mat8 &
00356                          - mat7*mat5*mat3 &
00357                          - mat8*mat6*mat1 &
00358                          - mat9*mat4*mat2
00359 
00360              deti =  dthp*mat5*mat9 &
00361                   +  mat2*mat6*dzp  &
00362                   +  mat3*dphp*mat8 &
00363                   -  dzp*mat5*mat3  &
00364                   -  mat8*mat6*dthp &
00365                   -  mat9*dphp*mat2
00366 
00367              detj =  mat1*dphp*mat9 &
00368                   +  dthp*mat6*mat7 &
00369                   +  mat3*mat4*dzp  &
00370                   -  mat7*dphp*mat3 &
00371                   -  dzp*mat6*mat1  &
00372                   -  mat9*mat4*dthp
00373 
00374              detk =  mat1*mat5*dzp  &
00375                   +  mat2*dphp*mat7 &
00376                   +  dthp*mat4*mat8 &
00377                   -  mat7*mat5*dthp &
00378                   -  mat8*dphp*mat1 &
00379                   -  dzp*mat4*mat2
00380 
00381              IF ( determinant == 0.0 ) THEN
00382 #ifdef DEBUG
00383                  PRINT *, 'determinant = 0.0, EXIT', determinant
00384 #endif
00385                  EXIT 
00386              ENDIF
00387 
00388              deli = deti/determinant
00389              delj = detj/determinant
00390              delk = detk/determinant
00391 
00392              IF (ABS(deli) < converge .and. &
00393                  ABS(delj) < converge .and. &
00394                  ABS(delk) < converge) EXIT
00395 
00396              iguess = iguess + deli
00397              jguess = jguess + delj
00398              kguess = kguess + delk
00399 
00400           END DO
00401 
00402           IF (iter <= PSMILe_trans_Max_iter) THEN
00403 
00404              !***
00405              !*** successfully found i,j - compute weights
00406              !***
00407 
00408              IF (iguess.lt.0.) iguess = 0.0
00409              IF (jguess.lt.0.) jguess = 0.0
00410              IF (kguess.lt.0.) kguess = 0.0
00411 
00412              IF (iguess.gt.1.) iguess = 1.0
00413              IF (jguess.gt.1.) jguess = 1.0
00414              IF (kguess.gt.1.) kguess = 1.0
00415 
00416              dda_weights(ib,1) = abs((one-iguess)*(one-jguess)*(one-kguess))
00417              dda_weights(ib,2) = abs(iguess*(one-jguess)*(one-kguess))
00418              dda_weights(ib,3) = abs(iguess*jguess*(one-kguess))
00419              dda_weights(ib,4) = abs((one-iguess)*jguess*(one-kguess))
00420              dda_weights(ib,5) = abs((one-iguess)*(one-jguess)*kguess)
00421              dda_weights(ib,6) = abs(iguess*(one-jguess)*kguess)
00422              dda_weights(ib,7) = abs(iguess*jguess*kguess)
00423              dda_weights(ib,8) = abs((one-iguess)*jguess*kguess)
00424 
00425              IF ( ABS(SUM(dda_weights(ib,1:8))-1.0D0) > tolerance ) THEN
00426                ! Something went wrong with the weights,
00427                ! reset and try again with nearest neighbor.
00428 #ifdef VERBOSE
00429                PRINT *, '| | | | | | | | We switch to nearest neighbor.', ABS(SUM(dda_weights(ib,1:8))-1.0D0)
00430 #endif
00431                dda_weights(ib,:) = 0.0
00432                nearest_neighbor = .true.
00433              ENDIF
00434 
00435           ELSE
00436 
00437 #ifdef VERBOSE
00438              PRINT *, '| | | | | | | | ib: ', ib
00439              PRINT *, '| | | | | | | | Point coords: ',plat,plon
00440              PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00441              PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00442              PRINT *, '| | | | | | | | Current i,j : ',iguess, jguess
00443              PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00444              PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00445 #endif
00446              nearest_neighbor = .true.
00447 
00448           ENDIF
00449 
00450        ENDIF ! switch to nearest neighbor
00451 
00452 ! We switch to nearest neighbor in case the iteration did not converge or produced
00453 ! wrong weights or in case the number of neighbors found was not sufficient to
00454 ! perform bilinear interpolation.
00455 
00456        IF ( nearest_neighbor ) THEN
00457 
00458           dist2v = 0.0D0
00459 
00460           DO ib_bis = 1, 8
00461              IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00462           ENDDO
00463 
00464           nb = ib_bis - 1
00465 
00466 #ifdef DEBUG
00467           PRINT *, 'nrst_neighbor, ib, lat, lon, z, nb= ',ib,plat,plon,pz,nb
00468 #endif
00469 
00470           DO ib_bis = 1, nb
00471 
00472              srch_add = ida_neighbor_index(ib, ib_bis)
00473 
00474              arg = SIN(plat)*SIN(dda_src_lat(srch_add)) +  &
00475                    COS(plat)*COS(dda_src_lat(srch_add)) *  &
00476                   (COS(plon)*COS(dda_src_lon(srch_add)) +  &
00477                    SIN(plon)*SIN(dda_src_lon(srch_add)))
00478 
00479              IF (arg < -1.0d0) THEN
00480                 arg = -1.0d0
00481              ELSE IF (arg > 1.0d0) THEN
00482                 arg = 1.0d0
00483              END IF
00484 
00485              distance = (a+pz)*ACOS(arg)
00486 
00487              dist2v(ib_bis) = SQRT(distance**2 + (pz - dda_src_z(srch_add))**2)
00488 
00489           END DO
00490        
00491           DO ib_bis = 1, nb
00492              IF (dist2v(ib_bis) == 0.0D0 ) THEN
00493 #ifdef DEBUG
00494                  PRINT *, 'dist2v(ib_bis) == 0.0, EXIT'
00495 #endif
00496                  EXIT
00497              ENDIF
00498           END DO
00499 
00500           IF (ib_bis <= nb) THEN
00501 !             ... Point to be interpolated is located on Point i
00502              dda_weights(ib,ib_bis) = 1.0d0
00503           ELSE
00504              dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00505              rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00506              dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00507           ENDIF
00508 
00509        END IF
00510        
00511        nearest_neighbor = .false.
00512 
00513     END IF
00514     
00515  END DO
00516 
00517 #else
00518 !
00519 ! 1. In the loop
00520 !
00521 ! 1.1. Set the lat and lon of the TARGET point
00522 
00523   plat = dda_tgt_lat
00524   plon = dda_tgt_lon
00525   pz   = dda_tgt_z
00526 
00527   DO ib = 1, id_tgt_size
00528      if ( ida_neighbor_index(ib, 8) <= 0  .AND. &
00529           ida_neighbor_index(ib, 1)  > 0  .AND. &
00530           ida_tgt_mask(ib) == 1 ) nearest_neighbor(ib) = .TRUE.
00531   END DO
00532 
00533   DO ib = 1, id_tgt_size
00534 
00535 ! 1.2. For the target point ib, check that the eight neighours were found
00536 
00537      IF ( .not. nearest_neighbor(ib) .AND. &
00538           ida_tgt_mask(ib) == 1      .AND. &
00539           ida_neighbor_index(ib, 1) > 0 ) THEN
00540 
00541 ! 1.2.2. Compute the weights for a trilinear method
00542 
00543           DO ib_bis = 1, 8
00544              src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00545           END DO
00546 
00547           DO ib_bis = 1, 8
00548             vec_lon(ib_bis) = src_lons (ib_bis) - plon(ib)
00549             IF (vec_lon(ib_bis) > dble_pi) THEN
00550                 src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00551             ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00552                 src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00553             ENDIF
00554           ENDDO
00555 
00556         dth1(ib) = dda_src_lat (ida_neighbor_index(ib, 2)) - &
00557                    dda_src_lat (ida_neighbor_index(ib, 1))  ! src_lats(2) - src_lats(1)
00558         dth2(ib) = dda_src_lat (ida_neighbor_index(ib, 4)) - &
00559                    dda_src_lat (ida_neighbor_index(ib, 1))  ! src_lats(4) - src_lats(1)
00560         dth3(ib) = dda_src_lat (ida_neighbor_index(ib, 5)) - &
00561                    dda_src_lat (ida_neighbor_index(ib, 1))  ! src_lats(5) - src_lats(1)
00562 
00563         dth4(ib) = dda_src_lat (ida_neighbor_index(ib, 3)) - &
00564                    dda_src_lat (ida_neighbor_index(ib, 4)) - dth1(ib) ! src_lats(3) - src_lats(4) - dth1
00565         dth5(ib) = dda_src_lat (ida_neighbor_index(ib, 6)) - &
00566                    dda_src_lat (ida_neighbor_index(ib, 5)) - dth1(ib) ! src_lats(6) - src_lats(5) - dth1
00567         dth6(ib) = dda_src_lat (ida_neighbor_index(ib, 8)) - &
00568                    dda_src_lat (ida_neighbor_index(ib, 5)) - dth2(ib) ! src_lats(8) - src_lats(5) - dth2
00569 
00570         dth7(ib) = dda_src_lat (ida_neighbor_index(ib, 2)) - &
00571                    dda_src_lat (ida_neighbor_index(ib, 1)) + & ! src_lats(2) - src_lats(1) +
00572                    dda_src_lat (ida_neighbor_index(ib, 4)) - &
00573                    dda_src_lat (ida_neighbor_index(ib, 3)) + & ! src_lats(4) - src_lats(3) +
00574                    dda_src_lat (ida_neighbor_index(ib, 5)) - &
00575                    dda_src_lat (ida_neighbor_index(ib, 6)) + & ! src_lats(5) - src_lats(6) +
00576                    dda_src_lat (ida_neighbor_index(ib, 7)) - &
00577                    dda_src_lat (ida_neighbor_index(ib, 8))     ! src_lats(7) - src_lats(8)
00578 
00579         dph1(ib) = src_lons(2) - src_lons(1)
00580         dph1(ib) = src_lons(2) - src_lons(1)
00581         dph2(ib) = src_lons(4) - src_lons(1)
00582         dph3(ib) = src_lons(5) - src_lons(1)
00583         
00584         IF (dph1(ib) >  three*dble_pih) dph1(ib) = dph1(ib) - dble_pi2
00585         IF (dph2(ib) >  three*dble_pih) dph2(ib) = dph2(ib) - dble_pi2
00586         IF (dph3(ib) >  three*dble_pih) dph3(ib) = dph3(ib) - dble_pi2
00587         IF (dph1(ib) < -three*dble_pih) dph1(ib) = dph1(ib) + dble_pi2
00588         IF (dph2(ib) < -three*dble_pih) dph2(ib) = dph2(ib) + dble_pi2
00589         IF (dph3(ib) < -three*dble_pih) dph3(ib) = dph3(ib) + dble_pi2
00590 
00591         dphtmp1 = src_lons(3) - src_lons(4)
00592         dphtmp2 = src_lons(6) - src_lons(5)
00593         dphtmp3 = src_lons(8) - src_lons(5)
00594 
00595         IF (dphtmp1 >  three*dble_pih) dphtmp1 = dphtmp1 - dble_pi2
00596         IF (dphtmp2 >  three*dble_pih) dphtmp2 = dphtmp2 - dble_pi2
00597         IF (dphtmp3 >  three*dble_pih) dphtmp3 = dphtmp3 - dble_pi2
00598         IF (dphtmp1 < -three*dble_pih) dphtmp1 = dphtmp1 + dble_pi2
00599         IF (dphtmp2 < -three*dble_pih) dphtmp2 = dphtmp2 + dble_pi2
00600         IF (dphtmp3 < -three*dble_pih) dphtmp3 = dphtmp3 + dble_pi2
00601 
00602         dph4(ib) = dphtmp1 - dph1(ib)
00603         dph5(ib) = dphtmp2 - dph1(ib)
00604         dph6(ib) = dphtmp3 - dph2(ib)
00605 
00606         dphtmp3 = src_lons(7) - src_lons(8)
00607 
00608         IF (dphtmp3 >  three*dble_pih) dphtmp3 = dphtmp3 - dble_pi2
00609         IF (dphtmp3 < -three*dble_pih) dphtmp3 = dphtmp3 + dble_pi2
00610 
00611         dph7(ib) = dph1(ib) - dphtmp1 - dphtmp2 + dphtmp3
00612 
00613         dz1(ib) = dda_src_z (ida_neighbor_index(ib, 2)) - &
00614                   dda_src_z (ida_neighbor_index(ib, 1))  ! src_zs(2) - src_zs(1)
00615         dz2(ib) = dda_src_z (ida_neighbor_index(ib, 4)) - &
00616                   dda_src_z (ida_neighbor_index(ib, 1))  ! src_zs(4) - src_zs(1)
00617         dz3(ib) = dda_src_z (ida_neighbor_index(ib, 5)) - &
00618                   dda_src_z (ida_neighbor_index(ib, 1))  ! src_zs(5) - src_zs(1)
00619 
00620         dz4(ib) = dda_src_z (ida_neighbor_index(ib, 3)) - &
00621                   dda_src_z (ida_neighbor_index(ib, 4)) - dz1(ib) ! src_zs(3) - src_zs(4) - dz1
00622         dz5(ib) = dda_src_z (ida_neighbor_index(ib, 6)) - &
00623                   dda_src_z (ida_neighbor_index(ib, 5)) - dz1(ib) ! src_zs(6) - src_zs(5) - dz1
00624         dz6(ib) = dda_src_z (ida_neighbor_index(ib, 8)) - &
00625                   dda_src_z (ida_neighbor_index(ib, 5)) - dz2(ib) ! src_zs(8) - src_zs(5) - dz2
00626 
00627         dz7(ib) = dda_src_z (ida_neighbor_index(ib, 2)) - &
00628                   dda_src_z (ida_neighbor_index(ib, 1)) + & ! src_zs(2) - src_zs(1) +
00629                   dda_src_z (ida_neighbor_index(ib, 4)) - &
00630                   dda_src_z (ida_neighbor_index(ib, 3)) + & ! src_zs(4) - src_zs(3) +
00631                   dda_src_z (ida_neighbor_index(ib, 5)) - &
00632                   dda_src_z (ida_neighbor_index(ib, 6)) + & ! src_zs(5) - src_zs(6) +
00633                   dda_src_z (ida_neighbor_index(ib, 7)) - &
00634                   dda_src_z (ida_neighbor_index(ib, 8))     ! src_zs(7) - src_zs(8)
00635      ENDIF
00636 
00637   ENDDO
00638 
00639   DO ib_part = 1, id_tgt_size, 5000
00640   DO ib = ib_part, min(id_tgt_size,ib_part+5000-1)
00641 
00642      IF ( .not. nearest_neighbor(ib) .AND. &
00643           ida_tgt_mask(ib) == 1      .AND. &
00644           ida_neighbor_index(ib, 1) > 0 ) THEN
00645 
00646         iguess = half
00647         jguess = half
00648         kguess = half
00649        
00650         ! Direct methods exist to calculate the weights. We need to
00651         ! investigate whether this is more performant. rr. 
00652         DO iter = 1, PSMILe_trans_Max_iter
00653 
00654            dthp = plat(ib) - dda_src_lat(ida_neighbor_index(ib, 1)) &
00655                            - dth1(ib)*iguess        &
00656                            - dth2(ib)*jguess        &
00657                            - dth3(ib)*kguess        &
00658                            - dth4(ib)*iguess*jguess &
00659                            - dth5(ib)*iguess*kguess &
00660                            - dth6(ib)*jguess*kguess &
00661                            - dth7(ib)*iguess*jguess*kguess
00662 
00663            dphp = plon(ib) - dda_src_lon (ida_neighbor_index(ib, 1))
00664 
00665            IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00666            IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00667 
00668            dphp = dphp - dph1(ib)*iguess        &
00669                        - dph2(ib)*jguess        &
00670                        - dph3(ib)*kguess        &
00671                        - dph4(ib)*iguess*jguess &
00672                        - dph5(ib)*iguess*kguess &
00673                        - dph6(ib)*jguess*kguess &
00674                        - dph7(ib)*iguess*jguess*kguess
00675 
00676 
00677            dzp = pz(ib) - dda_src_z (ida_neighbor_index(ib, 1))         &
00678                         - dz1(ib)*iguess        &
00679                         - dz2(ib)*jguess        &
00680                         - dz3(ib)*kguess        &
00681                         - dz4(ib)*iguess*jguess &
00682                         - dz5(ib)*iguess*kguess &
00683                         - dz6(ib)*jguess*kguess &
00684                         - dz7(ib)*iguess*jguess*kguess
00685 
00686 
00687            mat1 = dth1(ib) + dth4(ib)*jguess + dth5(ib)*kguess + dth7(ib)*jguess*kguess
00688            mat2 = dth2(ib) + dth4(ib)*iguess + dth6(ib)*kguess + dth7(ib)*iguess*kguess
00689            mat3 = dth3(ib) + dth5(ib)*iguess + dth6(ib)*jguess + dth7(ib)*iguess*jguess
00690 
00691            mat4 = dph1(ib) + dph4(ib)*jguess + dph5(ib)*kguess + dph7(ib)*jguess*kguess
00692            mat5 = dph2(ib) + dph4(ib)*iguess + dph6(ib)*kguess + dph7(ib)*iguess*kguess
00693            mat6 = dph3(ib) + dph5(ib)*iguess + dph6(ib)*jguess + dph7(ib)*iguess*jguess
00694 
00695            mat7 = dz1(ib) + dz4(ib)*jguess + dz5(ib)*kguess + dz7(ib)*jguess*kguess
00696            mat8 = dz2(ib) + dz4(ib)*iguess + dz6(ib)*kguess + dz7(ib)*iguess*kguess
00697            mat9 = dz3(ib) + dz5(ib)*iguess + dz6(ib)*jguess + dz7(ib)*iguess*jguess
00698 
00699 
00700            determinant = mat1*mat5*mat9 &
00701                        + mat2*mat6*mat7 &
00702                        + mat3*mat4*mat8 &
00703                        - mat7*mat5*mat3 &
00704                        - mat8*mat6*mat1 &
00705                        - mat9*mat4*mat2
00706 
00707            IF ( determinant == 0.0 ) EXIT 
00708 
00709            deti =  dthp*mat5*mat9 &
00710                 +  mat2*mat6*dzp  &
00711                 +  mat3*dphp*mat8 &
00712                 -  dzp*mat5*mat3  &
00713                 -  mat8*mat6*dthp &
00714                 -  mat9*dphp*mat2
00715 
00716            detj =  mat1*dphp*mat9 &
00717                 +  dthp*mat6*mat7 &
00718                 +  mat3*mat4*dzp  &
00719                 -  mat7*dphp*mat3 &
00720                 -  dzp*mat6*mat1  &
00721                 -  mat9*mat4*dthp
00722 
00723            detk =  mat1*mat5*dzp  &
00724                 +  mat2*dphp*mat7 &
00725                 +  dthp*mat4*mat8 &
00726                 -  mat7*mat5*dthp &
00727                 -  mat8*dphp*mat1 &
00728                 -  dzp*mat4*mat2
00729 
00730            deli = deti/determinant
00731            delj = detj/determinant
00732            delk = detk/determinant
00733 
00734            IF (ABS(deli) < converge .and. &
00735                ABS(delj) < converge .and. &
00736                ABS(delk) < converge) EXIT
00737 
00738            iguess = iguess + deli
00739            jguess = jguess + delj
00740            kguess = kguess + delk
00741 
00742         END DO
00743 
00744         IF (iter <= PSMILe_trans_Max_iter) THEN
00745 
00746            !***
00747            !*** successfully found i,j - compute weights
00748            !***
00749            
00750            IF (iguess < 0.) iguess = 0.0
00751            IF (jguess < 0.) jguess = 0.0
00752            IF (kguess < 0.) kguess = 0.0
00753 
00754            IF (iguess > 1.) iguess = 1.0
00755            IF (jguess > 1.) jguess = 1.0
00756            IF (kguess > 1.) kguess = 1.0
00757 
00758            dda_weights(ib,1) = abs((1.0D0-iguess)*(1.0D0-jguess)*(1.0D0-kguess))
00759            dda_weights(ib,2) = abs(       iguess *(1.0D0-jguess)*(1.0D0-kguess))
00760            dda_weights(ib,3) = abs(       iguess *       jguess *(1.0D0-kguess))
00761            dda_weights(ib,4) = abs((1.0D0-iguess)*       jguess *(1.0D0-kguess))
00762            dda_weights(ib,5) = abs((1.0D0-iguess)*(1.0D0-jguess)*       kguess)
00763            dda_weights(ib,6) = abs(       iguess *(1.0D0-jguess)*       kguess)
00764            dda_weights(ib,7) = abs(       iguess *       jguess *       kguess)
00765            dda_weights(ib,8) = abs((1.0D0-iguess)*       jguess *       kguess)
00766 
00767            IF ( ABS(SUM(dda_weights(ib,1:8))-1.0D0) > tolerance ) THEN
00768               ! Something went wrong with the weights,
00769               ! reset and try again with nearest neighbor.
00770 #ifdef VERBOSE
00771               PRINT *, '| | | | | | | | We switch to nearest neighbor.', ABS(SUM(dda_weights(ib,1:8))-1.0D0)
00772 #endif
00773               dda_weights(ib,:) = 0.0
00774               nearest_neighbor(ib) = .true.
00775            ENDIF
00776 
00777         ELSE
00778 
00779 #ifdef VERBOSE
00780            PRINT *, '| | | | | | | | ib: ', ib
00781            PRINT *, '| | | | | | | | Point coords: ',plat(ib),plon(ib)
00782            PRINT *, '| | | | | | | | Dest grid lats: ',dda_src_lat (ida_neighbor_index(ib, :))
00783            PRINT *, '| | | | | | | | Dest grid lons: ',dda_src_lon (ida_neighbor_index(ib, :))
00784            PRINT *, '| | | | | | | | Current i,j : ',iguess, jguess
00785            PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00786            PRINT *, '| | | | | | | | | We switch to nearest neighbor.'
00787 #endif
00788            nearest_neighbor(ib) = .true.
00789 
00790         ENDIF
00791 
00792      ENDIF
00793 
00794   END DO
00795   END DO
00796 
00797 ! We switch to nearest neighbor in case the iteration did not converge or produced
00798 ! wrong weights or in case the number of neighbors found was not sufficient to
00799 ! perform bilinear interpolation.
00800 
00801   DO ib = 1, id_tgt_size
00802 
00803      IF ( nearest_neighbor(ib) ) THEN
00804 
00805         dist2v = 0.0D0
00806 
00807         DO ib_bis = 1, 8
00808            IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00809         ENDDO
00810 
00811         nb = ib_bis - 1
00812 
00813         DO ib_bis = 1, nb
00814 
00815            srch_add = ida_neighbor_index(ib, ib_bis)
00816 
00817            arg = SIN(plat(ib))*SIN(dda_src_lat(srch_add)) +  &
00818                  COS(plat(ib))*COS(dda_src_lat(srch_add)) *  &
00819                 (COS(plon(ib))*COS(dda_src_lon(srch_add)) +  &
00820                  SIN(plon(ib))*SIN(dda_src_lon(srch_add)))
00821 
00822            IF (arg < -1.0d0) THEN
00823               arg = -1.0d0
00824            ELSE IF (arg > 1.0d0) THEN
00825               arg = 1.0d0
00826            END IF
00827 
00828            distance = (a+pz(ib))*ACOS(arg)
00829 
00830            dist2v(ib_bis) = ABS(distance**2 + (pz(ib) - dda_src_z(srch_add))**2)
00831 
00832         END DO
00833 
00834         DO ib_bis = 1, nb
00835            IF (dist2v(ib_bis) == 0.0D0 ) EXIT
00836         END DO
00837 
00838         IF (ib_bis <= nb) THEN
00839            !             ... Point to be interpolated is located on Point i
00840            dda_weights(ib,ib_bis) = 1.0d0
00841         ELSE
00842            dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00843            rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00844            dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00845         ENDIF
00846 
00847      END IF
00848 
00849   END DO
00850 
00851 #endif
00852 !
00853 #ifdef VERBOSE
00854  PRINT *, '| | | | | | | Quit PRISMTrs_Trilinear_weight'
00855  call psmile_flushstd
00856 #endif
00857 
00858 END SUBROUTINE PRISMTrs_Trilinear_weight
00859 
00860 
00861 
00862