prismtrs_bicubic_weight_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_weight_2d
00009 !
00010 ! !INTERFACE
00011 subroutine prismtrs_bicubic_weight_2d(id_src_size,          &
00012                                        dda_src_lat,         &
00013                                        dda_src_lon,         &
00014                                        id_tgt_size,         &
00015                                        dda_tgt_lat,         &
00016                                        dda_tgt_lon,         &
00017                                        ida_tgt_mask,        &
00018                                        ida_neighbor_index,  & 
00019                                        ida_same_lat,        &
00020                                        dda_weights,         &
00021                                        id_err)
00022 
00023 !
00024 ! !USES:
00025 !
00026   USE PRISMDrv, dummy_INTERFACE => PRISMTrs_Bicubic_weight_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  
00038 ! latitudes and longitudes of the source grid
00039   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lat
00040   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lon
00041 
00042 ! latitudes and longitudes 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 
00046   INTEGER, DIMENSION(id_tgt_size), INTENT (IN) :: ida_tgt_mask
00047 
00048 ! indices of the found neighbors on the source grid for each TARGET point
00049   INTEGER, DIMENSION(id_tgt_size,16), INTENT (IN) :: ida_neighbor_index
00050 
00051 ! integer 0 or 1 to know if a rotation has to be performed at the pole for
00052 ! gaussian reduced source grid
00053   INTEGER, DIMENSION(id_tgt_size), INTENT (IN) :: ida_same_lat
00054 
00055 ! ! RETURN VALUE
00056 !
00057 ! bicubic weights for four corners
00058   DOUBLE PRECISION, DIMENSION(id_tgt_size,16), INTENT (Out) :: dda_weights           
00059 
00060   INTEGER, INTENT (Out)                        :: id_err   ! error value
00061 
00062 ! !DESCRIPTION
00063 ! SUBROUTINE "PRISMTrs_Bicubic_weight_2d" computes the weights with 
00064 ! the bicubic interpolation method.
00065 ! The inputs are the neighboring adresses computed in the PSMILe.
00066 !
00067 ! !REVISED HISTORY
00068 !   Date      Programmer   Description
00069 ! ----------  ----------   -----------
00070 ! 06/12/2002  D. Declat    Creation
00071 ! 22/02/2005  R. Redler    Switch to nearest neighbor interpolation corrected.
00072 !
00073 ! EOP
00074 !----------------------------------------------------------------------
00075 ! $Id: prismtrs_bicubic_weight_2d.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_bicubic_weight_2d.F90 2963 2011-02-17 14:52:53Z coquart $'
00083 
00084 ! loop indices
00085   INTEGER                             :: ib, ib_bis, ib_ter, index
00086   DOUBLE PRECISION, PARAMETER         :: tolerance = 1.0D-6
00087 
00088 ! vectors of lon and lat of the source points used to compute the weights
00089   DOUBLE PRECISION, DIMENSION(4,4)    :: src_lats
00090   DOUBLE PRECISION, DIMENSION(4,4)    :: src_lons
00091   DOUBLE PRECISION, DIMENSION(4,4)    :: dla_wght_lat
00092   DOUBLE PRECISION, DIMENSION(4,4)    :: dla_wght_lon
00093   DOUBLE PRECISION, DIMENSION(4,4)    :: vec_lon
00094   DOUBLE PRECISION                    :: common_grid_range(2)
00095   DOUBLE PRECISION                    :: shiftSize
00096 
00097 ! lat/lon coords of destination point
00098   DOUBLE PRECISION                    :: plat, plon,plon_new
00099 ! for computing dist-weighted avg
00100   DOUBLE PRECISION                    :: dist2v(16)
00101   DOUBLE PRECISION                    :: rdistance, arg
00102 !
00103 #ifdef DEBUGX
00104   DOUBLE PRECISION                    :: dbl_rad2deg
00105 #endif
00106 !
00107   DOUBLE PRECISION                    :: vareps
00108   DOUBLE PRECISION, DIMENSION(4)      :: temp
00109   INTEGER                             :: nb,i,info
00110   INTEGER                             :: icase
00111   INTEGER                             :: srch_add
00112   LOGICAL                             :: nearest_neighbor
00113 
00114 ! ---------------------------------------------------------------------
00115 !
00116 #ifdef VERBOSE
00117   PRINT *, '| | | | | | | Enter PRISMTrs_Bicubic_weight_2d'
00118   call psmile_flushstd
00119 #endif
00120 
00121   id_err      = 0
00122   dda_weights = 0.0
00123   nearest_neighbor = .FALSE.
00124   common_grid_range(1) = -dble_pi
00125   common_grid_range(2) =  dble_pi
00126   shiftSize = common_grid_range(2) - common_grid_range(1)
00127 !
00128 ! 0. Treat the lat and lon arrays
00129 !
00130 ! 0.1. convert longitudes to 0,2pi interval
00131 !
00132   DO ib = 1, id_src_size
00133     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00134       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00135     ENDDO
00136     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00137       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00138     ENDDO
00139   ENDDO
00140   DO ib = 1, id_tgt_size
00141     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00142       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00143     ENDDO
00144     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00145       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00146     ENDDO
00147   ENDDO
00148 !
00149 ! 0.2. make sure input latitude range is within the machine values
00150 !      for +/- pi/2
00151 !
00152   DO ib = 1, id_src_size
00153     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00154         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00155         CALL psmile_abort
00156     ENDIF
00157   ENDDO
00158   DO ib = 1, id_tgt_size
00159     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00160         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00161         CALL psmile_abort
00162     ENDIF
00163   ENDDO
00164 
00165 !
00166 ! 1. In the loop
00167 
00168   vareps = EPSILON(ABS(dda_tgt_lat(1)))
00169 !
00170   DO ib = 1, id_tgt_size
00171     
00172     ! 1.1. Set the lat and lon of the TARGET point
00173     
00174     plat = dda_tgt_lat(ib)
00175     plon = dda_tgt_lon(ib)
00176 
00177     plon_new=plon+dble_pi
00178     
00179     ! 1.2. For the TARGET point ib, check that the 16 neighours were found
00180     
00181     ! if no neighbour at all, cycle
00182     
00183     IF ( ida_neighbor_index(ib, 1) <= 0 ) CYCLE
00184     
00185     IF ( ida_neighbor_index(ib, 16) > 0 ) THEN
00186         
00187             
00188             ! 1.2.1. For the TARGET point ib, set the lons of the sixteen neighours
00189             
00190             index = 0
00191             DO ib_ter = 1, 4
00192               DO ib_bis = 1, 4
00193                 index = index + 1
00194                 src_lons (ib_bis, ib_ter) = dda_src_lon (ida_neighbor_index(ib, index))
00195                 src_lats (ib_bis, ib_ter) = dda_src_lat (ida_neighbor_index(ib, index))
00196               END DO
00197             END DO
00198             DO ib_ter = 1, 4
00199               DO ib_bis = 1, 4
00200                 vec_lon(ib_bis,ib_ter) = src_lons (ib_bis, ib_ter) - plon
00201                 IF (vec_lon(ib_bis,ib_ter) > dble_pi) THEN
00202                   src_lons (ib_bis, ib_ter) = src_lons (ib_bis, ib_ter) - dble_pi2
00203                 ELSE IF (vec_lon(ib_bis,ib_ter) < -dble_pi) THEN
00204                    src_lons (ib_bis, ib_ter)  = src_lons (ib_bis, ib_ter) + dble_pi2
00205                 ENDIF
00206               ENDDO
00207             ENDDO
00208             !
00209 #ifdef DEBUGX
00210             dbl_rad2deg = 360.0/dble_pi2
00211             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00212                  ib==2312 ) THEN
00213                 PRINT* 
00214                 PRINT*, '+++++++++++++++++++++++++++++++'
00215                 PRINT*, 'EPIOT number :',ib
00216                 PRINT*, '+++++++++++++++++++++++++++++++'
00217                 PRINT*
00218                 PRINT*, 'lat and lon of TARGET point :',dbl_rad2deg*dda_tgt_lat(ib),dbl_rad2deg*dda_tgt_lon(ib)
00219                 DO ib_bis=1,16
00220                   PRINT*, 'EPIOS number, lat and lon of neighbours:',ida_neighbor_index(ib, ib_bis),&
00221                      dbl_rad2deg*dda_src_lat (ida_neighbor_index(ib, ib_bis)), &
00222                      dbl_rad2deg*dda_src_lon (ida_neighbor_index(ib, ib_bis))
00223                 ENDDO
00224                 PRINT*
00225             ENDIF
00226 #endif
00227             !
00228             DO ib_ter = 1, 4
00229               IF ( (ABS(src_lons(1,ib_ter) - src_lons(2,ib_ter)) .LE. vareps)  .OR. &
00230                  (ABS(src_lons(2,ib_ter) - src_lons(3,ib_ter)) .LE. vareps) .OR. &
00231                  (ABS(src_lons(3,ib_ter) - src_lons(4,ib_ter)) .LE. vareps) .OR. &
00232                  (ABS(src_lons(1,ib_ter) - src_lons(3,ib_ter)) .LE. vareps) .OR. &
00233                  (ABS(src_lons(1,ib_ter) - src_lons(4,ib_ter)) .LE. vareps) .OR. &
00234                  (ABS(src_lons(2,ib_ter) - src_lons(4,ib_ter)) .LE. vareps) ) THEN
00235                   PRINT*, 'WARNING : two points have same longitude '
00236                   PRINT*, 'Index and lons :',ib_ter,src_lats(1,ib_ter)*(180./3.141592653589793), &
00237                      src_lons(:,ib_ter)*(180./3.141592653589793)
00238                   PRINT*, 'We switch to nearest neighbor'
00239                   nearest_neighbor = .TRUE.
00240               ENDIF
00241             ENDDO
00242 
00243             IF (ida_same_lat(ib) .EQ. 1) THEN
00244                 !
00245                 IF ((ABS(src_lats(1,3) - src_lats(1,4)) .LE. vareps)) THEN
00246                     ! Upper source line is on the other side of the pole
00247                     IF (.NOT. nearest_neighbor) THEN
00248                         !
00249                         ! Rotate target point latitude by 90
00250                         CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00251                         ! Rotate first point of each line for the 3 bottom lines by 90; use of plon???
00252                         DO ib_bis = 1, 3
00253                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00254                         ENDDO
00255                         ! Rotate the first point of the line which is on the other side of the pole
00256                         DO ib_bis = 4, 4
00257                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon_new,dble_pih,vareps)
00258                         ENDDO
00259                         !
00260                         ! Apply rotation to the other points of each line
00261                         DO ib_bis = 2, 4
00262                           src_lats(ib_bis,:)=src_lats(1,:)
00263                         ENDDO
00264                         ! For the last row, inverse the longitude (maybe not needed?)
00265                         temp(1:4)=src_lons(1:4,3)
00266                         !
00267                         src_lons(1:4,4)=temp(4:1:-1)
00268                         !
00269                         ! NB 21/01/2011: a clean rotation of 90 degrees should be applied for the target point and the source
00270                         ! points instead of this weird partial rotation of the latitude only (in pole_transfo_rotation)
00271                     ENDIF
00272                     !
00273                 ELSE IF ((ABS(src_lats(1,2) - src_lats(1,3)) .LE. vareps) .AND. &
00274                    (ABS(src_lats(1,4) - src_lats(1,1)) .LE. vareps) ) THEN
00275                     ! Two upper source lines are on the other side of the pole
00276                     IF (.NOT. nearest_neighbor ) THEN
00277                         !
00278 
00279                         CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00280 
00281                         DO ib_bis = 1, 2
00282                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00283                         ENDDO
00284                         !
00285                         DO ib_bis = 3, 4
00286                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon_new,dble_pih,vareps)
00287                         ENDDO
00288                         !
00289 
00290                         !
00291                         DO ib_bis = 2, 4
00292                           src_lats(ib_bis,:)=src_lats(1,:)
00293                         ENDDO
00294                         !
00295                         temp(1:4)=src_lons(1:4,2)
00296                         !
00297                         src_lons(1:4,3)=temp(4:1:-1)
00298                         !
00299                         temp(1:4)=src_lons(1:4,1)
00300                         !
00301                         src_lons(1:4,4)=temp(4:1:-1)
00302                         !
00303                     ENDIF
00304                     !
00305                     !
00306                 ELSE IF ((ABS(src_lats(1,2) - src_lats(1,1)) .LE. vareps) ) THEN
00307                     !
00308                     IF (.NOT. nearest_neighbor) THEN
00309                         !
00310                         CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00311                         !
00312                         !
00313                         DO ib_bis = 1, 1
00314                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon_new,dble_pih,vareps)
00315                         ENDDO
00316                         !
00317                         DO ib_bis = 2, 4
00318                           CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00319                         ENDDO
00320                         !
00321                         DO ib_bis = 2, 4
00322                           src_lats(ib_bis,:)=src_lats(1,:)
00323                         ENDDO
00324                         !
00325                         temp(1:4)=src_lons(1:4,2)
00326                         !
00327                         src_lons(1:4,1)=temp(4:1:-1)
00328                         !
00329                     ENDIF
00330                     !
00331                 ENDIF
00332             ENDIF
00333             !
00334             !
00335             ! Correction of longitudes when crossing 0
00336             !
00337             DO ib_ter = 1, 4
00338               IF ( ABS(src_lons(1,ib_ter)-src_lons(4,ib_ter)) > dble_pi ) THEN
00339                   DO ib_bis = 1, 4
00340                     IF ( src_lons(ib_bis,ib_ter) < dble_pi ) &
00341                        src_lons(ib_bis,ib_ter) = src_lons(ib_bis,ib_ter) + dble_pi2
00342                   ENDDO
00343                   IF ( plon < dble_pi ) plon = plon + dble_pi2
00344               ENDIF
00345             ENDDO
00346             !
00347 #ifdef DEBUGX
00348             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00349                  ib==2312 .OR. &
00350                  ((ib .GE. 2150) .AND. (ib .LE. 2197)) .OR. ((ib .GE. 4872) .AND. (ib .LE. 4918))) THEN
00351                 PRINT*
00352                 PRINT*, 'After rotation :'
00353                 PRINT*, 'EPIOT number, tgt_plat, tgt_plon :',ib,dbl_rad2deg*plat,dbl_rad2deg*plon
00354                 DO ib_bis = 1, 4
00355                   PRINT*, 'src_lats by line:',dbl_rad2deg*src_lats(1:4,ib_bis)
00356                   PRINT*, 'src_lons by line:',dbl_rad2deg*src_lons(1:4,ib_bis)
00357                 ENDDO
00358                 PRINT*
00359             ENDIF
00360 #endif
00361             !
00362 #ifdef DEBUGX
00363             !
00364             ! test of the rotation for a target point not at the pole
00365             !
00366             !
00367             ! North pole
00368             !
00369             IF (ib==2131) THEN
00370                 !
00371                 CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00372                 !
00373                 DO ib_bis = 1, 4
00374                   CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00375                 ENDDO
00376                 !
00377                 DO ib_bis = 2, 4
00378                   src_lats(ib_bis,:)=src_lats(1,:)
00379                 ENDDO
00380                 !
00381             ENDIF
00382             !
00383             ! South pole
00384             !
00385             IF (ib==333) THEN
00386                 !
00387                 CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00388                 !
00389                 DO ib_bis = 1, 4
00390                   CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00391                 ENDDO
00392                 !
00393                 DO ib_bis = 2, 4
00394                   src_lats(ib_bis,:)=src_lats(1,:)
00395                 ENDDO
00396                 !
00397             ENDIF
00398             !
00399 #endif
00400             ! Calculation of weights for longitudes
00401             !
00402             IF ( .NOT. nearest_neighbor) THEN
00403                 !
00404                 DO ib_ter = 1, 4
00405                   CALL calcul_wght_irreg(src_lons(:,ib_ter), plon, dla_wght_lon(:,ib_ter))
00406                 ENDDO
00407                 !
00408               ENDIF
00409            !
00410             !
00411            ! Calculation of weights for latitudes.
00412            ! It is required that all four lat colums are identical
00413             !
00414               DO ib_bis = 1, 4
00415                 IF ( (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,2)) .LE. vareps) .OR. & 
00416                    (ABS(src_lats(ib_bis,2) - src_lats(ib_bis,3)) .LE. vareps) .OR. &
00417                    (ABS(src_lats(ib_bis,3) - src_lats(ib_bis,4)) .LE. vareps) .OR. &
00418                    (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,3)) .LE. vareps) .OR. &
00419                    (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,4)) .LE. vareps) .OR. &
00420                    (ABS(src_lats(ib_bis,2) - src_lats(ib_bis,4)) .LE. vareps) ) THEN
00421                     PRINT*, 'WARNING : two points have same latitude :'
00422                     PRINT*, 'Index and lats :',ib_bis,src_lons(ib_bis,1)*(180./3.141592653589793)&
00423                       ,src_lats(ib_bis,:)*(180./3.141592653589793)
00424                     PRINT*, 'We switch to nearest neighbor'
00425                     nearest_neighbor = .TRUE.
00426                 ENDIF
00427               ENDDO
00428             !
00429               IF ( .NOT. nearest_neighbor ) THEN
00430                   !
00431 !!$                  CALL calcul_wght_irreg(src_lats(1,:), plat, dla_wght_lat(1,:))
00432 !!$                  !
00433 !!$                  DO ib_bis=2,4
00434 !!$                    dla_wght_lat(ib_bis,:) = dla_wght_lat(1,:)
00435 !!$                  ENDDO
00436                   !
00437                   DO ib_bis = 1, 4
00438                     CALL calcul_wght_irreg(src_lats(ib_bis,:), plat, dla_wght_lat(ib_bis,:))
00439                   ENDDO
00440               ENDIF
00441             !
00442            ! 
00443            ! Calculation of total weight, elementwise multiplication
00444            !
00445            IF ( .NOT. nearest_neighbor ) THEN
00446                !
00447                index = 0
00448                DO ib_ter = 1, 4
00449                  DO ib_bis = 1, 4
00450                    index = index + 1
00451                    dda_weights(ib,index)=dla_wght_lat(ib_bis,ib_ter)*dla_wght_lon(ib_bis,ib_ter)
00452                  END DO
00453                END DO
00454                !
00455 #ifdef DEBUGX
00456             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00457                  ib==2312 ) THEN
00458                    PRINT*
00459                    PRINT*, 'EPIOT number, weights lon :',ib,dla_wght_lon(:,:)
00460                    PRINT*, 'EPIOT number, weights lat :',ib,dla_wght_lat(:,:)
00461                    PRINT*
00462                ENDIF
00463 #endif
00464                !
00465            ENDIF
00466 
00467            IF ( ABS(SUM(dda_weights(ib,1:16))-1.0D0) > tolerance .OR. nearest_neighbor) THEN
00468               
00469               ! Something went wrong with the weights,
00470               ! reset and try again with nearest neighbor.
00471               
00472 #ifdef VERBOSE
00473               PRINT *, '| | | | | | | | We switch to nearest neighbor for :',ib,ABS(SUM(dda_weights(ib,1:16)))
00474 #endif
00475               dda_weights(ib,:) = 0.0
00476               nearest_neighbor = .TRUE.
00477           ENDIF
00478       
00479       ENDIF ! switch to nearest neighbor
00480 
00481 ! We switch to nearest neighbor in case the iteration did not converge or produced
00482 ! wrong weights or in case the number of neighbors found was not sufficient to
00483 ! perform bicubic interpolation.
00484 
00485      IF ( ida_neighbor_index(ib, 16) <= 0 .OR. nearest_neighbor ) THEN
00486 
00487 #ifdef VERBOSE
00488         PRINT *, '| | | | | | | | We switch to nearest neighbor'
00489 #endif
00490 
00491         plat = dda_tgt_lat(ib)
00492         plon = dda_tgt_lon(ib)
00493 
00494         dist2v = 0.0d0
00495 
00496         DO ib_bis = 1, 16
00497            IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00498         ENDDO
00499 
00500         nb = ib_bis - 1
00501 
00502         DO ib_bis = 1, nb
00503 
00504            srch_add = ida_neighbor_index(ib, ib_bis)
00505            arg = SIN(dda_tgt_lat(ib))*SIN(dda_src_lat(srch_add)) +  &
00506                  COS(dda_tgt_lat(ib))*COS(dda_src_lat(srch_add)) *  &
00507                 (COS(dda_tgt_lon(ib))*COS(dda_src_lon(srch_add)) +  &
00508                  SIN(dda_tgt_lon(ib))*SIN(dda_src_lon(srch_add)))
00509 
00510            IF (arg < -1.0d0) THEN
00511               arg = -1.0d0
00512            ELSE IF (arg > 1.0d0) THEN
00513               arg = 1.0d0
00514            END IF
00515 
00516 
00517            dist2v(ib_bis) = ACOS(arg)
00518            dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00519 
00520         END DO
00521 
00522 
00523         ! Test if the target point coincide with the source point
00524         DO ib_bis = 1, nb
00525            IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00526         END DO
00527 
00528         IF (ib_bis <= nb) THEN
00529 
00530             !Point to be interpolated is located on Point i
00531             dda_weights(ib,ib_bis) = 1.0d0
00532 
00533         ELSE IF (nb == 1) THEN
00534 
00535             !Only one nearest neighbour
00536             dda_weights(ib,nb) = 1.0d0
00537 
00538         ELSE IF ( nb > 1 ) THEN
00539 
00540            dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00541            rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00542            dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00543 
00544         ENDIF
00545         !
00546         !
00547        END IF
00548        !
00549 #ifdef DEBUGX
00550             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00551                  ib==2312 ) THEN
00552            PRINT*
00553            PRINT*, 'EPIOT number, EPIOS numbers of neighbours:',ib,ida_neighbor_index(ib,:)
00554            PRINT*, 'dda_weights(ib,:) :',dda_weights(ib,:)
00555            PRINT*
00556        ENDIF
00557 #endif
00558 
00559      nearest_neighbor = .FALSE.
00560 
00561    END DO !ib=1,id_tgt_size
00562    !
00563 !
00564 #ifdef VERBOSE
00565   PRINT *, '| | | | | | | Quit PRISMTrs_Bicubic_weight_2d'
00566   call psmile_flushstd
00567 #endif
00568 
00569 END SUBROUTINE PRISMTrs_Bicubic_weight_2d
00570 
00571 
00572 
00573 !-----------------------------------------------------------------------
00574 ! BOP
00575 !
00576 ! !IROUTINE:  calcul_wght_irreg
00577 !
00578 ! !INTERFACE:
00579 ! 
00580   SUBROUTINE calcul_wght_irreg(rda_x, rd_pt, rda_wght)
00581 !  
00582 ! !USES:
00583 !              
00584 ! !RETURN VALUE:
00585 !
00586 ! array of weights for the points x
00587     DOUBLE PRECISION, DIMENSION(4), INTENT(out) :: rda_wght
00588 
00589 !      
00590 ! !PARAMETERS:
00591 ! 
00592 ! array of positions on source grid, lat or lon
00593     DOUBLE PRECISION, DIMENSION(4), INTENT(in)  :: rda_x 
00594 
00595 ! position of target point to interpolate
00596     DOUBLE PRECISION,INTENT(in)                 :: rd_pt
00597 
00598 !              
00599 ! !DESCRIPTION:
00600 !  Calculates a the weights of four points for a bicubic interpolation. 
00601 !  The distances between the points can be irregulier. 
00602 !
00603 ! !REVISION HISTORY:
00604 !  2002.10.21  J.Ghattas  created
00605 !  2005.05.04  R.Redler   adapted to OASIS4
00606 !
00607 ! EOP
00608 !
00609 ! Local declarations
00610 !
00611     DOUBLE PRECISION :: 
00612        rl_t1, rl_t2, rl_t3, rl_t4, rl_t5, rl_t6, rl_t7, rl_t8, rl_t9, 
00613        rl_u1, rl_u2, rl_u3, rl_u4, 
00614        rl_k1, rl_k2, rl_k3, 
00615        rl_d1, rl_d2, rl_d3, rl_d4, 
00616        rl_c1, rl_c2, rl_c3, rl_c4, 
00617        rl_b1, rl_b2, rl_b3, rl_b4, 
00618        rl_a1, rl_a2, rl_a3, rl_a4, 
00619        rl_y1, rl_y2, rl_y3, 
00620        rl_a1_y, rl_a2_y, rl_a3_y, rl_a4_y, 
00621        rl_b1_y, rl_b2_y, rl_b3_y, rl_b4_y, 
00622        rl_c1_y, rl_c2_y, rl_c3_y, rl_c4_y
00623 
00624     IF (rda_x(1)/=0 .and. rda_x(2)/=0 .and. rda_x(3)/=0 .and. rda_x(4)/=0) THEN
00625         
00626 
00627     rl_t1 = 1/rda_x(1) - 1/rda_x(2)
00628     rl_t2 = 1/rda_x(1)**2 - 1/rda_x(2)**2
00629     rl_t3 = 1/rda_x(1)**3 - 1/rda_x(2)**3
00630     rl_t4 = 1/rda_x(1) - 1/rda_x(3)
00631     rl_t5 = 1/rda_x(1)**2 - 1/rda_x(3)**2
00632     rl_t6 = 1/rda_x(1)**3 - 1/rda_x(3)**3
00633     rl_t7 = 1/rda_x(1) - 1/rda_x(4)
00634     rl_t8 = 1/rda_x(1)**2 - 1/rda_x(4)**2
00635     rl_t9 = 1/rda_x(1)**3 - 1/rda_x(4)**3
00636 
00637     rl_u1 = rl_t2/rl_t1 - rl_t5/rl_t4
00638     rl_u2 = rl_t3/rl_t1 - rl_t6/rl_t4
00639     rl_u3 = rl_t2/rl_t1 - rl_t8/rl_t7
00640     rl_u4 = rl_t3/rl_t1 - rl_t9/rl_t7
00641 
00642     rl_k1 = (1/(rl_t1*rl_u1)-1/(rl_t1*rl_u3)) / (rl_u2/rl_u1-rl_u4/rl_u3)
00643     rl_k2 = -1/(rl_t4*rl_u1) / (rl_u2/rl_u1-rl_u4/rl_u3)
00644     rl_k3 = 1/(rl_t7*rl_u3) / (rl_u2/rl_u1-rl_u4/rl_u3)
00645 
00646     rl_d1=(rl_k1+rl_k2+rl_k3)/rda_x(1)**3
00647     rl_d2 = -rl_k1/rda_x(2)**3
00648     rl_d3 = -rl_k2/rda_x(3)**3
00649     rl_d4 = -rl_k3/rda_x(4)**3
00650 
00651     rl_c1 = 1/rl_u1*(1/(rl_t1*rda_x(1)**3)-1/(rl_t4*rda_x(1)**3)- &
00652        rl_u2*rl_d1)
00653     rl_c2 = 1/rl_u1*(1/(-rl_t1*rda_x(2)**3)-rl_u2*rl_d2)
00654     rl_c3 = 1/rl_u1*(1/(rl_t4*rda_x(3)**3)-rl_u2*rl_d3)
00655     rl_c4 = 1/rl_u1*(-rl_u2*rl_d4)
00656     
00657     rl_b1 = 1/rl_t1/rda_x(1)**3-rl_t2/rl_t1*rl_c1-rl_t3/rl_t1*rl_d1
00658     rl_b2 = -1/rl_t1/rda_x(2)**3-rl_t2/rl_t1*rl_c2-rl_t3/rl_t1*rl_d2
00659     rl_b3 = -rl_t2/rl_t1*rl_c3-rl_t3/rl_t1*rl_d3
00660     rl_b4 = -rl_t2/rl_t1*rl_c4-rl_t3/rl_t1*rl_d4
00661 
00662     rl_a1 = 1/rda_x(1)**3-1/rda_x(1)*rl_b1-1/rda_x(1)**2*rl_c1- &
00663        1/rda_x(1)**3*rl_d1
00664     rl_a2 = -1/rda_x(1)*rl_b2-1/rda_x(1)**2*rl_c2-1/rda_x(1)**3*rl_d2
00665     rl_a3 = -1/rda_x(1)*rl_b3-1/rda_x(1)**2*rl_c3-1/rda_x(1)**3*rl_d3
00666     rl_a4 = -1/rda_x(1)*rl_b4-1/rda_x(1)**2*rl_c4-1/rda_x(1)**3*rl_d4
00667 
00668        ! Weights  
00669     rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt + rl_d1
00670     rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt + rl_d2
00671     rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt + rl_d3
00672     rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt + rl_d4
00673 
00674     ELSE ! there is one point = 0
00675       
00676     rl_d1=0; rl_d2=0; rl_d3=0; rl_d4=0
00677     
00678         ! Transformation for each case
00679     IF (rda_x(1)==0) THEN
00680         rl_y1=rda_x(2); rl_y2=rda_x(3); rl_y3=rda_x(4)
00681         rl_d1=1
00682     ELSE IF (rda_x(2)==0) THEN
00683         rl_y1=rda_x(1); rl_y2=rda_x(3); rl_y3=rda_x(4)
00684         rl_d2=1
00685     ELSE IF (rda_x(3)==0) THEN
00686         rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(4)
00687         rl_d3=1
00688     ELSE
00689         rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(3)
00690         rl_d4=1
00691     END IF
00692 
00693         ! Solving the system 
00694     rl_t1 = 1/rl_y1-1/rl_y2
00695     rl_t2 = 1/rl_y1**2-1/rl_y2**2
00696     rl_t3 = 1/rl_y1-1/rl_y3
00697     rl_t4 = 1/rl_y1**2-1/rl_y3**2
00698 
00699     rl_c1_y =(1/rl_y1**3/rl_t1-1/rl_y1**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
00700     rl_c2_y = -1/rl_y2**3/rl_t1/(rl_t2/rl_t1-rl_t4/rl_t3)
00701     rl_c3_y = 1/rl_y3**3/rl_t3/(rl_t2/rl_t1-rl_t4/rl_t3)
00702     rl_c4_y=(-1/rl_y1**3/rl_t1+1/rl_y2**3/rl_t1+ &
00703        1/rl_y1**3/rl_t3-1/rl_y3**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
00704 
00705     rl_b1_y = 1/rl_y1**3/rl_t1 - rl_c1_y*rl_t2/rl_t1
00706     rl_b2_y = -1/rl_y2**3/rl_t1 - rl_c2_y*rl_t2/rl_t1
00707     rl_b3_y = -rl_c3_y*rl_t2/rl_t1
00708     rl_b4_y = -1/rl_y1**3/rl_t1 + 1/rl_y2**3/rl_t1 - rl_c4_y*rl_t2/rl_t1
00709 
00710     rl_a1_y = 1/rl_y1**3 - rl_b1_y/rl_y1 - rl_c1_y/rl_y1**2
00711     rl_a2_y = -rl_b2_y/rl_y1 - rl_c2_y/rl_y1**2
00712     rl_a3_y = -rl_b3_y/rl_y1 - rl_c3_y/rl_y1**2
00713     rl_a4_y = -1/rl_y1**3 - rl_b4_y/rl_y1 - rl_c4_y/rl_y1**2
00714 
00715         ! Retransformation
00716     IF (rda_x(1)==0) THEN
00717         rl_a1=rl_a4_y; rl_a2=rl_a1_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
00718         rl_b1=rl_b4_y; rl_b2=rl_b1_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
00719         rl_c1=rl_c4_y; rl_c2=rl_c1_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
00720     ELSE IF (rda_x(2)==0) THEN
00721         rl_a1=rl_a1_y; rl_a2=rl_a4_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
00722         rl_b1=rl_b1_y; rl_b2=rl_b4_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
00723         rl_c1=rl_c1_y; rl_c2=rl_c4_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
00724     ELSE IF (rda_x(3)==0) THEN
00725         rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a4_y; rl_a4=rl_a3_y
00726         rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b4_y; rl_b4=rl_b3_y
00727         rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c4_y; rl_c4=rl_c3_y
00728     ELSE 
00729         rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a3_y; rl_a4=rl_a4_y
00730         rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b3_y; rl_b4=rl_b4_y
00731         rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c3_y; rl_c4=rl_c4_y
00732     END IF
00733 
00734         ! Weights  
00735     rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt +rl_d1
00736     rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt +rl_d2
00737     rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt +rl_d3
00738     rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt +rl_d4
00739 
00740     END IF
00741 
00742   END SUBROUTINE calcul_wght_irreg
00743 
00744   SUBROUTINE pole_transfo_rotation(pt_lat,pt_lon,db_pih,vareps)
00745 
00746 ! Routine given by Rene (psmile_transrot_dble) from C. Koeberle, R. Gerdes AWI Bremerhaven.
00747 ! The grid is rotated by 90 degree. The points are translated 
00748 ! by 90 deg towards the Equator. North and South pole are
00749 ! located on the Equator on the transformed grid.
00750 ! The data are in radian.
00751     !
00752   IMPLICIT NONE
00753 
00754     DOUBLE PRECISION, INTENT(INOUT)    :: pt_lat
00755     DOUBLE PRECISION, INTENT(IN)       :: pt_lon
00756     !
00757     DOUBLE PRECISION, INTENT(IN)       :: db_pih,vareps
00758     !
00759 ! Latitude
00760     !
00761     IF (ABS(pt_lon - 0.) .LE. vareps .OR. &
00762        ABS(pt_lon - 2*db_pih) .LE. vareps .OR. &
00763        ABS(pt_lon - 4*db_pih) .LE. vareps) THEN
00764         pt_lat = - COS(pt_lon) * (db_pih - pt_lat) 
00765     ELSE
00766         pt_lat = ASIN(COS( pt_lat ) * SIN( pt_lon ))
00767     ENDIF
00768     !
00769   END SUBROUTINE pole_transfo_rotation