|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE BILIN (KLUOUT,PX1,PY1,PFIELD1,PX2,PY2,PFIELD2,OINTERP) 00003 ! ######################################################################### 00004 ! 00005 !!**** *BILINEAR * - subroutine to interpolate surface FIELD 00006 !! 00007 !! PURPOSE 00008 !! ------- 00009 !! 00010 !! 00011 !!** METHOD 00012 !! ------ 00013 !! 00014 !! Interpolation is bilinear, and uses 9 grid points, located in the 00015 !! center of model 1 grid mesh, and at the boundaries of this grid 00016 !! mesh (2 X limits, 2 Y limits and 4 corners). 00017 !! This implies that the grid mesh values located around the model 1 00018 !! grid mesh are not used directly. The values at the boundaries of the 00019 !! grid mesh are defined by the average between the middle point 00020 !! (this grid mesh value), and the one in the considered direction. 00021 !! So the eight grid meshes around the considered grid mesh are used 00022 !! equally. 00023 !! This is important to note that these average values are erased 00024 !! and replaced by zero if they are at the limit of any grid 00025 !! mesh with the zero value. This allows to insure zero value in model 2 00026 !! grid meshes where there was not the considered class in corresponding 00027 !! model 1 grid mesh, and to insure continuity of the FIELD type 00028 !! at such boundaries. 00029 !! 00030 !! 00031 !! The arrays and array index are defined on the following (model1) grid: 00032 !! 00033 !! 00034 !! XFIELD XFIELD XFIELD 00035 !! * * * 00036 !! i-1,j+1 i,j+1 i+1,j+1 00037 !! 00038 !! 00039 !! 00040 !! ZFIELD_XY ZFIELD_Y ZFIELD_XY 00041 !! * * * 00042 !! i,j+1 i,j+1 i+1,j+1 00043 !! 00044 !! 00045 !! 00046 !! XFIELD ZFIELD_X XFIELD ZFIELD_X XFIELD 00047 !! * * * * * 00048 !! i-1,j i,j i,j i+1,j i+1,j 00049 !! 00050 !! 00051 !! 00052 !! ZFIELD_XY ZFIELD_Y ZFIELD_XY 00053 !! * * * 00054 !! i,j i,j i+1,j 00055 !! 00056 !! 00057 !! 00058 !! XFIELD XFIELD XFIELD 00059 !! * * * 00060 !! i-1,j-1 i,j-1 i+1,j-1 00061 !! 00062 !! 00063 !! 00064 !! 00065 !! 00066 !! AUTHOR 00067 !! ------ 00068 !! 00069 !! V. Masson * METEO-FRANCE * 00070 !! 00071 !! MODIFICATIONS 00072 !! ------------- 00073 !! 00074 !! Original 01/2004 00075 !------------------------------------------------------------------------------- 00076 ! 00077 !* 0. DECLARATIONS 00078 ! ------------ 00079 ! 00080 USE MODI_HOR_EXTRAPOL_SURF 00081 ! 00082 USE MODD_SURF_PAR, ONLY : XUNDEF 00083 ! 00084 ! 00085 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00086 USE PARKIND1 ,ONLY : JPRB 00087 ! 00088 IMPLICIT NONE 00089 ! 00090 !* 0.1 Declarations of dummy arguments : 00091 ! 00092 INTEGER, INTENT(IN) :: KLUOUT ! output listing logical unit 00093 REAL, DIMENSION(:), INTENT(IN) :: PX1 ! X coordinate of the regular input grid 00094 REAL, DIMENSION(:), INTENT(IN) :: PY1 ! Y coordinate of the regular input grid 00095 REAL, DIMENSION(:,:), INTENT(IN) :: PFIELD1 ! FIELD on regular input grid 00096 REAL, DIMENSION(:), INTENT(IN) :: PX2 ! X coordinate of all points of output grid 00097 REAL, DIMENSION(:), INTENT(IN) :: PY2 ! Y coordinate of all points of output grid 00098 REAL, DIMENSION(:), INTENT(OUT) :: PFIELD2 ! FIELD on model 2 00099 LOGICAL, DIMENSION(:),INTENT(IN) :: OINTERP ! .true. where physical value is needed 00100 ! 00101 ! 00102 !* 0.2 Declarations of local variables for print on FM file 00103 ! 00104 ! 00105 REAL, DIMENSION (:,:), ALLOCATABLE :: ZW ! weight. 1 if defined, 00106 ! 0 if FIELD=XUNDEF 00107 REAL, DIMENSION (:,:), ALLOCATABLE :: ZW_FIELD ! weight. * FIELD 00108 REAL, DIMENSION (:,:), ALLOCATABLE :: ZFIELD_X ! FIELD at mesh interface 00109 REAL, DIMENSION (:,:), ALLOCATABLE :: ZFIELD_Y ! FIELD at mesh interface 00110 REAL, DIMENSION (:,:), ALLOCATABLE :: ZFIELD_XY! FIELD at mesh corner 00111 REAL, DIMENSION (SIZE(PX1)+1) :: ZX ! X coordinate of left limit of input meshes 00112 REAL, DIMENSION (SIZE(PY1)+1) :: ZY ! Y coordinate of bottom limit of input meshes 00113 REAL, DIMENSION (:), ALLOCATABLE :: ZX1 ! input X coordinate of all points 00114 REAL, DIMENSION (:), ALLOCATABLE :: ZY1 ! input Y coordinate of all points 00115 REAL, DIMENSION (:), ALLOCATABLE :: ZFIELD1 ! input field of all points 00116 ! 00117 REAL :: ZC1_X ! coefficient for left points 00118 REAL :: ZC2_X ! coefficient for middle points 00119 REAL :: ZC3_X ! coefficient for right points 00120 REAL :: ZC1_Y ! coefficient for bottom points 00121 REAL :: ZC2_Y ! coefficient for middle points 00122 REAL :: ZC3_Y ! coefficient for top points 00123 ! 00124 INTEGER :: IIU ! model 1 X size 00125 INTEGER :: IJU ! model 1 Y size 00126 ! 00127 INTEGER :: JI ! grid 1 x index 00128 INTEGER :: JJ ! grid 1 y index 00129 INTEGER :: IIN ! loop counter on all input points 00130 ! 00131 INTEGER :: JL ! grid 2 index 00132 ! 00133 REAL :: ZEPS=1.E-3 00134 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00135 !------------------------------------------------------------------------------- 00136 ! 00137 IF (LHOOK) CALL DR_HOOK('BILIN',0,ZHOOK_HANDLE) 00138 IIU=SIZE(PFIELD1,1) 00139 IJU=SIZE(PFIELD1,2) 00140 ! 00141 !* weighting factor 00142 ! 00143 ALLOCATE(ZW(IIU,IJU)) 00144 WHERE (PFIELD1/=XUNDEF) 00145 ZW=1. 00146 ELSEWHERE 00147 ZW=0. 00148 END WHERE 00149 ! 00150 !* weighted FIELD 00151 ! 00152 ALLOCATE(ZW_FIELD(IIU,IJU)) 00153 ZW_FIELD=ZW*PFIELD1 00154 ! 00155 !------------------------------------------------------------------------------- 00156 ! 00157 !* 1. FIELD type at grid mesh interfaces (in X directions) 00158 ! ---------------------------------- 00159 ! 00160 ALLOCATE(ZFIELD_X(IIU+1,IJU)) 00161 ! 00162 !* 1.1 Standard case 00163 ! ------------- 00164 ! 00165 IF (IIU>1) & 00166 ZFIELD_X(2:IIU ,:) = (ZW_FIELD(1:IIU-1,:)+ZW_FIELD(2:IIU,:)) & 00167 / MAX(ZW (1:IIU-1,:)+ZW (2:IIU,:),ZEPS) 00168 ZFIELD_X(1 ,:) = ZW_FIELD(1 ,:) 00169 ZFIELD_X( IIU+1,:) = ZW_FIELD(IIU,:) 00170 ! 00171 ! 00172 !* 1.2 FIELD type value is 0 in the grid mesh 00173 ! -------------------------------------- 00174 ! 00175 WHERE ( PFIELD1(1:IIU,:)==0.) 00176 ZFIELD_X(1:IIU,:) = 0. 00177 END WHERE 00178 ! 00179 WHERE ( PFIELD1(1:IIU,:)==0.) 00180 ZFIELD_X(2:IIU+1,:) = 0. 00181 END WHERE 00182 ! 00183 !------------------------------------------------------------------------------- 00184 ! 00185 !* 2. FIELD type at grid mesh interfaces (in X directions) 00186 ! ---------------------------------- 00187 ! 00188 ALLOCATE(ZFIELD_Y(IIU,IJU+1)) 00189 ! 00190 !* 2.1 Standard case 00191 ! ------------- 00192 ! 00193 IF (IJU>1) & 00194 ZFIELD_Y(:,2:IJU ) = (ZW_FIELD(:,1:IJU-1)+ZW_FIELD(:,2:IJU)) & 00195 / MAX(ZW (:,1:IJU-1)+ZW (:,2:IJU),ZEPS) 00196 ZFIELD_Y(:,1 ) = ZW_FIELD(:,1 ) 00197 ZFIELD_Y(:, IJU+1) = ZW_FIELD(:, IJU ) 00198 ! 00199 ! 00200 !* 2.3 FIELD type value is 0 in the grid mesh 00201 ! -------------------------------------- 00202 ! 00203 WHERE ( PFIELD1(:,1:IJU)==0.) 00204 ZFIELD_Y(:,1:IJU) = 0. 00205 END WHERE 00206 ! 00207 WHERE ( PFIELD1(:,1:IJU)==0.) 00208 ZFIELD_Y(:,2:IJU+1) = 0. 00209 END WHERE 00210 ! 00211 !------------------------------------------------------------------------------- 00212 ! 00213 !* 3. FIELD type at grid mesh corners 00214 ! ------------------------------- 00215 ! 00216 ALLOCATE(ZFIELD_XY(IIU+1,IJU+1)) 00217 ! 00218 !* 3.1 Standard case 00219 ! ------------- 00220 ! 00221 00222 IF (IIU>1 .AND. IJU>1) & 00223 ZFIELD_XY(2:IIU ,2:IJU ) = ( ZW_FIELD(1:IIU-1,1:IJU-1) & 00224 + ZW_FIELD(1:IIU-1,2:IJU ) & 00225 + ZW_FIELD(2:IIU ,1:IJU-1) & 00226 + ZW_FIELD(2:IIU ,2:IJU ) ) & 00227 / MAX( ZW (1:IIU-1,1:IJU-1) & 00228 + ZW (1:IIU-1,2:IJU ) & 00229 + ZW (2:IIU ,1:IJU-1) & 00230 + ZW (2:IIU ,2:IJU ) , ZEPS) 00231 ! 00232 IF (IJU>1) & 00233 ZFIELD_XY(1 ,2:IJU ) = ( ZW_FIELD(1 ,1:IJU-1) & 00234 + ZW_FIELD(1 ,2:IJU ) ) & 00235 / MAX( ZW (1 ,1:IJU-1) & 00236 + ZW (1 ,2:IJU ) , ZEPS) 00237 ! 00238 IF (IJU>1) & 00239 ZFIELD_XY( IIU+1,2:IJU ) = ( ZW_FIELD(IIU ,1:IJU-1) & 00240 + ZW_FIELD(IIU ,2:IJU ) ) & 00241 / MAX( ZW (IIU ,1:IJU-1) & 00242 + ZW (IIU ,2:IJU ) , ZEPS) 00243 ! 00244 IF (IIU>1) & 00245 ZFIELD_XY(2:IIU ,1 ) = ( ZW_FIELD(1:IIU-1,1 ) & 00246 + ZW_FIELD(2:IIU ,1 ) ) & 00247 / MAX( ZW (1:IIU-1,1 ) & 00248 + ZW (2:IIU ,1 ) , ZEPS) 00249 ! 00250 IF (IIU>1) & 00251 ZFIELD_XY(2:IIU ,IJU+1 ) = ( ZW_FIELD(1:IIU-1,IJU ) & 00252 + ZW_FIELD(2:IIU ,IJU ) ) & 00253 / MAX( ZW (1:IIU-1,IJU ) & 00254 + ZW (2:IIU ,IJU ) , ZEPS) 00255 ! 00256 ZFIELD_XY(1 ,1 ) = ( PFIELD1(1 ,1 ) ) 00257 ZFIELD_XY(IIU+1 ,1 ) = ( PFIELD1(IIU ,1 ) ) 00258 ZFIELD_XY(1 ,IJU+1 ) = ( PFIELD1(1 ,IJU ) ) 00259 ZFIELD_XY(IIU+1 ,IJU+1 ) = ( PFIELD1(IIU ,IJU ) ) 00260 ! 00261 !* 3.2 FIELD type value is 0 in one grid mesh 00262 ! -------------------------------------- 00263 ! 00264 WHERE ( PFIELD1(1:IIU,1:IJU)==0.) 00265 ZFIELD_XY(1:IIU,1:IJU) = 0. 00266 END WHERE 00267 ! 00268 WHERE ( PFIELD1(1:IIU,1:IJU)==0.) 00269 ZFIELD_XY(1:IIU,2:IJU+1) = 0. 00270 END WHERE 00271 ! 00272 WHERE ( PFIELD1(1:IIU,1:IJU)==0.) 00273 ZFIELD_XY(2:IIU+1,1:IJU) = 0. 00274 END WHERE 00275 ! 00276 WHERE ( PFIELD1(1:IIU,1:IJU)==0.) 00277 ZFIELD_XY(2:IIU+1,2:IJU+1) = 0. 00278 END WHERE 00279 ! 00280 !------------------------------------------------------------------------------- 00281 ! 00282 !* 4. Coordinates of points between input grid points 00283 ! ----------------------------------------------- 00284 ! 00285 IF (IIU>1) THEN 00286 ZX(IIU+1) = 1.5*PX1(IIU) -0.5*PX1(IIU-1) 00287 ZX(2:IIU) = 0.5*PX1(1:IIU-1)+0.5*PX1(2:IIU) 00288 ZX(1) = 1.5*PX1(1) -0.5*PX1(2) 00289 ELSE 00290 ZX(1) = PX1(1) - 1.E6 ! uniform field in X direction if only 1 point is 00291 ZX(2) = PX1(1) + 1.E6 ! available. Arbitrary mesh length of 2000km assumed 00292 END IF 00293 ! 00294 IF (IJU>1) THEN 00295 ZY(IJU+1) = 1.5*PY1(IJU) -0.5*PY1(IJU-1) 00296 ZY(2:IJU) = 0.5*PY1(1:IJU-1)+0.5*PY1(2:IJU) 00297 ZY(1) = 1.5*PY1(1) -0.5*PY1(2) 00298 ELSE 00299 ZY(1) = PY1(1) - 1.E6 ! uniform field in Y direction if only 1 point is 00300 ZY(2) = PY1(1) + 1.E6 ! available. Arbitrary mesh length of 2000km assumed 00301 END IF 00302 ! 00303 !------------------------------------------------------------------------------- 00304 ! 00305 !* 5. Interpolation 00306 ! ------------- 00307 ! 00308 PFIELD2(:) = XUNDEF 00309 ! 00310 !* loop on all output grid points 00311 ! 00312 DO JL=1,SIZE(PFIELD2,1) 00313 ! 00314 !* interpolation weights in X direction 00315 ! 00316 JI=COUNT(ZX(:)<=PX2(JL)) 00317 JI=MAX(MIN(JI,IIU),0) 00318 IF (PX1(JI)<=PX2(JL)) THEN 00319 ZC2_X = (PX2(JL)-ZX(JI+1))/(PX1(JI)-ZX(JI+1)) 00320 ZC2_X = MAX(MIN(ZC2_X,1.),0.) 00321 ZC3_X = 1. - ZC2_X 00322 ZC1_X = 0. 00323 ELSE 00324 ZC2_X = (PX2(JL)-ZX(JI))/(PX1(JI)-ZX(JI)) 00325 ZC2_X = MAX(MIN(ZC2_X,1.),0.) 00326 ZC1_X = 1. - ZC2_X 00327 ZC3_X = 0. 00328 END IF 00329 ! 00330 !* interpolation weights in Y direction 00331 ! 00332 JJ=COUNT(ZY(:)<=PY2(JL)) 00333 JJ=MAX(MIN(JJ,IJU),0) 00334 IF (PY1(JJ)<=PY2(JL)) THEN 00335 ZC2_Y = (PY2(JL)-ZY(JJ+1))/(PY1(JJ)-ZY(JJ+1)) 00336 ZC2_Y = MAX(MIN(ZC2_Y,1.),0.) 00337 ZC3_Y = 1. - ZC2_Y 00338 ZC1_Y = 0. 00339 ELSE 00340 ZC2_Y = (PY2(JL)-ZY(JJ))/(PY1(JJ)-ZY(JJ)) 00341 ZC2_Y = MAX(MIN(ZC2_Y,1.),0.) 00342 ZC1_Y = 1. - ZC2_Y 00343 ZC3_Y = 0. 00344 END IF 00345 ! 00346 !* interpolation 00347 ! 00348 IF(PFIELD1(JI,JJ) /= XUNDEF) & 00349 PFIELD2(JL) = & 00350 ZC1_Y * ( ZC1_X * ZFIELD_XY(JI ,JJ ) & 00351 + ZC2_X * ZFIELD_Y (JI ,JJ ) & 00352 + ZC3_X * ZFIELD_XY(JI+1,JJ ) ) & 00353 + ZC2_Y * ( ZC1_X * ZFIELD_X (JI ,JJ ) & 00354 + ZC2_X * PFIELD1 (JI ,JJ ) & 00355 + ZC3_X * ZFIELD_X (JI+1,JJ ) ) & 00356 + ZC3_Y * ( ZC1_X * ZFIELD_XY(JI ,JJ+1) & 00357 + ZC2_X * ZFIELD_Y (JI ,JJ+1) & 00358 + ZC3_X * ZFIELD_XY(JI+1,JJ+1) ) 00359 00360 END DO 00361 ! 00362 !------------------------------------------------------------------------------- 00363 ! 00364 DEALLOCATE(ZW) 00365 DEALLOCATE(ZW_FIELD) 00366 DEALLOCATE(ZFIELD_X ) 00367 DEALLOCATE(ZFIELD_Y ) 00368 DEALLOCATE(ZFIELD_XY) 00369 ! 00370 !------------------------------------------------------------------------------- 00371 ! 00372 WHERE(ABS(PFIELD2-XUNDEF)<1.E-6) PFIELD2=XUNDEF 00373 ! 00374 !------------------------------------------------------------------------------ 00375 ! 00376 !* 6. EXTRAPOLATIONS IF SOME POINTS WERE NOT INTERPOLATED 00377 ! --------------------------------------------------- 00378 ! 00379 !* no missing point 00380 IF (COUNT(PFIELD2(:)==XUNDEF .AND. OINTERP(:))==0 .AND. LHOOK) CALL DR_HOOK('BILIN',1,ZHOOK_HANDLE) 00381 IF (COUNT(PFIELD2(:)==XUNDEF .AND. OINTERP(:))==0) RETURN 00382 ! 00383 !* no data point 00384 IF (COUNT(PFIELD1(:,:)/=XUNDEF)==0 .AND. LHOOK) CALL DR_HOOK('BILIN',1,ZHOOK_HANDLE) 00385 IF (COUNT(PFIELD1(:,:)/=XUNDEF)==0) RETURN 00386 ! 00387 WRITE(KLUOUT,*) ' Remaining horizontal extrapolations' 00388 WRITE(KLUOUT,*) ' Total number of input data : ',COUNT(PFIELD1(:,:)/=XUNDEF) 00389 WRITE(KLUOUT,*) ' Number of points to interpolate: ',COUNT(PFIELD2(:)==XUNDEF .AND. OINTERP(:)) 00390 ! 00391 !* input grid coordinates 00392 ! 00393 ALLOCATE(ZX1(IIU*IJU)) 00394 ALLOCATE(ZY1(IIU*IJU)) 00395 ALLOCATE(ZFIELD1(IIU*IJU)) 00396 IIN=0 00397 DO JJ=1,IJU 00398 DO JI=1,IIU 00399 IIN = IIN +1 00400 ZX1 (IIN) = PX1 (JI) 00401 ZY1 (IIN) = PY1 (JJ) 00402 ZFIELD1(IIN) = PFIELD1(JI,JJ) 00403 END DO 00404 END DO 00405 ! 00406 CALL HOR_EXTRAPOL_SURF(KLUOUT,'XY ',ZY1,ZX1,ZFIELD1,PY2,PX2,PFIELD2,OINTERP) 00407 ! 00408 DEALLOCATE(ZX1) 00409 DEALLOCATE(ZY1) 00410 DEALLOCATE(ZFIELD1) 00411 IF (LHOOK) CALL DR_HOOK('BILIN',1,ZHOOK_HANDLE) 00412 !------------------------------------------------------------------------------- 00413 ! 00414 END SUBROUTINE BILIN
1.8.0
