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 3241 2011-06-23 09:55:11Z 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 3241 2011-06-23 09:55:11Z 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) = 0.
00125   common_grid_range(2) =  dble_pi2
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 #ifdef DEBUGX
00336             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00337                  ib==2312 .OR. &
00338                  ((ib .GE. 2150) .AND. (ib .LE. 2197)) .OR. ((ib .GE. 4872) .AND. (ib .LE. 4918))) THEN
00339                 PRINT*
00340                 PRINT*, 'After rotation :'
00341                 PRINT*, 'EPIOT number, tgt_plat, tgt_plon :',ib,dbl_rad2deg*plat,dbl_rad2deg*plon
00342                 DO ib_bis = 1, 4
00343                   PRINT*, 'src_lats by line:',dbl_rad2deg*src_lats(1:4,ib_bis)
00344                   PRINT*, 'src_lons by line:',dbl_rad2deg*src_lons(1:4,ib_bis)
00345                 ENDDO
00346                 PRINT*
00347             ENDIF
00348 #endif
00349             !
00350 #ifdef DEBUGX
00351             !
00352             ! test of the rotation for a target point not at the pole
00353             !
00354             !
00355             ! North pole
00356             !
00357             IF (ib==2131) THEN
00358                 !
00359                 CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00360                 !
00361                 DO ib_bis = 1, 4
00362                   CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00363                 ENDDO
00364                 !
00365                 DO ib_bis = 2, 4
00366                   src_lats(ib_bis,:)=src_lats(1,:)
00367                 ENDDO
00368                 !
00369             ENDIF
00370             !
00371             ! South pole
00372             !
00373             IF (ib==333) THEN
00374                 !
00375                 CALL pole_transfo_rotation(plat,plon,dble_pih,vareps)
00376                 !
00377                 DO ib_bis = 1, 4
00378                   CALL pole_transfo_rotation(src_lats(1,ib_bis),plon,dble_pih,vareps)
00379                 ENDDO
00380                 !
00381                 DO ib_bis = 2, 4
00382                   src_lats(ib_bis,:)=src_lats(1,:)
00383                 ENDDO
00384                 !
00385             ENDIF
00386             !
00387 #endif
00388             ! Calculation of weights for longitudes
00389             !
00390             IF ( .NOT. nearest_neighbor) THEN
00391                 !
00392                 DO ib_ter = 1, 4
00393                   CALL calcul_wght_irreg(src_lons(:,ib_ter), plon, dla_wght_lon(:,ib_ter))
00394                 ENDDO
00395                 !
00396               ENDIF
00397            !
00398             !
00399            ! Calculation of weights for latitudes.
00400            ! It is required that all four lat colums are identical
00401             !
00402               DO ib_bis = 1, 4
00403                 IF ( (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,2)) .LE. vareps) .OR. & 
00404                    (ABS(src_lats(ib_bis,2) - src_lats(ib_bis,3)) .LE. vareps) .OR. &
00405                    (ABS(src_lats(ib_bis,3) - src_lats(ib_bis,4)) .LE. vareps) .OR. &
00406                    (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,3)) .LE. vareps) .OR. &
00407                    (ABS(src_lats(ib_bis,1) - src_lats(ib_bis,4)) .LE. vareps) .OR. &
00408                    (ABS(src_lats(ib_bis,2) - src_lats(ib_bis,4)) .LE. vareps) ) THEN
00409                     PRINT*, 'WARNING : two points have same latitude :'
00410                     PRINT*, 'Index and lats :',ib_bis,src_lons(ib_bis,1)*(180./3.141592653589793)&
00411                       ,src_lats(ib_bis,:)*(180./3.141592653589793)
00412                     PRINT*, 'We switch to nearest neighbor'
00413                     nearest_neighbor = .TRUE.
00414                 ENDIF
00415               ENDDO
00416             !
00417               IF ( .NOT. nearest_neighbor ) THEN
00418                   !
00419 !!$                  CALL calcul_wght_irreg(src_lats(1,:), plat, dla_wght_lat(1,:))
00420 !!$                  !
00421 !!$                  DO ib_bis=2,4
00422 !!$                    dla_wght_lat(ib_bis,:) = dla_wght_lat(1,:)
00423 !!$                  ENDDO
00424                   !
00425                   DO ib_bis = 1, 4
00426                     CALL calcul_wght_irreg(src_lats(ib_bis,:), plat, dla_wght_lat(ib_bis,:))
00427                   ENDDO
00428               ENDIF
00429             !
00430            ! 
00431            ! Calculation of total weight, elementwise multiplication
00432            !
00433            IF ( .NOT. nearest_neighbor ) THEN
00434                !
00435                index = 0
00436                DO ib_ter = 1, 4
00437                  DO ib_bis = 1, 4
00438                    index = index + 1
00439                    dda_weights(ib,index)=dla_wght_lat(ib_bis,ib_ter)*dla_wght_lon(ib_bis,ib_ter)
00440                  END DO
00441                END DO
00442                !
00443 #ifdef DEBUGX
00444             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00445                  ib==2312 ) THEN
00446                    PRINT*
00447                    PRINT*, 'EPIOT number, weights lon :',ib,dla_wght_lon(:,:)
00448                    PRINT*, 'EPIOT number, weights lat :',ib,dla_wght_lat(:,:)
00449                    PRINT*
00450                ENDIF
00451 #endif
00452                !
00453            ENDIF
00454 
00455            IF ( ABS(SUM(dda_weights(ib,1:16))-1.0D0) > tolerance .OR. nearest_neighbor) THEN
00456               
00457               ! Something went wrong with the weights,
00458               ! reset and try again with nearest neighbor.
00459               
00460 #ifdef VERBOSE
00461               PRINT *, '| | | | | | | | We switch to nearest neighbor for :',ib,ABS(SUM(dda_weights(ib,1:16)))
00462 #endif
00463               dda_weights(ib,:) = 0.0
00464               nearest_neighbor = .TRUE.
00465           ENDIF
00466       
00467       ENDIF ! switch to nearest neighbor
00468 
00469 ! We switch to nearest neighbor in case the iteration did not converge or produced
00470 ! wrong weights or in case the number of neighbors found was not sufficient to
00471 ! perform bicubic interpolation.
00472 
00473      IF ( ida_neighbor_index(ib, 16) <= 0 .OR. nearest_neighbor ) THEN
00474 
00475 #ifdef VERBOSE
00476         PRINT *, '| | | | | | | | We switch to nearest neighbor'
00477 #endif
00478 
00479         plat = dda_tgt_lat(ib)
00480         plon = dda_tgt_lon(ib)
00481 
00482         dist2v = 0.0d0
00483 
00484         DO ib_bis = 1, 16
00485            IF ( ida_neighbor_index(ib,ib_bis) <= 0 ) EXIT
00486         ENDDO
00487 
00488         nb = ib_bis - 1
00489 
00490         DO ib_bis = 1, nb
00491 
00492            srch_add = ida_neighbor_index(ib, ib_bis)
00493            arg = SIN(dda_tgt_lat(ib))*SIN(dda_src_lat(srch_add)) +  &
00494                  COS(dda_tgt_lat(ib))*COS(dda_src_lat(srch_add)) *  &
00495                 (COS(dda_tgt_lon(ib))*COS(dda_src_lon(srch_add)) +  &
00496                  SIN(dda_tgt_lon(ib))*SIN(dda_src_lon(srch_add)))
00497 
00498            IF (arg < -1.0d0) THEN
00499               arg = -1.0d0
00500            ELSE IF (arg > 1.0d0) THEN
00501               arg = 1.0d0
00502            END IF
00503 
00504 
00505            dist2v(ib_bis) = ACOS(arg)
00506            dist2v(ib_bis) = dist2v(ib_bis)*dist2v(ib_bis)
00507 
00508         END DO
00509 
00510 
00511         ! Test if the target point coincide with the source point
00512         DO ib_bis = 1, nb
00513            IF (dist2v(ib_bis) == 0.0d0 ) EXIT
00514         END DO
00515 
00516         IF (ib_bis <= nb) THEN
00517 
00518             !Point to be interpolated is located on Point i
00519             dda_weights(ib,ib_bis) = 1.0d0
00520 
00521         ELSE IF (nb == 1) THEN
00522 
00523             !Only one nearest neighbour
00524             dda_weights(ib,nb) = 1.0d0
00525 
00526         ELSE IF ( nb > 1 ) THEN
00527 
00528            dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00529            rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00530            dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00531 
00532         ENDIF
00533         !
00534         !
00535        END IF
00536        !
00537 #ifdef DEBUGX
00538             IF ( ((ib .GE. 2198) .AND. (ib .LE. 2245)) .OR. ((ib .GE. 4919) .AND. (ib .LE. 4965)) .OR. &
00539                  ib==2312 ) THEN
00540            PRINT*
00541            PRINT*, 'EPIOT number, EPIOS numbers of neighbours:',ib,ida_neighbor_index(ib,:)
00542            PRINT*, 'dda_weights(ib,:) :',dda_weights(ib,:)
00543            PRINT*
00544        ENDIF
00545 #endif
00546 
00547      nearest_neighbor = .FALSE.
00548 
00549    END DO !ib=1,id_tgt_size
00550    !
00551 !
00552 #ifdef VERBOSE
00553   PRINT *, '| | | | | | | Quit PRISMTrs_Bicubic_weight_2d'
00554   call psmile_flushstd
00555 #endif
00556 
00557 END SUBROUTINE PRISMTrs_Bicubic_weight_2d
00558 
00559 
00560 
00561 !-----------------------------------------------------------------------
00562 ! BOP
00563 !
00564 ! !IROUTINE:  calcul_wght_irreg
00565 !
00566 ! !INTERFACE:
00567 ! 
00568   SUBROUTINE calcul_wght_irreg(rda_x, rd_pt, rda_wght)
00569 !  
00570 ! !USES:
00571 !              
00572 ! !RETURN VALUE:
00573 !
00574 ! array of weights for the points x
00575     DOUBLE PRECISION, DIMENSION(4), INTENT(out) :: rda_wght
00576 
00577 !      
00578 ! !PARAMETERS:
00579 ! 
00580 ! array of positions on source grid, lat or lon
00581     DOUBLE PRECISION, DIMENSION(4), INTENT(in)  :: rda_x 
00582 
00583 ! position of target point to interpolate
00584     DOUBLE PRECISION,INTENT(in)                 :: rd_pt
00585 
00586 !              
00587 ! !DESCRIPTION:
00588 !  Calculates a the weights of four points for a bicubic interpolation. 
00589 !  The distances between the points can be irregulier. 
00590 !
00591 ! !REVISION HISTORY:
00592 !  2002.10.21  J.Ghattas  created
00593 !  2005.05.04  R.Redler   adapted to OASIS4
00594 !
00595 ! EOP
00596 !
00597 ! Local declarations
00598 !
00599     DOUBLE PRECISION :: 
00600        rl_t1, rl_t2, rl_t3, rl_t4, rl_t5, rl_t6, rl_t7, rl_t8, rl_t9, 
00601        rl_u1, rl_u2, rl_u3, rl_u4, 
00602        rl_k1, rl_k2, rl_k3, 
00603        rl_d1, rl_d2, rl_d3, rl_d4, 
00604        rl_c1, rl_c2, rl_c3, rl_c4, 
00605        rl_b1, rl_b2, rl_b3, rl_b4, 
00606        rl_a1, rl_a2, rl_a3, rl_a4, 
00607        rl_y1, rl_y2, rl_y3, 
00608        rl_a1_y, rl_a2_y, rl_a3_y, rl_a4_y, 
00609        rl_b1_y, rl_b2_y, rl_b3_y, rl_b4_y, 
00610        rl_c1_y, rl_c2_y, rl_c3_y, rl_c4_y
00611 
00612     IF (rda_x(1)/=0 .and. rda_x(2)/=0 .and. rda_x(3)/=0 .and. rda_x(4)/=0) THEN
00613         
00614 
00615     rl_t1 = 1/rda_x(1) - 1/rda_x(2)
00616     rl_t2 = 1/rda_x(1)**2 - 1/rda_x(2)**2
00617     rl_t3 = 1/rda_x(1)**3 - 1/rda_x(2)**3
00618     rl_t4 = 1/rda_x(1) - 1/rda_x(3)
00619     rl_t5 = 1/rda_x(1)**2 - 1/rda_x(3)**2
00620     rl_t6 = 1/rda_x(1)**3 - 1/rda_x(3)**3
00621     rl_t7 = 1/rda_x(1) - 1/rda_x(4)
00622     rl_t8 = 1/rda_x(1)**2 - 1/rda_x(4)**2
00623     rl_t9 = 1/rda_x(1)**3 - 1/rda_x(4)**3
00624 
00625     rl_u1 = rl_t2/rl_t1 - rl_t5/rl_t4
00626     rl_u2 = rl_t3/rl_t1 - rl_t6/rl_t4
00627     rl_u3 = rl_t2/rl_t1 - rl_t8/rl_t7
00628     rl_u4 = rl_t3/rl_t1 - rl_t9/rl_t7
00629 
00630     rl_k1 = (1/(rl_t1*rl_u1)-1/(rl_t1*rl_u3)) / (rl_u2/rl_u1-rl_u4/rl_u3)
00631     rl_k2 = -1/(rl_t4*rl_u1) / (rl_u2/rl_u1-rl_u4/rl_u3)
00632     rl_k3 = 1/(rl_t7*rl_u3) / (rl_u2/rl_u1-rl_u4/rl_u3)
00633 
00634     rl_d1=(rl_k1+rl_k2+rl_k3)/rda_x(1)**3
00635     rl_d2 = -rl_k1/rda_x(2)**3
00636     rl_d3 = -rl_k2/rda_x(3)**3
00637     rl_d4 = -rl_k3/rda_x(4)**3
00638 
00639     rl_c1 = 1/rl_u1*(1/(rl_t1*rda_x(1)**3)-1/(rl_t4*rda_x(1)**3)- &
00640        rl_u2*rl_d1)
00641     rl_c2 = 1/rl_u1*(1/(-rl_t1*rda_x(2)**3)-rl_u2*rl_d2)
00642     rl_c3 = 1/rl_u1*(1/(rl_t4*rda_x(3)**3)-rl_u2*rl_d3)
00643     rl_c4 = 1/rl_u1*(-rl_u2*rl_d4)
00644     
00645     rl_b1 = 1/rl_t1/rda_x(1)**3-rl_t2/rl_t1*rl_c1-rl_t3/rl_t1*rl_d1
00646     rl_b2 = -1/rl_t1/rda_x(2)**3-rl_t2/rl_t1*rl_c2-rl_t3/rl_t1*rl_d2
00647     rl_b3 = -rl_t2/rl_t1*rl_c3-rl_t3/rl_t1*rl_d3
00648     rl_b4 = -rl_t2/rl_t1*rl_c4-rl_t3/rl_t1*rl_d4
00649 
00650     rl_a1 = 1/rda_x(1)**3-1/rda_x(1)*rl_b1-1/rda_x(1)**2*rl_c1- &
00651        1/rda_x(1)**3*rl_d1
00652     rl_a2 = -1/rda_x(1)*rl_b2-1/rda_x(1)**2*rl_c2-1/rda_x(1)**3*rl_d2
00653     rl_a3 = -1/rda_x(1)*rl_b3-1/rda_x(1)**2*rl_c3-1/rda_x(1)**3*rl_d3
00654     rl_a4 = -1/rda_x(1)*rl_b4-1/rda_x(1)**2*rl_c4-1/rda_x(1)**3*rl_d4
00655 
00656        ! Weights  
00657     rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt + rl_d1
00658     rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt + rl_d2
00659     rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt + rl_d3
00660     rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt + rl_d4
00661 
00662     ELSE ! there is one point = 0
00663       
00664     rl_d1=0; rl_d2=0; rl_d3=0; rl_d4=0
00665     
00666         ! Transformation for each case
00667     IF (rda_x(1)==0) THEN
00668         rl_y1=rda_x(2); rl_y2=rda_x(3); rl_y3=rda_x(4)
00669         rl_d1=1
00670     ELSE IF (rda_x(2)==0) THEN
00671         rl_y1=rda_x(1); rl_y2=rda_x(3); rl_y3=rda_x(4)
00672         rl_d2=1
00673     ELSE IF (rda_x(3)==0) THEN
00674         rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(4)
00675         rl_d3=1
00676     ELSE
00677         rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(3)
00678         rl_d4=1
00679     END IF
00680 
00681         ! Solving the system 
00682     rl_t1 = 1/rl_y1-1/rl_y2
00683     rl_t2 = 1/rl_y1**2-1/rl_y2**2
00684     rl_t3 = 1/rl_y1-1/rl_y3
00685     rl_t4 = 1/rl_y1**2-1/rl_y3**2
00686 
00687     rl_c1_y =(1/rl_y1**3/rl_t1-1/rl_y1**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
00688     rl_c2_y = -1/rl_y2**3/rl_t1/(rl_t2/rl_t1-rl_t4/rl_t3)
00689     rl_c3_y = 1/rl_y3**3/rl_t3/(rl_t2/rl_t1-rl_t4/rl_t3)
00690     rl_c4_y=(-1/rl_y1**3/rl_t1+1/rl_y2**3/rl_t1+ &
00691        1/rl_y1**3/rl_t3-1/rl_y3**3/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
00692 
00693     rl_b1_y = 1/rl_y1**3/rl_t1 - rl_c1_y*rl_t2/rl_t1
00694     rl_b2_y = -1/rl_y2**3/rl_t1 - rl_c2_y*rl_t2/rl_t1
00695     rl_b3_y = -rl_c3_y*rl_t2/rl_t1
00696     rl_b4_y = -1/rl_y1**3/rl_t1 + 1/rl_y2**3/rl_t1 - rl_c4_y*rl_t2/rl_t1
00697 
00698     rl_a1_y = 1/rl_y1**3 - rl_b1_y/rl_y1 - rl_c1_y/rl_y1**2
00699     rl_a2_y = -rl_b2_y/rl_y1 - rl_c2_y/rl_y1**2
00700     rl_a3_y = -rl_b3_y/rl_y1 - rl_c3_y/rl_y1**2
00701     rl_a4_y = -1/rl_y1**3 - rl_b4_y/rl_y1 - rl_c4_y/rl_y1**2
00702 
00703         ! Retransformation
00704     IF (rda_x(1)==0) THEN
00705         rl_a1=rl_a4_y; rl_a2=rl_a1_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
00706         rl_b1=rl_b4_y; rl_b2=rl_b1_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
00707         rl_c1=rl_c4_y; rl_c2=rl_c1_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
00708     ELSE IF (rda_x(2)==0) THEN
00709         rl_a1=rl_a1_y; rl_a2=rl_a4_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
00710         rl_b1=rl_b1_y; rl_b2=rl_b4_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
00711         rl_c1=rl_c1_y; rl_c2=rl_c4_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
00712     ELSE IF (rda_x(3)==0) THEN
00713         rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a4_y; rl_a4=rl_a3_y
00714         rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b4_y; rl_b4=rl_b3_y
00715         rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c4_y; rl_c4=rl_c3_y
00716     ELSE 
00717         rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a3_y; rl_a4=rl_a4_y
00718         rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b3_y; rl_b4=rl_b4_y
00719         rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c3_y; rl_c4=rl_c4_y
00720     END IF
00721 
00722         ! Weights  
00723     rda_wght(1) = rl_a1*rd_pt**3 + rl_b1*rd_pt**2 + rl_c1*rd_pt +rl_d1
00724     rda_wght(2) = rl_a2*rd_pt**3 + rl_b2*rd_pt**2 + rl_c2*rd_pt +rl_d2
00725     rda_wght(3) = rl_a3*rd_pt**3 + rl_b3*rd_pt**2 + rl_c3*rd_pt +rl_d3
00726     rda_wght(4) = rl_a4*rd_pt**3 + rl_b4*rd_pt**2 + rl_c4*rd_pt +rl_d4
00727 
00728     END IF
00729 
00730   END SUBROUTINE calcul_wght_irreg
00731 
00732   SUBROUTINE pole_transfo_rotation(pt_lat,pt_lon,db_pih,vareps)
00733 
00734 ! Routine given by Rene (psmile_transrot_dble) from C. Koeberle, R. Gerdes AWI Bremerhaven.
00735 ! The grid is rotated by 90 degree. The points are translated 
00736 ! by 90 deg towards the Equator. North and South pole are
00737 ! located on the Equator on the transformed grid.
00738 ! The data are in radian.
00739     !
00740   IMPLICIT NONE
00741 
00742     DOUBLE PRECISION, INTENT(INOUT)    :: pt_lat
00743     DOUBLE PRECISION, INTENT(IN)       :: pt_lon
00744     !
00745     DOUBLE PRECISION, INTENT(IN)       :: db_pih,vareps
00746     !
00747 ! Latitude
00748     !
00749     IF (ABS(pt_lon - 0.) .LE. vareps .OR. &
00750        ABS(pt_lon - 2*db_pih) .LE. vareps .OR. &
00751        ABS(pt_lon - 4*db_pih) .LE. vareps) THEN
00752         pt_lat = - COS(pt_lon) * (db_pih - pt_lat) 
00753     ELSE
00754         pt_lat = ASIN(COS( pt_lat ) * SIN( pt_lon ))
00755     ENDIF
00756     !
00757   END SUBROUTINE pole_transfo_rotation
Generated on 1 Dec 2011 for Oasis4 by  doxygen 1.6.1