|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE HORIBL_SURF(PILA1,PILO1,PILA2,PILO2,KINLA,KINLO,KILEN,PARIN, & 00003 KOLEN,PXOUT,PYOUT,PAROUT,ODVECT,KLUOUT,OINTERP, & 00004 KLSMIN,KLSMOUT ) 00005 ! ########################################################################### 00006 ! 00007 !!**** *HORIBL_SURF* - horitontal bilinear interpolation 00008 !! 00009 !! PURPOSE 00010 !! ------- 00011 !! 00012 !! Interpolates a field, supports masks. 00013 !! 00014 !! METHOD 00015 !! ------ 00016 !! 00017 !! This routine performs a bilinear interpolation based on the 12 surrounding 00018 !! points. It begins with an interpolation along the latitudes (with third order 00019 !! polynoms interpolation with 4 points and linear interpolation for 2 points) 00020 !! and then a second along the longitude (third order polynoms interpolation). 00021 !! Two interpolations are performed : first along the parallels then between the 00022 !! four resulting points. 00023 !! 00024 !! The disposition of the points is the following : 00025 !! 00026 !! 00027 !! N 1 2 00028 !! 00029 !! ^ 3 4 5 6 00030 !! | x 00031 !! | 7 8 9 10 00032 !! | 00033 !! 11 12 00034 !! S 00035 !! W ---------------> E 00036 !! 00037 !! Note : the name 'south', 'north', may not be exact if the last data point is 00038 !! to the south of first (delta latitude < 0). This does not affect computations. 00039 !! 00040 !! The formula used to compute the weight is : 00041 !! (Lon - Lon.i) . (Lon - Lon.i) . (Lon - Lon.i) 00042 !! Wi = --------------------------------------------------- 00043 !! (Lon.i - Lon.j) . (Lon.i - Lon.k) . (Lon.i - Lon.l) 00044 !! Where j,k,l are the other points of the line. 00045 !! 00046 !! When masks are used, points with different types than the output points are 00047 !! not taken in account (in the formula, the corresponding coefficient is set 00048 !! to 1). If no points of the same nature are available, the interpolation is 00049 !! performed anyway with the 12 points. It is the task of the calling program 00050 !! to react to this situation. 00051 !! 00052 !! When the inputs parameters define a circular map (or global), the inputs data 00053 !! are extended. The value of the parameter ODVECT is used to know if the datas 00054 !! are vectorial or scalar (this affects the sign of extended values). 00055 !! 00056 !! EXTERNAL 00057 !! -------- 00058 !! 00059 !! subroutine FMLOOK_ll : to retrieve the logical unit number of the listing file 00060 !! 00061 !! IMPLICIT ARGUMENTS 00062 !! ------------------ 00063 !! 00064 !! REFERENCE 00065 !! --------- 00066 !! 00067 !! This routine is based on the one used by the software FULL-POS from Meteo France. 00068 !! More informations may be found in 'Book 1' 00069 !! 00070 !! AUTHOR 00071 !! ------ 00072 !! 00073 !! J.Pettre & V.Bousquet 00074 !! 00075 !! MODIFICATIONS 00076 !! ------------- 00077 !! 00078 !! Original 07/01/1999 00079 !! 21/04/1999 (V. Masson) set correct prefixes and bug in 00080 !! a logical definition 00081 !! 21/04/1999 (V. Masson) bug in north and south poles 00082 !! extension for input map land-sea mask 00083 !! 27/05/1999 (V. Masson) bug in 'grib south pole' 00084 !! extrapolation (number of point per parallel) 00085 !! 27/05/1999 (V. Masson) bug in 'grib pole' extrapolation 00086 !! extra latitudes are now computed symetrically 00087 !! to the poles. 00088 !! 17/03/2010 (P. LeMoigne) bug in weights computations 00089 !! 16/06/2010 (G. Tanguy) bug in 'grib north pole" 00090 !! extrapolation (tabular ZARIN not totaly filled) 00091 !! 00092 !------------------------------------------------------------------------------ 00093 ! 00094 ! 00095 !* 0. DECLARATIONS 00096 ! --------------- 00097 ! 00098 USE MODI_HOR_EXTRAPOL_SURF 00099 ! 00100 USE MODD_SURF_PAR, ONLY : XUNDEF 00101 ! 00102 ! 00103 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00104 USE PARKIND1 ,ONLY : JPRB 00105 ! 00106 USE MODI_ABOR1_SFX 00107 ! 00108 IMPLICIT NONE 00109 ! 00110 !* 0.1. Declaration of arguments 00111 ! 00112 REAL, INTENT(IN) :: PILA1 ! Lat. (y) of first input point 00113 REAL, INTENT(IN) :: PILO1 ! Lon. (x) of first input point 00114 REAL, INTENT(IN) :: PILA2 ! Lat. (y) of last input point 00115 REAL, INTENT(IN) :: PILO2 ! Lon. (x) of last input point 00116 INTEGER, INTENT(IN) :: KINLA ! Number of parallels 00117 INTEGER, DIMENSION(KINLA), INTENT(IN) :: KINLO ! Number of point along a parallel 00118 INTEGER, INTENT(IN) :: KILEN ! size of input arrays 00119 REAL, DIMENSION(KILEN), INTENT(IN) :: PARIN ! input array 00120 INTEGER, INTENT(IN) :: KOLEN ! size of output array 00121 REAL, DIMENSION(KOLEN), INTENT(IN) :: PXOUT ! X (lon.) of output points 00122 REAL, DIMENSION(KOLEN), INTENT(IN) :: PYOUT ! Y (lat.) of output points 00123 REAL, DIMENSION(KOLEN), INTENT(OUT) :: PAROUT ! output array 00124 LOGICAL, DIMENSION(KOLEN), INTENT(IN) :: OINTERP ! .true. where physical value is needed 00125 LOGICAL, INTENT(IN) :: ODVECT ! data is vectorial (True/False) 00126 INTEGER, INTENT(IN) :: KLUOUT ! output listing logical unit 00127 INTEGER, DIMENSION(KILEN), INTENT(IN), OPTIONAL :: KLSMIN ! input land/sea mask 00128 INTEGER, DIMENSION(KOLEN), INTENT(IN), OPTIONAL :: KLSMOUT ! output land/sea mask 00129 ! 00130 !* 0.2. Declaration of local variables 00131 ! 00132 ! Variables used to perform the interpolation 00133 REAL :: ZOLA ! Latitude of the output point 00134 REAL :: ZOLO ! Longitude of the output point 00135 REAL :: ZIDLA ! Delta latitude 00136 REAL :: ZIDLO ! Delta longitude 00137 INTEGER, DIMENSION(:), ALLOCATABLE :: IOFS ! Offset of each parallel in the array 00138 ! Number of the surrounding latitudes 00139 INTEGER :: IOSS,IOS,ION,IONN 00140 ! Posiiton in the array of the twelwe surrounding points 00141 INTEGER :: IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8,IP9,IP10, 00142 IP11,IP12 00143 ! Latitudes and longitudes of the surrounding points 00144 REAL :: ZLANN,ZLAN,ZLAS,ZLASS 00145 REAL :: ZLOP1,ZLOP2,ZLOP3,ZLOP4 ,ZLOP5 ,ZLOP6, 00146 ZLOP7,ZLOP8,ZLOP9,ZLOP10,ZLOP11,ZLOP12 00147 ! Weights of the latitudes and of the points 00148 REAL :: ZWNN,ZWN,ZWS,ZWSS 00149 REAL :: ZW1,ZW2,ZW3,ZW4,ZW5,ZW6,ZW7,ZW8,ZW9,ZW10, 00150 ZW11,ZW12 00151 ! Land/sea mask coefficient for each point : 0 -> point not taken in account, 00152 ! 1 -> point taken in account 00153 REAL :: ZLSM1,ZLSM2 ,ZLSM3 ,ZLSM4 ,ZLSM5 ,ZLSM6,ZLSM7,ZLSM8, 00154 ZLSM9,ZLSM10,ZLSM11,ZLSM12,ZLSMNN,ZLSMN,ZLSMS,ZLSMSS, 00155 ZLSMTOT 00156 ! Variables implied in the extension procedure 00157 REAL :: ZILO1 ! Longitude of the first data point 00158 REAL :: ZILO2 ! Longitude of the last data point 00159 LOGICAL :: GGLOBLON ! True if the map is circular 00160 LOGICAL :: GGLOBN ! True if the map has the north pole 00161 LOGICAL :: GGLOBS ! True if the map has the south pole 00162 INTEGER :: IBIGSIZE ! Size of the extended map 00163 INTEGER :: IMIDDLE ! Used for extensions around the poles 00164 INTEGER :: IOFFSET1 ! Offset in map 00165 INTEGER :: IOFFSET2 ! Offset in map 00166 REAL :: ZSOUTHPOLE! south pole latitude (-90 or 90) 00167 REAL :: ZNORTHPOLE! north pole latitude ( 90 or -90) 00168 REAL, DIMENSION(:), ALLOCATABLE :: ZLA ! input "latitude" coordinate 00169 REAL, DIMENSION(:), ALLOCATABLE :: ZLO ! input "longitude" coordinate 00170 REAL, DIMENSION(:), ALLOCATABLE :: ZARIN ! Extended input datas 00171 INTEGER, DIMENSION(:), ALLOCATABLE :: ILSMIN ! Extended land/sea mask 00172 INTEGER, DIMENSION(:), ALLOCATABLE :: IINLO ! Extended KINLO 00173 INTEGER :: IINLA ! Number of parallel 00174 REAL :: ZVECT ! -1 if input is vectorial 00175 LOGICAL :: LDLSM ! Specify if land/sea mask is present or not 00176 ! Loop counters 00177 INTEGER :: JOPOS ! Output position 00178 INTEGER :: JIPOS ! Input position 00179 INTEGER :: JLOOP1 ! Dummy counter 00180 ! 00181 ! 00182 !------------------------------------------------------------------------------ 00183 REAL :: ZMAX ! Max of 12 surrounding values 00184 REAL :: ZMIN ! Min of 12 surrounding values 00185 INTEGER :: JLOOP2 ! Dummy counter 00186 INTEGER, DIMENSION(12) :: IP ! Array for IPn 00187 INTEGER :: JLAT ! latitude loop counter 00188 INTEGER :: JLON ! longitude loop counter 00189 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00190 !------------------------------------------------------------------------------ 00191 ! 00192 !* 1. DETERMINATION of the latitude of the poles (depending of the latitude 00193 ! ------------- of the first data point) 00194 ! 00195 IF (LHOOK) CALL DR_HOOK('HORIBL_SURF',0,ZHOOK_HANDLE) 00196 IF (PILA1>0.) THEN 00197 ZSOUTHPOLE= 90. 00198 ZNORTHPOLE=-90. 00199 ELSE 00200 ZSOUTHPOLE=-90. 00201 ZNORTHPOLE= 90. 00202 END IF 00203 ! 00204 !------------------------------------------------------------------------------ 00205 ! 00206 !* 2. EXTEND DATA GRID 00207 ! ---------------- 00208 ! Land / Sea mask 00209 LDLSM = .FALSE. 00210 IF (PRESENT(KLSMIN) .AND. PRESENT(KLSMOUT)) LDLSM = .TRUE. 00211 ! 00212 !* 2.1 Alias input data 00213 ! 00214 ZILO1 = PILO1 00215 ZILO2 = PILO2 00216 ZVECT = 1. 00217 IF (ODVECT) ZVECT=-1. 00218 ! 00219 !* 2.2 Center input domain in order to have Lo1 < Lo 2 00220 ! 00221 IF (ZILO2 < 0.) ZILO2 = ZILO2 + 360. 00222 IF (ZILO1 < 0.) ZILO1 = ZILO1 + 360. 00223 IF (ZILO2 < ZILO1) ZILO1 = ZILO1 - 360. 00224 ! 00225 !* 2.3 Extend one point (needed for reduced grids) 00226 ! 00227 ! Longitude coordinate of points are found by : 00228 ! i 00229 ! Lon(i) = Lon1 + ------------- . (Lon2 - Lon1) 00230 ! Npts(Lat)-1 00231 ! Where i goes from 0 to Npts(Lat)-1. The result of this is that the last point of 00232 ! each parallel is located at Lon2. This is not the case for reduced grid where the 00233 ! position of the last point depends upon the number of points of the parallel. For 00234 ! reduced grid, the right formula to use is the following : 00235 ! i 00236 ! Lon(i) = Lon1 + ----------- . (Lon2' - Lon1) 00237 ! Npts(Lat) 00238 ! Where Lon2' = Lon1 + 2.PI. 00239 ! 00240 ! Lon2 - Lon1 00241 ! This can be generalized with Lon2' = Lon2 + ------------- 00242 ! Nptsmax - 1 00243 ! 00244 JOPOS = MAXVAL(KINLO(1:KINLA)) 00245 ZILO2 = ZILO1 + (ZILO2 - ZILO1) * JOPOS / (JOPOS - 1.) 00246 ! 00247 ! 00248 !* 2.4 Test if the input is global or partially global 00249 ! 00250 ! Note that we must have a global map to make extension around the poles 00251 GGLOBN = .FALSE. 00252 GGLOBS = .FALSE. 00253 GGLOBLON = .FALSE. 00254 IF (ZILO2-360.>ZILO1-1.E-3) GGLOBLON = .TRUE. 00255 ZIDLA = (PILA2 - PILA1) / (KINLA - 1) 00256 IF ((PILA1-ZIDLA>= 90.) .OR. (PILA1-ZIDLA<=-90.)) GGLOBS=GGLOBLON 00257 IF ((PILA2+ZIDLA>= 90.) .OR. (PILA2+ZIDLA<=-90.)) GGLOBN=GGLOBLON 00258 ! Aladin case (input PILA2, PILO2 are in meters) no extension 00259 IF ( PILA2 > 100. ) THEN 00260 GGLOBN = .FALSE. 00261 GGLOBS = .FALSE. 00262 GGLOBLON = .FALSE. 00263 END IF 00264 ! 00265 ! 00266 !* 2.5 Compute the size of the resulting map 00267 ! 00268 IBIGSIZE = KILEN 00269 IF (GGLOBS ) IBIGSIZE=IBIGSIZE+(4+KINLO( 1))+(4+KINLO( 2)) 00270 IF (GGLOBN ) IBIGSIZE=IBIGSIZE+(4+KINLO(KINLA))+(4+KINLO(KINLA-1)) 00271 IF (GGLOBLON) IBIGSIZE=IBIGSIZE+ 4*KINLA 00272 ! 00273 !* 2.6 Compute the resulting map 00274 ! 00275 ALLOCATE (ZARIN(IBIGSIZE)) 00276 ALLOCATE (ILSMIN(IBIGSIZE)) 00277 ! 00278 ! 2.6.1 Compute the longitude extension 00279 ! 00280 ! This is a basic copy of the data. If extension is possible, the first and last 00281 ! two lines are copied twice this way : 00282 ! 00283 ! /---------------\ 00284 ! | | 00285 ! [.] [.] [.... ...] [.] [.] 00286 ! | | 00287 ! \------------/ 00288 ! 00289 ! A point represent a data. 00290 ! 00291 JIPOS = 1 00292 JOPOS = 1 00293 IF (GGLOBS) JOPOS=JOPOS+(4+KINLO(1))+(4+KINLO(2)) 00294 IF (GGLOBLON) THEN 00295 DO JLOOP1 = 1, KINLA 00296 ZARIN(JOPOS ) = PARIN(JIPOS+KINLO(JLOOP1)-2) 00297 ZARIN(JOPOS+1) = PARIN(JIPOS+KINLO(JLOOP1)-1) 00298 ZARIN(JOPOS+2:JOPOS+2+KINLO(JLOOP1)-1) = PARIN(JIPOS:JIPOS+KINLO(JLOOP1)-1) 00299 ZARIN(JOPOS+2+KINLO(JLOOP1) ) = PARIN(JIPOS ) 00300 ZARIN(JOPOS+2+KINLO(JLOOP1)+1) = PARIN(JIPOS+1) 00301 IF (LDLSM) THEN 00302 ILSMIN(JOPOS ) = KLSMIN(JIPOS+KINLO(JLOOP1)-2) 00303 ILSMIN(JOPOS+1) = KLSMIN(JIPOS+KINLO(JLOOP1)-1) 00304 ILSMIN(JOPOS+2:JOPOS+2+KINLO(JLOOP1)-1) = KLSMIN(JIPOS:JIPOS+KINLO(JLOOP1)-1) 00305 ILSMIN(JOPOS+2+KINLO(JLOOP1) ) = KLSMIN(JIPOS ) 00306 ILSMIN(JOPOS+2+KINLO(JLOOP1)+1) = KLSMIN(JIPOS+1) 00307 END IF 00308 JIPOS = JIPOS + KINLO(JLOOP1) 00309 JOPOS = JOPOS + KINLO(JLOOP1) + 4 00310 END DO 00311 ELSE 00312 ZARIN(JOPOS:JOPOS+KILEN-1) = PARIN(JIPOS:JIPOS+KILEN-1) 00313 IF (LDLSM) THEN 00314 ILSMIN(JOPOS:JOPOS+KILEN-1) = KLSMIN(JIPOS:JIPOS+KILEN-1) 00315 END IF 00316 END IF 00317 ! 00318 ! 2.6.2 Compute the south pole extension 00319 ! 00320 ! Pole extension is performed by copying the first half datas to the last half 00321 ! datas of the extension parallel : 00322 ! 00323 ! [.] [.] [....] [....] [.] [.] 00324 ! |||| 00325 ! /-------/ 00326 ! |||| 00327 ! [.] [.] [....] [....] [.] [.] 00328 ! 00329 IF (GGLOBS) THEN ! South pole (south meaning begining of the grib) 00330 IOFFSET1 = 4 + KINLO(2) 00331 IOFFSET2 = IOFFSET1 + 4 + KINLO(1) 00332 IMIDDLE = (KINLO(1)+4) / 2 00333 ZARIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00334 ZVECT*ZARIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00335 ZARIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(1)+4) = & 00336 ZVECT*ZARIN(IOFFSET2+1+2:IOFFSET2+KINLO(1)+4-IMIDDLE+2) 00337 IF (LDLSM) THEN 00338 ILSMIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00339 ILSMIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00340 ILSMIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(1)+4) = & 00341 ILSMIN(IOFFSET2+1+2:IOFFSET2+KINLO(1)+4-IMIDDLE+2) 00342 END IF 00343 IOFFSET2 = IOFFSET2 + 4 + KINLO(1) 00344 IMIDDLE = (KINLO(2)+4) / 2 00345 ZARIN(1:IMIDDLE) = ZVECT*ZARIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00346 ZARIN(IMIDDLE+1:KINLO(2)+4) = & 00347 ZVECT*ZARIN(IOFFSET2+1+2:IOFFSET2+KINLO(2)+4-IMIDDLE+2) 00348 IF (LDLSM) THEN 00349 ILSMIN(1:IMIDDLE) = ILSMIN(IOFFSET2+1+IMIDDLE:IOFFSET2+2*IMIDDLE) 00350 ILSMIN(IMIDDLE+1:KINLO(2)+4) = ILSMIN(IOFFSET2+1+2:IOFFSET2+KINLO(2)+4-IMIDDLE+2) 00351 END IF 00352 END IF 00353 ! 00354 ! 2.6.3 Compute the north pole extension 00355 ! 00356 IF (GGLOBN) THEN ! North pole (north meaning end of the grib) 00357 IOFFSET1 = IBIGSIZE - (4+KINLO(KINLA-1)) - (4+KINLO(KINLA)) 00358 IOFFSET2 = IOFFSET1 - (4+KINLO(KINLA)) 00359 IMIDDLE = (KINLO(KINLA)+4) / 2 00360 ZARIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00361 ZVECT*ZARIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00362 ZARIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(KINLA)+4) = & 00363 ZVECT*ZARIN(IOFFSET2+1+2:IOFFSET2+KINLO(KINLA)+4-IMIDDLE+2) 00364 IF (LDLSM) THEN 00365 ILSMIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00366 ILSMIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00367 ILSMIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(KINLA)+4) = & 00368 ILSMIN(IOFFSET2+1+2:IOFFSET2+KINLO(KINLA)+4-IMIDDLE+2) 00369 END IF 00370 IOFFSET1 = IOFFSET1 + (4+KINLO(KINLA)) 00371 IOFFSET2 = IOFFSET2 - (4+KINLO(KINLA-1)) 00372 IMIDDLE = (KINLO(KINLA-1)+4) / 2 00373 ZARIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00374 ZVECT*ZARIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00375 ZARIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(KINLA-1)+4) = & 00376 ZVECT*ZARIN(IOFFSET2+1+2:IOFFSET2+KINLO(KINLA-1)+4-IMIDDLE+2) 00377 IF (LDLSM) THEN 00378 ILSMIN(IOFFSET1+1:IOFFSET1+IMIDDLE) = & 00379 ILSMIN(IOFFSET2+1+IMIDDLE-2:IOFFSET2+2*IMIDDLE-2) 00380 ILSMIN(IOFFSET1+IMIDDLE+1:IOFFSET1+KINLO(KINLA-1)+4) = & 00381 ILSMIN(IOFFSET2+1+2:IOFFSET2+KINLO(KINLA-1)+4-IMIDDLE+2) 00382 END IF 00383 END IF 00384 ! 00385 !* 2.7 Compute the resulting parameters of the map 00386 ! 00387 IINLA = KINLA 00388 IF (GGLOBS) IINLA = IINLA + 2 00389 IF (GGLOBN) IINLA = IINLA + 2 00390 ! 00391 ALLOCATE (IINLO(IINLA)) 00392 IOFFSET1 = 0 00393 IF (GGLOBS) THEN 00394 IINLO(IOFFSET1+1) = KINLO(2) 00395 IINLO(IOFFSET1+2) = KINLO(1) 00396 IOFFSET1 = IOFFSET1 + 2 00397 END IF 00398 IINLO(IOFFSET1+1:IOFFSET1+KINLA) = KINLO(1:KINLA) 00399 IOFFSET1 = IOFFSET1 + KINLA 00400 IF (GGLOBN) THEN 00401 IINLO(IOFFSET1+1) = KINLO(KINLA) 00402 IINLO(IOFFSET1+2) = KINLO(KINLA-1) 00403 IOFFSET1 = IOFFSET1 + 2 00404 END IF 00405 ! 00406 !* 2.8 Compute Offset array used to acces the lines 00407 ! 00408 ALLOCATE (IOFS(IINLA)) 00409 IOFS(1) = 1 00410 IF (GGLOBLON) IOFS(1)=IOFS(1)+2 00411 DO JLOOP1=2, IINLA 00412 IOFS(JLOOP1) = IOFS(JLOOP1-1) + IINLO(JLOOP1-1) 00413 IF (GGLOBLON) IOFS(JLOOP1) = IOFS(JLOOP1) + 4 00414 END DO 00415 ! 00416 !------------------------------------------------------------------------------ 00417 ! 00418 !* 3. LOOP OVER ALL THE POINTS OF THE OUTPUT GRID 00419 ! ------------------------------------------- 00420 ! 00421 PAROUT(:) = XUNDEF 00422 ! 00423 JOPOS = 0 00424 DO JLOOP1 = 1, KOLEN 00425 JOPOS = JOPOS + 1 00426 IF (.NOT. OINTERP(JOPOS)) CYCLE 00427 ZOLA = PYOUT(JOPOS) 00428 ZOLO = PXOUT(JOPOS) 00429 IF (ZOLO < ZILO1) ZOLO = ZOLO + 360. 00430 IF (ZOLO > ZILO2) ZOLO = ZOLO - 360. 00431 ! 00432 !* 3.1 Locate the 12 input points around the interpolated output point 00433 !* 00434 !* N 1 2 00435 !* 00436 !* ^ 3 4 5 6 00437 !* | x 00438 !* | 7 8 9 10 00439 !* | 00440 !* 11 12 00441 !* S 00442 !* W ---------------> E 00443 !* 00444 !* Note : the name 'south', 'north', may not be exact if the point 2 is 00445 !* to the south of point 1 (IDLA < 0). This does not affect computation. 00446 ! 00447 ! 3.1.1. find positions of latitudes 00448 IOS = NINT( (ZOLA-PILA1)/ZIDLA -0.5) ! because of the zero 00449 IF ( IOS<-1 ) THEN 00450 CALL ABOR1_SFX('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE') 00451 END IF 00452 ZLAS = PILA1 + IOS * ZIDLA 00453 IOS = IOS + 1 00454 IF (GGLOBS) THEN 00455 IOS = IOS + 2 00456 IOSS = IOS - 1 00457 ION = IOS + 1 00458 IONN = ION + 1 00459 ELSE 00460 IOSS = MAX(IOS - 1,1) 00461 ION = MIN(IOS + 1,IINLA) 00462 IONN = MIN(ION + 1,IINLA) 00463 IOS = MAX(IOS, 1) 00464 ENDIF 00465 ZLASS = ZLAS - ZIDLA 00466 ZLAN = ZLAS + ZIDLA 00467 ZLANN = ZLAN + ZIDLA 00468 ! 00469 ! extra latitudes are computed symetrically compared to the poles 00470 ! 00471 IF (GGLOBS .AND. IOS==2) THEN 00472 ZLASS = 2. * ZSOUTHPOLE - ZLANN 00473 ZLAS = 2. * ZSOUTHPOLE - ZLAN 00474 END IF 00475 IF (GGLOBS .AND. IOS==3) THEN 00476 ZLASS = 2. * ZSOUTHPOLE - ZLAS 00477 END IF 00478 IF (GGLOBN .AND. IOS==IINLA-2) THEN 00479 ZLANN = 2. * ZNORTHPOLE - ZLASS 00480 ZLAN = 2. * ZNORTHPOLE - ZLAS 00481 END IF 00482 IF (GGLOBN .AND. IOS==IINLA-3) THEN 00483 ZLANN = 2. * ZNORTHPOLE - ZLAN 00484 END IF 00485 ! 00486 IF ((IOSS<1).OR.(IONN>IINLA).OR. & 00487 (IOSS<1).OR.(IONN>IINLA)) THEN 00488 CALL ABOR1_SFX('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE') 00489 END IF 00490 ! 00491 ! 3.1.2. northern 00492 ZIDLO = (ZILO2 - ZILO1) / (IINLO(IONN)) 00493 IP1 = INT((ZOLO - ZILO1) / ZIDLO) 00494 IF (GGLOBLON) THEN 00495 IP2 = IP1 + 1 00496 ELSE 00497 IP2 = MIN(IINLO(IONN)-1, IP1+1) 00498 ENDIF 00499 ZLOP1 = ZILO1 + IP1 * ZIDLO 00500 ZLOP2 = ZLOP1 + ZIDLO 00501 ! 00502 ! 3.1.3. north 00503 ZIDLO = (ZILO2 - ZILO1) / (IINLO(ION )) 00504 IP4 = INT((ZOLO - ZILO1) / ZIDLO) 00505 IF (GGLOBLON) THEN 00506 IP3 = IP4 - 1 00507 IP5 = IP4 + 1 00508 IP6 = IP5 + 1 00509 ELSE 00510 IP3 = MAX(0,IP4-1) 00511 IP5 = MIN(IINLO(ION)-1,IP4+1) 00512 IP6 = MIN(IINLO(ION)-1,IP5+1) 00513 ENDIF 00514 ZLOP4 = ZILO1 + IP4 * ZIDLO 00515 ZLOP3 = ZLOP4 - ZIDLO 00516 ZLOP5 = ZLOP4 + ZIDLO 00517 ZLOP6 = ZLOP5 + ZIDLO 00518 ! 00519 ! 3.1.4. south 00520 ZIDLO = (ZILO2 - ZILO1) / (IINLO(IOS )) 00521 IP8 = INT((ZOLO - ZILO1) / ZIDLO) 00522 IF (GGLOBLON) THEN 00523 IP7 = IP8 - 1 00524 IP9 = IP8 + 1 00525 IP10 = IP9 + 1 00526 ELSE 00527 IP7 = MAX(0,IP8-1) 00528 IP9 = MIN(IINLO(IOS)-1,IP8+1) 00529 IP10 = MIN(IINLO(IOS)-1,IP9+1) 00530 ENDIF 00531 ZLOP8 = ZILO1 + IP8 * ZIDLO 00532 ZLOP7 = ZLOP8 - ZIDLO 00533 ZLOP9 = ZLOP8 + ZIDLO 00534 ZLOP10= ZLOP9 + ZIDLO 00535 ! 00536 ! 3.1.5. southern 00537 ZIDLO = (ZILO2 - ZILO1) / (IINLO(IOSS)) 00538 IP11 = INT((ZOLO - ZILO1) / ZIDLO) 00539 IF (GGLOBLON) THEN 00540 IP12 = IP11 + 1 00541 ELSE 00542 IP12 = MIN(IINLO(IOSS)-1,IP11+1) 00543 ENDIF 00544 ZLOP11= ZILO1 + IP11* ZIDLO 00545 ZLOP12= ZLOP11+ ZIDLO 00546 ! 00547 ! 3.1.6. check position of points 00548 IF (GGLOBLON ) THEN 00549 IF ((IP1 <-2) .OR. (IP2 >IINLO(IONN)+1) .OR. & 00550 (IP3 <-2) .OR. (IP6 >IINLO(ION )+1) .OR. & 00551 (IP7 <-2) .OR. (IP10>IINLO(IOS )+1) .OR. & 00552 (IP11<-2) .OR. (IP12>IINLO(IOSS)+1)) THEN 00553 CALL ABOR1_SFX('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LONGITUDE GLOBAL') 00554 END IF 00555 ELSE 00556 IF ((IP1 <0) .OR. (IP2 >IINLO(IONN)-1) .OR. & 00557 (IP3 <0) .OR. (IP6 >IINLO(ION )-1) .OR. & 00558 (IP7 <0) .OR. (IP10>IINLO(IOS )-1) .OR. & 00559 (IP11<0) .OR. (IP12>IINLO(IOSS)-1)) THEN 00560 CALL ABOR1_SFX('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LONGITUDE LOCAL') 00561 END IF 00562 END IF 00563 ! 00564 ! 3.1.7. add parallel offset 00565 IP1 =IP1 + IOFS(IONN) 00566 IP2 =IP2 + IOFS(IONN) 00567 IP3 =IP3 + IOFS(ION ) 00568 IP4 =IP4 + IOFS(ION ) 00569 IP5 =IP5 + IOFS(ION ) 00570 IP6 =IP6 + IOFS(ION ) 00571 IP7 =IP7 + IOFS(IOS ) 00572 IP8 =IP8 + IOFS(IOS ) 00573 IP9 =IP9 + IOFS(IOS ) 00574 IP10=IP10+ IOFS(IOS ) 00575 IP11=IP11+ IOFS(IOSS) 00576 IP12=IP12+ IOFS(IOSS) 00577 ! 00578 !* 3.2 Land / Sea mask 00579 ! 00580 ZLSM1 = 1. 00581 ZLSM2 = 1. 00582 ZLSM3 = 1. 00583 ZLSM4 = 1. 00584 ZLSM5 = 1. 00585 ZLSM6 = 1. 00586 ZLSM7 = 1. 00587 ZLSM8 = 1. 00588 ZLSM9 = 1. 00589 ZLSM10 = 1. 00590 ZLSM11 = 1. 00591 ZLSM12 = 1. 00592 ZLSMNN = 1. 00593 ZLSMN = 1. 00594 ZLSMS = 1. 00595 ZLSMSS = 1. 00596 IF (LDLSM) THEN 00597 IF (ILSMIN(IP1 ).NE.KLSMOUT(JOPOS)) ZLSM1 = 0. 00598 IF (ILSMIN(IP2 ).NE.KLSMOUT(JOPOS)) ZLSM2 = 0. 00599 IF (ILSMIN(IP3 ).NE.KLSMOUT(JOPOS)) ZLSM3 = 0. 00600 IF (ILSMIN(IP4 ).NE.KLSMOUT(JOPOS)) ZLSM4 = 0. 00601 IF (ILSMIN(IP5 ).NE.KLSMOUT(JOPOS)) ZLSM5 = 0. 00602 IF (ILSMIN(IP6 ).NE.KLSMOUT(JOPOS)) ZLSM6 = 0. 00603 IF (ILSMIN(IP7 ).NE.KLSMOUT(JOPOS)) ZLSM7 = 0. 00604 IF (ILSMIN(IP8 ).NE.KLSMOUT(JOPOS)) ZLSM8 = 0. 00605 IF (ILSMIN(IP9 ).NE.KLSMOUT(JOPOS)) ZLSM9 = 0. 00606 IF (ILSMIN(IP10).NE.KLSMOUT(JOPOS)) ZLSM10 = 0. 00607 IF (ILSMIN(IP11).NE.KLSMOUT(JOPOS)) ZLSM11 = 0. 00608 IF (ILSMIN(IP12).NE.KLSMOUT(JOPOS)) ZLSM12 = 0. 00609 ZLSMNN = MIN(ZLSM1 +ZLSM2,1.) 00610 ZLSMN = MIN(ZLSM3 +ZLSM4 +ZLSM5 +ZLSM6,1.) 00611 ZLSMS = MIN(ZLSM7 +ZLSM8 +ZLSM9 +ZLSM10,1.) 00612 ZLSMSS = MIN(ZLSM11+ZLSM12,1.) 00613 ZLSMTOT = MIN(ZLSMNN+ZLSMN+ZLSMS+ZLSMSS,1.) 00614 IF (ZLSMNN < 1.E-3) THEN 00615 ZLSM1 = 1. 00616 ZLSM2 = 1. 00617 END IF 00618 IF (ZLSMN < 1.E-3) THEN 00619 ZLSM3 = 1. 00620 ZLSM4 = 1. 00621 ZLSM5 = 1. 00622 ZLSM6 = 1. 00623 END IF 00624 IF (ZLSMS < 1.E-3) THEN 00625 ZLSM7 = 1. 00626 ZLSM8 = 1. 00627 ZLSM9 = 1. 00628 ZLSM10= 1. 00629 END IF 00630 IF (ZLSMSS < 1.E-3) THEN 00631 ZLSM11= 1. 00632 ZLSM12= 1. 00633 END IF 00634 IF (ZLSMTOT < 1.E-3) THEN 00635 ZLSMNN = 1. 00636 ZLSMN = 1. 00637 ZLSMS = 1. 00638 ZLSMSS = 1. 00639 END IF 00640 ENDIF 00641 ! 00642 !* 3.3 Weight of points 00643 ! 00644 ! 3.3.1 northern 00645 ZW1 = ZLSM1 * (1.+ZLSM2 *(ZOLO -ZLOP1 )/(ZLOP1 -ZLOP2 )) 00646 ZW2 = 1. - ZW1 00647 ZWNN = ZLSMNN* (1.+ZLSMN *(ZOLA -ZLANN)/(ZLANN-ZLAN )) & 00648 * (1.+ZLSMS *(ZOLA -ZLANN)/(ZLANN-ZLAS )) & 00649 * (1.+ZLSMSS*(ZOLA -ZLANN)/(ZLANN-ZLASS)) 00650 ! 00651 ! 3.3.2. north 00652 ZW3 = ZLSM3 * (1.+ZLSM4 *(ZOLO -ZLOP3 )/(ZLOP3 -ZLOP4 )) & 00653 * (1.+ZLSM5 *(ZOLO -ZLOP3 )/(ZLOP3 -ZLOP5 )) & 00654 * (1.+ZLSM6 *(ZOLO -ZLOP3 )/(ZLOP3 -ZLOP6 )) 00655 ZW4 = ZLSM4 * (1.+ZLSM3 *(ZOLO -ZLOP4 )/(ZLOP4 -ZLOP3 )) & 00656 * (1.+ZLSM5 *(ZOLO -ZLOP4 )/(ZLOP4 -ZLOP5 )) & 00657 * (1.+ZLSM6 *(ZOLO -ZLOP4 )/(ZLOP4 -ZLOP6 )) 00658 ZW5 = ZLSM5 * (1.+ZLSM3 *(ZOLO -ZLOP5 )/(ZLOP5 -ZLOP3 )) & 00659 * (1.+ZLSM4 *(ZOLO -ZLOP5 )/(ZLOP5 -ZLOP4 )) & 00660 * (1.+ZLSM6 *(ZOLO -ZLOP5 )/(ZLOP5 -ZLOP6 )) 00661 ZW6 = 1. - ZW3 - ZW4 - ZW5 00662 ZWN = ZLSMN * (1.+ZLSMNN*(ZOLA -ZLAN )/(ZLAN -ZLANN)) & 00663 * (1.+ZLSMS *(ZOLA -ZLAN )/(ZLAN -ZLAS )) & 00664 * (1.+ZLSMSS*(ZOLA -ZLAN )/(ZLAN -ZLASS)) 00665 ! 00666 ! 3.3.3. south 00667 ZW7 = ZLSM7 * (1.+ZLSM8 *(ZOLO -ZLOP7 )/(ZLOP7 -ZLOP8 )) & 00668 * (1.+ZLSM9 *(ZOLO -ZLOP7 )/(ZLOP7 -ZLOP9 )) & 00669 * (1.+ZLSM10*(ZOLO -ZLOP7 )/(ZLOP7 -ZLOP10)) 00670 ZW8 = ZLSM8 * (1.+ZLSM7 *(ZOLO -ZLOP8 )/(ZLOP8 -ZLOP7 )) & 00671 * (1.+ZLSM9 *(ZOLO -ZLOP8 )/(ZLOP8 -ZLOP9 )) & 00672 * (1.+ZLSM10*(ZOLO -ZLOP8 )/(ZLOP8 -ZLOP10)) 00673 ZW9 = ZLSM9 * (1.+ZLSM7 *(ZOLO -ZLOP9 )/(ZLOP9 -ZLOP7 )) & 00674 * (1.+ZLSM8 *(ZOLO -ZLOP9 )/(ZLOP9 -ZLOP8 )) & 00675 * (1.+ZLSM10*(ZOLO -ZLOP9 )/(ZLOP9 -ZLOP10)) 00676 ZW10 = 1. - ZW7 - ZW8 - ZW9 00677 ZWS = ZLSMS * (1.+ZLSMNN*(ZOLA -ZLAS )/(ZLAS -ZLANN)) & 00678 * (1.+ZLSMN *(ZOLA -ZLAS )/(ZLAS -ZLAN )) & 00679 * (1.+ZLSMSS*(ZOLA -ZLAS )/(ZLAS -ZLASS)) 00680 ! 00681 ! 3.3.4. southern 00682 ZW11 = ZLSM11* (1.+ZLSM12*(ZOLO -ZLOP11)/(ZLOP11-ZLOP12)) 00683 ZW12 = 1. - ZW11 00684 ZWSS = 1. - ZWNN - ZWN - ZWS 00685 ! 00686 ! In order to exclude undef values from computation of PAROUT, 00687 ! weights w2, w6, w10, w12 and wss which can be numerically very low 00688 ! because they are residual, are set to 0 00689 ! 00690 IF (ABS(ZW2 ) < 1.E-10) ZW2=0. 00691 IF (ABS(ZW6 ) < 1.E-10) ZW6=0. 00692 IF (ABS(ZW10) < 1.E-10) ZW10=0. 00693 IF (ABS(ZW12) < 1.E-10) ZW12=0. 00694 IF (ABS(ZWSS) < 1.E-10) ZWSS=0. 00695 ! 00696 ! 3.3.5. longitude weight x latitude weight 00697 ZW1 = ZW1 * ZWNN 00698 ZW2 = ZW2 * ZWNN 00699 ZW3 = ZW3 * ZWN 00700 ZW4 = ZW4 * ZWN 00701 ZW5 = ZW5 * ZWN 00702 ZW6 = ZW6 * ZWN 00703 ZW7 = ZW7 * ZWS 00704 ZW8 = ZW8 * ZWS 00705 ZW9 = ZW9 * ZWS 00706 ZW10 = ZW10 * ZWS 00707 ZW11 = ZW11 * ZWSS 00708 ZW12 = ZW12 * ZWSS 00709 ! 00710 PAROUT (JOPOS) = ZW1 * ZARIN(IP1 ) + & 00711 ZW2 * ZARIN(IP2 ) + & 00712 ZW3 * ZARIN(IP3 ) + & 00713 ZW4 * ZARIN(IP4 ) + & 00714 ZW5 * ZARIN(IP5 ) + & 00715 ZW6 * ZARIN(IP6 ) + & 00716 ZW7 * ZARIN(IP7 ) + & 00717 ZW8 * ZARIN(IP8 ) + & 00718 ZW9 * ZARIN(IP9 ) + & 00719 ZW10 * ZARIN(IP10) + & 00720 ZW11 * ZARIN(IP11) + & 00721 ZW12 * ZARIN(IP12) 00722 ! 00723 ! For surface fields, the interpoalted value is bounded 00724 ! by the min max values of the initial field 00725 00726 IF (PRESENT(KLSMIN)) THEN 00727 00728 IP(1)=IP1 00729 IP(2)=IP2 00730 IP(3)=IP3 00731 IP(4)=IP4 00732 IP(5)=IP5 00733 IP(6)=IP6 00734 IP(7)=IP7 00735 IP(8)=IP8 00736 IP(9)=IP9 00737 IP(10)=IP10 00738 IP(11)=IP11 00739 IP(12)=IP12 00740 00741 ZMIN=XUNDEF 00742 ZMAX=XUNDEF 00743 00744 DO JLOOP2=1,12 00745 IF (ZARIN(IP(JLOOP2))==XUNDEF) CYCLE 00746 00747 IF ((ZMAX==XUNDEF)) THEN 00748 ZMAX=ZARIN(IP(JLOOP2)) 00749 ZMIN=ZARIN(IP(JLOOP2)) 00750 ELSE 00751 ZMAX=MAX(ZMAX,ZARIN(IP(JLOOP2))) 00752 ZMIN=MIN(ZMIN,ZARIN(IP(JLOOP2))) 00753 ENDIF 00754 00755 END DO 00756 00757 PAROUT(JOPOS) = MAX(MIN(PAROUT(JOPOS),ZMAX),ZMIN) 00758 00759 ENDIF 00760 00761 END DO 00762 ! 00763 DEALLOCATE (IINLO) 00764 DEALLOCATE (ZARIN) 00765 DEALLOCATE (ILSMIN) 00766 DEALLOCATE (IOFS) 00767 ! 00768 WHERE(ABS(PAROUT-XUNDEF)<1.E-6) PAROUT=XUNDEF 00769 ! 00770 !------------------------------------------------------------------------------ 00771 ! 00772 !* 4. EXTRAPOLATIONS IF SOME POINTS WERE NOT INTERPOLATED 00773 ! --------------------------------------------------- 00774 ! 00775 !* no missing point 00776 IF (COUNT(PAROUT(:)==XUNDEF .AND. OINTERP(:))==0 .AND. LHOOK) CALL DR_HOOK('HORIBL_SURF',1,ZHOOK_HANDLE) 00777 IF (COUNT(PAROUT(:)==XUNDEF .AND. OINTERP(:))==0) RETURN 00778 ! 00779 !* no data point 00780 IF (COUNT(PARIN(:)/=XUNDEF)==0 .AND. LHOOK) CALL DR_HOOK('HORIBL_SURF',1,ZHOOK_HANDLE) 00781 IF (COUNT(PARIN(:)/=XUNDEF)==0) RETURN 00782 ! 00783 WRITE(KLUOUT,*) ' Remaining horizontal extrapolations' 00784 WRITE(KLUOUT,*) ' Total number of input data : ',COUNT(PARIN(:)/=XUNDEF) 00785 WRITE(KLUOUT,*) ' Number of points to interpolate: ',COUNT(PAROUT(:)==XUNDEF .AND. OINTERP(:)) 00786 ! 00787 !* input grid coordinates 00788 ! 00789 ALLOCATE(ZLA(KILEN)) 00790 ALLOCATE(ZLO(KILEN)) 00791 ! 00792 JIPOS = 0 00793 ZIDLA = (PILA2 - PILA1) / (KINLA - 1) 00794 DO JLAT=1,KINLA 00795 ZIDLO = (ZILO2-ZILO1) / KINLO(JLAT) 00796 DO JLON=1,KINLO(JLAT) 00797 JIPOS = JIPOS + 1 00798 ZLA(JIPOS) = PILA1 + (JLAT-1) * ZIDLA 00799 ZLO(JIPOS) = ZILO1 + (JLON-1) * ZIDLO 00800 END DO 00801 END DO 00802 ! 00803 CALL HOR_EXTRAPOL_SURF(KLUOUT,'LALO',ZLA,ZLO,PARIN,PYOUT,PXOUT,PAROUT,OINTERP) 00804 ! 00805 DEALLOCATE(ZLA) 00806 DEALLOCATE(ZLO) 00807 IF (LHOOK) CALL DR_HOOK('HORIBL_SURF',1,ZHOOK_HANDLE) 00808 ! 00809 !------------------------------------------------------------------------------ 00810 ! 00811 ! 00812 END SUBROUTINE HORIBL_SURF
1.8.0
