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