prismtrs_gauswght_weight_3d.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_3d
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_gauswght_weight_3d(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                        dd_gaus_var,         &
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_Gauswght_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 ! variance of the gaussian variance
00055   DOUBLE PRECISION, INTENT (IN)                :: dd_gaus_var
00056 
00057 ! indices of the foud neighbors on the source grid for each target point
00058   INTEGER, INTENT (IN)                          :: id_nb_neighbors
00059   INTEGER, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (IN) :: ida_neighbor_index
00060 
00061 !
00062 ! ! RETURN VALUE
00063 !
00064 ! gauswght weights for 8 corners
00065   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors), INTENT (Out) :: dda_weights
00066 
00067   INTEGER, INTENT (Out)                      :: id_err   ! error value
00068 
00069 ! !DESCRIPTION
00070 ! Subroutine "PRISMTrs_Gauswght_weight_comp" computes the weights with the 
00071 ! 3D gaussian nearest neighbors interpolation method, a 3D adaptation and
00072 ! extension of the SCRIP 1.4 library 2D method.
00073 ! The inputs are the neighboring adresses computed in the PSMILe.
00074 !
00075 ! !REVISED HISTORY
00076 !   Date      Programmer   Description
00077 ! ----------  ----------   -----------
00078 ! 26/06/2003  D. Declat     Creation
00079 !
00080 ! EOP
00081 !----------------------------------------------------------------------
00082 ! $Id: prismtrs_gauswght_weight_3d.F90 2963 2011-02-17 14:52:53Z coquart $
00083 ! $Author: coquart $
00084 !----------------------------------------------------------------------
00085 !
00086 ! Local declarations
00087 !
00088   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00089      '$Id: prismtrs_gauswght_weight_3d.F90 2963 2011-02-17 14:52:53Z coquart $'
00090 
00091 ! Earth radius
00092   DOUBLE PRECISION, PARAMETER :: a = 6400.0
00093 
00094 ! loop indices
00095   INTEGER                 :: ib, ib_bis
00096 
00097 ! vectors of lon and lat of the source points used to compute the weights
00098   DOUBLE PRECISION, DIMENSION(8)      :: src_lats
00099   DOUBLE PRECISION, DIMENSION(8)      :: src_lons
00100   DOUBLE PRECISION, DIMENSION(8)      :: src_zs
00101   DOUBLE PRECISION, DIMENSION(8)      :: vec_lon
00102   DOUBLE PRECISION                    :: common_grid_range(2) 
00103   DOUBLE PRECISION                    :: shiftSize
00104 
00105 ! lat/lon coords of destination point
00106   DOUBLE PRECISION                    :: plat, plon, pz, slat, slon, sz       
00107 
00108 ! for computing gaus-weighted avg
00109   DOUBLE PRECISION                                         :: arg
00110   DOUBLE PRECISION, DIMENSION(id_tgt_size)                 :: distance
00111   DOUBLE PRECISION, DIMENSION(id_tgt_size,id_nb_neighbors) :: dist
00112 
00113 ! ---------------------------------------------------------------------
00114 !
00115 #ifdef VERBOSE
00116   PRINT *, '| | | | | | | Enter PRISMTrs_Gauswght_weight_3d'
00117   call psmile_flushstd
00118 #endif
00119 
00120   id_err = 0
00121   dda_weights(:,:) = 0.0d0
00122   common_grid_range(1) = -dble_pi
00123   common_grid_range(2) =  dble_pi
00124   shiftSize = common_grid_range(2) - common_grid_range(1)
00125 !
00126 ! 0. Treat the lat and lon arrays
00127 !
00128 ! 0.1. convert longitudes to 0,2pi interval
00129 !
00130   DO ib = 1, id_src_size
00131     DO WHILE (dda_src_lon(ib) < common_grid_range(1))
00132       dda_src_lon(ib) = dda_src_lon(ib) + shiftSize
00133     ENDDO
00134     DO WHILE (dda_src_lon(ib) >= common_grid_range(2))
00135       dda_src_lon(ib) = dda_src_lon(ib) - shiftSize
00136     ENDDO
00137   ENDDO
00138   DO ib = 1, id_tgt_size
00139     DO WHILE (dda_tgt_lon(ib) < common_grid_range(1))
00140       dda_tgt_lon(ib) = dda_tgt_lon(ib) + shiftSize
00141     ENDDO
00142     DO WHILE (dda_tgt_lon(ib) >= common_grid_range(2))
00143       dda_tgt_lon(ib) = dda_tgt_lon(ib) - shiftSize
00144     ENDDO
00145   ENDDO
00146 !
00147 ! 0.2. make sure input latitude range is within the machine values
00148 !      for +/- pi/2
00149 !  
00150   DO ib = 1, id_src_size
00151     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih)) THEN
00152         WRITE(*,*) 'Problem :Latitudes of source are greater than +/-90 degrees'
00153         CALL psmile_abort
00154     ENDIF
00155   ENDDO
00156   DO ib = 1, id_tgt_size
00157     IF ( (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00158         WRITE(*,*) 'Problem :Latitudes of target are greater than +/-90 degrees'
00159         CALL psmile_abort
00160     ENDIF
00161   ENDDO
00162 !
00163 ! In the loop
00164 !
00165   DO ib = 1, id_tgt_size 
00166 
00167     IF (ida_tgt_mask(ib) .EQ. 1) THEN
00168 
00169 ! Set the lat and lon of the target point
00170         plat = dda_tgt_lat(ib)
00171         plon = dda_tgt_lon(ib)
00172         pz   = dda_tgt_z(ib)
00173 
00174 ! For the target point ib, check that at least one neighour was found
00175     
00176 ! For the target point ib, set the lat and lon of the 8 neighours
00177             DO ib_bis = 1, id_nb_neighbors
00178 
00179               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00180 
00181                   src_zs (ib_bis) = dda_src_z (ida_neighbor_index(ib, ib_bis))
00182                   src_lats (ib_bis) = &
00183                      dda_src_lat (ida_neighbor_index(ib, ib_bis))
00184                   src_lons (ib_bis) = &
00185                      dda_src_lon (ida_neighbor_index(ib, ib_bis))
00186                   
00187                   vec_lon(ib_bis) = src_lons (ib_bis) - plon
00188                   IF (vec_lon(ib_bis) > dble_pi) THEN
00189                       src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00190                   ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00191                       src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00192 
00193                   ENDIF
00194                   
00195               ENDIF
00196             END DO
00197 
00198 !
00199 ! Compute the weights for a gauswght method
00200 !
00201         ! In a 3d cartesian coordinates system, the coodrinates of the 
00202         ! target point is
00203         !xt = (a+pz)*COS(plat)*COS(plon)
00204         !yt = (a+pz)*COS(plat)*SIN(plon)
00205         !zt = (a+pz)*SIN(plat)
00206 
00207             distance = 0.
00208             DO ib_bis = 1, id_nb_neighbors
00209 
00210               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00211               
00212                   slat = src_lats(ib_bis)
00213                   slon = src_lons(ib_bis)
00214                   sz   = src_zs(ib_bis)
00215                   
00216                   ! In a 3d cartesian coordinates system, the coordinates of the
00217                   ! source point is
00218                   !xs = (a+sz)*COS(slat)*COS(slon)
00219                   !ys = (a+sz)*COS(slat)*SIN(slon)
00220                   !zs = (a+sz)*SIN(slat)
00221                   
00222                   
00223                   ! The linear distance between the two points is
00224                   !dist = SQRT((xt-xs)**2 + (yt-ys)**2 + (zt-zs)**2)
00225                   arg = SIN(plat)*SIN(slat) + &
00226                      COS(plat)*COS(slat) * &
00227                      (COS(plon)*COS(slon) + &
00228                      SIN(plon)*SIN(slon))
00229                   IF (arg < -1.0d0) THEN
00230                       arg = -1.0d0
00231                   ELSE IF (arg > 1.0d0) THEN
00232                       arg = 1.0d0
00233                   END IF
00234                   dist(ib,ib_bis) = (a + pz)*ACOS (arg)
00235                   
00236                   !             distance = SQRT(dist**2 + (pz-sz)**2)
00237                   dist(ib,ib_bis)  = SQRT(dist(ib,ib_bis)**2 + (pz-sz)**2)
00238               ENDIF
00239             ENDDO
00240 
00241             DO ib_bis = 1, id_nb_neighbors
00242               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00243                   distance(ib) = distance(ib) + dist(ib,ib_bis)
00244               ENDIF
00245             END DO
00246               
00247               distance(ib)=distance(ib)/id_nb_neighbors
00248               distance(ib)=dd_gaus_var*distance(ib)*distance(ib)
00249               
00250             DO ib_bis = 1, id_nb_neighbors
00251               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00252                   dda_weights(ib,ib_bis) = &
00253                      EXP(.5*dist(ib,ib_bis)*dist(ib,ib_bis)/distance(ib))
00254                   dda_weights(ib,ib_bis) = 1./dda_weights(ib,ib_bis)
00255               ENDIF
00256             ENDDO
00257 ! Normalisation
00258               distance(ib) = 0.
00259               DO ib_bis = 1, id_nb_neighbors
00260                 IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00261                     distance(ib) = distance(ib) + dda_weights(ib,ib_bis)
00262                 ENDIF
00263               END DO
00264 
00265             DO ib_bis = 1, id_nb_neighbors
00266               IF (ida_neighbor_index(ib, ib_bis) > 0) THEN
00267                   dda_weights(ib,ib_bis) = dda_weights(ib,ib_bis)/distance(ib)
00268               ENDIF
00269             ENDDO
00270      
00271           ENDIF
00272       
00273         ENDDO
00274  
00275 !
00276 #ifdef VERBOSE
00277   PRINT *, '| | | | | | | Quit PRISMTrs_Gauswght_weight_3d'
00278   call psmile_flushstd
00279 #endif
00280 
00281 END SUBROUTINE PRISMTrs_Gauswght_weight_3d
00282 
00283 
00284 
00285