prismtrs_distwght_weight_2d1d.F90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------
00002 ! Copyright 2006-2010, CERFACS, Toulouse, France.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PRISMTrs_Distwght_weight_2d1d
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_distwght_weight_2d1d(id_src_size,         &
00011                                          dda_src_lat,         &
00012                                          dda_src_lon,         &
00013                                          dda_src_z,           &
00014                                          ida_src_mask,        &
00015                                          id_tgt_size,         &
00016                                          dda_tgt_lat,         &
00017                                          dda_tgt_lon,         &
00018                                          dda_tgt_z,           &
00019                                          ida_tgt_mask,        &
00020                                          id_nb_neighbors,     &
00021                                          ida_neighbor_index,  & 
00022                                          dda_weights,         &
00023                                          id_err)
00024 
00025 !
00026 ! !USES:
00027 !
00028   USE PRISMDrv, dummy_interface => PRISMTrs_Distwght_weight_2d1d
00029 
00030   IMPLICIT NONE
00031 
00032 !
00033 ! !PARAMETERS:
00034 !
00035 ! size of the source and target grid to interpolate
00036   INTEGER, INTENT (IN)                         :: id_src_size
00037   INTEGER, INTENT (IN)                         :: id_tgt_size
00038  
00039 ! latitudes and longitudes and z coordinates of the source grid
00040   DOUBLE PRECISION, DIMENSION(id_src_size)              :: dda_src_lat
00041   DOUBLE PRECISION, DIMENSION(id_src_size)              :: dda_src_lon
00042   DOUBLE PRECISION, DIMENSION(id_src_size)              :: dda_src_z
00043   INTEGER, DIMENSION(id_src_size), INTENT (IN) :: ida_src_mask
00044 
00045 ! z coordinates of the source grid
00046 
00047 ! latitudes and longitudes and z coordinates 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   INTEGER, DIMENSION(id_tgt_size), INTENT (IN) :: ida_tgt_mask
00052 
00053 ! indices of the found neighbors on the source grid for each target point
00054   INTEGER, INTENT (IN) :: id_nb_neighbors
00055   INTEGER, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (InOut):: ida_neighbor_index
00056 
00057 !
00058 ! ! RETURN VALUE
00059 !
00060 ! distwght weights for the neighbors
00061   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (Out) :: dda_weights
00062 
00063 ! reference of the file that gathers the computed weights
00064 !  INTEGER, INTENT (Out)               :: file_weight
00065 
00066   INTEGER, INTENT (Out)                      :: id_err   ! error value
00067 
00068 ! !DESCRIPTION
00069 ! Subroutine "PRISMTrs_Distwght_weight_2d1d" computes the weights with the 
00070 ! nearest neighbors interpolation method, derived from the SCRIP 1.4
00071 ! library available from Los Alamos National Laboratory.. 
00072 ! The inputs are the results of the search neighboring done in the PSMILe.
00073 ! The weights are stored in files which reference is given to the 
00074 ! transformer as output of the routine.
00075 !
00076 ! !REVISED HISTORY
00077 !   Date      Programmer   Description
00078 ! ----------  ----------   -----------
00079 ! 06/12/2002  D. Declat    Creation
00080 !
00081 ! EOP
00082 !----------------------------------------------------------------------
00083 ! $Id: prismtrs_distwght_weight_2d1d.F90 2963 2011-02-17 14:52:53Z coquart $
00084 ! $Author: coquart $
00085 !----------------------------------------------------------------------
00086 !
00087 ! Local declarations
00088 !
00089   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00090      '$Id: prismtrs_distwght_weight_2d1d.F90 2963 2011-02-17 14:52:53Z coquart $'
00091 
00092 ! loop indices
00093   INTEGER                 :: ib, ib_bis, il_match
00094 
00095 ! vectors of lon and lat of the source points used to compute the weights
00096   DOUBLE PRECISION, DIMENSION(id_nb_neighbors/2) :: src_lats
00097   DOUBLE PRECISION, DIMENSION(id_nb_neighbors/2) :: src_lons
00098   DOUBLE PRECISION, DIMENSION(id_nb_neighbors/2) :: vec_lon
00099   DOUBLE PRECISION                               :: common_grid_range(2)
00100   DOUBLE PRECISION                               :: shiftSize
00101 
00102 ! lat/lon coords of destination point
00103   DOUBLE PRECISION :: plat, plon       
00104 
00105 ! for computing dist-weighted avg
00106   DOUBLE PRECISION :: distance, arg
00107 
00108   DOUBLE PRECISION, PARAMETER :: converge = 0.01
00109 
00110 ! ---------------------------------------------------------------------
00111 !
00112 #ifdef VERBOSE
00113   PRINT *, '| | | | | | | Enter PRISMTrs_Distwght_weight_2d1d'
00114   call psmile_flushstd
00115 #endif
00116 
00117   id_err = 0
00118   common_grid_range(1) = -dble_pi
00119   common_grid_range(2) =  dble_pi
00120   shiftSize = common_grid_range(2) - common_grid_range(1)
00121 
00122 ! 0. Treat the lat and lon arrays
00123 !
00124 ! 0.1. convert longitudes to 0,2pi interval
00125 !
00126   DO ib = 1, id_src_size
00127     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00128       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00129     ENDDO
00130     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00131       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00132     ENDDO
00133   ENDDO
00134   DO ib = 1, id_tgt_size
00135     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00136       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00137     ENDDO
00138     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00139       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00140     ENDDO
00141   ENDDO
00142 !
00143 ! 0.2. make sure input latitude range is within the machine values
00144 !      for +/- pi/2
00145 !  
00146   DO ib = 1, id_src_size
00147     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00148         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00149         CALL psmile_abort
00150     ENDIF
00151   ENDDO
00152   DO ib = 1, id_tgt_size
00153     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00154         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00155         CALL psmile_abort
00156     ENDIF
00157   ENDDO
00158 !
00159 ! 1. In the loop
00160 !
00161   DO ib = 1, id_tgt_size
00162 
00163     IF (ida_tgt_mask(ib) .eq. 1) THEN
00164 
00165 ! 1.1. Set the lat and lon of the target point
00166         plat = dda_tgt_lat(ib)
00167         plon = dda_tgt_lon(ib)
00168     
00169 ! 1.2. For the target point ib, check that the N neighours were found down side
00170 
00171         DO ib_bis = 1, id_nb_neighbors/2
00172 
00173           IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00174     
00175 ! 1.2.1. For the target point ib, set the lat and lon of the four neighours
00176               src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00177               src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00178 
00179               vec_lon(ib_bis) = src_lons (ib_bis) - plon
00180               IF (vec_lon(ib_bis) > dble_pi) THEN
00181                   src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00182               ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00183                   src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00184               ENDIF
00185 
00186           END IF
00187 
00188         END DO
00189 
00190 ! 1.2.2. Compute the weights for a distwght method for half of the neighbors
00191 
00192         distance = 0.0d0
00193         il_match = 0
00194         DO ib_bis = 1, id_nb_neighbors/2
00195 
00196           IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00197               arg = SIN(plat)*SIN(src_lats(ib_bis))     &
00198                  +  COS(plat)*COS(src_lats(ib_bis))   &
00199                  * (COS(plon)*COS(src_lons(ib_bis))  &
00200                  +  SIN(plon)*SIN(src_lons(ib_bis)))
00201               IF (arg < -1.0d0) THEN
00202                   arg = -1.0d0
00203               ELSE IF (arg > 1.0d0) THEN
00204                   arg = 1.0d0
00205               END IF
00206               dda_weights(ib,ib_bis) = ACOS(arg) 
00207 
00208               ! dd if the source point matches with the target point,
00209               ! only consider this point.
00210               IF (dda_weights(ib,ib_bis) .eq. 0) THEN 
00211                   il_match = ib_bis 
00212               ELSE
00213                   distance = distance + 1/dda_weights(ib,ib_bis)
00214 
00215               END IF
00216           ELSE
00217               dda_weights(ib,ib_bis) = 0
00218           END IF
00219 
00220         END DO
00221 
00222         DO ib_bis = 1, id_nb_neighbors/2
00223 
00224           IF (il_match .eq. 0) THEN
00225               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00226                   dda_weights(ib,ib_bis) = 1/dda_weights(ib,ib_bis)/distance
00227               ELSE
00228                   dda_weights(ib,ib_bis) = 0
00229               END IF
00230           ELSE
00231               IF (ib_bis .eq. il_match) THEN
00232                   dda_weights(ib,ib_bis) = 1 
00233               ELSE
00234                   dda_weights(ib,ib_bis) = 0
00235               END IF  
00236           END IF
00237 
00238         END DO
00239 
00240 ! 1.3. For the target point ib, set the wieights for the up side
00241 
00242         DO ib_bis = 1, id_nb_neighbors/2
00243 
00244       dda_weights(ib,ib_bis+ id_nb_neighbors/2) = dda_weights(ib,ib_bis)
00245 
00246     END DO
00247 
00248 
00249     ELSE
00250 
00251         dda_weights(ib,:) = 0
00252 
00253     END IF
00254 
00255   END DO
00256 !
00257 
00258 !
00259 #ifdef VERBOSE
00260   PRINT *, '| | | | | | | Quit PRISMTrs_Distwght_weight_2d1d'
00261   call psmile_flushstd
00262 #endif
00263 
00264 END SUBROUTINE PRISMTrs_Distwght_weight_2d1d
00265 
00266 
00267 
00268