SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
horibl_surf.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########
6  SUBROUTINE horibl_surf(PILA1,PILO1,PILA2,PILO2,KINLA,KINLO,KILEN,PARIN, &
7  kolen,pxout,pyout,parout,odvect,kluout,ointerp, &
8  klsmin,klsmout )
9 ! ###########################################################################
10 !
11 !!**** *HORIBL_SURF* - horitontal bilinear interpolation
12 !!
13 !! PURPOSE
14 !! -------
15 !!
16 !! Interpolates a field, supports masks.
17 !!
18 !! METHOD
19 !! ------
20 !!
21 !! This routine performs a bilinear interpolation based on the 12 surrounding
22 !! points. It begins with an interpolation along the latitudes (with third order
23 !! polynoms interpolation with 4 points and linear interpolation for 2 points)
24 !! and then a second along the longitude (third order polynoms interpolation).
25 !! Two interpolations are performed : first along the parallels then between the
26 !! four resulting points.
27 !!
28 !! The disposition of the points is the following :
29 !!
30 !!
31 !! N 1 2
32 !!
33 !! ^ 3 4 5 6
34 !! | x
35 !! | 7 8 9 10
36 !! |
37 !! 11 12
38 !! S
39 !! W ---------------> E
40 !!
41 !! Note : the name 'south', 'north', may not be exact if the last data point is
42 !! to the south of first (delta latitude < 0). This does not affect computations.
43 !!
44 !! The formula used to compute the weight is :
45 !! (Lon - Lon.i) . (Lon - Lon.i) . (Lon - Lon.i)
46 !! Wi = ---------------------------------------------------
47 !! (Lon.i - Lon.j) . (Lon.i - Lon.k) . (Lon.i - Lon.l)
48 !! Where j,k,l are the other points of the line.
49 !!
50 !! When masks are used, points with different types than the output points are
51 !! not taken in account (in the formula, the corresponding coefficient is set
52 !! to 1). If no points of the same nature are available, the interpolation is
53 !! performed anyway with the 12 points. It is the task of the calling program
54 !! to react to this situation.
55 !!
56 !! When the inputs parameters define a circular map (or global), the inputs data
57 !! are extended. The value of the parameter ODVECT is used to know if the datas
58 !! are vectorial or scalar (this affects the sign of extended values).
59 !!
60 !! EXTERNAL
61 !! --------
62 !!
63 !! subroutine FMLOOK_ll : to retrieve the logical unit number of the listing file
64 !!
65 !! IMPLICIT ARGUMENTS
66 !! ------------------
67 !!
68 !! REFERENCE
69 !! ---------
70 !!
71 !! This routine is based on the one used by the software FULL-POS from Meteo France.
72 !! More informations may be found in 'Book 1'
73 !!
74 !! AUTHOR
75 !! ------
76 !!
77 !! J.Pettre & V.Bousquet
78 !!
79 !! MODIFICATIONS
80 !! -------------
81 !!
82 !! Original 07/01/1999
83 !! 21/04/1999 (V. Masson) set correct prefixes and bug in
84 !! a logical definition
85 !! 21/04/1999 (V. Masson) bug in north and south poles
86 !! extension for input map land-sea mask
87 !! 27/05/1999 (V. Masson) bug in 'grib south pole'
88 !! extrapolation (number of point per parallel)
89 !! 27/05/1999 (V. Masson) bug in 'grib pole' extrapolation
90 !! extra latitudes are now computed symetrically
91 !! to the poles.
92 !! 17/03/2010 (P. LeMoigne) bug in weights computations
93 !! 16/06/2010 (G. Tanguy) bug in 'grib north pole"
94 !! extrapolation (tabular ZARIN not totaly filled)
95 !!
96 !------------------------------------------------------------------------------
97 !
98 !
99 !* 0. DECLARATIONS
100 ! ---------------
101 !
102 USE modi_hor_extrapol_surf
103 !
104 USE modd_surf_par, ONLY : xundef
105 !
106 !
107 USE yomhook ,ONLY : lhook, dr_hook
108 USE parkind1 ,ONLY : jprb
109 !
110 USE modi_abor1_sfx
111 !
112 IMPLICIT NONE
113 !
114 !* 0.1. Declaration of arguments
115 !
116 REAL, INTENT(IN) :: pila1 ! Lat. (y) of first input point KDGSA
117 REAL, INTENT(IN) :: pilo1 ! Lon. (x) of first input point
118 REAL, INTENT(IN) :: pila2 ! Lat. (y) of last input point KDGEN
119 REAL, INTENT(IN) :: pilo2 ! Lon. (x) of last input point
120 INTEGER, INTENT(IN) :: kinla ! Number of parallels
121 INTEGER, DIMENSION(KINLA), INTENT(IN) :: kinlo ! Number of point along a parallel
122 INTEGER, INTENT(IN) :: kilen ! size of input arrays
123 REAL, DIMENSION(KILEN), INTENT(IN) :: parin ! input array
124 INTEGER, INTENT(IN) :: kolen ! size of output array
125 REAL, DIMENSION(KOLEN), INTENT(IN) :: pxout ! X (lon.) of output points
126 REAL, DIMENSION(KOLEN), INTENT(IN) :: pyout ! Y (lat.) of output points
127 REAL, DIMENSION(KOLEN), INTENT(OUT) :: parout ! output array
128 LOGICAL, DIMENSION(KOLEN), INTENT(IN) :: ointerp ! .true. where physical value is needed
129 LOGICAL, INTENT(IN) :: odvect ! data is vectorial (True/False)
130 INTEGER, INTENT(IN) :: kluout ! output listing logical unit
131 INTEGER, DIMENSION(KILEN), INTENT(IN), OPTIONAL :: klsmin ! input land/sea mask
132 INTEGER, DIMENSION(KOLEN), INTENT(IN), OPTIONAL :: klsmout ! output land/sea mask
133 !
134 !* 0.2. Declaration of local variables
135 !
136  ! Variables used to perform the interpolation
137 REAL :: zola ! Latitude of the output point
138 REAL :: zolo ! Longitude of the output point
139 REAL :: zidla ! Delta latitude
140 REAL :: zidlo ! Delta longitude
141 INTEGER, DIMENSION(:), ALLOCATABLE :: iofs ! Offset of each parallel in the array
142  ! Number of the surrounding latitudes
143 INTEGER :: ioss,ios,ion,ionn
144  ! Posiiton in the array of the twelwe surrounding points
145 INTEGER :: ip1,ip2,ip3,ip4,ip5,ip6,ip7,ip8,ip9,ip10, &
146  ip11,ip12
147  ! Latitudes and longitudes of the surrounding points
148 REAL :: zlann,zlan,zlas,zlass
149 REAL :: zlop1,zlop2,zlop3,zlop4 ,zlop5 ,zlop6, &
150  zlop7,zlop8,zlop9,zlop10,zlop11,zlop12
151  ! Weights of the latitudes and of the points
152 REAL :: zwnn,zwn,zws,zwss
153 REAL :: zw1,zw2,zw3,zw4,zw5,zw6,zw7,zw8,zw9,zw10, &
154  zw11,zw12
155  ! Land/sea mask coefficient for each point : 0 -> point not taken in account,
156  ! 1 -> point taken in account
157 REAL :: zlsm1,zlsm2 ,zlsm3 ,zlsm4 ,zlsm5 ,zlsm6,zlsm7,zlsm8, &
158  zlsm9,zlsm10,zlsm11,zlsm12,zlsmnn,zlsmn,zlsms,zlsmss,&
159  zlsmtot
160  ! Variables implied in the extension procedure
161 REAL :: zilo1 ! Longitude of the first data point
162 REAL :: zilo2 ! Longitude of the last data point
163 LOGICAL :: ggloblon ! True if the map is circular
164 LOGICAL :: gglobn ! True if the map has the north pole
165 LOGICAL :: gglobs ! True if the map has the south pole
166 INTEGER :: ibigsize ! Size of the extended map
167 INTEGER :: imiddle ! Used for extensions around the poles
168 INTEGER :: ioffset1 ! Offset in map
169 INTEGER :: ioffset2 ! Offset in map
170 REAL :: zsouthpole! south pole latitude (-90 or 90)
171 REAL :: znorthpole! north pole latitude ( 90 or -90)
172 REAL, DIMENSION(:), ALLOCATABLE :: zla ! input "latitude" coordinate
173 REAL, DIMENSION(:), ALLOCATABLE :: zlo ! input "longitude" coordinate
174 REAL, DIMENSION(:), ALLOCATABLE :: zarin ! Extended input datas
175 INTEGER, DIMENSION(:), ALLOCATABLE :: ilsmin ! Extended land/sea mask
176 INTEGER, DIMENSION(:), ALLOCATABLE :: iinlo ! Extended KINLO
177 INTEGER :: iinla ! Number of parallel
178 REAL :: zvect ! -1 if input is vectorial
179 LOGICAL :: ldlsm ! Specify if land/sea mask is present or not
180  ! Loop counters
181 INTEGER :: jopos ! Output position
182 INTEGER :: jipos ! Input position
183 INTEGER :: jloop1 ! Dummy counter
184 !
185 !
186 !------------------------------------------------------------------------------
187 REAL :: zmax ! Max of 12 surrounding values
188 REAL :: zmin ! Min of 12 surrounding values
189 INTEGER :: jloop2 ! Dummy counter
190 INTEGER, DIMENSION(12) :: ip ! Array for IPn
191 INTEGER :: jlat ! latitude loop counter
192 INTEGER :: jlon ! longitude loop counter
193 REAL(KIND=JPRB) :: zhook_handle
194 !------------------------------------------------------------------------------
195 !
196 !* 1. DETERMINATION of the latitude of the poles (depending of the latitude
197 ! ------------- of the first data point)
198 !
199 IF (lhook) CALL dr_hook('HORIBL_SURF',0,zhook_handle)
200 IF (pila1>0.) THEN
201  zsouthpole= 90.
202  znorthpole=-90.
203 ELSE
204  zsouthpole=-90.
205  znorthpole= 90.
206 END IF
207 !
208 !------------------------------------------------------------------------------
209 !
210 !* 2. EXTEND DATA GRID
211 ! ----------------
212  ! Land / Sea mask
213 ldlsm = .false.
214 IF (present(klsmin) .AND. present(klsmout)) ldlsm = .true.
215 !
216 !* 2.1 Alias input data
217 !
218 zilo1 = pilo1
219 zilo2 = pilo2
220 zvect = 1.
221 IF (odvect) zvect=-1.
222 !
223 !* 2.2 Center input domain in order to have Lo1 < Lo 2
224 !
225 IF (zilo2 < 0.) zilo2 = zilo2 + 360.
226 IF (zilo1 < 0.) zilo1 = zilo1 + 360.
227 IF (zilo2 < zilo1) zilo1 = zilo1 - 360.
228 !
229 !* 2.3 Extend one point (needed for reduced grids)
230 !
231 ! Longitude coordinate of points are found by :
232 ! i
233 ! Lon(i) = Lon1 + ------------- . (Lon2 - Lon1)
234 ! Npts(Lat)-1
235 ! Where i goes from 0 to Npts(Lat)-1. The result of this is that the last point of
236 ! each parallel is located at Lon2. This is not the case for reduced grid where the
237 ! position of the last point depends upon the number of points of the parallel. For
238 ! reduced grid, the right formula to use is the following :
239 ! i
240 ! Lon(i) = Lon1 + ----------- . (Lon2' - Lon1)
241 ! Npts(Lat)
242 ! Where Lon2' = Lon1 + 2.PI.
243 !
244 ! Lon2 - Lon1
245 ! This can be generalized with Lon2' = Lon2 + -------------
246 ! Nptsmax - 1
247 !
248 jopos = maxval(kinlo(1:kinla))
249 zilo2 = zilo1 + (zilo2 - zilo1) * jopos / (jopos - 1.)
250 !
251 !
252 !* 2.4 Test if the input is global or partially global
253 !
254 ! Note that we must have a global map to make extension around the poles
255 gglobn = .false.
256 gglobs = .false.
257 ggloblon = .false.
258 IF (zilo2-360.>zilo1-1.e-3) ggloblon = .true.
259 zidla = (pila2 - pila1) / (kinla - 1)
260 IF ((pila1-zidla>= 90.) .OR. (pila1-zidla<=-90.)) gglobs=ggloblon
261 IF ((pila2+zidla>= 90.) .OR. (pila2+zidla<=-90.)) gglobn=ggloblon
262 ! Aladin case (input PILA2, PILO2 are in meters) no extension
263 IF ( pila2 > 100. ) THEN
264  gglobn = .false.
265  gglobs = .false.
266  ggloblon = .false.
267 END IF
268 !
269 !
270 !* 2.5 Compute the size of the resulting map
271 !
272 ibigsize = kilen
273 IF (gglobs ) ibigsize=ibigsize+(4+kinlo( 1))+(4+kinlo( 2))
274 IF (gglobn ) ibigsize=ibigsize+(4+kinlo(kinla))+(4+kinlo(kinla-1))
275 IF (ggloblon) ibigsize=ibigsize+ 4*kinla
276 !
277 !* 2.6 Compute the resulting map
278 !
279 ALLOCATE (zarin(ibigsize))
280 ALLOCATE (ilsmin(ibigsize))
281 !
282 ! 2.6.1 Compute the longitude extension
283 !
284 ! This is a basic copy of the data. If extension is possible, the first and last
285 ! two lines are copied twice this way :
286 !
287 ! /---------------\
288 ! | |
289 ! [.] [.] [.... ...] [.] [.]
290 ! | |
291 ! \------------/
292 !
293 ! A point represent a data.
294 !
295 jipos = 1
296 jopos = 1
297 IF (gglobs) jopos=jopos+(4+kinlo(1))+(4+kinlo(2))
298 IF (ggloblon) THEN
299  DO jloop1 = 1, kinla
300  zarin(jopos ) = parin(jipos+kinlo(jloop1)-2)
301  zarin(jopos+1) = parin(jipos+kinlo(jloop1)-1)
302  zarin(jopos+2:jopos+2+kinlo(jloop1)-1) = parin(jipos:jipos+kinlo(jloop1)-1)
303  zarin(jopos+2+kinlo(jloop1) ) = parin(jipos )
304  zarin(jopos+2+kinlo(jloop1)+1) = parin(jipos+1)
305  IF (ldlsm) THEN
306  ilsmin(jopos ) = klsmin(jipos+kinlo(jloop1)-2)
307  ilsmin(jopos+1) = klsmin(jipos+kinlo(jloop1)-1)
308  ilsmin(jopos+2:jopos+2+kinlo(jloop1)-1) = klsmin(jipos:jipos+kinlo(jloop1)-1)
309  ilsmin(jopos+2+kinlo(jloop1) ) = klsmin(jipos )
310  ilsmin(jopos+2+kinlo(jloop1)+1) = klsmin(jipos+1)
311  END IF
312  jipos = jipos + kinlo(jloop1)
313  jopos = jopos + kinlo(jloop1) + 4
314  END DO
315 ELSE
316  zarin(jopos:jopos+kilen-1) = parin(jipos:jipos+kilen-1)
317  IF (ldlsm) THEN
318  ilsmin(jopos:jopos+kilen-1) = klsmin(jipos:jipos+kilen-1)
319  END IF
320 END IF
321 !
322 ! 2.6.2 Compute the south pole extension
323 !
324 ! Pole extension is performed by copying the first half datas to the last half
325 ! datas of the extension parallel :
326 !
327 ! [.] [.] [....] [....] [.] [.]
328 ! ||||
329 ! /-------/
330 ! ||||
331 ! [.] [.] [....] [....] [.] [.]
332 !
333 IF (gglobs) THEN ! South pole (south meaning begining of the grib)
334  ioffset1 = 4 + kinlo(2)
335  ioffset2 = ioffset1 + 4 + kinlo(1)
336  imiddle = (kinlo(1)+4) / 2
337  zarin(ioffset1+1:ioffset1+imiddle) = &
338  zvect*zarin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
339  zarin(ioffset1+imiddle+1:ioffset1+kinlo(1)+4) = &
340  zvect*zarin(ioffset2+1+2:ioffset2+kinlo(1)+4-imiddle+2)
341  IF (ldlsm) THEN
342  ilsmin(ioffset1+1:ioffset1+imiddle) = &
343  ilsmin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
344  ilsmin(ioffset1+imiddle+1:ioffset1+kinlo(1)+4) = &
345  ilsmin(ioffset2+1+2:ioffset2+kinlo(1)+4-imiddle+2)
346  END IF
347  ioffset2 = ioffset2 + 4 + kinlo(1)
348  imiddle = (kinlo(2)+4) / 2
349  zarin(1:imiddle) = zvect*zarin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
350  zarin(imiddle+1:kinlo(2)+4) = &
351  zvect*zarin(ioffset2+1+2:ioffset2+kinlo(2)+4-imiddle+2)
352  IF (ldlsm) THEN
353  ilsmin(1:imiddle) = ilsmin(ioffset2+1+imiddle:ioffset2+2*imiddle)
354  ilsmin(imiddle+1:kinlo(2)+4) = ilsmin(ioffset2+1+2:ioffset2+kinlo(2)+4-imiddle+2)
355  END IF
356 END IF
357 !
358 ! 2.6.3 Compute the north pole extension
359 !
360 IF (gglobn) THEN ! North pole (north meaning end of the grib)
361  ioffset1 = ibigsize - (4+kinlo(kinla-1)) - (4+kinlo(kinla))
362  ioffset2 = ioffset1 - (4+kinlo(kinla))
363  imiddle = (kinlo(kinla)+4) / 2
364  zarin(ioffset1+1:ioffset1+imiddle) = &
365  zvect*zarin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
366  zarin(ioffset1+imiddle+1:ioffset1+kinlo(kinla)+4) = &
367  zvect*zarin(ioffset2+1+2:ioffset2+kinlo(kinla)+4-imiddle+2)
368  IF (ldlsm) THEN
369  ilsmin(ioffset1+1:ioffset1+imiddle) = &
370  ilsmin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
371  ilsmin(ioffset1+imiddle+1:ioffset1+kinlo(kinla)+4) = &
372  ilsmin(ioffset2+1+2:ioffset2+kinlo(kinla)+4-imiddle+2)
373  END IF
374  ioffset1 = ioffset1 + (4+kinlo(kinla))
375  ioffset2 = ioffset2 - (4+kinlo(kinla-1))
376  imiddle = (kinlo(kinla-1)+4) / 2
377  zarin(ioffset1+1:ioffset1+imiddle) = &
378  zvect*zarin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
379  zarin(ioffset1+imiddle+1:ioffset1+kinlo(kinla-1)+4) = &
380  zvect*zarin(ioffset2+1+2:ioffset2+kinlo(kinla-1)+4-imiddle+2)
381  IF (ldlsm) THEN
382  ilsmin(ioffset1+1:ioffset1+imiddle) = &
383  ilsmin(ioffset2+1+imiddle-2:ioffset2+2*imiddle-2)
384  ilsmin(ioffset1+imiddle+1:ioffset1+kinlo(kinla-1)+4) = &
385  ilsmin(ioffset2+1+2:ioffset2+kinlo(kinla-1)+4-imiddle+2)
386  END IF
387 END IF
388 !
389 !* 2.7 Compute the resulting parameters of the map
390 !
391 iinla = kinla
392 IF (gglobs) iinla = iinla + 2
393 IF (gglobn) iinla = iinla + 2
394 !
395 ALLOCATE (iinlo(iinla))
396 ioffset1 = 0
397 IF (gglobs) THEN
398  iinlo(ioffset1+1) = kinlo(2)
399  iinlo(ioffset1+2) = kinlo(1)
400  ioffset1 = ioffset1 + 2
401 END IF
402 iinlo(ioffset1+1:ioffset1+kinla) = kinlo(1:kinla)
403 ioffset1 = ioffset1 + kinla
404 IF (gglobn) THEN
405  iinlo(ioffset1+1) = kinlo(kinla)
406  iinlo(ioffset1+2) = kinlo(kinla-1)
407  ioffset1 = ioffset1 + 2
408 END IF
409 !
410 !* 2.8 Compute Offset array used to acces the lines
411 !
412 ALLOCATE (iofs(iinla))
413 iofs(1) = 1
414 IF (ggloblon) iofs(1)=iofs(1)+2
415 DO jloop1=2, iinla
416  iofs(jloop1) = iofs(jloop1-1) + iinlo(jloop1-1)
417  IF (ggloblon) iofs(jloop1) = iofs(jloop1) + 4
418 END DO
419 !
420 !------------------------------------------------------------------------------
421 !
422 !* 3. LOOP OVER ALL THE POINTS OF THE OUTPUT GRID
423 ! -------------------------------------------
424 !
425 parout(:) = xundef
426 !
427 jopos = 0
428 DO jloop1 = 1, kolen
429  jopos = jopos + 1
430  IF (.NOT. ointerp(jopos)) cycle
431  zola = pyout(jopos)
432  zolo = pxout(jopos)
433  IF (zolo < zilo1) zolo = zolo + 360.
434  IF (zolo > zilo2) zolo = zolo - 360.
435 !
436 !* 3.1 Locate the 12 input points around the interpolated output point
437 !*
438 !* N 1 2
439 !*
440 !* ^ 3 4 5 6
441 !* | x
442 !* | 7 8 9 10
443 !* |
444 !* 11 12
445 !* S
446 !* W ---------------> E
447 !*
448 !* Note : the name 'south', 'north', may not be exact if the point 2 is
449 !* to the south of point 1 (IDLA < 0). This does not affect computation.
450 !
451  ! 3.1.1. find positions of latitudes
452  ios = nint( (zola-pila1)/zidla -0.5) ! because of the zero
453  IF ( ios<-1 ) THEN
454  CALL abor1_sfx('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE')
455  END IF
456  zlas = pila1 + ios * zidla
457  ios = ios + 1
458  IF (gglobs) THEN
459  ios = ios + 2
460  ioss = ios - 1
461  ion = ios + 1
462  ionn = ion + 1
463  ELSE
464  ioss = max(ios - 1,1)
465  ion = min(ios + 1,iinla)
466  ionn = min(ion + 1,iinla)
467  ios = max(ios, 1)
468  ENDIF
469  zlass = zlas - zidla
470  zlan = zlas + zidla
471  zlann = zlan + zidla
472  !
473  ! extra latitudes are computed symetrically compared to the poles
474  !
475  IF (gglobs .AND. ios==2) THEN
476  zlass = 2. * zsouthpole - zlann
477  zlas = 2. * zsouthpole - zlan
478  END IF
479  IF (gglobs .AND. ios==3) THEN
480  zlass = 2. * zsouthpole - zlas
481  END IF
482  IF (gglobn .AND. ios==iinla-2) THEN
483  zlann = 2. * znorthpole - zlass
484  zlan = 2. * znorthpole - zlas
485  END IF
486  IF (gglobn .AND. ios==iinla-3) THEN
487  zlann = 2. * znorthpole - zlan
488  END IF
489 !
490  IF ((ioss<1).OR.(ionn>iinla).OR. &
491  (ioss<1).OR.(ionn>iinla)) THEN
492  CALL abor1_sfx('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE')
493  END IF
494 !
495  ! 3.1.2. northern
496  zidlo = (zilo2 - zilo1) / (iinlo(ionn))
497  ip1 = int((zolo - zilo1) / zidlo)
498  IF (ggloblon) THEN
499  ip2 = ip1 + 1
500  ELSE
501  ip2 = min(iinlo(ionn)-1, ip1+1)
502  ENDIF
503  zlop1 = zilo1 + ip1 * zidlo
504  zlop2 = zlop1 + zidlo
505 !
506  ! 3.1.3. north
507  zidlo = (zilo2 - zilo1) / (iinlo(ion ))
508  ip4 = int((zolo - zilo1) / zidlo)
509  IF (ggloblon) THEN
510  ip3 = ip4 - 1
511  ip5 = ip4 + 1
512  ip6 = ip5 + 1
513  ELSE
514  ip3 = max(0,ip4-1)
515  ip5 = min(iinlo(ion)-1,ip4+1)
516  ip6 = min(iinlo(ion)-1,ip5+1)
517  ENDIF
518  zlop4 = zilo1 + ip4 * zidlo
519  zlop3 = zlop4 - zidlo
520  zlop5 = zlop4 + zidlo
521  zlop6 = zlop5 + zidlo
522 !
523  ! 3.1.4. south
524  zidlo = (zilo2 - zilo1) / (iinlo(ios ))
525  ip8 = int((zolo - zilo1) / zidlo)
526  IF (ggloblon) THEN
527  ip7 = ip8 - 1
528  ip9 = ip8 + 1
529  ip10 = ip9 + 1
530  ELSE
531  ip7 = max(0,ip8-1)
532  ip9 = min(iinlo(ios)-1,ip8+1)
533  ip10 = min(iinlo(ios)-1,ip9+1)
534  ENDIF
535  zlop8 = zilo1 + ip8 * zidlo
536  zlop7 = zlop8 - zidlo
537  zlop9 = zlop8 + zidlo
538  zlop10= zlop9 + zidlo
539 !
540  ! 3.1.5. southern
541  zidlo = (zilo2 - zilo1) / (iinlo(ioss))
542  ip11 = int((zolo - zilo1) / zidlo)
543  IF (ggloblon) THEN
544  ip12 = ip11 + 1
545  ELSE
546  ip12 = min(iinlo(ioss)-1,ip11+1)
547  ENDIF
548  zlop11= zilo1 + ip11* zidlo
549  zlop12= zlop11+ zidlo
550 !
551  ! 3.1.6. check position of points
552  IF (ggloblon ) THEN
553  IF ((ip1 <-2) .OR. (ip2 >iinlo(ionn)+1) .OR. &
554  (ip3 <-2) .OR. (ip6 >iinlo(ion )+1) .OR. &
555  (ip7 <-2) .OR. (ip10>iinlo(ios )+1) .OR. &
556  (ip11<-2) .OR. (ip12>iinlo(ioss)+1)) THEN
557  CALL abor1_sfx('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LONGITUDE GLOBAL')
558  END IF
559  ELSE
560  IF ((ip1 <0) .OR. (ip2 >iinlo(ionn)-1) .OR. &
561  (ip3 <0) .OR. (ip6 >iinlo(ion )-1) .OR. &
562  (ip7 <0) .OR. (ip10>iinlo(ios )-1) .OR. &
563  (ip11<0) .OR. (ip12>iinlo(ioss)-1)) THEN
564  CALL abor1_sfx('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LONGITUDE LOCAL')
565  END IF
566  END IF
567 !
568  ! 3.1.7. add parallel offset
569  ip1 =ip1 + iofs(ionn)
570  ip2 =ip2 + iofs(ionn)
571  ip3 =ip3 + iofs(ion )
572  ip4 =ip4 + iofs(ion )
573  ip5 =ip5 + iofs(ion )
574  ip6 =ip6 + iofs(ion )
575  ip7 =ip7 + iofs(ios )
576  ip8 =ip8 + iofs(ios )
577  ip9 =ip9 + iofs(ios )
578  ip10=ip10+ iofs(ios )
579  ip11=ip11+ iofs(ioss)
580  ip12=ip12+ iofs(ioss)
581 !
582 !* 3.2 Land / Sea mask
583 !
584  zlsm1 = 1.
585  zlsm2 = 1.
586  zlsm3 = 1.
587  zlsm4 = 1.
588  zlsm5 = 1.
589  zlsm6 = 1.
590  zlsm7 = 1.
591  zlsm8 = 1.
592  zlsm9 = 1.
593  zlsm10 = 1.
594  zlsm11 = 1.
595  zlsm12 = 1.
596  zlsmnn = 1.
597  zlsmn = 1.
598  zlsms = 1.
599  zlsmss = 1.
600  IF (ldlsm) THEN
601  IF (ilsmin(ip1 ).NE.klsmout(jopos)) zlsm1 = 0.
602  IF (ilsmin(ip2 ).NE.klsmout(jopos)) zlsm2 = 0.
603  IF (ilsmin(ip3 ).NE.klsmout(jopos)) zlsm3 = 0.
604  IF (ilsmin(ip4 ).NE.klsmout(jopos)) zlsm4 = 0.
605  IF (ilsmin(ip5 ).NE.klsmout(jopos)) zlsm5 = 0.
606  IF (ilsmin(ip6 ).NE.klsmout(jopos)) zlsm6 = 0.
607  IF (ilsmin(ip7 ).NE.klsmout(jopos)) zlsm7 = 0.
608  IF (ilsmin(ip8 ).NE.klsmout(jopos)) zlsm8 = 0.
609  IF (ilsmin(ip9 ).NE.klsmout(jopos)) zlsm9 = 0.
610  IF (ilsmin(ip10).NE.klsmout(jopos)) zlsm10 = 0.
611  IF (ilsmin(ip11).NE.klsmout(jopos)) zlsm11 = 0.
612  IF (ilsmin(ip12).NE.klsmout(jopos)) zlsm12 = 0.
613  zlsmnn = min(zlsm1 +zlsm2,1.)
614  zlsmn = min(zlsm3 +zlsm4 +zlsm5 +zlsm6,1.)
615  zlsms = min(zlsm7 +zlsm8 +zlsm9 +zlsm10,1.)
616  zlsmss = min(zlsm11+zlsm12,1.)
617  zlsmtot = min(zlsmnn+zlsmn+zlsms+zlsmss,1.)
618  IF (zlsmnn < 1.e-3) THEN
619  zlsm1 = 1.
620  zlsm2 = 1.
621  END IF
622  IF (zlsmn < 1.e-3) THEN
623  zlsm3 = 1.
624  zlsm4 = 1.
625  zlsm5 = 1.
626  zlsm6 = 1.
627  END IF
628  IF (zlsms < 1.e-3) THEN
629  zlsm7 = 1.
630  zlsm8 = 1.
631  zlsm9 = 1.
632  zlsm10= 1.
633  END IF
634  IF (zlsmss < 1.e-3) THEN
635  zlsm11= 1.
636  zlsm12= 1.
637  END IF
638  IF (zlsmtot < 1.e-3) THEN
639  zlsmnn = 1.
640  zlsmn = 1.
641  zlsms = 1.
642  zlsmss = 1.
643  END IF
644  ENDIF
645 !
646 !* 3.3 Weight of points
647 !
648  ! 3.3.1 northern
649  zw1 = zlsm1 * (1.+zlsm2 *(zolo -zlop1 )/(zlop1 -zlop2 ))
650  zw2 = 1. - zw1
651  zwnn = zlsmnn* (1.+zlsmn *(zola -zlann)/(zlann-zlan )) &
652  * (1.+zlsms *(zola -zlann)/(zlann-zlas )) &
653  * (1.+zlsmss*(zola -zlann)/(zlann-zlass))
654 !
655  ! 3.3.2. north
656  zw3 = zlsm3 * (1.+zlsm4 *(zolo -zlop3 )/(zlop3 -zlop4 )) &
657  * (1.+zlsm5 *(zolo -zlop3 )/(zlop3 -zlop5 )) &
658  * (1.+zlsm6 *(zolo -zlop3 )/(zlop3 -zlop6 ))
659  zw4 = zlsm4 * (1.+zlsm3 *(zolo -zlop4 )/(zlop4 -zlop3 )) &
660  * (1.+zlsm5 *(zolo -zlop4 )/(zlop4 -zlop5 )) &
661  * (1.+zlsm6 *(zolo -zlop4 )/(zlop4 -zlop6 ))
662  zw5 = zlsm5 * (1.+zlsm3 *(zolo -zlop5 )/(zlop5 -zlop3 )) &
663  * (1.+zlsm4 *(zolo -zlop5 )/(zlop5 -zlop4 )) &
664  * (1.+zlsm6 *(zolo -zlop5 )/(zlop5 -zlop6 ))
665  zw6 = 1. - zw3 - zw4 - zw5
666  zwn = zlsmn * (1.+zlsmnn*(zola -zlan )/(zlan -zlann)) &
667  * (1.+zlsms *(zola -zlan )/(zlan -zlas )) &
668  * (1.+zlsmss*(zola -zlan )/(zlan -zlass))
669 !
670  ! 3.3.3. south
671  zw7 = zlsm7 * (1.+zlsm8 *(zolo -zlop7 )/(zlop7 -zlop8 )) &
672  * (1.+zlsm9 *(zolo -zlop7 )/(zlop7 -zlop9 )) &
673  * (1.+zlsm10*(zolo -zlop7 )/(zlop7 -zlop10))
674  zw8 = zlsm8 * (1.+zlsm7 *(zolo -zlop8 )/(zlop8 -zlop7 )) &
675  * (1.+zlsm9 *(zolo -zlop8 )/(zlop8 -zlop9 )) &
676  * (1.+zlsm10*(zolo -zlop8 )/(zlop8 -zlop10))
677  zw9 = zlsm9 * (1.+zlsm7 *(zolo -zlop9 )/(zlop9 -zlop7 )) &
678  * (1.+zlsm8 *(zolo -zlop9 )/(zlop9 -zlop8 )) &
679  * (1.+zlsm10*(zolo -zlop9 )/(zlop9 -zlop10))
680  zw10 = 1. - zw7 - zw8 - zw9
681  zws = zlsms * (1.+zlsmnn*(zola -zlas )/(zlas -zlann)) &
682  * (1.+zlsmn *(zola -zlas )/(zlas -zlan )) &
683  * (1.+zlsmss*(zola -zlas )/(zlas -zlass))
684 !
685  ! 3.3.4. southern
686  zw11 = zlsm11* (1.+zlsm12*(zolo -zlop11)/(zlop11-zlop12))
687  zw12 = 1. - zw11
688  zwss = 1. - zwnn - zwn - zws
689 !
690 ! In order to exclude undef values from computation of PAROUT,
691 ! weights w2, w6, w10, w12 and wss which can be numerically very low
692 ! because they are residual, are set to 0
693 !
694 IF (abs(zw2 ) < 1.e-10) zw2=0.
695 IF (abs(zw6 ) < 1.e-10) zw6=0.
696 IF (abs(zw10) < 1.e-10) zw10=0.
697 IF (abs(zw12) < 1.e-10) zw12=0.
698 IF (abs(zwss) < 1.e-10) zwss=0.
699 !
700  ! 3.3.5. longitude weight x latitude weight
701  zw1 = zw1 * zwnn
702  zw2 = zw2 * zwnn
703  zw3 = zw3 * zwn
704  zw4 = zw4 * zwn
705  zw5 = zw5 * zwn
706  zw6 = zw6 * zwn
707  zw7 = zw7 * zws
708  zw8 = zw8 * zws
709  zw9 = zw9 * zws
710  zw10 = zw10 * zws
711  zw11 = zw11 * zwss
712  zw12 = zw12 * zwss
713 !
714  parout(jopos) = zw1 * zarin(ip1 ) + &
715  zw2 * zarin(ip2 ) + &
716  zw3 * zarin(ip3 ) + &
717  zw4 * zarin(ip4 ) + &
718  zw5 * zarin(ip5 ) + &
719  zw6 * zarin(ip6 ) + &
720  zw7 * zarin(ip7 ) + &
721  zw8 * zarin(ip8 ) + &
722  zw9 * zarin(ip9 ) + &
723  zw10 * zarin(ip10) + &
724  zw11 * zarin(ip11) + &
725  zw12 * zarin(ip12)
726 !
727 ! For surface fields, the interpoalted value is bounded
728 ! by the min max values of the initial field
729 
730  IF (present(klsmin)) THEN
731 
732  ip(1)=ip1
733  ip(2)=ip2
734  ip(3)=ip3
735  ip(4)=ip4
736  ip(5)=ip5
737  ip(6)=ip6
738  ip(7)=ip7
739  ip(8)=ip8
740  ip(9)=ip9
741  ip(10)=ip10
742  ip(11)=ip11
743  ip(12)=ip12
744 
745  zmin=xundef
746  zmax=xundef
747 
748  DO jloop2=1,12
749  IF (zarin(ip(jloop2))==xundef) cycle
750 
751  IF ((zmax==xundef)) THEN
752  zmax=zarin(ip(jloop2))
753  zmin=zarin(ip(jloop2))
754  ELSE
755  zmax=max(zmax,zarin(ip(jloop2)))
756  zmin=min(zmin,zarin(ip(jloop2)))
757  ENDIF
758 
759  END DO
760 
761  parout(jopos) = max(min(parout(jopos),zmax),zmin)
762 
763  ENDIF
764 
765 END DO
766 !
767 DEALLOCATE (iinlo)
768 DEALLOCATE (zarin)
769 DEALLOCATE (ilsmin)
770 DEALLOCATE (iofs)
771 !
772 WHERE(abs(parout-xundef)<1.e-6) parout=xundef
773 !
774 !------------------------------------------------------------------------------
775 !
776 !* 4. EXTRAPOLATIONS IF SOME POINTS WERE NOT INTERPOLATED
777 ! ---------------------------------------------------
778 !
779 !* no missing point
780 IF (count(parout(:)==xundef .AND. ointerp(:))==0 .AND. lhook) CALL dr_hook('HORIBL_SURF',1,zhook_handle)
781 IF (count(parout(:)==xundef .AND. ointerp(:))==0) RETURN
782 !
783 !* no data point
784 IF (count(parin(:)/=xundef)==0 .AND. lhook) CALL dr_hook('HORIBL_SURF',1,zhook_handle)
785 IF (count(parin(:)/=xundef)==0) RETURN
786 !
787 WRITE(kluout,*) ' Remaining horizontal extrapolations'
788 WRITE(kluout,*) ' Total number of input data : ',count(parin(:)/=xundef)
789 WRITE(kluout,*) ' Number of points to interpolate: ',count(parout(:)==xundef .AND. ointerp(:))
790 !
791 !* input grid coordinates
792 !
793 ALLOCATE(zla(kilen))
794 ALLOCATE(zlo(kilen))
795 !
796 jipos = 0
797 zidla = (pila2 - pila1) / (kinla - 1)
798 DO jlat=1,kinla
799  zidlo = (zilo2-zilo1) / kinlo(jlat)
800  DO jlon=1,kinlo(jlat)
801  jipos = jipos + 1
802  zla(jipos) = pila1 + (jlat-1) * zidla
803  zlo(jipos) = zilo1 + (jlon-1) * zidlo
804  END DO
805 END DO
806 !
807  CALL hor_extrapol_surf(kluout,'LALO',zla,zlo,parin,pyout,pxout,parout,ointerp)
808 !
809 DEALLOCATE(zla)
810 DEALLOCATE(zlo)
811 IF (lhook) CALL dr_hook('HORIBL_SURF',1,zhook_handle)
812 !
813 !------------------------------------------------------------------------------
814 !
815 !
816 END SUBROUTINE horibl_surf
horibl_surf
subroutine horibl_surf(PILA1, PILO1, PILA2, PILO2, KINLA, KINLO, KILEN, PARIN, KOLEN, PXOUT, PYOUT, PAROUT, ODVECT, KLUOUT, OINTERP, KLSMIN, KLSMOUT)
Definition: horibl_surf.F90:6
abor1_sfx
subroutine abor1_sfx(YTEXT)
Definition: abor1_sfx.F90:6
hor_extrapol_surf
subroutine hor_extrapol_surf(KLUOUT, HCOORTYPE, PLAT_IN, PLON_IN, PFIELD_IN, PLAT, PLON, PFIELD, OINTERP)
Definition: hor_extrapol_surf.F90:6
modd_surf_par
Definition: modd_surf_par.F90:6
Generated on Fri May 13 2016 09:31:04 for SURFEX by   doxygen 1.8.5