prismtrs_distwght_weight_2d.F90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------
00002 ! Copyright 2006-2010, CERFACS, Toulouse, France.
00003 ! Copyright 2006-2010, NEC High Performance Computing, Duesseldorf, Germany.
00004 ! All rights reserved. Use is subject to OASIS4 license terms.
00005 !-----------------------------------------------------------------------
00006 !BOP
00007 !
00008 ! !ROUTINE: PRISMTrs_Distwght_weight_2d
00009 !
00010 ! !INTERFACE
00011 subroutine prismtrs_distwght_weight_2d(id_src_size,         &
00012                                        dda_src_lat,         &
00013                                        dda_src_lon,         &
00014                                        ida_src_mask,        &
00015                                        id_tgt_size,         &
00016                                        dda_tgt_lat,         &
00017                                        dda_tgt_lon,         &
00018                                        ida_tgt_mask,        &
00019                                        id_nb_neighbors,     &
00020                                        ida_neighbor_index,  & 
00021                                        dda_weights,         &
00022                                        id_err)
00023 
00024 !
00025 ! !USES:
00026 !
00027   USE PRISMDrv, dummy_interface => PRISMTrs_Distwght_weight_2d
00028 
00029   IMPLICIT NONE
00030 
00031 !
00032 ! !PARAMETERS:
00033 !
00034 ! size of the source and target grid to interpolate
00035   INTEGER, INTENT (IN) :: id_src_size
00036   INTEGER, INTENT (IN) :: id_tgt_size
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   INTEGER, DIMENSION(id_src_size), INTENT (IN)          :: ida_src_mask
00042 
00043 ! latitudes and longitudes of the target grid
00044   DOUBLE PRECISION, DIMENSION(id_tgt_size)              :: dda_tgt_lat
00045   DOUBLE PRECISION, DIMENSION(id_tgt_size)              :: dda_tgt_lon
00046   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, INTENT (IN) :: id_nb_neighbors
00050   INTEGER, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (IN) :: ida_neighbor_index
00051 
00052 !
00053 ! ! RETURN VALUE
00054 !
00055 ! distwght weights for the neighbors
00056   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (Out) :: dda_weights
00057 
00058   INTEGER, INTENT (Out)                      :: id_err   ! error value
00059 
00060 ! !DESCRIPTION
00061 ! Subroutine "PRISMTrs_Distwght_weight_2d" computes the weights with the 
00062 ! nearest neighbors interpolation method, derived from the SCRIP 1.4
00063 ! library available from Los Alamos National Laboratory. 
00064 ! The inputs are the neighboring adresses computed in the PSMILe.
00065 !
00066 ! !REVISED HISTORY
00067 !   Date      Programmer   Description
00068 ! ----------  ----------   -----------
00069 ! 06/12/2002  D. Declat    Creation
00070 !
00071 ! EOP
00072 !----------------------------------------------------------------------
00073 ! $Id: prismtrs_distwght_weight_2d.F90 3241 2011-06-23 09:55:11Z coquart $
00074 ! $Author: coquart $
00075 !----------------------------------------------------------------------
00076 !
00077 ! Local declarations
00078 !
00079   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00080      '$Id: prismtrs_distwght_weight_2d.F90 3241 2011-06-23 09:55:11Z coquart $'
00081 
00082 ! loop indices
00083   INTEGER                 :: i, ib, nb
00084 
00085 ! for computing dist-weighted avg
00086   DOUBLE PRECISION                                :: rdistance, dist, arg
00087 
00088   DOUBLE PRECISION                  :: common_grid_range(2)
00089   DOUBLE PRECISION                  :: shiftSize
00090 
00091 #ifndef SX
00092   DOUBLE PRECISION, DIMENSION(id_nb_neighbors)    :: dist2v
00093   DOUBLE PRECISION                                :: plat, plon
00094   DOUBLE PRECISION                                :: slat, slon, vec_lon
00095 #else
00096   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors) :: dist2v 
00097   DOUBLE PRECISION, DIMENSION(id_tgt_size)        :: plon, src_lons, vec_lon
00098   DOUBLE PRECISION, DIMENSION(id_tgt_size)        :: splat, splon, cplon, cplat
00099   DOUBLE PRECISION, DIMENSION(id_src_size)        :: sslat, sslon, cslon, cslat
00100   INTEGER                                         :: flag(id_tgt_size)
00101 #endif
00102 
00103 ! ---------------------------------------------------------------------
00104 !
00105 #ifdef VERBOSE
00106   PRINT *, '| | | | | | | Enter PRISMTrs_Distwght_weight_2d'
00107   call psmile_flushstd
00108 #endif
00109 
00110   id_err = 0
00111   dda_weights(:,:) = 0.0d0
00112   common_grid_range(1) = 0.
00113   common_grid_range(2) =  dble_pi2
00114   shiftSize = common_grid_range(2) - common_grid_range(1)
00115 
00116 ! 0. Treat the lat and lon arrays
00117 !
00118 ! 0.1. convert longitudes to 0,2pi interval
00119 !
00120   DO ib = 1, id_src_size
00121     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00122       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00123     ENDDO
00124     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00125       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00126     ENDDO
00127   ENDDO
00128   DO ib = 1, id_tgt_size
00129     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00130       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00131     ENDDO
00132     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00133       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00134     ENDDO
00135   ENDDO
00136 !
00137 ! 0.2. make sure input latitude range is within the machine values
00138 !      for +/- pi/2
00139 !  
00140   DO ib = 1, id_src_size
00141     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00142         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00143         CALL psmile_abort
00144     ENDIF
00145   ENDDO
00146   DO ib = 1, id_tgt_size
00147     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00148         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00149         CALL psmile_abort
00150     ENDIF
00151   ENDDO
00152 
00153 #ifndef SX
00154 ! Hubert: Sin and Cos values of all values should be computed once
00155 !         and used within this loop instead of computing the sin
00156 !         cos values multiple times.
00157 
00158 ! Do ib = 1, id_tgt_size
00159 !    IF (ida_tgt_mask(ib) .eq. 1) dda_weights(ib,:) = 0
00160 ! end do
00161 
00162   DO ib = 1, id_tgt_size
00163 
00164      IF (ida_tgt_mask(ib) .eq. 1) THEN
00165 
00166         ! For the target point ib, check that at least one neighour was found
00167 
00168            do nb = 1, id_nb_neighbors
00169            if (ida_neighbor_index(ib, nb) <= 0) exit
00170 
00171            end do
00172         nb = nb - 1
00173 
00174 ! Hubert: Does it make sense to compute this weights if nb < 8
00175 
00176         IF (nb > 0) THEN
00177 
00178            ! Set the lat and lon of the target point
00179 
00180            plat = dda_tgt_lat(ib)
00181            plon = dda_tgt_lon(ib)
00182 
00183            !
00184            ! Compute the weights for a distwght method in case we have a neighbor point.
00185            !
00186            ! In a 3d cartesian coordinates system, the coodrinates of the
00187            ! target point is
00188            !xt = COS(plat)*COS(plon)
00189            !yt = COS(plat)*SIN(plon)
00190            !
00191 
00192            DO i = 1, nb
00193               slat = dda_src_lat (ida_neighbor_index(ib, i))
00194               slon = dda_src_lon (ida_neighbor_index(ib, i))
00195 
00196               vec_lon = slon - plon
00197               IF (vec_lon > dble_pi) THEN
00198                   slon = slon - dble_pi2
00199               ELSE IF (vec_lon < -dble_pi) THEN
00200                   slon = slon + dble_pi2
00201               ENDIF
00202 
00203               ! In a 3d cartesian coordinates system, the coordinates of the
00204               ! source point is
00205               ! xs = COS(slat)*COS(slon)
00206               ! ys = COS(slat)*SIN(slon)
00207               ! zs = SIN(slat)
00208 
00209               ! The linear distance between the two points is
00210               ! dist = SQRT((xt-xs)**2 + (yt-ys)**2)
00211 
00212               arg = SIN(plat)*SIN(slat) + &
00213                     COS(plat)*COS(slat) * &
00214                    (COS(plon)*COS(slon) + &
00215                     SIN(plon)*SIN(slon))
00216               IF (arg < -1.0d0) THEN
00217                   arg = -1.0d0
00218               ELSE IF (arg > 1.0d0) THEN
00219                   arg = 1.0d0
00220               END IF
00221               dist = ACOS(arg)
00222 
00223 !             distance = SQRT(dist**2)
00224               dist2v(i) = SQRT(dist**2)
00225            END DO
00226 
00227 !          Hubert: A tolerance instead of 0.0d0 would be better
00228 
00229 #ifdef DEBUG
00230         if (ib == 1) print *, 'dist2v', dist2v (1:nb)
00231 #endif
00232            do i = 1, nb
00233            if (dist2v (i) .eq. 0.0d0) exit
00234            end do
00235 
00236            if (i <= nb) then
00237 !             ... Point to be interpolated is located on Point i
00238               dda_weights(ib,i) = 1.0d0
00239            else
00240            
00241               dda_weights(ib,1:nb) = 1.0d0 / dist2v(1:nb)
00242               rdistance = 1.0d0 / SUM (dda_weights(ib,1:nb))
00243               dda_weights(ib,1:nb) = dda_weights(ib,1:nb) * rdistance
00244            endif
00245 
00246         END IF
00247 
00248      END IF
00249 
00250   END DO
00251 
00252 #else
00253 
00254   flag(:) = 0
00255 
00256   DO ib = 1, id_tgt_size 
00257 
00258      IF (ida_tgt_mask(ib) .EQ. 1) THEN
00259         ! For the target point ib, check that at least one neighour was found
00260         DO nb = 1, id_nb_neighbors
00261            IF (ida_neighbor_index(ib, nb) > 0) flag(ib)= flag(ib)+1
00262         END DO
00263      END IF
00264 
00265   END DO
00266 
00267   !
00268   ! Set the lat and lon of the source points
00269   !
00270   DO ib = 1, id_src_size
00271     src_lons(ib) = dda_src_lon(ib)
00272     cslat(ib) = cos(dda_src_lat(ib))
00273     sslat(ib) = sin(dda_src_lat(ib))
00274   END DO
00275   !
00276   ! Set the lat and lon of the target points
00277   !
00278   DO ib = 1, id_tgt_size
00279     plon(ib) = dda_tgt_lon(ib)
00280     cplat(ib) = cos(dda_tgt_lat(ib))
00281     cplon(ib) = cos(dda_tgt_lon(ib))
00282     splat(ib) = sin(dda_tgt_lat(ib))
00283     splon(ib) = sin(dda_tgt_lon(ib))
00284   END DO
00285 
00286   DO nb = 1, id_nb_neighbors
00287      DO ib = 1, id_tgt_size
00288         !
00289         ! Set the lat and lon of the target point
00290         !
00291         !
00292         ! Compute the weights for a distwght method in case we have a neighbor point.
00293         !
00294         ! In a cartesian coordinates system, the coordinates of the 
00295         ! target point are 
00296         !
00297         !   xt =  COS(plat)*COS(plon)
00298         !   yt =  COS(plat)*SIN(plon)
00299         !
00300         ! We compute something here even if we have not found any neighbors
00301         ! just for the purpose of better vectorisation. Those non-valid points
00302         ! will be irgnored for the calculation of weights later.
00303         !
00304         IF ( ida_neighbor_index(ib, nb) > 0 ) THEN
00305 
00306             !
00307             vec_lon(ib) = src_lons (ida_neighbor_index(ib, nb)) - plon(ib)
00308             IF (vec_lon(ib) > dble_pi) THEN
00309                 src_lons (ida_neighbor_index(ib, nb))= src_lons (ida_neighbor_index(ib, nb)) - dble_pi2
00310              ELSE IF (vec_lon(ib) < -dble_pi) THEN
00311                  src_lons (ida_neighbor_index(ib, nb)) = src_lons (ida_neighbor_index(ib, nb)) + dble_pi2
00312              ENDIF
00313 
00314              cslon(ida_neighbor_index(ib, nb)) = COS(src_lons (ida_neighbor_index(ib, nb)))
00315              sslon(ida_neighbor_index(ib, nb)) = SIN(src_lons (ida_neighbor_index(ib, nb)))
00316 
00317            ! In a 3d cartesian coordinates system, the coordinates of the
00318            ! source point is
00319            !
00320            !   xs = COS(slat)*COS(slon)
00321            !   ys = COS(slat)*SIN(slon)
00322 
00323            ! The linear distance between the two points is
00324            !   dist = SQRT((xt-xs)**2 + (yt-ys)**2)
00325 
00326            arg = splat(ib)*sslat(ida_neighbor_index(ib, nb)) + &
00327                  cplat(ib)*cslat(ida_neighbor_index(ib, nb)) * &
00328                 (cplon(ib)*cslon(ida_neighbor_index(ib, nb)) + &
00329                  splon(ib)*sslon(ida_neighbor_index(ib, nb)))
00330            IF (arg < -1.0d0) THEN
00331                arg = -1.0d0
00332            ELSE IF (arg > 1.0d0) THEN
00333                arg = 1.0d0
00334            END IF
00335            dist = ACOS(arg)
00336 
00337            dist2v(ib,nb) = ABS(dist**2)
00338         END IF
00339 
00340      ENDDO
00341   ENDDO
00342 
00343   DO ib = 1, id_tgt_size
00344      !
00345      ! Hubert: Does it make sense to compute this weights if flag(ib) < id_nb_neighbors
00346      ! Hubert: A tolerance instead of 0.0d0 would be better
00347      !
00348      DO nb = 1, flag(ib)
00349         IF ( dist2v(ib,nb) == 0.0d0 ) exit
00350      end do
00351 
00352      if (nb <= flag(ib)) then
00353         dda_weights(ib,nb) = 1.0d0
00354      else
00355         dda_weights(ib,1:flag(ib)) = 1.0d0 / dist2v(ib,1:flag(ib))
00356      endif
00357   END DO
00358         
00359   DO ib = 1, id_tgt_size
00360      IF ( flag(ib) > 0 ) THEN
00361         rdistance = 1.0d0 / SUM (dda_weights(ib,1:id_nb_neighbors))
00362         dda_weights(ib,1:id_nb_neighbors) = dda_weights(ib,1:id_nb_neighbors) * rdistance
00363      ENDIF
00364   END DO
00365 #endif
00366 !
00367 #ifdef VERBOSE
00368   PRINT *, '| | | | | | | Quit PRISMTrs_Distwght_weight_2d'
00369   call psmile_flushstd
00370 #endif
00371 
00372 END SUBROUTINE PRISMTrs_Distwght_weight_2d
00373 
00374 
00375 
00376