prismtrs_bilinear_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_Bilinear_weight_2d1d
00008 !
00009 ! !INTERFACE
00010 subroutine prismtrs_bilinear_weight_2d1d(id_src_size,         &
00011                                          dda_src_lat,         &
00012                                          dda_src_lon,         &
00013                                          dda_src_z,           &
00014                                          id_tgt_size,         &
00015                                          dda_tgt_lat,         &
00016                                          dda_tgt_lon,         &
00017                                          dda_tgt_z,           &
00018                                          ida_neighbor_index,  & 
00019                                          dda_weights,         &
00020                                          id_err)
00021 
00022 !
00023 ! !USES:
00024 !
00025   USE PRISMDrv, dummy_interface => PRISMTrs_Bilinear_weight_2d1d
00026 
00027   IMPLICIT NONE
00028 
00029 !
00030 ! !PARAMETERS:
00031 !
00032 ! size of the source and target grid to interpolate
00033   INTEGER, INTENT (IN)                         :: id_src_size
00034   INTEGER, INTENT (IN)                         :: id_tgt_size
00035   
00036 ! latitudes and longitudes and z coordinates of the source grid
00037   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lat
00038   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_lon
00039   DOUBLE PRECISION, DIMENSION(id_src_size)     :: dda_src_z
00040 
00041 ! latitudes and longitudes and z coordinates of the target grid
00042   DOUBLE PRECISION, DIMENSION(id_tgt_size)     :: dda_tgt_lat
00043   DOUBLE PRECISION, DIMENSION(id_tgt_size)     :: dda_tgt_lon
00044   DOUBLE PRECISION, DIMENSION(id_tgt_size)     :: dda_tgt_z
00045 
00046 ! indices of the foud neighbors on the source grid for each target point
00047   INTEGER, DIMENSION(id_tgt_size,8), INTENT(INOUT) :: ida_neighbor_index
00048 
00049 !
00050 ! ! RETURN VALUE
00051 !
00052 ! bilinear dda_weights for four corners
00053   DOUBLE PRECISION, DIMENSION(id_tgt_size,8), INTENT (Out) :: dda_weights
00054 
00055   INTEGER, INTENT (Out)                        :: id_err   ! error value
00056 
00057 ! !DESCRIPTION
00058 ! Subroutine "PRISMTrs_Bilinear_weight_2d1d" computes the weights with 
00059 ! the bilinear interpolation method, derived from the SCRIP 1.4
00060 ! library available from Los Alamos National Laboratory.
00061 ! The inputs are the neighboring adresses computed in the PSMILe.
00062 !
00063 ! !REVISED HISTORY
00064 !   Date      Programmer   Description
00065 ! ----------  ----------   -----------
00066 ! 06/12/2002  D. Declat    Creation
00067 !
00068 ! EOP
00069 !----------------------------------------------------------------------
00070 ! $Id: prismtrs_bilinear_weight_2d1d.F90 2963 2011-02-17 14:52:53Z coquart $
00071 ! $Author: coquart $
00072 !----------------------------------------------------------------------
00073 !
00074 ! Local declarations
00075 !
00076   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00077      '$Id: prismtrs_bilinear_weight_2d1d.F90 2963 2011-02-17 14:52:53Z coquart $'
00078 
00079 ! loop indices
00080   INTEGER                 :: ib, ib_bis, iter, srch_add
00081   DOUBLE PRECISION, PARAMETER         :: converge = 0.01
00082 
00083 ! vectors of lon and lat of the source points used to compute the weights
00084   DOUBLE PRECISION, DIMENSION(4)      :: src_lats
00085   DOUBLE PRECISION, DIMENSION(4)      :: src_lons
00086   DOUBLE PRECISION, DIMENSION(4)      :: vec_lon
00087   DOUBLE PRECISION                    :: common_grid_range(2)
00088   DOUBLE PRECISION                    :: shiftSize
00089 
00090 ! lat/lon coords of destination point
00091   DOUBLE PRECISION :: plat, plon       
00092 ! current guess for bilinear coordinate
00093   DOUBLE PRECISION :: iguess, jguess
00094 ! corrections to i,j
00095   DOUBLE PRECISION :: deli, delj
00096 ! some latitude  differences
00097   DOUBLE PRECISION :: dth1, dth2, dth3
00098 ! some longitude differences
00099   DOUBLE PRECISION :: dph1, dph2, dph3
00100  ! difference between point and sw corner
00101   DOUBLE PRECISION :: dthp, dphp
00102 ! matrix elements
00103   DOUBLE PRECISION :: mat1, mat2, mat3, mat4
00104 ! matrix determinant
00105   DOUBLE PRECISION :: determinant
00106 
00107 ! for computing dist-weighted avg
00108   DOUBLE PRECISION :: distance, dist_tot, arg
00109 !
00110 ! ---------------------------------------------------------------------
00111 !
00112 #ifdef VERBOSE
00113   PRINT *, '| | | | | | | Enter PRISMTrs_Bilinear_weight_2d1d'
00114   call psmile_flushstd
00115 #endif
00116 
00117 !
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_tgt_size
00147     IF ( (dda_src_lat(ib) >  dble_pih) .OR. (dda_src_lat(ib) < -dble_pih) .OR. &
00148        (dda_tgt_lat(ib) >  dble_pih) .OR. (dda_tgt_lat(ib) < -dble_pih) ) THEN
00149         WRITE(*,*) 'Problem :Latitudes are greater than +/-90 degrees'
00150         CALL psmile_abort
00151     ENDIF
00152   ENDDO
00153 ! 1. In the loop
00154 !
00155   DO ib = 1, id_tgt_size
00156 
00157 ! 1.1. Set the lat and lon of the target point
00158     plat = dda_tgt_lat(ib)
00159     plon = dda_tgt_lon(ib)
00160     
00161 ! 1.2. For the target point ib, check that the 4 neighours were found down side
00162     IF (ida_neighbor_index(ib, 4) > 0) THEN
00163     
00164 ! 1.2.1. For the target point ib, set the lat and lon of the four neighours
00165         src_lats(1:4) = 0
00166         src_lons(1:4) = 0
00167         DO ib_bis = 1, 4
00168           src_lats (ib_bis) = dda_src_lat (ida_neighbor_index(ib, ib_bis))
00169           src_lons (ib_bis) = dda_src_lon (ida_neighbor_index(ib, ib_bis))
00170         END DO
00171 
00172         DO ib_bis = 1, 4
00173           vec_lon(ib_bis) = src_lons (ib_bis) - plon
00174           IF (vec_lon(ib_bis) > dble_pi) THEN
00175               src_lons (ib_bis)= src_lons (ib_bis) - dble_pi2
00176           ELSE IF (vec_lon(ib_bis) < -dble_pi) THEN
00177               src_lons (ib_bis) = src_lons (ib_bis) + dble_pi2
00178           ENDIF
00179         ENDDO
00180 
00181 ! 1.2.2. Compute the weights for a bilinear method
00182 
00183         dth1 = src_lats(2) - src_lats(1)
00184         dth2 = src_lats(4) - src_lats(1)
00185         dth3 = src_lats(3) - src_lats(2) - dth2
00186 
00187         dph1 = src_lons(2) - src_lons(1)
00188         dph2 = src_lons(4) - src_lons(1)
00189         dph3 = src_lons(3) - src_lons(2)
00190 
00191         IF (dph1 >  three*dble_pih) dph1 = dph1 - dble_pi2
00192         IF (dph2 >  three*dble_pih) dph2 = dph2 - dble_pi2
00193         IF (dph3 >  three*dble_pih) dph3 = dph3 - dble_pi2
00194         IF (dph1 < -three*dble_pih) dph1 = dph1 + dble_pi2
00195         IF (dph2 < -three*dble_pih) dph2 = dph2 + dble_pi2
00196         IF (dph3 < -three*dble_pih) dph3 = dph3 + dble_pi2
00197 
00198         dph3 = dph3 - dph2
00199 
00200         iguess = half
00201         jguess = half
00202 
00203         iter_loop1: DO iter= 1, PSMILe_trans_Max_iter
00204 
00205           dthp = plat - src_lats(1) - dth1*iguess -   &
00206              dth2*jguess - dth3*iguess*jguess
00207           dphp = plon - src_lons(1)
00208 
00209           IF (dphp >  three*dble_pih) dphp = dphp - dble_pi2
00210           IF (dphp < -three*dble_pih) dphp = dphp + dble_pi2
00211 
00212           dphp = dphp - dph1*iguess - dph2*jguess -   &
00213              dph3*iguess*jguess
00214 
00215           mat1 = dth1 + dth3*jguess
00216           mat2 = dth2 + dth3*iguess
00217           mat3 = dph1 + dph3*jguess
00218           mat4 = dph2 + dph3*iguess
00219 
00220           determinant = mat1*mat4 - mat2*mat3
00221 
00222           deli = (dthp*mat4 - mat2*dphp)/determinant
00223           delj = (mat1*dphp - dthp*mat3)/determinant
00224 
00225           IF (ABS(deli) < converge .and. &
00226              ABS(delj) < converge) EXIT iter_loop1
00227 
00228           iguess = iguess + deli
00229           jguess = jguess + delj
00230 
00231         END DO iter_loop1
00232 
00233         IF (iter <= PSMILe_trans_Max_iter) THEN
00234 
00235             !***
00236             !*** successfully found i,j - compute weights
00237             !***
00238 
00239             dda_weights(ib,1) = abs((one-iguess)*(one-jguess))
00240             dda_weights(ib,2) = abs(iguess*(one-jguess))
00241             dda_weights(ib,3) = abs(iguess*jguess)
00242             dda_weights(ib,4) = abs((one-iguess)*jguess)
00243 
00244         ELSE
00245 
00246             PRINT *, '| | | | | | | | Point coords: ',plat,plon
00247             PRINT *, '| | | | | | | | Dest grid lats: ',src_lats
00248             PRINT *, '| | | | | | | | Dest grid lons: ',src_lons
00249             PRINT *, '| | | | | | | | Current i,j : ',iguess, jguess
00250             call psmile_flushstd
00251             id_err = 28
00252             PRINT *, '| | | | | | | | | Iteration for i,j exceed max iteration count'
00253             call psmile_abort
00254         ENDIF
00255 
00256     ELSE IF ((ida_neighbor_index(ib,1)>0) .or. (ida_neighbor_index(ib,2)>0) &
00257        .or. (ida_neighbor_index(ib,3)>0)) THEN
00258 
00259         dist_tot = 0
00260         DO ib_bis = 1, 3
00261           
00262           IF (ida_neighbor_index(ib, ib_bis) .eq. 0) THEN
00263               dda_weights(ib,ib_bis) = 0
00264           ELSE
00265               srch_add = ida_neighbor_index(ib, ib_bis)
00266               arg = COS(plat)*COS(dda_src_lat(srch_add))  *  &
00267                    (COS(plon)*COS(dda_src_lon(srch_add))  +  &
00268                     SIN(plon)*SIN(dda_src_lon(srch_add))) +  &
00269                     SIN(plat)*SIN(dda_src_lat(srch_add))
00270               IF (arg < -1.0d0) THEN
00271                   arg = -1.0d0
00272               ELSE IF (arg > 1.0d0) THEN
00273                   arg = 1.0d0
00274               END IF
00275               distance=ACOS(arg)
00276               dda_weights(ib,ib_bis) = 1/distance
00277               dist_tot = dist_tot + 1/distance
00278           END IF
00279 
00280         END DO
00281 
00282         DO ib_bis = 1, 3
00283           IF (dda_weights(ib,ib_bis) .ne. 0) THEN
00284               dda_weights(ib,ib_bis) =  dda_weights(ib,ib_bis)/dist_tot
00285           END IF
00286         END DO
00287 
00288     END IF
00289 
00290 
00291 ! 1.3. For the target point ib, set the wieights for the up side
00292     DO ib_bis = 1, 4
00293       dda_weights(ib,ib_bis+4) = dda_weights(ib,ib_bis)
00294     END DO
00295   
00296   END DO
00297 !
00298 #ifdef VERBOSE
00299   PRINT *, '| | | | | | | Quit PRISMTrs_Bilinear_weight_2d1d'
00300   call psmile_flushstd
00301 #endif
00302   
00303 END SUBROUTINE PRISMTrs_Bilinear_weight_2d1d
00304 
00305 
00306 
00307 
00308