Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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, &
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) :: &
97      REAL(KIND=wp), DIMENSION(2,2,kpk2), INTENT(OUT) ::  &
98         & pweig                ! Weights for interpolation
99      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
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      !------------------------------------------------------------------------
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
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.       &
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 ) ) )    &
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) :: &
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) :: &
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)  :: &
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) :: &
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) :: &
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 )  &
915         p2dmp(jk) = ABS(   zmati(2,1)        + zmati(2,2) * plam    &
916            &             + zmati(2,3) * pphi + zmati(2,4) * zplp )  &
918         p2dpm(jk) = ABS(   zmati(3,1)        + zmati(3,2) * plam    &
919            &             + zmati(3,3) * pphi + zmati(3,4) * zplp )  &
921         p2dpp(jk) = ABS(   zmati(4,1)        + zmati(4,2) * plam    &
922            &             + zmati(4,3) * pphi + zmati(4,4) * zplp )  &
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
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      ! Main computation
1244      pflt = 1.0_wp
1245      DO ji = 1, kdim1
1246         zaamax = 0.0_wp
1247         DO jj = 1, kdim1
1248            IF ( ABS( pmatin(ji,jj) ) > zaamax ) zaamax = ABS( pmatin(ji,jj) )
1249         END DO
1250         IF ( zaamax == 0.0_wp ) THEN
1251            CALL ctl_stop( 'singular matrix' )
1252         ENDIF
1253         zvv(ji) = 1.0_wp / zaamax
1254      END DO
1255      DO jj = 1, kdim1
1256         DO ji = 1, jj-1
1257            zsum = pmatin(ji,jj)
1258            DO jk = 1, ji-1
1259          zsum = zsum - pmatin(ji,jk) * pmatin(jk,jj)
1260       END DO
1261            pmatin(ji,jj) = zsum
1262         END DO
1263         zaamax = 0.0_wp
1264         DO ji = jj, kdim1
1265            zsum = pmatin(ji,jj)
1266            DO jk = 1, jj-1
1267               zsum = zsum - pmatin(ji,jk) * pmatin(jk,jj)
1268            END DO
1269            pmatin(ji,jj) = zsum
1270            zdum = zvv(ji) * ABS( zsum )
1271            IF ( zdum >= zaamax ) THEN
1272               imax = ji
1273               zaamax = zdum
1274            ENDIF
1275    END DO
1276         IF ( jj /= imax ) THEN
1277            DO jk = 1, kdim1
1278               zdum = pmatin(imax,jk)
1279               pmatin(imax,jk) = pmatin(jj,jk)
1280               pmatin(jj,jk) = zdum
1281       END DO
1282            pflt = -pflt
1283            zvv(imax) = zvv(jj)
1284         ENDIF
1285         kindex(jj) = imax
1286         IF ( pmatin(jj,jj) == 0.0_wp ) pmatin(jj,jj) = pptiny
1287         IF ( jj /= kdim1 ) THEN
1288            zdum = 1.0_wp / pmatin(jj,jj)
1289            DO ji = jj+1, kdim1
1290               pmatin(ji,jj) = pmatin(ji,jj) * zdum
1291       END DO
1292         ENDIF
1293      END DO
1294
1295   END SUBROUTINE lu_decomp
1296
1297   SUBROUTINE lu_backsb( pmat, kdim1, kdim2, kindex, pvect )
1298      !!-----------------------------------------------------------------------
1299      !!
1300      !!                   ***  ROUTINE lu_backsb  ***
1301      !!
1302      !! ** Purpose : Back substitution
1303      !!
1304      !! ** Method  :
1305      !!
1306      !! ** Action  :
1307      !!
1308      !! References :
1309      !!
1310      !! History
1311      !!      ! 02-11 (A. Weaver, N. Daget)
1312      !!      ! 06-03 (A. Vidard)
1313      !!      ! 06-10 (A. Weaver) Cleanup
1314      !!-----------------------------------------------------------------------
1315
1316      !! * Arguments
1317      INTEGER, INTENT(IN) :: &
1318         & kdim1, &     ! Array dimensions
1319         & kdim2
1320      INTEGER, DIMENSION(kdim1), INTENT(IN) :: &
1321         & kindex
1322      REAL(KIND=wp), DIMENSION(kdim1), INTENT(INOUT) :: &
1323         & pvect
1324      REAL(KIND=wp), DIMENSION(kdim2,kdim2), INTENT(IN) :: &
1325         & pmat
1326
1327      !! * Local declarations
1328      INTEGER :: &
1329         & ji,  &
1330         & jii, &
1331         & jj,  &
1332         & jll
1333      REAL(KIND=wp) :: &
1334         & zsum
1335
1336      ! Main computation
1337      jii = 0
1338      DO ji = 1, kdim1
1339         jll = kindex(ji)
1340         zsum = pvect(jll)
1341         pvect(jll) = pvect(ji)
1342         IF ( jii /= 0 ) THEN
1343            DO jj = jii, ji-1
1344               zsum = zsum - pmat(ji,jj) * pvect(jj)
1345            END DO
1346         ELSEIF ( zsum /= 0.0_wp ) THEN
1347            jii = ji
1348         ENDIF
1349         pvect(ji) = zsum
1350      END DO
1351      DO ji = kdim1, 1, -1
1352         zsum = pvect(ji)
1353         DO jj = ji+1, kdim1
1354            zsum = zsum - pmat(ji,jj) * pvect(jj)
1355         END DO
1356         pvect(ji) = zsum / pmat(ji,ji)
1357      END DO
1358
1359   END SUBROUTINE lu_backsb
Note: See TracBrowser for help on using the repository browser.