|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE INTERPOL_NPTS(HPROGRAM,KLUOUT,KNPTS,KCODE,PX,PY,PFIELD) 00003 ! ######################################################### 00004 ! 00005 !!**** *INTERPOL_NPTS* interpolates with ###ine f77 programs a 2D field 00006 !! from all grid points valid values 00007 !! 00008 !! PURPOSE 00009 !! ------- 00010 !! 00011 !! The points are all on only one grid (defined with the coordinates 00012 !! of all the points). The code to apply for each point is: 00013 !! 00014 !! KCODE>0 : data point (with field valid for interpolation) 00015 !! KCODE=-1: point to ignore 00016 !! KCODE=0 : point to interpolate 00017 !! 00018 !! 00019 !! 00020 !! METHOD 00021 !! ------ 00022 !! 00023 !! EXTERNAL 00024 !! -------- 00025 !! 00026 !! IMPLICIT ARGUMENTS 00027 !! ------------------ 00028 !! 00029 !! 00030 !! 00031 !! REFERENCE 00032 !! --------- 00033 !! 00034 !! AUTHOR 00035 !! ------ 00036 !! 00037 !! V. Masson Meteo-France 00038 !! 00039 !! MODIFICATION 00040 !! ------------ 00041 !! 00042 !! Original 03/2004 00043 !! Modification 00044 !---------------------------------------------------------------------------- 00045 ! 00046 !* 0. DECLARATION 00047 ! ----------- 00048 ! 00049 USE MODD_SURF_PAR, ONLY : XUNDEF 00050 USE MODD_SURF_ATM_GRID_n, ONLY : CGRID, XGRID_PAR, NGRID_PAR, NNEAR 00051 USE MODD_SURF_ATM_n, ONLY : NSIZE_FULL, NDIM_FULL 00052 ! 00053 USE MODI_GET_INTERP_HALO 00054 USE MODI_GET_NEAR_MESHES 00055 ! 00056 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00057 USE PARKIND1 ,ONLY : JPRB 00058 ! 00059 IMPLICIT NONE 00060 ! 00061 !* 0.1 Declaration of arguments 00062 ! ------------------------ 00063 ! 00064 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! host program 00065 INTEGER, INTENT(IN) :: KLUOUT ! output listing 00066 INTEGER, INTENT(IN) :: KNPTS ! number of points to interpolate with 00067 INTEGER,DIMENSION(:), INTENT(INOUT) :: KCODE ! code for each point 00068 ! >0 point used for interpolation 00069 ! 0 point to interpolate 00070 ! -1 point not used 00071 ! -2 point not used 00072 ! ! -3 if spline is no computed 00073 ! ! for this point 00074 REAL, DIMENSION(:), INTENT(IN) :: PX ! x of each grid mesh. 00075 REAL, DIMENSION(:), INTENT(IN) :: PY ! y of each grid mesh. 00076 REAL, DIMENSION(:,:),INTENT(INOUT) :: PFIELD ! pgd field on grid mesh. 00077 ! 00078 !* 0.2 Declaration of local variables 00079 ! ------------------------------ 00080 ! 00081 INTEGER :: IL ! number of points 00082 INTEGER :: JD ! data point index 00083 INTEGER :: JS ! loop counter on data points 00084 INTEGER :: JL ! loop counter on points to initialize 00085 INTEGER :: JP, JPP ! loops counter on KNPTS 00086 REAL :: ZDIST ! square distance between two interpolating and interpolated points 00087 REAL, DIMENSION(0:KNPTS) :: ZNDIST ! 3 nearest square distances 00088 REAL, DIMENSION(0:KNPTS,SIZE(PFIELD,2)) :: ZNVAL ! 3 corresponding field values 00089 REAL, DIMENSION(SIZE(PFIELD,2)) :: ZSUM 00090 ! 00091 INTEGER :: INEAR_NBR ! number of points to scan 00092 INTEGER :: JLIST ! loop counter on points to interpolate 00093 INTEGER :: ICOUNT ! counter 00094 INTEGER :: INPTS 00095 INTEGER :: ISCAN ! number of points to scan 00096 INTEGER, DIMENSION(:), ALLOCATABLE :: IINDEX ! list of index to scan 00097 INTEGER :: IHALO ! halo available 00098 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00099 !------------------------------------------------------------------------------- 00100 ! 00101 IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS',0,ZHOOK_HANDLE) 00102 IL = SIZE(PFIELD,1) 00103 ! 00104 CALL GET_INTERP_HALO(HPROGRAM,CGRID,IHALO) 00105 ! 00106 INEAR_NBR = (2*IHALO+1)**2 00107 ! 00108 ! 00109 ALLOCATE(IINDEX(IL)) 00110 IINDEX(:) = 0 00111 ! 00112 ! 00113 IF (.NOT.ASSOCIATED(NNEAR)) THEN 00114 ALLOCATE(NNEAR(IL,INEAR_NBR)) 00115 NNEAR(:,:) = 0 00116 CALL GET_NEAR_MESHES(CGRID,NGRID_PAR,NSIZE_FULL,XGRID_PAR,INEAR_NBR,NNEAR) 00117 ENDIF 00118 ! 00119 DO JL=1,IL 00120 00121 IF (KCODE(JL)/=0) CYCLE 00122 00123 ZNDIST (1:KNPTS) = 1.E20 00124 ZNDIST (0) = 0. 00125 ZNVAL(0:KNPTS,:) = 0. 00126 ! 00127 ICOUNT = 0 00128 DO JD=1,INEAR_NBR 00129 IF (NNEAR(JL,JD)>0) THEN 00130 IF (KCODE(NNEAR(JL,JD))>0) THEN 00131 ICOUNT = ICOUNT+1 00132 IINDEX(ICOUNT) = NNEAR(JL,JD) 00133 END IF 00134 END IF 00135 END DO 00136 ! 00137 IF (ICOUNT>=1) THEN 00138 ISCAN = ICOUNT 00139 INPTS = MIN(ICOUNT,KNPTS) 00140 ELSE 00141 KCODE(JL) = -4 00142 CYCLE 00143 END IF 00144 ! 00145 ! 00146 DO JS=1,ISCAN 00147 ! 00148 JD = IINDEX(JS) 00149 ! 00150 ZDIST= ( ( PX(JD)-PX(JL) ) ** 2 ) + ( ( PY(JD)-PY(JL) ) ** 2 ) 00151 ! 00152 IF ( ZDIST>ZNDIST(INPTS) ) CYCLE 00153 ! 00154 DO JP = INPTS,1,-1 00155 ! 00156 IF ( ZDIST>ZNDIST(JP-1) ) THEN 00157 ! 00158 IF ( JP<INPTS ) THEN 00159 DO JPP = INPTS,JP+1,-1 00160 ZNDIST(JPP) = ZNDIST(JPP-1) 00161 ZNVAL(JPP,:) = ZNVAL(JPP-1,:) 00162 ENDDO 00163 ENDIF 00164 ! 00165 ZNDIST(JP) = ZDIST 00166 ZNVAL(JP,:) = PFIELD(JD,:) 00167 ! 00168 EXIT 00169 ! 00170 ENDIF 00171 ! 00172 ENDDO 00173 ! 00174 ENDDO 00175 ! 00176 ZNDIST(:) = SQRT(ZNDIST(:)) 00177 ! 00178 PFIELD(JL,:) = 0. 00179 ZSUM(:) = 0. 00180 DO JP = 1, INPTS 00181 PFIELD(JL,:) = PFIELD(JL,:) + ZNVAL(JP,:)/ZNDIST(JP) 00182 ZSUM(:) = ZSUM(:) + 1./ZNDIST(JP) 00183 ENDDO 00184 PFIELD(JL,:) = PFIELD(JL,:) / ZSUM(:) 00185 ! 00186 END DO 00187 ! 00188 DEALLOCATE(IINDEX ) 00189 IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS',1,ZHOOK_HANDLE) 00190 !------------------------------------------------------------------------------- 00191 ! 00192 END SUBROUTINE INTERPOL_NPTS
1.8.0
