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