Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obsinter_h2d.h90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago
GMED 462 add flush
Property svn:mime-type set to `text/x-fortran`
File size: 48.4 KB

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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/OBS/obsinter_h2d.h90 @ 10986

Download in other formats: