prismtrs_gauswght_weight_2d.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_Gauswght_weight_2d
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_gauswght_weight_2d(id_src_size,         &
00011                                        dda_src_lat,         &
00012                                        dda_src_lon,         &
00013                                        ida_src_mask,        &
00014                                        id_tgt_size,         &
00015                                        dda_tgt_lat,         &
00016                                        dda_tgt_lon,         &
00017                                        ida_tgt_mask,        &
00018                        dd_gaus_var,         &
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_Gauswght_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 ! variance of the gaussian variance
00049   DOUBLE PRECISION, INTENT (IN)                :: dd_gaus_var
00050 
00051 ! indices of the found neighbors on the source grid for each target point
00052   INTEGER, INTENT (IN) :: id_nb_neighbors
00053   INTEGER, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (INOUT) :: ida_neighbor_index
00054 
00055 !
00056 ! ! RETURN VALUE
00057 !
00058 ! gauswght weights for the neighbors
00059   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (Out) :: dda_weights
00060 
00061   INTEGER, INTENT (Out)                      :: id_err   ! error value
00062 
00063 ! !DESCRIPTION
00064 ! Subroutine "PRISMTrs_Gauswght_weight_2d" computes the weights with the 
00065 ! gaussian nearest neighbors interpolation method, adapted from the SCRIP 1.4
00066 ! library available from Los Alamos National Laboratory. 
00067 ! The inputs are the neighboring adresses computed in the PSMILe.
00068 !
00069 ! !REVISED HISTORY
00070 !   Date      Programmer   Description
00071 ! ----------  ----------   -----------
00072 ! 06/12/2002  D. Declat    Creation
00073 !
00074 ! EOP
00075 !----------------------------------------------------------------------
00076 ! $Id: prismtrs_gauswght_weight_2d.F90 2963 2011-02-17 14:52:53Z coquart $
00077 ! $Author: coquart $
00078 !----------------------------------------------------------------------
00079 !
00080 ! Local declarations
00081 !
00082   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00083      '$Id: prismtrs_gauswght_weight_2d.F90 2963 2011-02-17 14:52:53Z coquart $'
00084 
00085 ! loop indices
00086   INTEGER                 :: ib, ib_bis, nb_match
00087 
00088   INTEGER, DIMENSION(id_nb_neighbors)   :: ila_match
00089 
00090 ! vectors of lon and lat of the source points used to compute the weights
00091   DOUBLE PRECISION, DIMENSION(id_nb_neighbors)      :: dla_src_lats
00092   DOUBLE PRECISION, DIMENSION(id_nb_neighbors)      :: dla_src_lons
00093   DOUBLE PRECISION, DIMENSION(id_nb_neighbors)      :: 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       
00099 
00100 ! for computing dist-weighted avg
00101   DOUBLE PRECISION                                         :: arg
00102   DOUBLE PRECISION, PARAMETER                              :: converge = 0.01
00103   DOUBLE PRECISION, DIMENSION(id_tgt_size)                 :: distance
00104   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors) :: dist
00105 
00106 ! ---------------------------------------------------------------------
00107 !
00108 #ifdef VERBOSE
00109   PRINT *, '| | | | | | | Enter PRISMTrs_Gauswght_weight_2d'
00110   call psmile_flushstd
00111 #endif
00112 !
00113   id_err = 0
00114   dda_weights(:,:) = 0.0d0
00115   common_grid_range(1) = -dble_pi
00116   common_grid_range(2) =  dble_pi
00117   shiftSize = common_grid_range(2) - common_grid_range(1)
00118 !
00119 ! 0. Treat the lat and lon arrays
00120 !
00121 ! 0.1. convert longitudes to 0,2pi interval
00122 !
00123   DO ib = 1, id_src_size
00124     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00125       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00126     ENDDO
00127     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00128       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00129     ENDDO
00130   ENDDO
00131   DO ib = 1, id_tgt_size
00132     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00133       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00134     ENDDO
00135     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00136       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00137     ENDDO
00138   ENDDO
00139 !
00140 ! 0.2. make sure input latitude range is within the machine values
00141 !      for +/- pi/2
00142 !  
00143   DO ib = 1, id_src_size
00144     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00145         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00146         CALL psmile_abort
00147     ENDIF
00148   ENDDO
00149   DO ib = 1, id_tgt_size
00150     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00151         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00152         CALL psmile_abort
00153     ENDIF
00154   ENDDO
00155 !
00156 ! 1. In the loop
00157 !
00158   DO ib = 1, id_tgt_size
00159 
00160     IF (ida_tgt_mask(ib) .EQ. 1) THEN
00161 
00162 ! 1.1. Set the lat and lon of the target point
00163         plat = dda_tgt_lat(ib)
00164         plon = dda_tgt_lon(ib)
00165     
00166 ! 1.2. For the target point ib, check that the 4 neighours were found
00167 
00168         DO ib_bis = 1, id_nb_neighbors
00169 
00170           IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00171     
00172 ! 1.2.1. For the target point ib, set the lat and lon of the four neighours
00173               dla_src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00174               dla_src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00175 
00176               vec_lon(ib_bis) = dla_src_lons(ib_bis) - plon
00177               IF (vec_lon(ib_bis) > dble_pi) THEN
00178                   dla_src_lons (ib_bis)= dla_src_lons (ib_bis) - dble_pi2
00179               ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00180                   dla_src_lons (ib_bis) = dla_src_lons (ib_bis) + dble_pi2
00181               ENDIF
00182           END IF
00183 
00184         END DO
00185 
00186 ! 1.2.2. Compute the weights for a gauswght method
00187 
00188         nb_match = 0
00189         ila_match = 0
00190         distance = 0.
00191         DO ib_bis = 1, id_nb_neighbors
00192           IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00193 
00194               arg = SIN(plat)*SIN(dla_src_lats(ib_bis))     &
00195                  + COS(plat)*COS(dla_src_lats(ib_bis))     &
00196                  *(COS(plon)*COS(dla_src_lons(ib_bis))     &
00197                  + SIN(plon)*SIN(dla_src_lons(ib_bis)))
00198               IF (arg < -1.0d0) THEN
00199                   arg = -1.0d0
00200               ELSE IF (arg > 1.0d0) THEN
00201                   arg = 1.0d0
00202               END IF
00203               
00204               dist(ib,ib_bis) = ACOS(arg)
00205           ENDIF
00206 
00207         END DO
00208           
00209             DO ib_bis = 1, id_nb_neighbors
00210               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00211                   distance(ib) = distance(ib) + dist(ib,ib_bis)
00212               ENDIF
00213             END DO
00214 
00215             distance(ib)=distance(ib)/id_nb_neighbors
00216             distance(ib)=dd_gaus_var*distance(ib)*distance(ib)
00217 
00218             DO ib_bis = 1, id_nb_neighbors
00219               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00220                   dda_weights(ib,ib_bis) = &
00221                      EXP(.5*dist(ib,ib_bis)*dist(ib,ib_bis)/distance(ib))
00222                   dda_weights(ib,ib_bis) = 1./dda_weights(ib,ib_bis)
00223               ENDIF
00224             ENDDO
00225 ! Normalisation
00226             distance(ib) = 0.
00227             DO ib_bis = 1, id_nb_neighbors
00228               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00229                   distance(ib) = distance(ib) + dda_weights(ib,ib_bis)
00230               ENDIF
00231             END DO
00232             
00233             DO ib_bis = 1, id_nb_neighbors
00234               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00235                   dda_weights(ib,ib_bis) = dda_weights(ib,ib_bis)/distance(ib)
00236               ENDIF
00237             ENDDO
00238 
00239         END IF
00240         
00241       END DO
00242 
00243 !
00244 #ifdef VERBOSE
00245   PRINT *, '| | | | | | | Quit PRISMTrs_Gauswght_weight_2d'
00246   call psmile_flushstd
00247 #endif
00248 
00249 END SUBROUTINE PRISMTrs_Gauswght_weight_2d
00250 
00251 
00252 
00253