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