source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90 @ 11202

Last change on this file since 11202 was 11202, checked in by jcastill, 15 months ago

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

  • Property svn:keywords set to Id
File size: 48.4 KB
Line 
1   !!----------------------------------------------------------------------
2   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
3   !! $Id$
4   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
5   !!----------------------------------------------------------------------
6
7   SUBROUTINE obs_int_h2d_init( kpk,   kpk2,  k2dint, plam,  pphi, &
8      &                         pglam, pgphi, pmask,  pweig, pobsmask, &
9      &                         iminpoints )
10      !!-----------------------------------------------------------------------
11      !!
12      !!                     ***  ROUTINE obs_int_h2d  ***
13      !!
14      !! ** Purpose : Computes weights for horizontal interpolation to the
15      !!              observation point.
16      !!
17      !! ** Method  : Horizontal interpolation to the observation point using
18      !!              model values at the corners of the surrounding grid
19      !!              points.
20      !!
21      !!    Interpolation Schemes :
22      !!
23      !!    1) k2dint = 0: Distance-weighted interpolation scheme 1
24      !!
25      !!       The interpolation weights are computed as a weighted
26      !!       sum of the distance between the model grid points (A)
27      !!       and the observation point (B). Distance (s) is computed
28      !!       using the great-circle distance formula:
29      !!
30      !!       s(AB) = arcos(   sin( phiA ) x sin( phiB )
31      !!                      + cos( phiA ) x cos( phiB )
32      !!                                    x cos( lamB - lamA ) )
33      !!
34      !!    2) k2dint = 1: Distance-weighted interpolation scheme 2
35      !!
36      !!       As k2dint = 0 but with distance (ds) computed using
37      !!       a small-angle approximation to the great-circle formula:
38      !!
39      !!       ds(AB) = sqrt(   ( phiB - phiA )^{2}
40      !!                    + ( ( lamB - lamA ) * cos( phiB ) )^{2} )
41      !!
42      !!    3) k2dint = 2: Bilinear interpolation on a geographical grid
43      !!
44      !!       The interpolation is split into two 1D interpolations in
45      !!       the longitude and latitude directions, respectively.
46      !!
47      !!    4) k2dint = 3: General bilinear remapping interpolation
48      !!
49      !!       An iterative scheme that involves first mapping a
50      !!       quadrilateral cell into a cell with coordinates
51      !!       (0,0), (1,0), (0,1) and (1,1).
52      !!
53      !!    5) k2dint = 4: Polynomial interpolation
54      !!
55      !!       The interpolation weights are computed by fitting a
56      !!       polynomial function of the form
57      !!             
58      !!       P(i) = a1(i) + a2(i) * phi + a3(i) * plam
59      !!                                  + a4(i) * phi * plam
60      !!   
61      !!       through the model values at the four surrounding grid points.
62      !!
63      !! ** Action  :
64      !!                   
65      !! References : Jones, P.: A users guide for SCRIP: A Spherical
66      !!                        Coordinate Remapping and Interpolation Package.
67      !!                        Version 1.4. Los Alomos.
68      !!                   
69      !!       http://www.acl.lanl.gov/climate/software/SCRIP/SCRIPmain.html
70      !!
71      !! History :
72      !!        ! 97-11 (A. Weaver, N. Daget)
73      !!        ! 06-03 (A. Vidard) NEMOVAR migration
74      !!        ! 06-10 (A. Weaver) Cleanup
75      !!        ! 07-08 (K. Mogensen) Split in two routines for easier adj.
76      !!-----------------------------------------------------------------------
77      !! * Modules used
78      !! * Arguments
79      INTEGER, INTENT(IN) :: &
80         & kpk,   &             ! Parameter values for automatic arrays
81         & kpk2,  &       
82         & k2dint               ! Interpolation scheme options
83                                ! = 0 distance-weighted (great circle)
84                                ! = 1 distance-weighted (small angle)
85                                ! = 2 bilinear (geographical grid)
86                                ! = 3 bilinear (quadrilateral grid)
87                                ! = 4 polynomial (quadrilateral grid)
88      REAL(KIND=wp), INTENT(INOUT) :: &
89         & plam, &
90         & pphi                 ! Geographical (lat,lon) coordinates of
91                                ! observation
92      REAL(KIND=wp), DIMENSION(2,2), INTENT(IN) :: &
93         & pglam, &             ! Model variable lat
94         & pgphi                ! Model variable lon
95      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
96         & pmask                ! Model variable mask
97      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(OUT) ::  &
98         & pweig                ! Weights for interpolation
99      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
100         & pobsmask             ! Vertical mask for observations
101      INTEGER, INTENT(IN), OPTIONAL :: &
102         & iminpoints           ! Reject point which is not surrounded
103                                ! by at least iminpoints sea points
104
105      !! * Local declarations
106      INTEGER :: &
107         & jk
108      INTEGER :: &
109         & ikmax, &
110         & iamb1, &
111         & iamb2
112      REAL(KIND=wp) :: &
113         & zphimm,  &
114         & zphimp,  &
115         & zphipm,  &
116         & zphipp,  &
117         & zlammm,  &
118         & zlammp,  &
119         & zlampm,  &
120         & zlampp,  &
121         & zphimin, &
122         & zphimax, &
123         & zlammin, &
124         & zlammax
125      REAL(KIND=wp), DIMENSION(kpk2) :: &
126         & z2dmm,  &
127         & z2dmp,  &
128         & z2dpm,  &
129         & z2dpp,  &
130         & z2dmmt, &
131         & z2dmpt, &
132         & z2dpmt, &
133         & z2dppt, &
134         & zsum
135      LOGICAL :: &
136         & ll_ds1, &
137         & ll_skip, &
138         & ll_fail
139       
140      !------------------------------------------------------------------------
141      ! Constants for the 360 degrees ambiguity
142      !------------------------------------------------------------------------
143      iamb1 = 10   ! dlam < iamb1 * dphi
144      iamb2 = 3    ! Special treatment if iamb2 * lam < max(lam)
145       
146      !------------------------------------------------------------------------
147      ! Initialize number of levels
148      !------------------------------------------------------------------------
149      IF ( kpk2 == 1 ) THEN
150         ikmax = 1
151      ELSEIF ( kpk2 == kpk) THEN
152         ikmax = kpk-1
153      ENDIF
154      !------------------------------------------------------------------------
155      ! Initialize the cell corners
156      !------------------------------------------------------------------------
157      zphimm = pgphi(1,1)
158      zphimp = pgphi(1,2)
159      zphipm = pgphi(2,1)
160      zphipp = pgphi(2,2)
161      zlammm = pglam(1,1)
162      zlammp = pglam(1,2)
163      zlampm = pglam(2,1)
164      zlampp = pglam(2,2)
165       
166      !------------------------------------------------------------------------
167      ! Treat the 360 degrees ambiguity
168      !------------------------------------------------------------------------
169      DO WHILE ( ( zlammm < 0.0_wp ).OR.( zlammm > 360.0_wp )  &
170         &   .OR.( zlampm < 0.0_wp ).OR.( zlampm > 360.0_wp )  &
171         &   .OR.( zlampp < 0.0_wp ).OR.( zlampp > 360.0_wp )  &
172         &   .OR.( zlammp < 0.0_wp ).OR.( zlammp > 360.0_wp ) )
173       
174         IF ( zlammm < 0.0_wp   ) zlammm = zlammm + 360.0_wp
175         IF ( zlammm > 360.0_wp ) zlammm = zlammm - 360.0_wp
176         IF ( zlammp < 0.0_wp   ) zlammp = zlammp + 360.0_wp
177         IF ( zlammp > 360.0_wp ) zlammp = zlammp - 360.0_wp
178         IF ( zlampm < 0.0_wp   ) zlampm = zlampm + 360.0_wp
179         IF ( zlampm > 360.0_wp ) zlampm = zlampm - 360.0_wp
180         IF ( zlampp < 0.0_wp   ) zlampp = zlampp + 360.0_wp
181         IF ( zlampp > 360.0_wp ) zlampp = zlampp - 360.0_wp
182       
183      END DO
184
185      DO WHILE ( ( plam < 0.0_wp ) .OR. ( plam > 360.0_wp ) )
186         IF ( plam < 0.0_wp   ) plam = plam + 360.0_wp
187         IF ( plam > 360.0_wp ) plam = plam - 360.0_wp
188      END DO
189
190      !------------------------------------------------------------------------
191      ! Special case for observation on grid points
192      !------------------------------------------------------------------------
193      ll_skip = .FALSE.
194      IF ( ( ABS( zphimm - pphi ) < 1.0e-6_wp ) .AND. &
195         & ( ABS( zlammm - plam ) < 1.0e-6_wp ) ) THEN
196         z2dmm(:) = 1.0_wp
197         z2dpm(:) = 0.0_wp
198         z2dmp(:) = 0.0_wp
199         z2dpp(:) = 0.0_wp
200         ll_skip = .TRUE.
201      ENDIF
202      IF ( ( ABS( zphipm - pphi ) < 1.0e-6_wp ) .AND. &
203         & ( ABS( zlampm - plam ) < 1.0e-6_wp ) ) THEN
204         z2dmm(:) = 0.0_wp
205         z2dpm(:) = 1.0_wp
206         z2dmp(:) = 0.0_wp
207         z2dpp(:) = 0.0_wp
208         ll_skip = .TRUE.
209      ENDIF
210      IF ( ( ABS( zphimp - pphi ) < 1.0e-6_wp ) .AND. &
211         & ( ABS( zlammp - plam ) < 1.0e-6_wp ) ) THEN
212         z2dmm(:) = 0.0_wp
213         z2dpm(:) = 0.0_wp
214         z2dmp(:) = 1.0_wp
215         z2dpp(:) = 0.0_wp
216         ll_skip = .TRUE.
217      ENDIF
218      IF ( ( ABS( zphipp - pphi ) < 1.0e-6_wp ) .AND. &
219         & ( ABS( zlampp - plam ) < 1.0e-6_wp ) ) THEN
220         z2dmm(:) = 0.0_wp
221         z2dpm(:) = 0.0_wp
222         z2dmp(:) = 0.0_wp
223         z2dpp(:) = 1.0_wp
224         ll_skip = .TRUE.
225      ENDIF
226
227      IF ( .NOT.ll_skip ) THEN
228
229         zphimin = MIN( zphimm, zphipm, zphipp, zphimp )
230         zphimax = MAX( zphimm, zphipm, zphipp, zphimp )
231         zlammin = MIN( zlammm, zlampm, zlampp, zlammp )
232         zlammax = MAX( zlammm, zlampm, zlampp, zlammp )
233
234         IF ( ( ( zlammax - zlammin ) / ( zphimax - zphimin ) ) > iamb1 ) THEN
235            IF ( iamb2 * zlammm < zlammax ) zlammm = zlammm + 360.0_wp
236            IF ( iamb2 * zlammp < zlammax ) zlammp = zlammp + 360.0_wp
237            IF ( iamb2 * zlampm < zlammax ) zlampm = zlampm + 360.0_wp
238            IF ( iamb2 * zlampp < zlammax ) zlampp = zlampp + 360.0_wp
239         ENDIF
240
241         zlammin = MIN( zlammm, zlampm, zlampp, zlammp )
242         IF ( zlammm > ( zlammin + 180.0_wp ) ) zlammm = zlammm - 360.0_wp
243         IF ( zlammp > ( zlammin + 180.0_wp ) ) zlammp = zlammp - 360.0_wp
244         IF ( zlampm > ( zlammin + 180.0_wp ) ) zlampm = zlampm - 360.0_wp
245         IF ( zlampp > ( zlammin + 180.0_wp ) ) zlampp = zlampp - 360.0_wp
246         
247         IF ( plam < zlammin ) plam = plam + 360.0_wp
248      z2dmm = 0.0_wp
249      z2dmp = 0.0_wp
250      z2dpm = 0.0_wp
251      z2dpp = 0.0_wp
252         SELECT CASE (k2dint)
253           
254         CASE(0)
255            CALL obs_int_h2d_ds1( kpk2, ikmax,                        &
256               &                  pphi, plam, pmask,                  &
257               &                  zphimm, zlammm, zphimp, zlammp,     &
258               &                  zphipm, zlampm, zphipp, zlampp,     &
259               &                  z2dmm, z2dmp, z2dpm, z2dpp )
260         CASE(1)
261            CALL obs_int_h2d_ds2( kpk2, ikmax,                        &
262               &                  pphi, plam, pmask,                  &
263               &                  zphimm, zlammm, zphimp, zlammp,     &
264               &                  zphipm, zlampm, zphipp, zlampp,     &
265               &                  z2dmm, z2dmp, z2dpm, z2dpp )
266         CASE(2)
267            CALL obs_int_h2d_bil( kpk2, ikmax,                        &
268               &                  pphi, plam, pmask,                  &
269               &                                          zlammp,     &
270               &                  zphipm,         zphipp, zlampp,     &
271               &                  z2dmm, z2dmp, z2dpm, z2dpp )
272         CASE(3)
273            CALL obs_int_h2d_bir( kpk2, ikmax,                        &
274               &                  pphi, plam, pmask,                  &
275               &                  zphimm, zlammm, zphimp, zlammp,     &
276               &                  zphipm, zlampm, zphipp, zlampp,     &
277               &                  z2dmm, z2dmp, z2dpm, z2dpp, ll_fail )
278            IF (ll_fail) THEN
279               IF(lwp) THEN
280                  WRITE(numout,*)'Bilinear weight computation failed'
281                  WRITE(numout,*)'Switching to great circle distance'
282                  WRITE(numout,*)
283               ENDIF
284               CALL obs_int_h2d_ds1( kpk2, ikmax,                        &
285                  &                  pphi, plam, pmask,                  &
286                  &                  zphimm, zlammm, zphimp, zlammp,     &
287                  &                  zphipm, zlampm, zphipp, zlampp,     &
288                  &                  z2dmm, z2dmp, z2dpm, z2dpp )
289            ENDIF
290         CASE(4) 
291            CALL obs_int_h2d_pol( kpk2, ikmax,                        &
292               &                  pphi, plam, pmask,                  &
293               &                  zphimm, zlammm, zphimp, zlammp,     &
294               &                  zphipm, zlampm, zphipp, zlampp,     &
295               &                  z2dmm, z2dmp, z2dpm, z2dpp )
296         END SELECT
297
298      ENDIF
299      !------------------------------------------------------------------------
300      ! Compute weights for interpolation to the observation point
301      !------------------------------------------------------------------------
302      pobsmask(:) = 0.0_wp
303      pweig(:,:,:) = 0.0_wp
304      ! ll_ds1 is used for failed interpolations
305      ll_ds1 = .FALSE.
306      DO jk = 1, ikmax
307         IF (PRESENT(iminpoints)) THEN
308            IF (NINT(SUM(pmask(:,:,jk)))<iminpoints) CYCLE
309         ENDIF
310         zsum(jk) = z2dmm(jk) + z2dmp(jk) + z2dpm(jk) + z2dpp(jk)
311         IF ( zsum(jk) /= 0.0_wp ) THEN
312            pweig(1,1,jk) = z2dmm(jk)
313            pweig(1,2,jk) = z2dmp(jk)
314            pweig(2,1,jk) = z2dpm(jk)
315            pweig(2,2,jk) = z2dpp(jk)
316            ! Set the vertical mask
317            IF ( ( ( z2dmm(jk) > 0.0_wp ) .AND.       &
318               &   ( pmask(1,1,jk) == 1.0_wp ) ) .OR. &
319               & ( ( z2dmp(jk) > 0.0_wp ) .AND.       &
320               &   ( pmask(1,2,jk) == 1.0_wp ) ) .OR. &
321               & ( ( z2dpm(jk) > 0.0_wp ) .AND.       &
322               &   ( pmask(2,1,jk) == 1.0_wp ) ) .OR. &
323               & ( ( z2dpp(jk) > 0.0_wp ) .AND.       &
324               &   ( pmask(2,2,jk) == 1.0_wp ) ) ) pobsmask(jk)=1.0_wp
325         ELSE
326            ! If the interpolation has failed due to the point
327            ! being on the intersect of two land points retry with
328            ! k2dint = 0
329            IF ( ( pmask(1,1,jk) /= 0.0_wp ).OR. &
330               & ( pmask(1,2,jk) /= 0.0_wp ).OR. &
331               & ( pmask(2,1,jk) /= 0.0_wp ).OR. &
332               & ( pmask(2,2,jk) /= 0.0_wp ) ) THEN
333               ! If ll_ds1 is false compute k2dint = 0 weights
334               IF ( .NOT.ll_ds1 ) THEN
335                  CALL obs_int_h2d_ds1( kpk2, ikmax,                        &
336                     &                  pphi, plam, pmask,                  &
337                     &                  zphimm, zlammm, zphimp, zlammp,     &
338                     &                  zphipm, zlampm, zphipp, zlampp,     &
339                     &                  z2dmmt, z2dmpt, z2dpmt, z2dppt )
340                  ll_ds1 = .TRUE.
341               ENDIF
342               zsum(jk) = z2dmmt(jk) + z2dmpt(jk) + z2dpmt(jk) + z2dppt(jk)
343               IF ( zsum(jk) /= 0.0_wp ) THEN
344                  pweig(1,1,jk) = z2dmmt(jk)
345                  pweig(1,2,jk) = z2dmpt(jk)
346                  pweig(2,1,jk) = z2dpmt(jk)
347                  pweig(2,2,jk) = z2dppt(jk)
348                  ! Set the vertical mask
349                  IF ( ( ( z2dmmt(jk) > 0.0_wp ) .AND.      &
350                     &   ( pmask(1,1,jk) == 1.0_wp ) ) .OR. &
351                     & ( ( z2dmpt(jk) > 0.0_wp ) .AND.      &
352                     &   ( pmask(1,2,jk) == 1.0_wp ) ) .OR. &
353                     & ( ( z2dpmt(jk) > 0.0_wp) .AND.      &
354                     &   ( pmask(2,1,jk) == 1.0_wp ) ) .OR. &
355                     & ( ( z2dppt(jk) > 0.0_wp ) .AND.      &
356                     &   ( pmask(2,2,jk) == 1.0_wp ) ) )    &
357                     & pobsmask(jk)=1.0_wp
358               ENDIF
359            ENDIF
360         ENDIF
361      END DO
362 
363   END SUBROUTINE obs_int_h2d_init
364
365   SUBROUTINE obs_int_h2d( kpk, kpk2, &
366      &                    pweig, pmod, pobsk )
367      !!-----------------------------------------------------------------------
368      !!
369      !!                     ***  ROUTINE obs_int_h2d  ***
370      !!
371      !! ** Purpose : Horizontal interpolation to the observation point.
372      !!
373      !! ** Method  : Horizontal interpolation to the observation point using
374      !!              model values at the corners of the surrounding grid
375      !!              points.
376      !!
377      !! ** Action  :
378      !!                   
379      !! References :
380      !!
381      !! History :
382      !!        ! 97-11 (A. Weaver, N. Daget)
383      !!        ! 06-03 (A. Vidard) NEMOVAR migration
384      !!        ! 06-10 (A. Weaver) Cleanup
385      !!        ! 07-08 (K. Mogensen) Split in two routines for easier adj.
386      !!-----------------------------------------------------------------------
387      !! * Modules used
388      !! * Arguments
389      INTEGER, INTENT(IN) :: &
390         & kpk,   &             ! Parameter values for automatic arrays
391         & kpk2
392      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
393         & pweig                ! Interpolation weights
394      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
395         & pmod                 ! Model variable to interpolate
396      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
397         & pobsk                ! Model profile interpolated to obs (i,j) pt
398
399      !! * Local declarations
400      INTEGER :: &
401         & jk
402      INTEGER :: &
403         & ikmax
404      REAL(KIND=wp) :: &
405         & zsum
406      !------------------------------------------------------------------------
407      ! Initialize number of levels
408      !------------------------------------------------------------------------
409      IF ( kpk2 == 1 ) THEN
410         ikmax = 1
411      ELSEIF ( kpk2 == kpk) THEN
412         ikmax = kpk-1
413      ENDIF
414      !------------------------------------------------------------------------
415      ! Interpolate to the observation point
416      !------------------------------------------------------------------------
417      pobsk(:) = obfillflt
418      DO jk = 1, ikmax
419         zsum = pweig(1,1,jk) + pweig(1,2,jk) + pweig(2,1,jk) +  pweig(2,2,jk)
420         IF ( zsum /= 0.0_wp ) THEN
421            pobsk(jk) = (  pweig(1,1,jk) * pmod(1,1,jk)         &
422               &         + pweig(1,2,jk) * pmod(1,2,jk)         &
423               &         + pweig(2,1,jk) * pmod(2,1,jk)         &
424               &         + pweig(2,2,jk) * pmod(2,2,jk)         &
425               &                                             ) / zsum
426         ENDIF
427      END DO
428
429   END SUBROUTINE obs_int_h2d
430 
431   SUBROUTINE obs_int_h2d_ds1( kpk2, kmax,                        &
432      &                        pphi, plam, pmask,                 &
433      &                        pphimm, plammm, pphimp, plammp,    &
434      &                        pphipm, plampm, pphipp, plampp,    &
435      &                        p2dmm, p2dmp, p2dpm, p2dpp )
436      !!-----------------------------------------------------------------------
437      !!
438      !!                     ***  ROUTINE obs_int_h2d_ds1  ***
439      !!
440      !! ** Purpose : Distance-weighted interpolation scheme (k2dint = 0)
441      !!
442      !! ** Method  : The interpolation weights are computed as a weighted
443      !!              sum of the distance between the model grid points (A)
444      !!              and the observation point (B).
445      !!
446      !!    Distance (s) is computed using the great-circle distance formula:
447      !!
448      !!    s(AB) = arcos(   sin( phiA ) x sin( phiB )
449      !!                   + cos( phiA ) x cos( phiB ) x cos( lamB - lamA )
450      !!
451      !! ** Action  :
452      !!
453      !! History :
454      !!        ! 97-11 (A. Weaver, N. Daget)
455      !!        ! 06-10 (A. Weaver) Cleanup
456      !!-----------------------------------------------------------------------
457
458      !! * Modules used
459
460      !! * Arguments
461      INTEGER, INTENT(IN) :: &
462         & kpk2, &             ! Parameter values for automatic arrays
463         & kmax
464      REAL(KIND=wp), INTENT(IN) :: &
465         & pphi,   &           ! Geographical location of observation
466         & plam,   &
467         & pphimm, &           ! Geographical location of surrounding
468         & pphimp, &           ! model grid points
469         & pphipm, &
470         & pphipp, &
471         & plammm, &
472         & plammp, &
473         & plampm, &
474         & plampp
475      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
476         & pmask               ! Model variable mask
477      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) :: &
478         & p2dmm, &            ! Interpolation weights
479         & p2dmp, &
480         & p2dpm, &
481         & p2dpp
482
483      !! * Local declarations
484      INTEGER :: &
485         & jk
486      REAL(KIND=wp) :: &
487         & zphi2,   &
488         & zlam2,   &
489         & zcola,   &
490         & za2,     &
491         & zb2,     &
492         & zc2,     &
493         & zphimm2, &
494         & zphimp2, &
495         & zphipm2, &
496         & zphipp2, &
497         & zlammm2, &
498         & zlammp2, &
499         & zlampm2, &
500         & zlampp2, &
501         & za1mm,   &
502         & za1mp,   &
503         & za1pm,   &
504         & za1pp,   &
505         & zcomm,   &
506         & zcomp,   &
507         & zcopm,   &
508         & zcopp,   &
509         & zb1mm,   &
510         & zb1mp,   &
511         & zb1pm,   &
512         & zb1pp,   &
513         & zc1mm,   &
514         & zc1mp,   &
515         & zc1pm,   &
516         & zc1pp,   &
517         & zsopmpp, &
518         & zsommmp, &
519         & zsomm,   &
520         & zsomp,   &
521         & zsopm,   &
522         & zsopp
523       
524      !------------------------------------------------------------------------
525      ! Distance-weighted interpolation using the great circle formula
526      !------------------------------------------------------------------------
527      zphi2 = pphi * rad
528      zlam2 = plam * rad
529      zcola = COS( zphi2 )
530      za2   = SIN( zphi2 )
531      zb2   = zcola * COS( zlam2 )
532      zc2   = zcola * SIN( zlam2 )
533       
534      zphimm2 = pphimm * rad
535      zphimp2 = pphimp * rad
536      zphipm2 = pphipm * rad
537      zphipp2 = pphipp * rad
538       
539      zlammm2 = plammm * rad
540      zlammp2 = plammp * rad
541      zlampm2 = plampm * rad
542      zlampp2 = plampp * rad
543       
544      za1mm = SIN( zphimm2 )
545      za1mp = SIN( zphimp2 )
546      za1pm = SIN( zphipm2 )
547      za1pp = SIN( zphipp2 )
548       
549      zcomm = COS( zphimm2 )
550      zcomp = COS( zphimp2 )
551      zcopm = COS( zphipm2 )
552      zcopp = COS( zphipp2 )
553       
554      zb1mm = zcomm * COS( zlammm2 )
555      zb1mp = zcomp * COS( zlammp2 )
556      zb1pm = zcopm * COS( zlampm2 )
557      zb1pp = zcopp * COS( zlampp2 )
558       
559      zc1mm = zcomm * SIN( zlammm2 )
560      zc1mp = zcomp * SIN( zlammp2 )
561      zc1pm = zcopm * SIN( zlampm2 )
562      zc1pp = zcopp * SIN( zlampp2 )
563       
564      ! Function for arcsin(sqrt(1-x^2) version of great-circle formula
565      zsomm = grt_cir_dis( za1mm, za2, zb1mm, zb2, zc1mm, zc2 )
566      zsomp = grt_cir_dis( za1mp, za2, zb1mp, zb2, zc1mp, zc2 )
567      zsopm = grt_cir_dis( za1pm, za2, zb1pm, zb2, zc1pm, zc2 )
568      zsopp = grt_cir_dis( za1pp, za2, zb1pp, zb2, zc1pp, zc2 )
569       
570      zsopmpp = zsopm * zsopp
571      zsommmp = zsomm * zsomp
572      DO jk = 1, kmax
573         p2dmm(jk) = zsomp * zsopmpp * pmask(1,1,jk)
574         p2dmp(jk) = zsomm * zsopmpp * pmask(1,2,jk)
575         p2dpm(jk) = zsopp * zsommmp * pmask(2,1,jk)
576         p2dpp(jk) = zsopm * zsommmp * pmask(2,2,jk)
577      END DO
578       
579   END SUBROUTINE obs_int_h2d_ds1
580 
581   SUBROUTINE obs_int_h2d_ds2( kpk2, kmax,                        &
582      &                        pphi, plam, pmask,                 &
583      &                        pphimm, plammm, pphimp, plammp,    &
584      &                        pphipm, plampm, pphipp, plampp,    &
585      &                        p2dmm, p2dmp, p2dpm, p2dpp )
586      !!-----------------------------------------------------------------------
587      !!
588      !!                     ***  ROUTINE obs_int_h2d_ds2  ***
589      !!
590      !! ** Purpose : Distance-weighted interpolation scheme (k2dint = 1)
591      !!
592      !! ** Method  : As k2dint = 0 but with distance (ds) computed using a
593      !!              small-angle approximation to the great-circle distance
594      !!              formula:
595      !!
596      !!    ds(AB) = sqrt(   ( phiB - phiA )^{2}
597      !!                   + ( ( lamB - lamA ) * cos( phiB ) )^{2} )
598      !!
599      !! ** Action  :
600      !!
601      !! History :
602      !!        ! 97-11 (A. Weaver, N. Daget)
603      !!        ! 06-10 (A. Weaver) Cleanup
604      !!-----------------------------------------------------------------------
605
606      !!-----------------------------------------------------------------------
607      !! * Modules used
608      !!-----------------------------------------------------------------------
609      !! * Arguments
610      INTEGER, INTENT(IN) :: &
611         & kpk2, &             ! Parameter values for automatic arrays
612         & kmax
613      REAL(KIND=wp), INTENT(IN) ::        &
614         & pphi,   &           ! Geographical location of observation
615         & plam,   &
616         & pphimm, &           ! Geographical location of surrounding
617         & pphimp, &           ! model grid points
618         & pphipm, &
619         & pphipp, &
620         & plammm, &
621         & plammp, &
622         & plampm, &
623         & plampp
624      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
625         & pmask               ! Model variable mask
626      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) :: &
627         & p2dmm, &            ! Interpolation weights
628         & p2dmp, &
629         & p2dpm, &
630         & p2dpp
631
632      !! * Local declarations
633      INTEGER :: &
634         & jk
635      REAL(KIND=wp) :: &
636         & zcosp,   &
637         & zdlmm,   &
638         & zdlmp,   &
639         & zdlpm,   &
640         & zdlpp,   &
641         & zdpmm,   &
642         & zdpmp,   &
643         & zdppm,   &
644         & zdppp,   &
645         & zsomm,   &
646         & zsomp,   &
647         & zsopm,   &
648         & zsopp,   &
649         & zsopmpp, &
650         & zsommmp
651       
652      !------------------------------------------------------------------------
653      ! Distance-weighted interpolation with a small angle approximation
654      !------------------------------------------------------------------------
655      zcosp = COS( pphi * rad )
656       
657      zdlmm = plammm - plam
658      zdlmp = plammp - plam
659      zdlpm = plampm - plam
660      zdlpp = plampp - plam
661       
662      zdpmm = pphimm - pphi
663      zdpmp = pphimp - pphi
664      zdppm = pphipm - pphi
665      zdppp = pphipp - pphi
666       
667      zsomm = grt_cir_dis_saa( zdlmm, zdpmm, zcosp )
668      zsomp = grt_cir_dis_saa( zdlmp, zdpmp, zcosp )
669      zsopm = grt_cir_dis_saa( zdlpm, zdppm, zcosp )
670      zsopp = grt_cir_dis_saa( zdlpp, zdppp, zcosp )
671       
672      zsopmpp = zsopm * zsopp
673      zsommmp = zsomm * zsomp
674       
675      DO jk = 1, kmax
676         p2dmm(jk) = zsomp * zsopmpp * pmask(1,1,jk)
677         p2dmp(jk) = zsomm * zsopmpp * pmask(1,2,jk)
678         p2dpm(jk) = zsopp * zsommmp * pmask(2,1,jk)
679         p2dpp(jk) = zsopm * zsommmp * pmask(2,2,jk)
680      END DO
681       
682   END SUBROUTINE obs_int_h2d_ds2
683 
684   SUBROUTINE obs_int_h2d_bil( kpk2, kmax,                        &
685      &                        pphi, plam, pmask,                 &
686      &                        plammp, pphipm, pphipp, plampp,    &
687      &                        p2dmm, p2dmp, p2dpm, p2dpp)
688      !!-----------------------------------------------------------------------
689      !!
690      !!                     ***  ROUTINE obs_int_h2d_bil  ***
691      !!
692      !! ** Purpose : Bilinear interpolation on a geographical grid (k2dint = 2)
693      !!
694      !! ** Method  : The interpolation is split into two 1D interpolations in
695      !!              the longitude and latitude directions, respectively.
696      !!
697      !!    An iterative scheme that involves first mapping a quadrilateral
698      !!    cell into a cell with coordinates (0,0), (1,0), (0,1) and (1,1).
699      !!
700      !! ** Action  :
701      !!
702      !! History :
703      !!        ! 97-11 (A. Weaver, N. Daget)
704      !!        ! 06-10 (A. Weaver) Cleanup
705      !!-----------------------------------------------------------------------
706
707      !! * Arguments
708      INTEGER, INTENT(IN) :: &
709         & kpk2, &             ! Parameter values for automatic arrays
710         & kmax
711      REAL(KIND=wp), INTENT(IN) :: &
712         & pphi,   &           ! Geographical location of observation
713         & plam,   &
714         & pphipm, &           ! Geographical location of surrounding
715         & pphipp, &           ! model grid points
716         & plammp, &
717         & plampp
718      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN)  :: &
719         & pmask               ! Model variable mask
720      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) :: &
721         & p2dmm, &            ! Interpolation weights
722         & p2dmp, &
723         & p2dpm, &
724         & p2dpp
725
726      !! * Local declarations
727      INTEGER :: &
728         & jk
729      REAL(KIND=wp) :: &
730         & zdlmp, &
731         & zdppm, &
732         & zdlpp, &
733         & zdppp
734       
735      !----------------------------------------------------------------------
736      ! Bilinear interpolation for geographical grid
737      !----------------------------------------------------------------------
738      zdlmp = ABS(plam   - plammp)
739      zdppm = ABS(pphi   - pphipm)
740      zdlpp = ABS(plampp - plam)
741      zdppp = ABS(pphipp - pphi)
742       
743      DO jk = 1, kmax
744         p2dmm(jk) = zdlpp * zdppp * pmask(1,1,jk)
745         p2dmp(jk) = zdlpp * zdppm * pmask(1,2,jk)
746         p2dpm(jk) = zdlmp * zdppp * pmask(2,1,jk)
747         p2dpp(jk) = zdlmp * zdppm * pmask(2,2,jk)
748      END DO
749       
750   END SUBROUTINE obs_int_h2d_bil
751 
752   SUBROUTINE obs_int_h2d_bir( kpk2, kmax,                        &
753      &                        pphi, plam, pmask,                 &
754      &                        pphimm, plammm, pphimp, plammp,    &
755      &                        pphipm, plampm, pphipp, plampp,    &
756      &                        p2dmm, p2dmp, p2dpm, p2dpp, ldfail )
757      !!-----------------------------------------------------------------------
758      !!
759      !!                     ***  ROUTINE obs_int_h2d_bir  ***
760      !!
761      !! ** Purpose : General bilinear remapping interpolation (k2dint = 3)
762      !!
763      !! ** Method  : An iterative scheme that involves first mapping a
764      !!              quadrilateral cell into a cell with coordinates
765      !!              (0,0), (1,0), (0,1) and (1,1).
766      !!
767      !! ** Action  :
768      !!
769      !! History :
770      !!        ! 97-11 (A. Weaver, N. Daget)
771      !!        ! 06-10 (A. Weaver) Cleanup
772      !!-----------------------------------------------------------------------
773
774      !! * Arguments
775      INTEGER, INTENT(IN) :: &
776         & kpk2, &             ! Parameter values for automatic arrays
777         & kmax
778      REAL(KIND=wp), INTENT(IN) ::        &
779         & pphi,   &           ! Geographical location of observation
780         & plam,   &
781         & pphimm, &           ! Geographical location of surrounding
782         & pphimp, &           ! model grid points
783         & pphipm, &
784         & pphipp, &
785         & plammm, &
786         & plammp, &
787         & plampm, &
788         & plampp
789      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
790         & pmask               ! Model variable mask
791      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) :: &
792         & p2dmm, &            ! Interpolation weights
793         & p2dmp, &
794         & p2dpm, &
795         & p2dpp
796      LOGICAL, INTENT(OUT) :: &
797         & ldfail
798      !! * Local declarations
799      INTEGER :: &
800         & jk
801      REAL(KIND=wp) :: &
802         & zbiwmm, &
803         & zbiwmp, &
804         & zbiwpm, &
805         & zbiwpp
806 
807      !----------------------------------------------------------------------
808      ! Bilinear remapping interpolation for general quadrilateral grid
809      !----------------------------------------------------------------------
810      CALL bil_wgt( pphimm, pphimp, pphipm, pphipp,  &
811         &          plammm, plammp, plampm, plampp,  &
812         &          zbiwmm, zbiwmp, zbiwpm, zbiwpp,  &
813         &          pphi  , plam,   ldfail           )
814
815      IF ( .NOT.ldfail ) THEN
816         DO jk = 1, kmax
817            p2dmm(jk) = zbiwmm * pmask(1,1,jk)
818            p2dmp(jk) = zbiwmp * pmask(1,2,jk)
819            p2dpm(jk) = zbiwpm * pmask(2,1,jk)
820            p2dpp(jk) = zbiwpp * pmask(2,2,jk)
821         END DO
822      ENDIF
823
824   END SUBROUTINE obs_int_h2d_bir
825
826   SUBROUTINE obs_int_h2d_pol( kpk2, kmax,                         &
827      &                        pphi, plam, pmask,                  &
828      &                        pphimm, plammm, pphimp, plammp,     &
829      &                        pphipm, plampm, pphipp, plampp,     &
830      &                        p2dmm, p2dmp, p2dpm, p2dpp )
831      !!-----------------------------------------------------------------------
832      !!
833      !!                     ***  ROUTINE obs_int_h2d_pol  ***
834      !!
835      !! ** Purpose : Polynomial interpolation (k2dint = 4)
836      !!
837      !! ** Method  : The interpolation weights are computed by fitting a
838      !!              polynomial function of the form
839      !!
840      !!    P(i) = a1(i) + a2(i) * phi + a3(i) * plam + a4(i) * phi * plam
841      !!
842      !!    through the model values at four surrounding grid pts (i=1,4).
843      !!    As k2dint = 0 but with distance (ds) computed using a small-
844      !!    angle approximation to the great-circle distance formula:
845      !!
846      !!    ds(AB) = sqrt(   ( phiB - phiA )^{2}
847      !!                   + ( ( lamB - lamA ) * cos( phiB ) )^{2} )
848      !!
849      !! ** Action  :
850      !!
851      !! History :
852      !!        ! 97-11 (A. Weaver, N. Daget)
853      !!        ! 06-10 (A. Weaver) Cleanup
854      !!-----------------------------------------------------------------------
855
856      !! * Arguments
857      INTEGER, INTENT(IN) :: & 
858         & kpk2,   &           ! Parameter values for automatic arrays
859         & kmax
860      REAL(KIND=wp), INTENT(IN) :: &
861         & pphi,   &           ! Geographical location of observation
862         & plam,   &
863         & pphimm, &           ! Geographical location of surrounding
864         & pphimp, &           ! model grid points
865         & pphipm, &
866         & pphipp, &
867         & plammm, &
868         & plammp, &
869         & plampm, &
870         & plampp
871      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(IN) :: &
872         & pmask               ! Model variable mask
873      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) :: &
874         & p2dmm,  &           ! Interpolation weights
875         & p2dmp,  &
876         & p2dpm,  &
877         & p2dpp
878 
879      !! * Local declarations
880      INTEGER :: &
881         & jk
882      REAL(KIND=wp) :: &
883         & zplp
884      REAL(KIND=wp), DIMENSION(4,4) :: &
885         & zmat,  &
886         & zmati
887       
888      !------------------------------------------------------------------------
889      ! Polynomial interpolation
890      !------------------------------------------------------------------------       
891      zmat(1,1) = 1.0_wp
892      zmat(1,2) = 1.0_wp
893      zmat(1,3) = 1.0_wp
894      zmat(1,4) = 1.0_wp
895      zmat(2,1) = plammm
896      zmat(2,2) = plammp
897      zmat(2,3) = plampm
898      zmat(2,4) = plampp
899      zmat(3,1) = pphimm
900      zmat(3,2) = pphimp
901      zmat(3,3) = pphipm
902      zmat(3,4) = pphipp
903      zmat(4,1) = plammm * pphimm
904      zmat(4,2) = plammp * pphimp
905      zmat(4,3) = plampm * pphipm
906      zmat(4,4) = plampp * pphipp
907       
908      CALL lu_invmat( zmat, 4, zmati )
909       
910      zplp = plam * pphi           
911      DO jk = 1, kmax
912         p2dmm(jk) = ABS(   zmati(1,1)        + zmati(1,2) * plam    &
913            &             + zmati(1,3) * pphi + zmati(1,4) * zplp )  &
914            &      * pmask(1,1,jk)
915         p2dmp(jk) = ABS(   zmati(2,1)        + zmati(2,2) * plam    &
916            &             + zmati(2,3) * pphi + zmati(2,4) * zplp )  &
917            &      * pmask(1,2,jk)
918         p2dpm(jk) = ABS(   zmati(3,1)        + zmati(3,2) * plam    &
919            &             + zmati(3,3) * pphi + zmati(3,4) * zplp )  &
920            &      * pmask(2,1,jk)
921         p2dpp(jk) = ABS(   zmati(4,1)        + zmati(4,2) * plam    &
922            &             + zmati(4,3) * pphi + zmati(4,4) * zplp )  &
923            &      * pmask(2,2,jk)
924      END DO
925 
926   END SUBROUTINE obs_int_h2d_pol
927
928   SUBROUTINE bil_wgt( pphimm, pphimp, pphipm, pphipp,   &
929      &               plammm, plammp, plampm, plampp,   &
930      &               pbiwmm, pbiwmp, pbiwpm, pbiwpp,   &
931      &               pphi  , plam,   ldfail            )
932      !!-------------------------------------------------------------------
933      !!
934      !!                   ***  ROUTINE bil_wgt  ***
935      !!
936      !! ** Purpose : Compute the weights for a bilinear remapping
937      !!              interpolation scheme.
938      !!
939      !! ** Method  : This scheme is appropriate for bilinear interpolation
940      !!              on a general quadrilateral grid.   
941      !!              This scheme is also used in OASIS.
942      !!
943      !!              This routine is a derivative of the SCRIP software.
944      !!              Copyright 1997, 1998 the Regents of the University
945      !!              of California. See SCRIP_Copyright.txt.
946      !!     
947      !! ** Action  :
948      !!
949      !! References : Jones, P.: A user's guide for SCRIP: A Spherical
950      !!                Coordinate Remapping and Interpolation Package.
951      !!                Version 1.4. Los Alamos.
952      !!                 
953      !!         http://www.acl.lanl.gov/climate/software/SCRIP/SCRIPmain.html
954      !!
955      !! History
956      !!      ! 97-11 (A. Weaver, N. Daget)
957      !!      ! 06-03 (A. Vidard)
958      !!      ! 06-10 (A. Weaver) Cleanup
959      !!-----------------------------------------------------------------------
960 
961      !! * Arguments
962      REAL(KIND=wp), INTENT(IN) :: &
963         & pphi,   &           ! Geographical location of observation
964         & plam,   &
965         & pphimm, &           ! Geographical location of surrounding
966         & pphimp, &           ! model grid points
967         & pphipm, &
968         & pphipp, &
969         & plammm, &
970         & plammp, &
971         & plampm, &
972         & plampp
973      REAL(KIND=wp), INTENT(OUT) :: &
974         & pbiwmm, &           ! Interpolation weights
975         & pbiwmp, &
976         & pbiwpm, &
977         & pbiwpp
978      LOGICAL, INTENT(out) :: &
979         & ldfail
980
981      !! * Local declarations
982      INTEGER :: &
983         & jiter
984      INTEGER :: &
985         & itermax
986      REAL(KIND=wp) :: &
987         & zphi,    &           ! Geographical location of observation
988         & zlam,    &
989         & zphimm,  &           ! Geographical location of surrounding
990         & zphimp,  &           ! model grid points
991         & zphipm,  &
992         & zphipp,  &
993         & zlammm,  &
994         & zlammp,  &
995         & zlampm,  &
996         & zlampp,  &
997         & zdth1,   &
998         & zdth2,   &
999         & zdth3,   &
1000         & zdthp,   &
1001         & zdph1,   &
1002         & zdph2,   &
1003         & zdph3,   &
1004         & zdphp,   &
1005         & zmat1,   &
1006         & zmat2,   &
1007         & zmat3,   &
1008         & zmat4,   &
1009         & zdeli,   &
1010         & zdelj,   &
1011         & ziguess, &
1012         & zjguess, &
1013         & zeps,    &
1014         & zdeterm, &
1015         & z2pi,    &
1016         & zhpi
1017     
1018      ! Initialization
1019
1020      ! Conversion to radians
1021
1022      zphi   = pphi   * rad
1023      zlam   = plam   * rad
1024      zphimm = pphimm * rad
1025      zphimp = pphimp * rad
1026      zphipm = pphipm * rad
1027      zphipp = pphipp * rad
1028      zlammm = plammm * rad
1029      zlammp = plammp * rad
1030      zlampm = plampm * rad
1031      zlampp = plampp * rad
1032
1033      ldfail = .FALSE.
1034
1035      zdth1 = zphipm - zphimm
1036      zdth2 = zphimp - zphimm
1037      zdth3 = zphipp - zphipm - zdth2
1038       
1039      zdph1 = zlampm - zlammm
1040      zdph2 = zlammp - zlammm
1041      zdph3 = zlampp - zlampm
1042       
1043      z2pi = 2.0_wp * rpi
1044       
1045      IF ( zdph1 >  3.0_wp * rpi ) zdph1 = zdph1 - z2pi
1046      IF ( zdph2 >  3.0_wp * rpi ) zdph2 = zdph2 - z2pi
1047      IF ( zdph3 >  3.0_wp * rpi ) zdph3 = zdph3 - z2pi
1048      IF ( zdph1 < -3.0_wp * rpi ) zdph1 = zdph1 + z2pi
1049      IF ( zdph2 < -3.0_wp * rpi ) zdph2 = zdph2 + z2pi
1050      IF ( zdph3 < -3.0_wp * rpi ) zdph3 = zdph3 + z2pi
1051       
1052      zdph3 = zdph3 - zdph2
1053       
1054      ziguess = 0.5_wp
1055      zjguess = 0.5_wp
1056       
1057      itermax = 100
1058
1059      IF ( wp == sp ) THEN
1060         zeps = 1.0e-6_wp  ! Single precision
1061      ELSE
1062         zeps = 1.0e-10_wp      ! Double precision
1063      ENDIF
1064 
1065      !------------------------------------------------------------------------
1066      ! Iterate to determine (i,j) in new coordinate system
1067      !------------------------------------------------------------------------
1068      jiter_loop: DO jiter = 1, itermax
1069         
1070         zdthp = zphi - zphimm - zdth1 * ziguess - zdth2 * zjguess &
1071            &  - zdth3 * ziguess * zjguess
1072         zdphp = zlam - zlammm
1073         
1074         zhpi = 0.5_wp * rpi
1075         IF ( zdphp >  3.0_wp * zhpi ) zdphp = zdphp - z2pi
1076         IF ( zdphp < -3.0_wp * zhpi ) zdphp = zdphp + z2pi
1077         
1078         zdphp = zdphp - zdph1 * ziguess - zdph2 * zjguess &
1079            &  - zdph3 * ziguess * zjguess
1080         
1081         zmat1 = zdth1 + zdth3 * zjguess
1082         zmat2 = zdth2 + zdth3 * ziguess
1083         zmat3 = zdph1 + zdph3 * zjguess
1084         zmat4 = zdph2 + zdph3 * ziguess
1085         
1086         ! Matrix determinant
1087         zdeterm = zmat1 * zmat4 - zmat2 * zmat3
1088         
1089         zdeli = ( zdthp * zmat4 - zmat2 * zdphp) / zdeterm
1090         zdelj = ( zmat1 * zdphp - zdthp * zmat3) / zdeterm
1091         
1092         IF ( ABS( zdeli ) < zeps .AND. ABS( zdelj ) < zeps ) EXIT jiter_loop
1093         
1094         ziguess = ziguess + zdeli
1095         zjguess = zjguess + zdelj
1096
1097         ! DJL prevent ziguess and zjguess from going outside the range
1098         ! 0 to 1
1099         ! prevents interpolated value going wrong
1100         ! for example sea ice concentration gt 1
1101         
1102         IF ( ziguess < 0 ) ziguess = 0.0_wp
1103         IF ( zjguess < 0 ) zjguess = 0.0_wp
1104         IF ( ziguess > 1 ) ziguess = 1.0_wp
1105         IF ( zjguess > 1 ) zjguess = 1.0_wp
1106         
1107      END DO jiter_loop
1108       
1109      IF ( jiter <= itermax ) THEN
1110           
1111         ! Successfully found i,j, now compute the weights
1112
1113         pbiwmm = ( 1.0_wp - ziguess ) * ( 1.0_wp - zjguess )
1114         pbiwmp = ( 1.0_wp - ziguess ) *            zjguess
1115         pbiwpm =            ziguess   * ( 1.0_wp - zjguess )
1116         pbiwpp =            ziguess   *            zjguess
1117         
1118      ELSEIF ( jiter > itermax ) THEN
1119           
1120         IF(lwp) THEN
1121           
1122            WRITE(numout,*)'Obs lat/lon  : ',pphi, plam
1123            WRITE(numout,*)'Grid lats    : ',pphimm, pphimp, pphipm, pphipp
1124            WRITE(numout,*)'Grid lons    : ',plammm, plammp, plampm, plampp
1125            WRITE(numout,*)'Current i,j  : ',ziguess, zjguess
1126            WRITE(numout,*)'Weights      : ',pbiwmm, pbiwmp, pbiwpm, pbiwpp
1127            WRITE(numout,*)'jiter        = ',jiter
1128            WRITE(numout,*)'zeps         = ',zeps
1129            WRITE(numout,*)'zdeli, zdelj = ',zdeli, zdelj
1130            WRITE(numout,*)' Iterations for i,j exceed max iteration count!'
1131            WRITE(numout,*)
1132           
1133            ldfail = .TRUE.
1134
1135         ENDIF
1136         
1137      ENDIF
1138   
1139   END SUBROUTINE bil_wgt
1140
1141   SUBROUTINE lu_invmat( pmatin, kdim, pmatou )
1142      !!-----------------------------------------------------------------------
1143      !!
1144      !!                   ***  ROUTINE lu_invmat  ***
1145      !!
1146      !! ** Purpose : Invert a matrix using LU decomposition.
1147      !!
1148      !! ** Method  :
1149      !!
1150      !! ** Action  :
1151      !!
1152      !! References :
1153      !!
1154      !! History
1155      !!      ! 02-11 (A. Weaver, N. Daget)
1156      !!      ! 06-03 (A. Vidard)
1157      !!      ! 06-10 (A. Weaver) Cleanup
1158      !!      ! 06-11 (NEMOVAR task force) Fix declaration of zd.
1159      !!-----------------------------------------------------------------------
1160
1161      !! * Arguments
1162      INTEGER, INTENT(IN) :: &
1163         & kdim             ! Array dimension
1164      REAL(KIND=wp), DIMENSION(kdim,kdim), INTENT(IN) :: &
1165         & pmatin
1166      REAL(KIND=wp), DIMENSION(kdim,kdim), INTENT(OUT) :: &
1167         & pmatou
1168 
1169      !! * Local declarations
1170      INTEGER :: &
1171         & ji, &
1172         & jj
1173      INTEGER, DIMENSION(kdim) :: &
1174         & indx
1175      REAL(KIND=wp), DIMENSION(kdim,kdim) :: &
1176         & zmat
1177      REAL(KIND=wp) :: &
1178         & zd
1179
1180      ! Invert the matrix       
1181      DO jj = 1, kdim
1182         DO ji = 1, kdim
1183            pmatou(ji,jj) = 0.0_wp
1184            zmat(ji,jj) = pmatin(ji,jj)
1185         END DO
1186         pmatou(jj,jj) = 1.0_wp
1187      END DO
1188      CALL lu_decomp( zmat, kdim, kdim, indx, zd )
1189      DO jj = 1, kdim
1190         CALL lu_backsb( zmat, kdim, kdim, indx, pmatou(1,jj) )
1191      END DO
1192     
1193   END SUBROUTINE lu_invmat
1194
1195   SUBROUTINE lu_decomp( pmatin, kdim1, kdim2, kindex, pflt )
1196      !!-----------------------------------------------------------------------
1197      !!
1198      !!                   ***  ROUTINE lu_decomp  ***
1199      !!
1200      !! ** Purpose : Compute the LU decomposition of a matrix
1201      !!
1202      !! ** Method  :
1203      !!
1204      !! ** Action  :
1205      !!
1206      !! References :
1207      !!
1208      !! History
1209      !!      ! 02-11 (A. Weaver, N. Daget)
1210      !!      ! 06-03 (A. Vidard)
1211      !!      ! 06-10 (A. Weaver) Cleanup
1212      !!-----------------------------------------------------------------------
1213       
1214      !! * Arguments
1215      INTEGER, INTENT(IN) :: &
1216         & kdim1, &   ! Array dimensions
1217         & kdim2
1218      INTEGER, DIMENSION(kdim1), INTENT(OUT) :: &
1219         & kindex
1220      REAL(KIND=wp), INTENT(OUT) :: &
1221         & pflt
1222      REAL(KIND=wp), DIMENSION(kdim2,kdim2), INTENT(INOUT) :: &
1223         & pmatin
1224 
1225      !! * Local declarations
1226      INTEGER, PARAMETER :: &
1227         & jpmax = 100
1228      REAL(KIND=wp), PARAMETER :: &
1229         & pptiny = 1.0e-20_wp
1230      REAL(KIND=wp), DIMENSION(jpmax) :: &
1231         & zvv
1232      INTEGER :: &
1233         & ji, &
1234         & jj, &
1235         & jk
1236      INTEGER :: &
1237         & imax
1238      REAL(KIND=wp) :: &
1239         & zsum,  &
1240         & zdum,  &
1241         & zaamax
1242     
1243      imax = -1
1244      ! Main computation
1245      pflt = 1.0_wp
1246      DO ji = 1, kdim1
1247         zaamax = 0.0_wp
1248         DO jj = 1, kdim1
1249            IF ( ABS( pmatin(ji,jj) ) > zaamax ) zaamax = ABS( pmatin(ji,jj) )
1250         END DO
1251         IF ( zaamax == 0.0_wp ) THEN
1252            CALL ctl_stop( 'singular matrix' )
1253         ENDIF
1254         zvv(ji) = 1.0_wp / zaamax
1255      END DO
1256      DO jj = 1, kdim1
1257         DO ji = 1, jj-1
1258            zsum = pmatin(ji,jj)
1259            DO jk = 1, ji-1
1260          zsum = zsum - pmatin(ji,jk) * pmatin(jk,jj)
1261       END DO
1262            pmatin(ji,jj) = zsum
1263         END DO
1264         zaamax = 0.0_wp
1265         DO ji = jj, kdim1
1266            zsum = pmatin(ji,jj)
1267            DO jk = 1, jj-1
1268               zsum = zsum - pmatin(ji,jk) * pmatin(jk,jj)
1269            END DO
1270            pmatin(ji,jj) = zsum
1271            zdum = zvv(ji) * ABS( zsum )
1272            IF ( zdum >= zaamax ) THEN
1273               imax = ji
1274               zaamax = zdum
1275            ENDIF
1276    END DO
1277         IF ( jj /= imax ) THEN
1278            DO jk = 1, kdim1
1279               zdum = pmatin(imax,jk)
1280               pmatin(imax,jk) = pmatin(jj,jk)
1281               pmatin(jj,jk) = zdum
1282       END DO
1283            pflt = -pflt
1284            zvv(imax) = zvv(jj)
1285         ENDIF
1286         kindex(jj) = imax
1287         IF ( pmatin(jj,jj) == 0.0_wp ) pmatin(jj,jj) = pptiny
1288         IF ( jj /= kdim1 ) THEN
1289            zdum = 1.0_wp / pmatin(jj,jj)
1290            DO ji = jj+1, kdim1
1291               pmatin(ji,jj) = pmatin(ji,jj) * zdum
1292       END DO
1293         ENDIF
1294      END DO
1295     
1296   END SUBROUTINE lu_decomp
1297     
1298   SUBROUTINE lu_backsb( pmat, kdim1, kdim2, kindex, pvect )
1299      !!-----------------------------------------------------------------------
1300      !!
1301      !!                   ***  ROUTINE lu_backsb  ***
1302      !!
1303      !! ** Purpose : Back substitution
1304      !!
1305      !! ** Method  :
1306      !!
1307      !! ** Action  :
1308      !!
1309      !! References :
1310      !!
1311      !! History
1312      !!      ! 02-11 (A. Weaver, N. Daget)
1313      !!      ! 06-03 (A. Vidard)
1314      !!      ! 06-10 (A. Weaver) Cleanup
1315      !!-----------------------------------------------------------------------
1316     
1317      !! * Arguments
1318      INTEGER, INTENT(IN) :: &
1319         & kdim1, &     ! Array dimensions
1320         & kdim2
1321      INTEGER, DIMENSION(kdim1), INTENT(IN) :: &
1322         & kindex
1323      REAL(KIND=wp), DIMENSION(kdim1), INTENT(INOUT) :: &
1324         & pvect
1325      REAL(KIND=wp), DIMENSION(kdim2,kdim2), INTENT(IN) :: &
1326         & pmat
1327 
1328      !! * Local declarations
1329      INTEGER :: &
1330         & ji,  &
1331         & jii, &
1332         & jj,  &
1333         & jll
1334      REAL(KIND=wp) :: &
1335         & zsum
1336
1337      ! Main computation
1338      jii = 0
1339      DO ji = 1, kdim1
1340         jll = kindex(ji)
1341         zsum = pvect(jll)
1342         pvect(jll) = pvect(ji)
1343         IF ( jii /= 0 ) THEN
1344            DO jj = jii, ji-1
1345               zsum = zsum - pmat(ji,jj) * pvect(jj)
1346            END DO
1347         ELSEIF ( zsum /= 0.0_wp ) THEN
1348            jii = ji
1349         ENDIF
1350         pvect(ji) = zsum
1351      END DO
1352      DO ji = kdim1, 1, -1
1353         zsum = pvect(ji)
1354         DO jj = ji+1, kdim1
1355            zsum = zsum - pmat(ji,jj) * pvect(jj)
1356         END DO
1357         pvect(ji) = zsum / pmat(ji,ji)
1358      END DO
1359     
1360   END SUBROUTINE lu_backsb
Note: See TracBrowser for help on using the repository browser.