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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_inter_h2d.h90 in branches/devukmo2010/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/devukmo2010/NEMO/OPA_SRC/OBS/obs_inter_h2d.h90 @ 2128

Last change on this file since 2128 was 2128, checked in by rfurner, 14 years ago
merged branches OBS, ASM, Rivers, BDY & mixed_dynldf ready for vn3.3 merge
File size: 48.3 KB

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

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

Download in other formats:

Original Format