Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

discendo.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago
Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).
File size: 123.0 KB

Line
1	C* ----------------------------------------------------------------------
2	C
3	C FSCINT
4	C ------
5	C *******************************
6	C * OASIS INTERPOLATION PACKAGE *
7	C * ----- ------------- ------- *
8	C *******************************
9	C
10	C* This is the Fast Scalar INTerpolator package initially written by Yves
11	C* Chartier and colleagues at RPN Canada. This software has been adapted
12	C* to the OASIS structure and provides several interpolation techniques
13	C* (see oasis documentation for details). This package has been split in
14	C* 2 files: discendo.f and semper.f (for people not fluent in latin,
15	C* "semper discendo" means "always learning" ... ;-)...). discendo.f
16	C* contains all the interpolation related routines while semper.f has
17	C* the memory allocation using fortran 90 features.
18	C
19	C History:
20	C -------
21	C Version Programmer Date Description
22	C ------- ---------- ---- -----------
23	C 1.1 L. Terray 94/06/01 Introduced
24	C 1.1 P. Braconnot 94/06/15 Extend longitude range for Z grids
25	C 2.0 L. Terray 96/09/25
26	C 2.2 G. Risari 97/12/01 new routine names
27	C 2.2 A.P, S.V, L.T 97/12/14 Change memory allocation
28	C 2.3 L. Terray 99/09/15 Introduction of Y grid and cleaning
29	C 2.4 E. Maisonnave 00/07/03 REAL NEWAXEY(I1:I2) (line 2806)
30	C
31	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32	C
33	FUNCTION FINDLON(VAL, TABLEAU, NBELEM, TMP)
34	C* ----------------------------------------------------------------------
35	C
36	C* New function to account for strange AGCM Z-grids (longitude<0 or >360)
37	C
38	C* ----------------------------------------------------------------------
39	USE mod_kinds_oasis
40	INTEGER (kind=ip_intwp_p) FINDLON
41	INTEGER (kind=ip_intwp_p) NBELEM
42	REAL(kind=ip_realwp_p) VAL
43	REAL(kind=ip_realwp_p) TABLEAU(NBELEM)
44	REAL(kind=ip_realwp_p) TMP
45	INTEGER (kind=ip_intwp_p) DEBUT, MILIEU, FIN
46	C
47	C* Get longitudes in the right interval
48	C
49	IF( VAL .GT. TABLEAU(NBELEM)) THEN
50	TMP = VAL - 360.
51	ELSEIF( VAL .LT. TABLEAU(1)) THEN
52	TMP = VAL + 360.
53	ELSE
54	TMP = VAL
55	ENDIF
56	C
57	C* Find grid point coordinates
58	C
59	DEBUT = 1
60	FIN = NBELEM
61	MILIEU = (DEBUT+FIN)*0.5
62	23000 IF( (MILIEU.NE. DEBUT))THEN
63	IF( TMP .LE. TABLEAU(MILIEU) ) THEN
64	FIN = MILIEU
65	ELSE
66	DEBUT = MILIEU
67	ENDIF
68	MILIEU = (DEBUT+FIN)*0.5
69	GOTO 23000
70	ENDIF
71	FINDLON = MILIEU
72	RETURN
73	END
74	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75	C
76	FUNCTION CHERCHE(VAL, TABLEAU, NBELEM)
77	USE mod_kinds_oasis
78	INTEGER (kind=ip_intwp_p) CHERCHE
79	INTEGER (kind=ip_intwp_p) NBELEM
80	REAL(kind=ip_realwp_p) VAL
81	REAL(kind=ip_realwp_p) TABLEAU(NBELEM)
82	INTEGER (kind=ip_intwp_p) DEBUT, MILIEU, FIN
83	DEBUT = 1
84	FIN = NBELEM
85	MILIEU = (DEBUT+FIN)*0.5
86	23000 IF( (MILIEU.NE. DEBUT))THEN
87	IF( (VAL.LE. TABLEAU(MILIEU)))THEN
88	FIN = MILIEU
89	ELSE
90	DEBUT = MILIEU
91	ENDIF
92	MILIEU = (DEBUT+FIN)*0.5
93	GOTO 23000
94	ENDIF
95	CHERCHE = MILIEU
96	RETURN
97	END
98	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99	C
100	SUBROUTINE CHKXTRAP(PX, PY, NPTS, NI, NJ)
101	USE mod_kinds_oasis
102	IMPLICIT NONE
103	INTEGER (kind=ip_intwp_p) NPTS, NI, NJ, OFFL, OFFR
104	REAL(kind=ip_realwp_p) PX(NPTS), PY(NPTS)
105	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
106	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
107	PARAMETER (VOISIN = 0)
108	PARAMETER (LINEAIR = 1)
109	PARAMETER (CUBIQUE = 3)
110	PARAMETER (OUI = 1)
111	PARAMETER (MINIMUM = 2)
112	PARAMETER (MAXIMUM = 3)
113	PARAMETER (VALEUR = 4)
114	PARAMETER (ABORT = 13)
115	LOGICAL FLGXTRAP
116	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
117	REAL(kind=ip_realwp_p) VALXTRAP
118	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
119	INTEGER (kind=ip_intwp_p) N, I, J
120	FLGXTRAP = .FALSE.
121	IF( (ORDINT.EQ. 3))THEN
122	OFFR = 3
123	OFFL = 2
124	ELSE
125	OFFR = 1
126	OFFL = 0
127	ENDIF
128	DO 23002 N=1, NPTS
129	I = INT(PX(N))
130	J = INT(PY(N))
131	IF( (I.LT. OFFL.OR. J.LT. OFFL.OR. I.GT. (NI-OFFR).OR. J
132	$ .GT. (NJ-OFFR)))THEN
133	FLGXTRAP = .TRUE.
134	ENDIF
135	23002 CONTINUE
136	IF( (FLGXTRAP.AND. CODXTRAP.EQ. ABORT))THEN
137	WRITE(6, *)
138	$'*****************************************************************
139	$****'
140	WRITE(6, *)
141	$ 'target grid contains points outside of source grid'
142	WRITE(6, *) 'We stop the program'
143	WRITE(6, *)
144	$'*****************************************************************
145	$****'
146	STOP
147	ENDIF
148	RETURN
149	END
150	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151	C
152	SUBROUTINE DGAUSS (N,ROOTS,KASE)
153	USE mod_kinds_oasis
154	REAL(kind=ip_realwp_p) ROOTS(*)
155	C
156	C AUTEUR- D. ROBERTSON
157	C
158	C OBJET(DGAUSS)
159	C - CALCULATES THE ZEROES OF THE ORDINARY LEGENDRE
160	C POLYNOMIAL OF ORDER N, I.E. DEFINE GAUSSIAN GRID
161	C
162	C ALGORITHME
163	C - THE POSITIVE ROOTS ARE APPROXIMATED BY THE BEST
164	C ASYMPTOTIC FORMULA AVAILABLE TO THE AUTHOR, FOUND IN
165	C ABRAMOWITZ AND STEGUN "HANDBOOK OF MATHEMATICAL FUNCTIONS".
166	C CHAPTER 22 FORMULA 22.16.6.
167	C NEWTON'S METHOD IS USED TO REFINE THE GUESS TO PRECISION
168	C DEFINED BY THE CONSTANT TOL. SINCE THE ROOTS ARE OF ORDER
169	C OF MAGNITUDE UNITY, ABSOLUTE PRECISION IS ADEQUATE, RATHER
170	C THAN A RELATIVE TEST.
171	C A STANDARD IDENTITY IS USED TO DETERMINE THE DERIVATIVE OF
172	C THE POLYNOMIAL IN TERMS OF THE VALUES OF P(N;X),P(N-1;X).
173	C (X*2-1.0)(DP/DX)=N(XP(N;,X)-P(N-1;X)).
174	C SEE ABRAMOWITZ AND STEGUN FORMULA 22.8.5
175	C NOTE THAT IN CONTRAST TO OTHER FORMULAS THIS REQUIRES ONLY
176	C 2 EVALUATIONS OF A LEGENDRE POLYNOMIAL PER ITERATION.
177	C NOTE THAT THE COORDINATE USED IS CONVENTIONALLY REFERRED TO
178	C AS MU=COS(THETA), RUNNING FROM +1 TO -1, FOR THETA FROM 0 TO
179	C PI. THE NEGATIVE ROOTS ARE FILLED BY SYMMETRY.
180	C FOR KASE=GLOBAL, ALL N ROOTS ARE FOUND, WHILE FOR
181	C DASE=NORTH/SOUTH ONLY THE +VE/-VE ROOTS ARE FOUND,
182	C (INCLUDING 0 IF N IS ODD) I.E. N/2+MOD(N,2) ROOTS.
183	C
184	C
185	C APPEL - CALL DGAUSS(N,ROOTS,KASE)
186	C
187	C ARGUMENTS
188	C IN - N - ORDER OF THE POLYNOMIALS
189	C OUT - ROOTS - ARRAY CONTAINING THE ZEROES OF THE
190	C ORDINARY LEGENDRE POLYNOMIALS
191	C IN - KASE - =0, GLOBAL
192	C =1, NORTH
193	C =2, SOUTH
194	C
195	C MODULES APPELES
196	C - ORDLEG
197	C
198	C ----------------------------------------------------------------------
199	C
200	C THE ANSWERS ARE RETURNED IN ROOTS.
201
202	REAL(kind=ip_realwp_p) NORMN,NORMNM
203	INTEGER (kind=ip_intwp_p) GLOBAL,NORTH,SOUTH,NORD,SUD
204	PARAMETER(GLOBAL=0,NORTH=1,NORD=1,SOUTH=2,SUD=2)
205	DATA TOL /1.0E-06/
206	C
207	C RDTODG = 180/PI, DGTORD = PI/180
208	C
209	DATA PI /3.14159265358979323846/
210	DATA RDTODG /57.295779513082/
211	DATA DGTORD /1.7453292519943E-2/
212	C
213	C ORDLEG RETURNS POLYNOMIALS NORMALIZED TO UNIT INTEGRAL.
214	C NORMN,NORMNMN RESTORE THE CONVENTION NORMALIZATION, P(N;1.0)=1.0.
215	C* ----------------------------------------------------------------------
216
217	NORMN =SQRT(2.0/(2.0*N+1.0))
218	NORMNM=SQRT(2.0/(2.0*N-1.0))
219	L = N/2
220	C
221	C CALCULATE ASYMPTOTIC APPROXIMATION
222	C
223	DO 23000 I=1,L
224	IF((KASE.NE.SOUTH))THEN
225	J = I
226	ENDIF
227	IF((KASE.EQ.SOUTH))THEN
228	J=I+L+MOD(N,2)
229	ENDIF
230	T = (4J-1)PI/FLOAT(4*N+2)
231	IF((KASE.NE.SOUTH))THEN
232	IRT = I
233	ENDIF
234	IF((KASE.EQ.SOUTH))THEN
235	IRT = I + MOD(N,2)
236	ENDIF
237	ROOTS(IRT)=COS(T+1.0/(8.0FLOAT(N2)TAN(T)))
238
239	C
240	23000 CONTINUE
241	DO 23010 I=1,L
242
243	C
244	C REPEAT 1 NEWTON ITERATION
245	C ** BEGIN
246	C
247	6 CALL ORDLEG(G,ROOTS(I),N)
248	CALL ORDLEG(GM,ROOTS(I),N-1)
249	PN = NORMN*G
250	PNM= NORMNM*GM
251	C
252	C ** GUESS(K+1)=GUESS(K)-P/(DP/DX)
253	C
254
255	RDPDX = (ROOTS(I)*2-1.0)/(N(ROOTS(I)*PN-PNM))
256	DELTA = -PN*RDPDX
257	ROOTS(I) = ROOTS(I)+DELTA
258	IF((ABS(DELTA).GT.TOL))THEN
259	GO TO 6
260	C
261	C ** END UNTIL ABS(DELTA).LE.TOL
262	C
263
264	ENDIF
265	ROOTS(N+1-I) = -ROOTS(I)
266
267	C
268	23010 CONTINUE
269	IF((MOD(N,2).EQ.0))THEN
270	RETURN
271	ENDIF
272	IF((KASE.NE.SOUTH))THEN
273	IRT = L+1
274	ENDIF
275	IF((KASE.EQ.SOUTH))THEN
276	IRT = 1
277	ENDIF
278	ROOTS(IRT) = 0.0
279	RETURN
280	END
281	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282	C
283	SUBROUTINE GAUSS(NRACP,RACP,PG,SIA,RAD,PGSSIN2,SINM1,SINM2,
284	$SIN2)
285	C* -----------------------------------------------------------------------
286	C CALCULE LES NRACP RACINES POSITIVES DU POLYNOME DE LEGENDRE DE
287	C DEGRE 2*NRACP (ICI-APRES NOTE PN) DEFINI SUR L INTERVALLE DES
288	C COLATITUDES ALLANT DE 0 (POLE NORD) A PI (POLE SUD). ON SAIT QUE
289	C LES 2*NRACP RACINES SONT ANTI-SYMETRIQUES P/R A L EQUATEUR PI/2,
290	C ETANT POSITIVES ENTRE COLAT=0 ET COLAT =PI/2.
291	C ON CALCULE ENSUITE LES POIDS DE GAUSS ASSOCIES AUX COLATITUDES
292	C GAUSSIENNES (ICI APRES NOTEES CG), AINSI QU UN CERTAIN NOMBRE DE
293	C FONCTIONS DE CG DEFINIES PLUS LOIN. ON RAPPELLE ENFIN QUE LA LATI-
294	C TUDE LAT=COLAT-PI/2, ET DONC QUE SIN(LAT)=COS(COLAT).
295	C NRACP : NOMBRE DE RACINES POSITIVES DU POLYNOME DE LEGENDRE
296	C : DE DEGRE 2*NRACP.
297	C RACP(I) : RACINES DE PN, =SIN(LG)=COS(CG).
298	C PG(I) : POIDS DE GAUSS CORRESPONDANTS.
299	C SIA(I) : SIN(CG)=COS(LG).
300	C RAD(I) : COLATITUDE CG EN RADIANS.
301	C PGSSIN2(I) : POIDS DE GAUSS / (SIN(CG))**2.
302	C SINM1(I) : (SIN(CG))**-1.
303	C SINM2(I) : (SIN(CG))**-2.
304	C VOIR NST 8, CHAP. A, PP.1-7, ET APPENDICE D12, PP. 26-27.
305	C VERSION REVISEE PAR MICHEL BELAND, 9 DECEMBRE 1980.
306	C *****************************************************************
307	C
308	C
309	USE mod_kinds_oasis
310	C -----------------------------------------------------------------
311
312	DIMENSION RACP(1),PG(1),SIA(1),RAD(1),PGSSIN2(1),SINM1(1),
313	$SINM2(1),SIN2(1)
314	C --------------------------------------------------------------
315	C
316	C ON DEMANDE UNE PRECISION DE 1.E-13 POUR LES RACINES DE PN.
317	C
318
319	XLIM=1.E-6
320	PI = 3.14159265358979323846
321	IR = 2*NRACP
322	FI=FLOAT(IR)
323	FI1=FI+1.
324	FN=FLOAT(NRACP)
325	C
326	C ON UTILISE UNE FORMULE ASYMPTOTIQUE POUR OBTENIR LES VALEURS
327	C APPROXIMATIVES DES COLATITUDES GAUSSIENNES
328	C CG(I) = (PI/2) * (2I-1)/(2NRACP).
329	C VOIR ABRAMOWITZ AND STEGUN, P. 787, EQU. 22.16.6 .
330	C
331
332	DO 23000 I=1,NRACP
333	DOT=FLOAT(I-1)
334	RACP(I)=-PI.5(DOT+.5)/FN + PI*.5
335	RACP(I) = SIN(RACP(I))
336
337	C
338	C ON CALCULE ENSUITE LES CONSTANTES FACTEURS DE P(N+1) ET P(N-1)
339	C DANS L EXPRESSION DE LA PSEUDO-DERIVEE DE PN.
340	C
341	23000 CONTINUE
342	DN = FI/SQRT(4.FIFI-1.)
343	DN1=FI1/SQRT(4.FI1FI1-1.)
344	A = DN1*FI
345	B = DN*FI1
346	IRP = IR + 1
347	IRM = IR -1
348	C
349	C ON EMPLOIE ENSUITE UNE METHODE DE NEWTON POUR AUGMENTER LA PREC.
350	C SI RACTEMP EST UNE SOL. APPROXIMATIVE DE PN(RACP)=0., ALORS LA
351	C SEQUENCE RACTEMP(I+1)=RACTEMP(I)-PN(RACTEMP(I))/DER.PN(RACTEMP(I))
352	C CONVERGE VERS RACP DE FACON QUADRATIQUE.
353	C VOIR ABRAMOWITZ AND STEGUN, P.18, EQU. 3.9.5.
354	C ORDLEG CALCULE LA VALEUR DE PN (RACP) , NORMALISE.
355	C
356
357	DO 23002 I=1,NRACP
358	42 CALL ORDLEG(G,RACP(I),IR)
359	CALL ORDLEG(GM,RACP(I),IRM)
360	CALL ORDLEG(GP,RACP(I),IRP)
361	GT = (AGP-BGM)/(RACP(I)*RACP(I)-1.)
362	RACTEMP = RACP(I) - G/GT
363	GTEMP = RACP(I) - RACTEMP
364	RACP(I) = RACTEMP
365	IF(( ABS(GTEMP).GT.XLIM))THEN
366	GO TO 42
367	ENDIF
368
369	C
370	C ON CALCULE ENSUITE LES POIDS DE GAUSS SELON L ALGORITHME
371	C PG(I) = 2./[(1.-RACP(I)*2)(DER.PN(RACP(I)))**2].
372	C VOIR ABRAMOWITZ AND STEGUN, P.887, EQU. 25.4.29.
373	C NOTE: ON DOIT MULTIPLIER LA PRECEDENTE FORMULE PAR UN FACTEUR
374	C DE DENORMALISATION, LES PN DONNES PAR ORDLEG ETANT NORMALISES.
375	C ON SE SERT D UNE FORMULE DE RECURRENCE POUR LA DERIVEE DE PN.
376	C
377	23002 CONTINUE
378	DO 23006 I=1,NRACP
379	A=2.(1.-RACP(I)*2)
380	CALL ORDLEG(B,RACP(I),IRM)
381	B = BBFI*FI
382	PG(I)=A*(FI-.5)/B
383	RAD(I) = ACOS(RACP(I))
384	SIA(I) = SIN(RAD(I))
385	C=(SIA(I))**2
386	SINM1(I) = 1./SIA(I)
387	SINM2(I) = 1./C
388	PGSSIN2(I) =PG(I)/C
389	SIN2(I)=C
390	23006 CONTINUE
391	RETURN
392	END
393	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
394	C
395	SUBROUTINE GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
396	USE mod_kinds_oasis
397	IMPLICIT NONE
398	INTEGER (kind=ip_intwp_p) LI,LJ
399	REAL(kind=ip_realwp_p) SPDO(LI,LJ), PSIO(LI,LJ), XLON(LI,LJ)
400	CHARACTER*1 GRTYP
401	INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
402	EXTERNAL CIGAXG
403	C
404	C AUTEUR - Y. CHARTIER - AVRIL 91
405	C
406	C OBJET(GDW2LLW)
407	C - PASSE DE VENT DE GRILLE (COMPOSANTES U ET V)
408	C - A VITESSE ET DIRECTION (SPEED, PSI)
409	C APPEL - CALL GDW2LLW(SPD,PSI,LI,LJ,IYP,XG1,XG2,XG3,XG4)
410	C
411	C MODULES - XGAIG
412	C
413	C ARGUMENTS
414	C IN/OUT - SPD - A L'ENTREE CONTIENT LA PSIOTESSE DU VENT ET
415	C A LA SORTIE LA COMPOSANTE U.
416	C IN/OUT - PSI - A L'ENTREE CONTIENT LA DIRECTION DU VENT ET
417	C A LA SORTIE LA COMPOSANTE V.
418	C IN - LI - PREMIERE DIMENSION DES CHAMPS SPD ET PSI
419	C IN - LJ - DEUXIEME DIMENSION DES CHAMPS SPD ET PSI
420	C IN - IGTYP - TYPE DE GRILLE (VOIR OUVRIR)
421	C IN - XG1 - ** DESCRIPTEUR DE GRILLE (REEL),
422	C IN - XG2 - IGTYP = 'N', PI, PJ, D60, DGRW
423	C IN - XG3 - IGTYP = 'L', LAT0, LON0, DLAT, DLON,
424	C IN - XG4 - IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
425	C = 1. NORD
426	C = 2. SUD **
427	C
428	C MESSAGES - "ERREUR MAUVAISE GRILLE (GDW2LLW)"
429	C
430	C-------------------------------------------------------------
431	C
432	C
433	C
434	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
435	C
436	C RDTODG = 180/PIE, DGTORD = PIE/180
437
438	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
439	DATA PIE /3.14159265358979323846/
440	DATA RDTODG /57.295779513082/
441	DATA DGTORD /1.7453292519943E-2/
442	C
443
444	INTEGER (kind=ip_intwp_p) I,J
445	REAL(kind=ip_realwp_p) SPD, PSI
446	REAL(kind=ip_realwp_p) XG1,XG2,XG3,XG4
447	IF( (GRTYP.EQ. 'N'))THEN
448	CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
449	DO 23002 I=1,LI
450	DO 23004 J=1,LJ
451	SPD=SQRT(SPDO(I,J)SPDO(I,J)+PSIO(I,J)PSIO(I,J))
452	IF( (SPD.EQ. 0.0))THEN
453	PSI= 0.0
454	ELSE
455	IF( (SPDO(I,J).EQ. 0.0))THEN
456	IF( (PSIO(I,J).GE. 0.0))THEN
457	PSI= XLON(I,J)+XG4-90.0
458	ELSE
459	PSI= XLON(I,J)+XG4+90.0
460	ENDIF
461	ELSE
462	PSI=XLON(I,J)+XG4-RDTODG*ATAN2(PSIO(I,J),SPDO(I
463	$ ,J))
464	ENDIF
465	ENDIF
466	PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
467	$ +360.0_ip_realwp_p,360.0_ip_realwp_p)
468	SPDO(I,J)=SPD
469	PSIO(I,J)=PSI
470	23004 CONTINUE
471	23002 CONTINUE
472	RETURN
473	ENDIF
474	IF( (GRTYP.EQ. 'S'))THEN
475	CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
476	DO 23014 I=1,LI
477	DO 23016 J=1,LJ
478	SPD=SQRT(SPDO(I,J)SPDO(I,J)+PSIO(I,J)PSIO(I,J))
479	IF( (SPD.EQ. 0.0))THEN
480	PSI= 0.0
481	ELSE
482	IF( (SPDO(I,J).EQ. 0.0))THEN
483	IF( (PSIO(I,J).GE. 0.0))THEN
484	PSI= 90.0 - XLON(I,J)+XG4
485	ELSE
486	PSI= 270.0 - XLON(I,J)+XG4
487	ENDIF
488	ELSE
489	PSI= 180.0 - XLON(I,J)+XG4-RDTODG*ATAN2(PSIO(I,
490	$ J),SPDO(I,J))
491	ENDIF
492	ENDIF
493	PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
494	$ +360.0_ip_realwp_p,360.0_ip_realwp_p)
495	SPDO(I,J)=SPD
496	PSIO(I,J)=PSI
497	23016 CONTINUE
498	23014 CONTINUE
499	RETURN
500	ENDIF
501	IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'G'.OR.GRTYP.EQ.
502	$'L'))THEN
503	DO 23026 I=1,LI
504	DO 23028 J=1,LJ
505	SPD=SQRT(SPDO(I,J)SPDO(I,J)+PSIO(I,J)PSIO(I,J))
506	IF( (SPD.EQ. 0.0))THEN
507	PSI= 0.0
508	ELSE
509	IF( (SPDO(I,J).EQ. 0.0))THEN
510	IF( (PSIO(I,J).GE. 0.0))THEN
511	PSI= 180.0
512	ELSE
513	PSI= 0.0
514	ENDIF
515	ELSE
516	PSI=270.0 - RDTODG*ATAN2(PSIO(I,J),SPDO(I,J))
517	ENDIF
518	ENDIF
519	PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
520	$ +360.0_ip_realwp_p,360.0_ip_realwp_p)
521	SPDO(I,J)=SPD
522	PSIO(I,J)=PSI
523	23028 CONTINUE
524	23026 CONTINUE
525	RETURN
526	ENDIF
527	WRITE(6, 600) GRTYP
528	600 FORMAT('0',' Error bad grid (GDW2LLW) - GRTYP = ', A1)
529	RETURN
530	END
531	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
532	C
533	SUBROUTINE GGGDINT(ZO,PX,PY,NPTS,AY,CY,Z,I1,I2,J1,J2)
534	USE mod_kinds_oasis
535	IMPLICIT NONE
536	C*******
537	C AUTEUR: Y.CHARTIER, DRPN
538	C FEVRIER 1991
539	C
540	C OBJET: INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR
541	C D'UNE GRILLE GAUSSIENNE.
542	C*******
543	INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
544	REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
545	REAL(kind=ip_realwp_p) AY(J1:J2),CY(J1:J2,6)
546	REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
547	C
548	C NPTS : NOMBRE DE POINTS A INTERPOLER
549	C I1:I2 : DIMENSION DE LA GRILLE SOURCE SELON X
550	C J1:J2 : DIMENSION DE LA GRILLE SOURCE SELON Y
551	C ZO : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
552	C PX : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
553	C VEUT INTERPOLER
554	C PY : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
555	C VEUT INTERPOLER
556	C AY : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES Y.
557	C CY : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON Y.
558	C Z : VALEURS DE LA GRILLE SOURCE.
559	C
560	C***************************************************************************
561	C
562	C * * * *
563	C
564	C * * * *
565	C # ==> PT (X,Y)
566	C * (=) * * ==> = PT (I, J)
567	C
568	C * * * *
569	C
570	C
571	C
572	C CY(I,1) = 1.0 / (X2-X1)
573	C CY(I,2) = 1.0 / (X3-X1)
574	C CY(I,3) = 1.0 / (X3-X2)
575	C CY(I,4) = 1.0 / (X4-X1)
576	C CY(I,5) = 1.0 / (X4-X2)
577	C CY(I,6) = 1.0 / (X4-X3)
578	C
579	C STRUCTURE IDENTIQUE POUR CY(J,1..6)
580
581	INTEGER (kind=ip_intwp_p) I, J, M, N,STRIDE
582	REAL(kind=ip_realwp_p) X, X1, X2, X3, X4
583	REAL(kind=ip_realwp_p) B1, B2, B3, B4
584	REAL(kind=ip_realwp_p) B11, B12, B13, B14
585	REAL(kind=ip_realwp_p) Y,Y1,Y2,Y3,Y4
586	REAL(kind=ip_realwp_p) Y11, Y12, Y13, Y14
587	REAL(kind=ip_realwp_p) AY1, AY2, AY3, AY4
588	REAL(kind=ip_realwp_p) FA, FA2, FA3, FA4
589	REAL(kind=ip_realwp_p) A1,A2,A3,A4,C1,C2,C3,C4,C5,C6
590	C DEFINITION DES FONCTIONS IN-LINE
591
592	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
593	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
594	PARAMETER (VOISIN = 0)
595	PARAMETER (LINEAIR = 1)
596	PARAMETER (CUBIQUE = 3)
597	PARAMETER (OUI = 1)
598	PARAMETER (MINIMUM = 2)
599	PARAMETER (MAXIMUM = 3)
600	PARAMETER (VALEUR = 4)
601	PARAMETER (ABORT = 13)
602	LOGICAL FLGXTRAP
603	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
604	REAL(kind=ip_realwp_p) VALXTRAP
605	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
606	REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
607	REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
608	CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)0.1666666666666 + 0.5(Z2-Z3)
609	$)DX + 0.5(Z1+Z3)-Z2)DX + Z3-0.1666666666666Z4-0.5*Z2-0.
610	$3333333333333Z1)DX+Z2
611	ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
612	FA(A1,A2,A3,A4,X,X1,X2,X3)=A1+(X-X1)(A2+(X-X2)(A3+A4*(X-X3))
613	$)
614	FA2(C1,A1,A2)=C1*(A2-A1)
615	FA3(C1,C2,C3,A1,A2,A3)=C2(C3(A3-A2)-C1*(A2-A1))
616	FA4(C1,C2,C3,C4,C5,C6,A1,A2,A3,A4)=C4(C5(C6(A4-A3)-C3(A3-
617	$A2)) - C2(C3(A3-A2)-C1*(A2-A1)))
618	STRIDE = 1
619	IF( (ORDINT.EQ. CUBIQUE))THEN
620	DO 23002 N=1,NPTS
621	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
622	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
623	DX = PX(N) - I
624	Y1=CUBIC(Z(I-1,J-1),Z(I ,J-1),Z(I+1,J-1),Z(I+2,J-1),DX)
625	Y2=CUBIC(Z(I-1,J ),Z(I ,J ),Z(I+1,J ),Z(I+2,J ),DX)
626	Y3=CUBIC(Z(I-1,J+1),Z(I ,J+1),Z(I+1,J+1),Z(I+2,J+1),DX)
627	Y4=CUBIC(Z(I-1,J+2),Z(I ,J+2),Z(I+1,J+2),Z(I+2,J+2),DX)
628	Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
629	C INTERPOLATION FINALE SELON Y
630
631	AY1=AY(J-1)
632	AY2=AY(J)
633	AY3=AY(J+1)
634	AY4=AY(J+2)
635	Y11 = Y1
636	Y12 = FA2(CY(J,1),Y1,Y2)
637	Y13 = FA3(CY(J,1),CY(J,2),CY(J,3),Y1,Y2,Y3)
638	Y14 = FA4(CY(J,1),CY(J,2),CY(J,3),CY(J,4),CY(J,5),CY(J,6
639	$ ),Y1,Y2,Y3,Y4)
640	ZO(N) = FA(Y11,Y12,Y13,Y14,Y,AY1,AY2,AY3)
641	23002 CONTINUE
642	ENDIF
643	IF( (ORDINT.EQ. LINEAIR))THEN
644	DO 23006 N=1,NPTS
645	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
646	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
647	DX = PX(N) - I
648	Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
649	DY = (Y -AY(J))/(AY(J+1)-AY(J))
650	Y2=ZLIN(Z(I,J ),Z(I+1,J ),DX)
651	Y3=ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
652	ZO(N)=ZLIN(Y2,Y3,DY)
653	23006 CONTINUE
654	ENDIF
655	IF( (ORDINT.EQ. VOISIN))THEN
656	DO 23010 N=1,NPTS
657	I = MIN(I2,MAX(I1,NINT(PX(N))))
658	J = MIN(J2,MAX(J1,NINT(PY(N))))
659	ZO(N) = Z(I,J)
660	23010 CONTINUE
661	ENDIF
662	RETURN
663	END
664	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
665	C
666	SUBROUTINE GGLL2GD(X,Y,XLAT,XLON,NPTS,NI,NJ,HEM,LROOTS)
667	C* -------------------------------------------------------------------
668	C* S/R GGLL2GD - COMPUTES THE GRID CO-ORDINATES OF A POINT ON
669	C A GAUSSIAN GRID
670	C* -------------------------------------------------------------------
671	USE mod_kinds_oasis
672	C* -------------------------------------------------------------------
673	IMPLICIT NONE
674	C* --------------------------------------------------------------------
675	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
676	PARAMETER (GLOBAL = 0)
677	PARAMETER (NORD = 1)
678	PARAMETER (SUD = 2)
679	PARAMETER (SUDNORD= 0)
680	PARAMETER (NORDSUD= 1)
681	INTEGER (kind=ip_intwp_p) NPTS, NI, NJ
682	REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), XLAT(NPTS), XLON(NPTS)
683	REAL(kind=ip_realwp_p) LROOTS(NJ)
684	REAL(kind=ip_realwp_p) ABS,DEL,EPSPHI
685	INTEGER (kind=ip_intwp_p) HEM
686	INTEGER (kind=ip_intwp_p) I,J
687	INTEGER (kind=ip_intwp_p) GUESS
688	INTEGER (kind=ip_intwp_p) NJHEM
689	REAL(kind=ip_realwp_p) TMPLAT, DELLAT, DELLON, XLAT0, XLON0
690	IF( (HEM.EQ. GLOBAL))THEN
691	NJHEM = NJ / 2
692	ELSE
693	NJHEM = NJ
694	ENDIF
695	DELLON = 360.0 / REAL(NI,kind=ip_realwp_p)
696	XLON0 = 0.0
697	DELLAT = 90.0 / REAL(NJHEM,kind=ip_realwp_p)
698	XLAT0 = DELLAT * 0.5
699	CVDIR NOVECTOR
700
701	DO 23002 I = 1, NPTS
702	X(I) = (XLON(I) - XLON0)/DELLON + 1.0
703	TMPLAT = XLAT(I)
704	IF( (TMPLAT.LT. 0))THEN
705	TMPLAT = -TMPLAT
706	ENDIF
707	IF( (TMPLAT.GT. LROOTS(NJHEM)))THEN
708	Y(I)=NJHEM+0.5*(TMPLAT-LROOTS(NJHEM))/(90.0-LROOTS(NJHEM
709	$ ))
710	ELSE
711	IF( (TMPLAT.LT. LROOTS(1)))THEN
712	Y(I)=0.5+(0.5*TMPLAT/LROOTS(1))
713	ELSE
714	GUESS=INT((TMPLAT-XLAT0)/DELLAT + 1.0)
715	23010 IF((LROOTS(GUESS+1).LT.TMPLAT))THEN
716	GUESS = GUESS+1
717	GOTO 23010
718	ENDIF
719	23012 IF((LROOTS(GUESS).GE.TMPLAT))THEN
720	GUESS = GUESS-1
721	GOTO 23012
722	ENDIF
723	Y(I)=GUESS+(TMPLAT-LROOTS(GUESS))/(LROOTS(GUESS+1)-
724	$ LROOTS(GUESS))
725	ENDIF
726	ENDIF
727	IF( (HEM.EQ. GLOBAL))THEN
728	IF( (0.EQ. MOD(NJ, 2)))THEN
729	IF( (XLAT(I).GE. 0.0))THEN
730	Y(I) = Y(I) + NJHEM
731	ELSE
732	Y(I) = NJHEM - Y(I) + 1
733	ENDIF
734	ELSE
735	IF( (XLAT(I).GE. 0.0))THEN
736	Y(I) = Y(I) + NJHEM + 1
737	ELSE
738	Y(I) = NJHEM - Y(I) + 1
739	ENDIF
740	ENDIF
741	ENDIF
742	IF( (HEM.EQ. NORD))THEN
743	IF( (XLAT(I).LT. 0.0))THEN
744	Y(I) = -Y(I) + 1
745	ENDIF
746	ENDIF
747	IF( (HEM.EQ. SUD))THEN
748	IF( (XLAT(I).GE. 0.0))THEN
749	Y(I) = Y(I) + NJ
750	ELSE
751	Y(I) = NJ - Y(I) + 1.0
752	ENDIF
753	ENDIF
754	23002 CONTINUE
755	RETURN
756	END
757	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
758	C
759	BLOCK DATA QQQCOMBLK
760	USE mod_kinds_oasis
761	IMPLICIT NONE
762	C* ----------------------------------------------------------------------
763	C* COMMON QQQCOM1
764	C
765	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
766	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
767	$ AYTMP,CXTMP,CYTMP
768	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
769	$ SPOLPTS
770	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
771	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
772	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
773	$ NPOLMAX,SPOLMAX
774	LOGICAL PTPOLN, PTPOLS
775	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
776	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
777	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
778	DATA XGDPTR, YGDPTR, ZPTR / 0, 0, 0 /
779	DATA AXPTR, AYPTR, CXPTR, CYPTR / 0, 0, 0, 0 /
780	DATA NPOLPTR,SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX / 0,
781	$ 0, 0, 0, 0, 0 /
782	DATA PTPOLN, PTPOLS / .FALSE., .FALSE. /
783	C
784	C-- COMMON QQQCOM2, QQQCOM3
785	C-
786
787	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
788	INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
789	INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
790	COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
791	$,LSTIG3,LSTIG4,LSTXOR
792	CHARACTER*1 LSTGTYP
793	COMMON /QQQCOM3/ LSTGTYP
794	DATA LSTLILJ, LSTNI, LSTNJ, LSTIG1, LSTIG2, LSTIG3, LSTIG4 /
795	$ 0, 0, 0, 0, 0, 0, 0 /
796	DATA LSTXOR, LSTGTYP / 0, ' ' /
797	END
798	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
799	C
800	SUBROUTINE IGSCINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,GRREF,IG1
801	$,IG2,IG3,IG4,SYM,AX,AY,VECFLG)
802	USE mod_kinds_oasis
803	IMPLICIT NONE
804	INTEGER (kind=ip_intwp_p) LI,LJ,NI,NJ,IG1,IG2,IG3,IG4,BIDON
805	REAL(kind=ip_realwp_p) ZO(LI,LJ),XLAT(LI,LJ),XLON(LI,LJ)
806	REAL(kind=ip_realwp_p) ZI(NI,NJ),AX(NI),AY(NJ)
807	LOGICAL SYM
808	CHARACTER*1 GRTYP, GRREF
809	LOGICAL QQQCMP
810	EXTERNAL QQQCMP, RGOPTI
811	EXTERNAL RGFREE, QQQCHK, QQQALOC, QQQKEEP, LL2RGD, LL2IGD
812	EXTERNAL XPNCOF, XPNGD, NWTNCOF, IRGDINT, RGDINT
813	EXTERNAL XPNAXEZ,XPNAXEG,GGGDINT,QQQGLAT,QQQXTRAP,POLRINT
814	EXTERNAL XPNAXEY
815	C
816	CDEFINITION DES VARIABLES LOCALES
817	C
818	C-- COMMON QQQCOM1
819	C--
820	C
821	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
822	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
823	$ AYTMP,CXTMP,CYTMP
824	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
825	$ SPOLPTS
826	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
827	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
828	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
829	$ NPOLMAX,SPOLMAX
830	LOGICAL PTPOLN, PTPOLS
831	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
832	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
833	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
834	C----------------------
835	C-- COMMON QQQCOM2
836	C-
837
838	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
839	INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
840	INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
841	COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
842	$,LSTIG3,LSTIG4,LSTXOR
843	C-
844	C------------
845	C* COMMON QQQCOM3
846	C
847
848	CHARACTER*1 LSTGTYP
849	COMMON /QQQCOM3/ LSTGTYP
850	C------------
851
852	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
853	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
854	PARAMETER (VOISIN = 0)
855	PARAMETER (LINEAIR = 1)
856	PARAMETER (CUBIQUE = 3)
857	PARAMETER (OUI = 1)
858	PARAMETER (MINIMUM = 2)
859	PARAMETER (MAXIMUM = 3)
860	PARAMETER (VALEUR = 4)
861	PARAMETER (ABORT = 13)
862	LOGICAL FLGXTRAP
863	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
864	REAL(kind=ip_realwp_p) VALXTRAP
865	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
866	EXTERNAL QQQCOMBLK
867	LOGICAL OK
868	INTEGER (kind=ip_intwp_p) XORSUM
869	LOGICAL VECTEUR
870	C
871	C* Test if field is scalar or vector
872	C* ---------------------------------
873	C
874	CHARACTER*6 VECFLG
875	IF (VECFLG .EQ. 'SCALAR') THEN
876	VECTEUR = .FALSE.
877	ELSEIF (VECFLG .EQ. 'VECTOR') THEN
878	VECTEUR = .TRUE.
879	ELSE
880	WRITE(6,*) ' WARNING !!! Wrong value for VECFLG '
881	WRITE(6,*) ' WE STOP IN igscint '
882	WRITE(6,*) ' '
883	STOP
884	ENDIF
885	C*----------------------
886	C
887	CVERIFICATION DU TYPE DE GRILLE
888	C
889
890	IF( (GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y')
891	$ .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
892	$ GRREF.EQ. 'L')) THEN
893	GOTO 10
894	ELSE
895	GOTO 101
896	ENDIF
897	ENTRY IGVINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,GRREF,IG1,IG2,
898	$IG3,IG4,SYM,AX,AY)
899	VECTEUR = .TRUE.
900	IF( (GRTYP.EQ.'Z' .OR. GRTYP .EQ. 'Y')
901	$ .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
902	$ GRREF.EQ. 'L')) THEN
903	GOTO 10
904	ELSE
905	GOTO 101
906	ENDIF
907	101 CONTINUE
908	IF( GRTYP .NE. 'Z' .OR. GRTYP .NE. 'Y')THEN
909	WRITE(6,*) 'Error, bad grid type (GRTYP)(IGSCINT)'
910	WRITE(6,*) '''GRTYP'' should be equal to ''Z'' or ''Y'''
911	ELSE
912	WRITE(6,*)
913	$ 'Error, bad grid type (GRREF)(IGSCINT)'
914	WRITE(6,*)
915	$ '''GRREF'' should be equal to ''N'', ''S'' or ''L'''
916	ENDIF
917	RETURN
918	ENTRY RGVINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,IG4
919	$,SYM)
920	VECTEUR = .TRUE.
921	GOTO 3
922	ENTRY RGSCINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
923	$IG4,SYM,VECFLG)
924	C
925	C* Test if field is scalar or vector
926	C* ---------------------------------
927	C
928	IF (VECFLG .EQ. 'SCALAR') THEN
929	VECTEUR = .FALSE.
930	ELSEIF (VECFLG .EQ. 'VECTOR') THEN
931	VECTEUR = .TRUE.
932	ELSE
933	WRITE(6,*) ' WARNING !!! Wrong value for VECFLG '
934	WRITE(6,*) ' WE STOP IN igscint '
935	WRITE(6,*) ' '
936	STOP
937	ENDIF
938	C*
939	3 CONTINUE
940	IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'L'.OR.GRTYP.EQ.
941	$'N'.OR.GRTYP.EQ.'S' .OR.GRTYP.EQ.'G'))THEN
942	GOTO 10
943	ENDIF
944	WRITE (6,200)
945	200 FORMAT('0','Error bad grid type (GRTYP)(RGSCINT)')
946	RETURN
947	10 CALL RGOPTI('INTERP', BIDON, .FALSE.)
948	OK = QQQCMP(XLAT,XLON,LI*LJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM
949	$)
950	IF( (.NOT.OK))THEN
951	CALL XPNCOF(I1,I2,J1,J2,NI,NJ,GRTYP,GRREF,IG1,IG1,IG3,IG4,
952	$ SYM,AX,AY)
953	CALL RGFREE
954	CALL QQQALOC(XLAT, XLON, LI, LJ, I1, I2, J1, J2)
955	CALL QQQKEEP(LI*LJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM)
956	CALL QQQCHK(XLAT, XLON, LI*LJ)
957	IF( GRTYP.EQ. 'Z' .OR. GRTYP.EQ. 'Y')THEN
958	IF( GRTYP.EQ. 'Z') CALL XPNAXEZ(AXTMP(AXPTR),AYTMP(AYPTR),
959	$ I1,I2,J1,J2,AX,AY,NI,NJ)
960	IF( GRTYP.EQ. 'Y') CALL XPNAXEY(AXTMP(AXPTR),AYTMP(AYPTR),
961	$ I1,I2,J1,J2,AX,AY,NI,NJ)
962	CALL NWTNCOF(CXTMP(CXPTR),CYTMP(CYPTR),AXTMP(AXPTR),
963	$ AYTMP(AYPTR),I1,I2,J1,J2)
964	CALL LL2IGD(XGD(XGDPTR),YGD(YGDPTR),XLAT,XLON,LI*LJ,NI,
965	$ NJ,GRTYP,GRREF,IG1,IG2,IG3,IG4,SYM,AX,AY)
966	ELSE
967	IF( (GRTYP.EQ. 'G'))THEN
968	CALL QQQGLAT(NJ,IG1)
969	CALL XPNAXEG(AXTMP(AXPTR),AYTMP(AYPTR),I1,I2,J1,J2,NI
970	$ ,NJ,IG1)
971	CALL NWTNCOF(CXTMP(CXPTR),CYTMP(CYPTR),AXTMP(AXPTR),
972	$ AYTMP(AYPTR),I1,I2,J1,J2)
973	ENDIF
974	CALL LL2RGD(XGD(XGDPTR),YGD(YGDPTR),XLAT,XLON,LI*LJ,NI,
975	$ NJ,GRTYP,IG1,IG2,IG3,IG4,SYM)
976	ENDIF
977	ENDIF
978	CALL XPNGD(ZTMP(ZPTR),I1,I2,J1,J2,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
979	$IG4,SYM,VECTEUR)
980	IF( (GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
981	CALL IRGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,AXTMP(AXPTR),
982	$ AYTMP(AYPTR),CXTMP(CXPTR),CYTMP(CYPTR),ZTMP(ZPTR),I1,I2,J1,
983	$ J2)
984	ENDIF
985	IF( (GRTYP.EQ. 'G'))THEN
986	CALL GGGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,AYTMP(AYPTR),
987	$ CYTMP(CYPTR),ZTMP(ZPTR),I1,I2,J1,J2)
988	ENDIF
989	IF( (GRTYP.NE. 'G'.AND. GRTYP.NE. 'Z'
990	$ .AND. GRTYP .NE. 'Y'))THEN
991	CALL RGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,ZTMP(ZPTR),I1,
992	$ I2,J1,J2)
993	ENDIF
994	IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
995	$ GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
996	CALL QQQXTRAP(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,ZTMP(ZPTR),
997	$ I1,I2,J1,J2)
998	ENDIF
999	IF( (PTPOLN.OR. PTPOLS))THEN
1000	CALL POLRINT(ZO,XLAT,XLON,LI*LJ,ZTMP(ZPTR),NI,NJ,I1,I2,J1,
1001	$ J2,GRTYP,IG1,IG2,IG3,IG4,VECTEUR)
1002	ENDIF
1003	RETURN
1004	END
1005	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1006	C
1007	SUBROUTINE IGUVINT(SPDO,PSIO,LI,LJ,XLAT,XLON,UI,VI,NI,NJ,GRTYP
1008	$,GRREF,IG1,IG2,IG3,IG4,SWS,AX,AY)
1009	USE mod_kinds_oasis
1010	IMPLICIT NONE
1011	INTEGER (kind=ip_intwp_p) LI,LJ,NI,NJ,IG1,IG2,IG3,IG4,BIDON
1012	REAL(kind=ip_realwp_p) SPDO(LI,LJ),PSIO(LI,LJ),XLAT(LI,LJ),
1013	$ XLON(LI,LJ)
1014	REAL(kind=ip_realwp_p) UI(NI,NJ),VI(NI,NJ),AX(NI),AY(NJ)
1015	LOGICAL SWS
1016	CHARACTER*1 GRTYP, GRREF
1017	EXTERNAL GDW2LLW,IGSCINT,RGVINT,MODULE,RGOPTI,IGVINT
1018	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
1019	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
1020	PARAMETER (VOISIN = 0)
1021	PARAMETER (LINEAIR = 1)
1022	PARAMETER (CUBIQUE = 3)
1023	PARAMETER (OUI = 1)
1024	PARAMETER (MINIMUM = 2)
1025	PARAMETER (MAXIMUM = 3)
1026	PARAMETER (VALEUR = 4)
1027	PARAMETER (ABORT = 13)
1028	LOGICAL FLGXTRAP
1029	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
1030	REAL(kind=ip_realwp_p) VALXTRAP
1031	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
1032	C* -------------------------------------------------------------------
1033	C
1034	CVERIFICATION DU TYPE DE GRILLE
1035	C
1036
1037	CALL RGOPTI('INTERP', BIDON, .FALSE.)
1038	IF( (GRTYP.EQ.'Z' .OR. GRTYP .EQ. 'Y')
1039	$ .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
1040	$ GRREF.EQ. 'L')) THEN
1041	CALL IGVINT(SPDO,LI,LJ,XLAT,XLON,UI,NI,NJ,GRTYP,GRREF,IG1,
1042	$ IG2,IG3,IG4,.TRUE.,AX,AY)
1043	CALL IGVINT(PSIO,LI,LJ,XLAT,XLON,VI,NI,NJ,GRTYP,GRREF,IG1,
1044	$ IG2,IG3,IG4,.FALSE.,AX,AY)
1045	IF( (SWS))THEN
1046	CALL GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRREF,IG1,IG2,IG3,IG4)
1047	ELSE
1048	CALL MODULE(SPDO,PSIO,LI,LJ)
1049	ENDIF
1050	ELSE
1051	IF( (GRTYP.NE. 'Z' .AND. GRTYP .NE. 'Y'))THEN
1052	WRITE(6,*)
1053	$ 'Error, Bad grid type (GRTYP)(IGSCINT)'
1054	WRITE(6,*)
1055	$ '''GRTYP'' should be equal to ''Z'' or ''Y'''
1056	ELSE
1057	WRITE(6,*)
1058	$ 'Error, Bad reference grid type(GRREF)(IGSCINT)'
1059	WRITE(6,*)
1060	$ '''GRREF'' should be equal to ''N'', ''S'' or ''L'''
1061	ENDIF
1062	ENDIF
1063	RETURN
1064	ENTRY RGUVINT(SPDO,PSIO,LI,LJ,XLAT,XLON,UI,VI,NI,NJ,GRTYP,IG1,
1065	$IG2,IG3,IG4,SWS)
1066	CALL RGOPTI('INTERP', BIDON, .FALSE.)
1067	IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'L'.OR.GRTYP.EQ.
1068	$'N'.OR.GRTYP.EQ.'S' .OR.GRTYP.EQ.'G'))THEN
1069	CALL RGVINT(SPDO,LI,LJ,XLAT,XLON,UI,NI,NJ,GRTYP,IG1,IG2,IG3
1070	$ ,IG4,.TRUE.)
1071	CALL RGVINT(PSIO,LI,LJ,XLAT,XLON,VI,NI,NJ,GRTYP,IG1,IG2,IG3
1072	$ ,IG4,.FALSE.)
1073	IF( (SWS))THEN
1074	CALL GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
1075	ELSE
1076	CALL MODULE(SPDO,PSIO,LI,LJ)
1077	ENDIF
1078	ELSE
1079	WRITE (6,200)
1080	200 FORMAT('0','Error bad grid type (GRTYP)(RGSCINT)')
1081	ENDIF
1082	RETURN
1083	END
1084	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1085	C
1086	SUBROUTINE IRGDINT(ZO,PX,PY,NPTS,AX,AY,CX,CY,Z,I1,I2,J1,J2)
1087	USE mod_kinds_oasis
1088	IMPLICIT NONE
1089	C*******
1090	C AUTEUR: Y.CHARTIER, DRPN
1091	C FEVRIER 1991
1092	C
1093	C OBJET: INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR
1094	C D'UNE GRILLE SOURCE IRREGULIERE.
1095	C*******
1096	INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
1097	REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
1098	REAL(kind=ip_realwp_p) FA, FA2, FA3, FA4
1099	REAL(kind=ip_realwp_p) AX(I1:I2),AY(J1:J2),CX(I1:I2,6)
1100	REAL(kind=ip_realwp_p) CY(J1:J2,6)
1101	REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
1102	C
1103	C NPTS : NOMBRE DE POINTS A INTERPOLER
1104	C I1:I2 : DIMENSION DE LA GRILLE SOURCE SELON X
1105	C J1:J2 : DIMENSION DE LA GRILLE SOURCE SELON Y
1106	C ZO : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
1107	C PX : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
1108	C VEUT INTERPOLER
1109	C PY : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
1110	C VEUT INTERPOLER
1111	C AX : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES X.
1112	C AY : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES Y.
1113	C CX : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON X.
1114	C CY : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON Y.
1115	C Z : VALEURS DE LA GRILLE SOURCE.
1116	C
1117	C***************************************************************************
1118	C
1119	C * * * *
1120	C
1121	C * * * *
1122	C # ==> PT (X,Y)
1123	C * (=) * * ==> = PT (I, J)
1124	C
1125	C * * * *
1126	C
1127	C
1128	C
1129	C CX(I,1) = 1.0 / (X2-X1)
1130	C CX(I,2) = 1.0 / (X3-X1)
1131	C CX(I,3) = 1.0 / (X3-X2)
1132	C CX(I,4) = 1.0 / (X4-X1)
1133	C CX(I,5) = 1.0 / (X4-X2)
1134	C CX(I,6) = 1.0 / (X4-X3)
1135	C
1136	C STRUCTURE IDENTIQUE POUR CY(J,1..6)
1137
1138	REAL(kind=ip_realwp_p) A11,A12,A13,A14,A21,A22,A23,A24
1139	REAL(kind=ip_realwp_p) A31,A32,A33,A34,A41,A42,A43,A44
1140	REAL(kind=ip_realwp_p) B1,B2,B3,B4,B11,B12,B13,B14
1141	REAL(kind=ip_realwp_p) X1,X2,X3,X4,Y1,Y2,Y3,Y4
1142	INTEGER (kind=ip_intwp_p) I, J, N
1143	REAL(kind=ip_realwp_p) A1,A2,A3,A4,X,Y,C1,C2,C3,C4,C5,C6
1144	C DEFINITION DES FONCTIONS IN-LINE
1145
1146	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
1147	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
1148	PARAMETER (VOISIN = 0)
1149	PARAMETER (LINEAIR = 1)
1150	PARAMETER (CUBIQUE = 3)
1151	PARAMETER (OUI = 1)
1152	PARAMETER (MINIMUM = 2)
1153	PARAMETER (MAXIMUM = 3)
1154	PARAMETER (VALEUR = 4)
1155	PARAMETER (ABORT = 13)
1156	LOGICAL FLGXTRAP
1157	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
1158	REAL(kind=ip_realwp_p) VALXTRAP
1159	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
1160	REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
1161	REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
1162	CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)0.1666666666666 + 0.5(Z2-Z3)
1163	$)DX + 0.5(Z1+Z3)-Z2)DX + Z3-0.1666666666666Z4-0.5*Z2-0.
1164	$3333333333333Z1)DX+Z2
1165	ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
1166	FA(A1,A2,A3,A4,X,X1,X2,X3)=A1+(X-X1)(A2+(X-X2)(A3+A4*(X-X3))
1167	$)
1168	FA2(C1,A1,A2)=C1*(A2-A1)
1169	FA3(C1,C2,C3,A1,A2,A3)=C2(C3(A3-A2)-C1*(A2-A1))
1170	FA4(C1,C2,C3,C4,C5,C6,A1,A2,A3,A4)=C4(C5(C6(A4-A3)-C3(A3-
1171	$A2)) - C2(C3(A3-A2)-C1*(A2-A1)))
1172	IF( (ORDINT.EQ. CUBIQUE))THEN
1173	DO 23002 N=1,NPTS
1174	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
1175	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
1176	X = AX(I) + (AX(I+1)-AX(I))*(PX(N)-I)
1177	Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
1178	X1=AX(I-1)
1179	X2=AX(I)
1180	X3=AX(I+1)
1181	X4=AX(I+2)
1182	Y1=AY(J-1)
1183	Y2=AY(J)
1184	Y3=AY(J+1)
1185	Y4=AY(J+2)
1186	C INTERPOLATION 1ERE RANGEE SELON X
1187
1188	Z1=Z(I-1,J-1)
1189	Z2=Z(I ,J-1)
1190	Z3=Z(I+1,J-1)
1191	Z4=Z(I+2,J-1)
1192	A11 = Z1
1193	A12 = FA2(CX(I,1),Z1,Z2)
1194	A13 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
1195	A14 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
1196	$ ),Z1,Z2,Z3,Z4)
1197	B1 = FA(A11,A12,A13,A14,X,X1,X2,X3)
1198	C INTERPOLATION 2EME RANGEE SELON X
1199
1200	Z1=Z(I-1,J)
1201	Z2=Z(I ,J)
1202	Z3=Z(I+1,J)
1203	Z4=Z(I+2,J)
1204	A21 = Z1
1205	A22 = FA2(CX(I,1),Z1,Z2)
1206	A23 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
1207	A24 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
1208	$ ),Z1,Z2,Z3,Z4)
1209	B2 = FA(A21,A22,A23,A24,X,X1,X2,X3)
1210	C INTERPOLATION 3EME RANGEE SELON X
1211
1212	Z1=Z(I-1,J+1)
1213	Z2=Z(I ,J+1)
1214	Z3=Z(I+1,J+1)
1215	Z4=Z(I+2,J+1)
1216	A31 = Z1
1217	A32 = FA2(CX(I,1),Z1,Z2)
1218	A33 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
1219	A34 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
1220	$ ),Z1,Z2,Z3,Z4)
1221	B3 = FA(A31,A32,A33,A34,X,X1,X2,X3)
1222	C INTERPOLATION 4EME RANGEE SELON X
1223
1224	Z1=Z(I-1,J+2)
1225	Z2=Z(I ,J+2)
1226	Z3=Z(I+1,J+2)
1227	Z4=Z(I+2,J+2)
1228	A41 = Z1
1229	A42 = FA2(CX(I,1),Z1,Z2)
1230	A43 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
1231	A44 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
1232	$ ),Z1,Z2,Z3,Z4)
1233	B4 = FA(A41,A42,A43,A44,X,X1,X2,X3)
1234	C INTERPOLATION FINALE SELON Y
1235
1236	B11 = B1
1237	B12 = FA2(CY(J,1),B1,B2)
1238	B13 = FA3(CY(J,1),CY(J,2),CY(J,3),B1,B2,B3)
1239	B14 = FA4(CY(J,1),CY(J,2),CY(J,3),CY(J,4),CY(J,5),CY(J,6
1240	$ ),B1,B2,B3,B4)
1241	ZO(N) = FA(B11,B12,B13,B14,Y,Y1,Y2,Y3)
1242	23002 CONTINUE
1243	ENDIF
1244	IF( (ORDINT.EQ. LINEAIR))THEN
1245	DO 23006 N=1,NPTS
1246	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
1247	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
1248	X = AX(I) + (AX(I+1)-AX(I))*(PX(N)-I)
1249	Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
1250	DX = (X - AX(I))/(AX(I+1)-AX(I))
1251	DY = (Y - AY(J))/(AY(J+1)-AY(J))
1252	Y1 = ZLIN(Z(I,J),Z(I+1,J),DX)
1253	Y2 = ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
1254	ZO(N) = ZLIN(Y1,Y2,DY)
1255	23006 CONTINUE
1256	ENDIF
1257	IF( (ORDINT.EQ. VOISIN))THEN
1258	DO 23010 N=1,NPTS
1259	I = MIN(I2,MAX(I1,NINT(PX(N))))
1260	J = MIN(J2,MAX(J1,NINT(PY(N))))
1261	ZO(N) = Z(I,J)
1262	23010 CONTINUE
1263	ENDIF
1264	RETURN
1265	END
1266	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1267	C
1268	SUBROUTINE LL2IGD(PX,PY,XLAT,XLON,NPTS,NI,NJ,GRTYP,GRREF,IG1,
1269	$IG2,IG3,IG4,SYM,AX,AY)
1270	C* ----------------------------------------------------------------------
1271	C**S/R LL2IGD - CONVERSION DE COORDONNEES LAT-LON A PTS DE GRILLE
1272	C* ----------------------------------------------------------------------
1273	USE mod_kinds_oasis
1274	IMPLICIT NONE
1275	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
1276	PARAMETER (GLOBAL = 0)
1277	PARAMETER (NORD = 1)
1278	PARAMETER (SUD = 2)
1279	PARAMETER (SUDNORD= 0)
1280	PARAMETER (NORDSUD= 1)
1281	EXTERNAL CIGAXG,VXYFLL,LLLL2GD,PERMUT,CHKXTRAP
1282	INTEGER (kind=ip_intwp_p) I, NPTS, NI, NJ
1283	REAL(kind=ip_realwp_p) PX(NPTS),PY(NPTS),XLAT(NPTS),XLON(NPTS)
1284	CHARACTER*1 GRTYP, GRREF
1285	INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
1286	LOGICAL SYM
1287	REAL(kind=ip_realwp_p) AX(NI), AY(NJ)
1288	REAL(kind=ip_realwp_p) PI,PJ,DGRW,D60,TMP
1289	REAL(kind=ip_realwp_p) DLAT, DLON, XLAT0, XLON0
1290	INTEGER (kind=ip_intwp_p) INDX, INDY
1291	INTEGER (kind=ip_intwp_p) CHERCHE, FINDLON
1292	EXTERNAL CHERCHE, FINDLON
1293	IF( (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR. GRREF.EQ. 'L')
1294	$)THEN
1295	IF( (GRREF.EQ. 'N'))THEN
1296	CALL CIGAXG(GRREF, PI, PJ, D60, DGRW, IG1, IG2, IG3,
1297	$ IG4)
1298	CALL VXYFLL(PX, PY, XLAT, XLON, NPTS, D60,DGRW,PI,PJ,1)
1299	ELSE
1300	IF( (GRREF.EQ. 'S'))THEN
1301	CALL CIGAXG(GRREF, PI, PJ, D60, DGRW, IG1, IG2, IG3,
1302	$ IG4)
1303	CALL VXYFLL(PX, PY, XLAT, XLON, NPTS, D60,DGRW,PI,PJ,
1304	$ 2)
1305	ELSE
1306	CALL CIGAXG(GRREF, XLAT0, XLON0, DLAT, DLON, IG1, IG2
1307	$ , IG3, IG4)
1308	CALL LLLL2GD(PX, PY, XLAT, XLON, NPTS, XLAT0, XLON0,
1309	$ DLAT, DLON)
1310	ENDIF
1311	ENDIF
1312	DO 23006 I=1,NPTS
1313	INDX = FINDLON(PX(I), AX, NI, TMP)
1314	INDY = CHERCHE(PY(I), AY, NJ)
1315	IF( (INDX.GE. NI))THEN
1316	INDX = NI - 1
1317	ENDIF
1318	IF( (INDY.GE. NJ))THEN
1319	INDY = NJ - 1
1320	ENDIF
1321	PX(I)=REAL(INDX,kind=ip_realwp_p)+
1322	$ (TMP-AX(INDX))/(AX(INDX+1)-AX(INDX))
1323	PY(I) = REAL(INDY,kind=ip_realwp_p)
1324	$ + (PY(I) - AY(INDY))/(AY(INDY+1)-AY(INDY))
1325	23006 CONTINUE
1326	CALL CHKXTRAP(PX, PY, NPTS, NI, NJ)
1327	ENDIF
1328	RETURN
1329	END
1330	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1331	C
1332	SUBROUTINE LL2RGD(PX,PY,XLAT,XLON,NPTS,NI,NJ,GRTYP,IG1,IG2,IG3
1333	$,IG4,SYM)
1334	C* -----------------------------------------------------------------
1335	C* S/R LL2RGD - CONVERSION DE COORDONNEES LAT-LON A PTS DE GRILLE
1336	C* ----------------------------------------------------------------------
1337	USE mod_kinds_oasis
1338	IMPLICIT NONE
1339	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
1340	PARAMETER (GLOBAL = 0)
1341	PARAMETER (NORD = 1)
1342	PARAMETER (SUD = 2)
1343	PARAMETER (SUDNORD= 0)
1344	PARAMETER (NORDSUD= 1)
1345	C
1346	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
1347	C
1348	C RDTODG = 180/PIE, DGTORD = PIE/180
1349
1350	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
1351	DATA PIE /3.14159265358979323846/
1352	DATA RDTODG /57.295779513082/
1353	DATA DGTORD /1.7453292519943E-2/
1354	C
1355	C-- COMMON GAUSSGD
1356	C REAL(kind=ip_realwp_p) ROOTS(1),LROOTS(1)
1357	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
1358	INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
1359	COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
1360	C----------------------
1361
1362	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
1363	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
1364	PARAMETER (VOISIN = 0)
1365	PARAMETER (LINEAIR = 1)
1366	PARAMETER (CUBIQUE = 3)
1367	PARAMETER (OUI = 1)
1368	PARAMETER (MINIMUM = 2)
1369	PARAMETER (MAXIMUM = 3)
1370	PARAMETER (VALEUR = 4)
1371	PARAMETER (ABORT = 13)
1372	LOGICAL FLGXTRAP
1373	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
1374	REAL(kind=ip_realwp_p) VALXTRAP
1375	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
1376	EXTERNAL CIGAXG,VXYFLL,LLLL2GD,PERMUT
1377	EXTERNAL QQQGLAT,GGLL2GD,CHKXTRAP
1378	INTEGER (kind=ip_intwp_p) NPTS, NI, NJ
1379	REAL(kind=ip_realwp_p) PX(NPTS),PY(NPTS),XLAT(NPTS),XLON(NPTS)
1380	CHARACTER*1 GRTYP
1381	INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
1382	LOGICAL SYM
1383	REAL(kind=ip_realwp_p) PI,PJ,DGRW,D60
1384	REAL(kind=ip_realwp_p) DELLAT, DELLON, XLAT0, XLON0
1385	INTEGER (kind=ip_intwp_p) I,J
1386	IF( (GRTYP.EQ.'N'))THEN
1387	CALL CIGAXG(GRTYP,PI,PJ,D60,DGRW,IG1,IG2,IG3,IG4)
1388	CALL VXYFLL(PX,PY,XLAT,XLON,NPTS,D60,DGRW,PI,PJ,NORD)
1389	CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
1390	ENDIF
1391	IF( (GRTYP.EQ.'S'))THEN
1392	CALL CIGAXG(GRTYP,PI,PJ,D60,DGRW,IG1,IG2,IG3,IG4)
1393	CALL VXYFLL(PX,PY,XLAT,XLON,NPTS,D60,DGRW,PI,PJ,SUD)
1394	CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
1395	ENDIF
1396	IF( (GRTYP.EQ.'A'))THEN
1397	DELLON = 360.0 / REAL(NI,ip_realwp_p)
1398	XLON0 = 0.0
1399	IF( (IG1.EQ. GLOBAL))THEN
1400	DELLAT = 180.0 / REAL(NJ,ip_realwp_p)
1401	XLAT0 = -90.0 + DELLAT * 0.5
1402	ENDIF
1403	IF( (IG1.EQ. NORD))THEN
1404	DELLAT = 90.0 / REAL(NJ,ip_realwp_p)
1405	XLAT0 = DELLAT * 0.5
1406	ENDIF
1407	IF( (IG1.EQ. SUD))THEN
1408	DELLAT = 90.0 / REAL(NJ,ip_realwp_p)
1409	XLAT0 = -90.0 + DELLAT * 0.5
1410	ENDIF
1411	FLGXTRAP = .FALSE.
1412	CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
1413	$ DELLON)
1414	ENDIF
1415	IF( (GRTYP.EQ.'B'))THEN
1416	DELLON = 360.0 / REAL(NI-1,ip_realwp_p)
1417	XLON0 = 0.0
1418	IF( (IG1.EQ. GLOBAL))THEN
1419	DELLAT = 180.0 / REAL(NJ-1,ip_realwp_p)
1420	XLAT0 = -90.0
1421	ENDIF
1422	IF( (IG1.EQ. NORD))THEN
1423	DELLAT = 90.0 / REAL(NJ-1,ip_realwp_p)
1424	XLAT0 = 0.0
1425	ENDIF
1426	IF( (IG1.EQ. SUD))THEN
1427	DELLAT = 90.0 / REAL(NJ-1,ip_realwp_p)
1428	XLAT0 = -90.0
1429	ENDIF
1430	FLGXTRAP = .FALSE.
1431	CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
1432	$ DELLON)
1433	ENDIF
1434	IF( (GRTYP.EQ. 'G'))THEN
1435	FLGXTRAP = .FALSE.
1436	CALL GGLL2GD(PX,PY,XLAT,XLON,NPTS,NI,NJ,IG1,LROOTS(ILROOTS)
1437	$ )
1438	ENDIF
1439	IF( (GRTYP.EQ. 'L'))THEN
1440	CALL CIGAXG(GRTYP,XLAT0,XLON0,DELLAT,DELLON,IG1,IG2,IG3,IG4
1441	$ )
1442	DO 23024 I=1,NPTS
1443	IF( (XLON(I).LT. XLON0))THEN
1444	XLON(I) = XLON(I) + 360.0
1445	ENDIF
1446	IF( (XLON(I).GT. (XLON0 + NI*DELLON)))THEN
1447	XLON(I) = XLON(I) - 360.0
1448	ENDIF
1449	23024 CONTINUE
1450	CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
1451	$ DELLON)
1452	CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
1453	ENDIF
1454	RETURN
1455	END
1456	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1457	C
1458	SUBROUTINE LLLL2GD(X,Y,DLAT,DLON,NPTS,XLAT0,XLON0, DELLAT,
1459	$ DELLON)
1460	C* ----------------------------------------------------------------------
1461	C* S/R LLLL2GD - COMPUTES THE GRID CO-ORDINATES OF A POINT
1462	C* ----------------------------------------------------------------------
1463	USE mod_kinds_oasis
1464	IMPLICIT NONE
1465	INTEGER (kind=ip_intwp_p) NPTS
1466	REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), DLAT(NPTS), DLON(NPTS)
1467	REAL(kind=ip_realwp_p) XLAT0, XLON0, DELLAT, DELLON
1468	INTEGER (kind=ip_intwp_p) I
1469	DO 23000 I=1,NPTS
1470	X(I) = (DLON(I) - XLON0)/DELLON + 1.0
1471	Y(I) = (DLAT(I) - XLAT0)/DELLAT + 1.0
1472	23000 CONTINUE
1473	RETURN
1474	END
1475	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1476	C
1477	SUBROUTINE MODULE(U,V,LI,LJ)
1478	USE mod_kinds_oasis
1479	INTEGER (kind=ip_intwp_p) LI,LJ
1480	REAL(kind=ip_realwp_p) U(LI,LJ), V(LI,LJ)
1481	INTEGER (kind=ip_intwp_p) I,J
1482	DO 23000 I=1,LI
1483	DO 23002 J=1,LJ
1484	U(I,J)=SQRT(U(I,J)U(I,J)+V(I,J)V(I,J))
1485	23002 CONTINUE
1486	23000 CONTINUE
1487	RETURN
1488	END
1489	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1490	C
1491	SUBROUTINE NWTNCOF(CX,CY,AX,AY,I1,I2,J1,J2)
1492	C*******
1493	C AUTEUR: Y. CHARTIER, DRPN
1494	C FEVRIER 1991
1495	C
1496	C OBJET: CALCUL DE COEFFICIENTS UTI1ISES DANS LA FORME NEWTONIENNE
1497	C DE L'INTERPOLATION DE LAGRANGE.
1498	C
1499	C===========================================
1500	C
1501	C ------------------------#----------------------->
1502	C X1 X2 X X3 X4
1503	C
1504	C===========================================
1505	C CX(I,1) = 1.0 / (X2-X1)
1506	C CX(I,2) = 1.0 / (X3-X1)
1507	C CX(I,3) = 1.0 / (X3-X2)
1508	C CX(I,4) = 1.0 / (X4-X1)
1509	C CX(I,5) = 1.0 / (X4-X2)
1510	C CX(I,6) = 1.0 / (X4-X3)
1511	C
1512	C STRUCTURE IDENTIQUE POUR CY(J,1..6)
1513	C*******
1514	USE mod_kinds_oasis
1515	IMPLICIT NONE
1516	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
1517	REAL(kind=ip_realwp_p) CX(I1:I2,6),CY(J1:J2,6),AX(I1:I2),AY(J1:J2)
1518	INTEGER (kind=ip_intwp_p) I,J
1519	DO 23000 I=I1+1,I2-2
1520	CX(I,1) = 1. / (AX(I ) - AX(I-1))
1521	CX(I,2) = 1. / (AX(I+1) - AX(I-1))
1522	CX(I,3) = 1. / (AX(I+1) - AX(I ))
1523	CX(I,4) = 1. / (AX(I+2) - AX(I-1))
1524	CX(I,5) = 1. / (AX(I+2) - AX(I ))
1525	CX(I,6) = 1. / (AX(I+2) - AX(I+1))
1526	23000 CONTINUE
1527	DO 23002 J=J1+1,J2-2
1528	CY(J,1) = 1. / (AY(J ) - AY(J-1))
1529	CY(J,2) = 1. / (AY(J+1) - AY(J-1))
1530	CY(J,3) = 1. / (AY(J+1) - AY(J ))
1531	CY(J,4) = 1. / (AY(J+2) - AY(J-1))
1532	CY(J,5) = 1. / (AY(J+2) - AY(J ))
1533	CY(J,6) = 1. / (AY(J+2) - AY(J+1))
1534	23002 CONTINUE
1535	RETURN
1536	END
1537	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1538	C
1539	SUBROUTINE POLRINT(ZO,XLAT,XLON,LILJ,ZI,NI,NJ,I1,I2,J1,J2,
1540	$GRTYP,IG1,IG2,IG3,IG4,VECTEUR)
1541	USE mod_kinds_oasis
1542	IMPLICIT NONE
1543	INTEGER (kind=ip_intwp_p) LILJ,NI,NJ,IG1,IG2,IG3,IG4
1544	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
1545	REAL(kind=ip_realwp_p) ZO(LILJ),XLAT(LILJ),XLON(LILJ)
1546	REAL(kind=ip_realwp_p) ZI(I1:I2,J1:J2)
1547	CHARACTER*1 GRTYP
1548	LOGICAL VECTEUR
1549	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
1550	PARAMETER (GLOBAL = 0)
1551	PARAMETER (NORD = 1)
1552	PARAMETER (SUD = 2)
1553	PARAMETER (SUDNORD= 0)
1554	PARAMETER (NORDSUD= 1)
1555	C-- COMMON QQQCOM1
1556	C--
1557	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
1558	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
1559	$ AYTMP,CXTMP,CYTMP
1560	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
1561	$ SPOLPTS
1562	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
1563	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
1564	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
1565	$ NPOLMAX,SPOLMAX
1566	LOGICAL PTPOLN, PTPOLS
1567	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
1568	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
1569	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
1570	C----------------------
1571	C
1572	CDEFINITION DES VARIABLES LOCALES
1573	C
1574
1575	INTEGER (kind=ip_intwp_p) I,J
1576	INTEGER (kind=ip_intwp_p) N1, N2, S1, S2
1577	REAL(kind=ip_realwp_p) VPOLNOR, VPOLSUD
1578	IF( (VECTEUR))THEN
1579	RETURN
1580	ENDIF
1581	IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
1582	$ GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
1583	RETURN
1584	ENDIF
1585	IF( (GRTYP.EQ. 'B'))THEN
1586	IF( (IG1.EQ. GLOBAL))THEN
1587	VPOLNOR = ZI(0, NJ)
1588	VPOLSUD = ZI(0, 1)
1589	ENDIF
1590	IF( (IG1.EQ. NORD))THEN
1591	VPOLNOR = ZI(0, NJ)
1592	VPOLSUD = ZI(0, -NJ+1)
1593	ENDIF
1594	IF( (IG1.EQ. SUD))THEN
1595	VPOLNOR = ZI(0, 2*NJ-1)
1596	VPOLSUD = ZI(0, 1)
1597	ENDIF
1598	ENDIF
1599	IF( (GRTYP.EQ. 'A'.OR. GRTYP.EQ. 'G'))THEN
1600	VPOLNOR = 0.0
1601	VPOLSUD = 0.0
1602	IF( (IG1.EQ. GLOBAL))THEN
1603	N1 = NJ
1604	N2 = NJ - 1
1605	S1 = 1
1606	S2 = 2
1607	ENDIF
1608	IF( (IG1.EQ. NORD))THEN
1609	N1 = NJ
1610	N2 = NJ - 1
1611	S1 = -NJ+1
1612	S2 = -NJ+2
1613	ENDIF
1614	IF( (IG1.EQ. SUD))THEN
1615	N1 = 2*NJ
1616	N2 = 2*NJ - 1
1617	S1 = 1
1618	S2 = 2
1619	ENDIF
1620	IF( (PTPOLN))THEN
1621	DO 23022 I=1,NI
1622	VPOLNOR = VPOLNOR + 9.0 * ZI(I, N1) + ZI(I,N2)
1623	23022 CONTINUE
1624	VPOLNOR = 0.1 * VPOLNOR / REAL(NI,ip_realwp_p)
1625	ENDIF
1626	IF( (PTPOLS))THEN
1627	DO 23026 I=1,NI
1628	VPOLSUD = VPOLSUD + 9.0 * ZI(I, S1) + ZI(I,S2)
1629	23026 CONTINUE
1630	VPOLSUD = 0.1 * VPOLSUD / REAL(NI,ip_realwp_p)
1631	ENDIF
1632	ENDIF
1633	DO 23028 I=1,NPOLNUM
1634	ZO(NPOLPTS(NPOLPTR+I-1)) = VPOLNOR
1635	23028 CONTINUE
1636	DO 23030 I=1,SPOLNUM
1637	ZO(SPOLPTS(SPOLPTR+I-1)) = VPOLSUD
1638	23030 CONTINUE
1639	RETURN
1640	END
1641	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1642	C
1643	SUBROUTINE QQQALOC(XLAT, XLON, LI, LJ, I1, I2, J1, J2)
1644	USE mod_kinds_oasis
1645	USE memoir
1646	IMPLICIT NONE
1647	C EXTERNAL MEMOIRH
1648	INTEGER (kind=ip_intwp_p) LI,LJ,I1,I2,J1,J2
1649	REAL(kind=ip_realwp_p) XLAT(LI,LJ), XLON(LI,LJ)
1650	C-- COMMON QQQCOM1
1651	C--
1652	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
1653	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
1654	$ AYTMP,CXTMP,CYTMP
1655	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
1656	$ SPOLPTS
1657	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
1658	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
1659	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
1660	$ NPOLMAX,SPOLMAX
1661	LOGICAL PTPOLN, PTPOLS
1662	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
1663	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
1664	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
1665	C----------------------
1666	CALL MEMOIRH(XGD, XGDPTR, LI*LJ, 0)
1667	CALL MEMOIRH(YGD, YGDPTR, LI*LJ, 0)
1668	CALL MEMOIRH(ZTMP,ZPTR,(I2-I1+1)*(J2-J1+1), 0)
1669	CALL MEMOIRH(AXTMP,AXPTR,(I2-I1+1), 0)
1670	CALL MEMOIRH(AYTMP,AYPTR,(J2-J1+1), 0)
1671	CALL MEMOIRH(CXTMP,CXPTR, (I2-I1+1) * 6, 0)
1672	CALL MEMOIRH(CYTMP,CYPTR, (J2-J1+1) * 6, 0)
1673	CALL MEMOIRH(NPOLPTS, NPOLPTR, 128, 0)
1674	CALL MEMOIRH(SPOLPTS, SPOLPTR, 128, 0)
1675	NPOLMAX = 128
1676	SPOLMAX = 128
1677	RETURN
1678	END
1679	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1680	C
1681	SUBROUTINE QQQCHK(XLAT, XLON, LILJ)
1682	USE mod_kinds_oasis
1683	USE memoir
1684	IMPLICIT NONE
1685
1686	C-- COMMON QQQCOM1
1687	C--
1688	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
1689	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
1690	$ AYTMP,CXTMP,CYTMP
1691	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
1692	$ SPOLPTS
1693	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
1694	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
1695	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
1696	$ NPOLMAX,SPOLMAX
1697	LOGICAL PTPOLN, PTPOLS
1698	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
1699	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
1700	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
1701	C----------------------
1702
1703	INTEGER (kind=ip_intwp_p) LILJ
1704	REAL(kind=ip_realwp_p) XLAT(LILJ), XLON(LILJ)
1705	REAL(kind=ip_realwp_p) TMP
1706	LOGICAL BADLAT
1707	INTEGER (kind=ip_intwp_p) NPOLPT2,SPOLPT2
1708	INTEGER (kind=ip_intwp_p) I,J
1709	C VERIF POUR POINT AU POLE
1710
1711	REAL(kind=ip_realwp_p) EPSILON
1712	DATA EPSILON /1.0E-6/
1713	BADLAT = .FALSE.
1714	DO 23000 I=1,LILJ
1715	TMP = XLON(I)
1716	XLON(I)=MOD(MOD(XLON(I),360.0_ip_realwp_p)
1717	$ +360.0_ip_realwp_p,360.0_ip_realwp_p)
1718	TMP = XLAT(I)
1719	XLAT(I)=MAX(-90.0_ip_realwp_p,
1720	$ MIN(90.0_ip_realwp_p,XLAT(I)))
1721	IF( (TMP.NE. XLAT(I)))THEN
1722	BADLAT = .TRUE.
1723	ENDIF
1724	23000 CONTINUE
1725	IF( (BADLAT))THEN
1726	WRITE(6,*) '<QQQGLAT> Latitudes gotten > 90.0 or < -90.0'
1727	WRITE(6,*)
1728	$ 'Put latitudes back into -90.0 and +90.0'
1729	WRITE(6,*) 'LAT. > 90.0 = 90.0, LAT < -90.0 = -90.0'
1730	ENDIF
1731	PTPOLN = .FALSE.
1732	NPOLNUM = 0
1733	DO 23006 I=1,LILJ
1734	IF( (EPSILON.GT. ABS(90.0 - XLAT(I))))THEN
1735	IF( (NPOLNUM.EQ. NPOLMAX))THEN
1736	CALL MEMOIRH(NPOLPTS, NPOLPT2, NPOLMAX+128, NPOLMAX)
1737	NPOLPTR = NPOLPT2
1738	NPOLMAX = NPOLMAX + 128
1739	ENDIF
1740	PTPOLN = .TRUE.
1741	NPOLPTS(NPOLPTR+NPOLNUM) = I
1742	NPOLNUM = NPOLNUM + 1
1743	ENDIF
1744	23006 CONTINUE
1745	PTPOLS = .FALSE.
1746	SPOLNUM = 0
1747	DO 23016 I=1,LILJ
1748	IF( (EPSILON.GT. ABS(-90.0 - XLAT(I))))THEN
1749	IF( (SPOLNUM.EQ. SPOLMAX))THEN
1750	CALL MEMOIRH(SPOLPTS, SPOLPT2, SPOLMAX+128, SPOLMAX)
1751	SPOLPTR = SPOLPT2
1752	SPOLMAX = SPOLMAX + 128
1753	ENDIF
1754	PTPOLS = .TRUE.
1755	SPOLPTS(SPOLPTR+SPOLNUM) = I
1756	SPOLNUM = SPOLNUM + 1
1757	ENDIF
1758	23016 CONTINUE
1759	RETURN
1760	END
1761	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1762	C
1763	LOGICAL FUNCTION QQQCMP(XLAT, XLON, LILJ, NI, NJ,GRTYP, IG1,
1764	$ IG2, IG3, IG4, XORSUM)
1765	USE mod_kinds_oasis
1766	INTEGER (kind=ip_intwp_p) LILJ
1767	REAL(kind=ip_realwp_p) XLAT(LILJ), XLON(LILJ)
1768	INTEGER (kind=ip_intwp_p) NI, NJ, IG1, IG2, IG3, IG4, XORSUM
1769	CHARACTER*1 GRTYP
1770	INTEGER (kind=ip_intwp_p) XORCALC
1771	CHARACTER*4 GRTMP
1772	C-- COMMON QQQCOM1
1773	C--
1774	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
1775	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
1776	$ AYTMP,CXTMP,CYTMP
1777	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
1778	$ SPOLPTS
1779	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
1780	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
1781	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
1782	$ NPOLMAX,SPOLMAX
1783	LOGICAL PTPOLN, PTPOLS
1784	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
1785	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
1786	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
1787	C----------------------
1788	C-- COMMON QQQCOM2
1789	C-
1790	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
1791	INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
1792	INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
1793	COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
1794	$,LSTIG3,LSTIG4,LSTXOR
1795	C-
1796	C------------
1797	C* COMMON QQQCOM3
1798	C
1799
1800	CHARACTER*1 LSTGTYP
1801	COMMON /QQQCOM3/ LSTGTYP
1802	C------------
1803
1804	QQQCMP = .TRUE.
1805	IF( (LILJ.NE.LSTLILJ.OR.NI.NE.LSTNI.OR.NJ.NE.LSTNJ))THEN
1806	QQQCMP = .FALSE.
1807	ENDIF
1808	GRTMP = GRTYP
1809	IF( (IG1.NE.LSTIG1.OR.IG2.NE.LSTIG2.OR.IG3.NE.LSTIG3.OR.IG4
1810	$.NE.LSTIG4.OR.GRTMP.NE.LSTGTYP))THEN
1811	QQQCMP = .FALSE.
1812	ENDIF
1813	XORSUM = XORCALC(XLAT, XLON, LILJ)
1814	IF( (XORSUM.NE.LSTXOR))THEN
1815	QQQCMP = .FALSE.
1816	ENDIF
1817	RETURN
1818	END
1819	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1820	C
1821	BLOCK DATA GAUSSGDBLK
1822	USE mod_kinds_oasis
1823	IMPLICIT NONE
1824	C
1825	C-- COMMON GAUSSGD
1826	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
1827	INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
1828	COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
1829	C----------------------
1830	DATA IROOTS, ILROOTS /0, 0/
1831	END
1832	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1833	C
1834	SUBROUTINE QQQGLAT(NJ,HEM)
1835	C* ----------------------------------------------------------------------
1836	C* S/R QQQGLAT - CALCUL DES LATITUDES D'UNE GRILLE GAUSSIENNE
1837	C*
1838	C AUTEUR: YVES CHARTIER. MARS 1991.
1839	C* ----------------------------------------------------------------------
1840	USE mod_kinds_oasis
1841	USE memoir
1842	IMPLICIT NONE
1843	INTEGER (kind=ip_intwp_p) NJ, HEM
1844	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
1845	PARAMETER (GLOBAL = 0)
1846	PARAMETER (NORD = 1)
1847	PARAMETER (SUD = 2)
1848	PARAMETER (SUDNORD= 0)
1849	PARAMETER (NORDSUD= 1)
1850	C
1851	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
1852	C
1853	C RDTODG = 180/PIE, DGTORD = PIE/180
1854
1855	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
1856	DATA PIE /3.14159265358979323846/
1857	DATA RDTODG /57.295779513082/
1858	DATA DGTORD /1.7453292519943E-2/
1859	C
1860	C-- COMMON GAUSSGD
1861	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
1862	INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
1863	COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
1864	C----------------------
1865	EXTERNAL DGAUSS
1866	EXTERNAL GAUSSGDBLK
1867	INTEGER (kind=ip_intwp_p) J,NPOLY
1868	INTEGER (kind=ip_intwp_p) NJ2
1869	IF( (IROOTS.NE. 0))THEN
1870	CALL MEMOIRH(ROOTS, IROOTS, 0, 0)
1871	ENDIF
1872	IF( (ILROOTS.NE. 0))THEN
1873	CALL MEMOIRH(LROOTS,ILROOTS, 0, 0)
1874	ENDIF
1875	IF( (HEM.NE. GLOBAL))THEN
1876	NPOLY = NJ * 2
1877	ELSE
1878	NPOLY = NJ
1879	ENDIF
1880	CALL MEMOIRH(ROOTS, IROOTS, NPOLY, 0)
1881	CALL MEMOIRH(LROOTS,ILROOTS,NPOLY, 0)
1882	CALL DGAUSS(NPOLY, ROOTS(IROOTS), HEM)
1883	IF( (HEM.EQ. GLOBAL))THEN
1884	NJ2 = NJ / 2
1885	DO 23008 J=1,NJ2
1886	LROOTS(J-1+ILROOTS)=90.0-RDTODG*ACOS(ROOTS(IROOTS+NJ2-J)
1887	$ )
1888	23008 CONTINUE
1889	ENDIF
1890	IF( (HEM.EQ. NORD))THEN
1891	DO 23012 J=1,NJ
1892	LROOTS(J-1+ILROOTS)=90.0-RDTODG*ACOS(ROOTS(IROOTS+NJ-J))
1893	23012 CONTINUE
1894	ENDIF
1895	IF( (HEM.EQ. SUD))THEN
1896	DO 23016 J=1,NJ
1897	LROOTS(ILROOTS-1+J)=-(90.0-RDTODG*ACOS(ROOTS(IROOTS-1+J)
1898	$ ))
1899	23016 CONTINUE
1900	ENDIF
1901	RETURN
1902	END
1903	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1904	C
1905	SUBROUTINE QQQKEEP(LILJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM)
1906	USE mod_kinds_oasis
1907	IMPLICIT NONE
1908	INTEGER (kind=ip_intwp_p) LILJ, NI, NJ, IG1, IG2, IG3, IG4, XORSUM
1909	CHARACTER*1 GRTYP
1910	C-- COMMON QQQCOM2
1911	C-
1912
1913	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
1914	INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
1915	INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
1916	COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
1917	$,LSTIG3,LSTIG4,LSTXOR
1918	C-
1919	C------------
1920	C* COMMON QQQCOM3
1921	C
1922
1923	CHARACTER*1 LSTGTYP
1924	COMMON /QQQCOM3/ LSTGTYP
1925	C------------
1926
1927	LSTLILJ = LILJ
1928	LSTNI = NI
1929	LSTNJ = NJ
1930	LSTGTYP = GRTYP
1931	LSTIG1 = IG1
1932	LSTIG2 = IG2
1933	LSTIG3 = IG3
1934	LSTIG4 = IG4
1935	LSTXOR = XORSUM
1936	RETURN
1937	END
1938	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1939	C
1940	SUBROUTINE QQQXTRAP(ZO, PX, PY, NPTS, Z, I1, I2, J1, J2)
1941	USE mod_kinds_oasis
1942	IMPLICIT NONE
1943	INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
1944	REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
1945	REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
1946	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
1947	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
1948	PARAMETER (VOISIN = 0)
1949	PARAMETER (LINEAIR = 1)
1950	PARAMETER (CUBIQUE = 3)
1951	PARAMETER (OUI = 1)
1952	PARAMETER (MINIMUM = 2)
1953	PARAMETER (MAXIMUM = 3)
1954	PARAMETER (VALEUR = 4)
1955	PARAMETER (ABORT = 13)
1956	LOGICAL FLGXTRAP
1957	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
1958	REAL(kind=ip_realwp_p) VALXTRAP
1959	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
1960	INTEGER (kind=ip_intwp_p) N, I, J, OFFL, OFFR
1961	REAL(kind=ip_realwp_p) RMIN, RMAX, TEMPMIN, TEMPMAX
1962	IF( (.NOT.FLGXTRAP))THEN
1963	RETURN
1964	ENDIF
1965	IF( (CODXTRAP.EQ. OUI))THEN
1966	RETURN
1967	ENDIF
1968	IF( (ORDINT.EQ. 3))THEN
1969	OFFR = 2
1970	OFFL = 1
1971	ELSE
1972	OFFR = 0
1973	OFFL = 0
1974	ENDIF
1975	RMIN = Z(I1, J1)
1976	RMAX = Z(I1, J1)
1977	DO 23006 J=J1, J2
1978	DO 23008 I=I1, I2
1979	IF( (Z(I,J).LT. RMIN))THEN
1980	RMIN = Z(I,J)
1981	ENDIF
1982	IF( (Z(I,J).GT. RMAX))THEN
1983	RMAX = Z(I,J)
1984	ENDIF
1985	23008 CONTINUE
1986	23006 CONTINUE
1987	TEMPMIN = RMIN - 0.25*(RMAX - RMIN)
1988	TEMPMAX = RMAX + 0.25*(RMAX - RMIN)
1989	RMIN = TEMPMIN
1990	RMAX = TEMPMAX
1991	DO 23014 N=1, NPTS
1992	I = INT(PX(N))
1993	J = INT(PY(N))
1994	IF( (I.LT.(I1+OFFL).OR. J.LT.(J1+OFFL).OR. I.GT. (I2-OFFR)
1995	$ .OR. J.GT. (J2-OFFR)))THEN
1996	IF( (CODXTRAP.EQ. VOISIN))THEN
1997	I = MIN(I2, MAX(I1, NINT(PX(N))))
1998	J = MIN(J2, MAX(J1, NINT(PY(N))))
1999	ZO(N) = Z(I,J)
2000	ENDIF
2001	IF( (CODXTRAP.EQ. MINIMUM))THEN
2002	ZO(N) = RMIN
2003	ENDIF
2004	IF( (CODXTRAP.EQ. MAXIMUM))THEN
2005	ZO(N) = RMAX
2006	ENDIF
2007	IF( (CODXTRAP.EQ. VALEUR))THEN
2008	ZO(N) = VALXTRAP
2009	ENDIF
2010	ENDIF
2011	23014 CONTINUE
2012	RETURN
2013	END
2014	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2015	C
2016	SUBROUTINE RGDINT(ZO,PX,PY,NPTS,Z,I1,I2,J1,J2)
2017	C*******
2018	C AUTEUR: Y.CHARTIER, DRPN
2019	C FEVRIER 1991
2020	C
2021	C OBJET: INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR D'UNE GRILLE
2022	C SOURCE REGULIERE.
2023	C
2024	C*******
2025	USE mod_kinds_oasis
2026	IMPLICIT NONE
2027	INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
2028	REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
2029	REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
2030	C
2031	C NPTS : NOMBRE DE POINTS A INTERPOLER
2032	C I1:I2 : DIMENSION DE LA GRILLE SOURCE SELON X
2033	C J1:J2 : DIMENSION DE LA GRILLE SOURCE SELON Y
2034	C ZO : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
2035	C PX : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
2036	C : VEUT INTERPOLER
2037	C PY : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
2038	C : VEUT INTERPOLER
2039	C Z : VALEURS DE LA GRILLE SOURCE.
2040	C
2041	C===========================================
2042	C
2043	C * * * *
2044	C
2045	C * * * *
2046	C # ==> PT (X,Y)
2047	C * (=) * * ==> = PT (IIND, JIND)
2048	C
2049	C * * * *
2050	C
2051	C===========================================
2052
2053	REAL(kind=ip_realwp_p) Y1,Y2,Y3,Y4
2054	INTEGER (kind=ip_intwp_p) M,N,I,J,STRIDE
2055	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
2056	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
2057	PARAMETER (VOISIN = 0)
2058	PARAMETER (LINEAIR = 1)
2059	PARAMETER (CUBIQUE = 3)
2060	PARAMETER (OUI = 1)
2061	PARAMETER (MINIMUM = 2)
2062	PARAMETER (MAXIMUM = 3)
2063	PARAMETER (VALEUR = 4)
2064	PARAMETER (ABORT = 13)
2065	LOGICAL FLGXTRAP
2066	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
2067	REAL(kind=ip_realwp_p) VALXTRAP
2068	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
2069	REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
2070	REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
2071	CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)0.1666666666666 + 0.5(Z2-Z3)
2072	$)DX + 0.5(Z1+Z3)-Z2)DX + Z3-0.1666666666666Z4-0.5*Z2-0.
2073	$3333333333333Z1)DX+Z2
2074	ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
2075	STRIDE = 1
2076	IF( (ORDINT.EQ. CUBIQUE))THEN
2077	DO 23002 N=1,NPTS
2078	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
2079	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
2080	DX = PX(N) - I
2081	DY = PY(N) - J
2082	Y1=CUBIC(Z(I-1,J-1),Z(I ,J-1),Z(I+1,J-1),Z(I+2,J-1),DX)
2083	Y2=CUBIC(Z(I-1,J ),Z(I ,J ),Z(I+1,J ),Z(I+2,J ),DX)
2084	Y3=CUBIC(Z(I-1,J+1),Z(I ,J+1),Z(I+1,J+1),Z(I+2,J+1),DX)
2085	Y4=CUBIC(Z(I-1,J+2),Z(I ,J+2),Z(I+1,J+2),Z(I+2,J+2),DX)
2086	ZO(N)=CUBIC(Y1,Y2,Y3,Y4,DY)
2087	23002 CONTINUE
2088	ENDIF
2089	IF( (ORDINT.EQ. LINEAIR))THEN
2090	DO 23006 N=1,NPTS
2091	I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
2092	J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
2093	DX = PX(N) - I
2094	DY = PY(N) - J
2095	Y2=ZLIN(Z(I,J ),Z(I+1,J ),DX)
2096	Y3=ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
2097	ZO(N)=ZLIN(Y2,Y3,DY)
2098	23006 CONTINUE
2099	ENDIF
2100	IF( (ORDINT.EQ. VOISIN))THEN
2101	DO 23010 N=1,NPTS
2102	I = MIN(I2,MAX(I1,NINT(PX(N))))
2103	J = MIN(J2,MAX(J1,NINT(PY(N))))
2104	ZO(N)=Z(I,J)
2105	23010 CONTINUE
2106	ENDIF
2107	RETURN
2108	END
2109	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2110	C
2111	SUBROUTINE RGFREE()
2112	USE mod_kinds_oasis
2113	USE memoir
2114	IMPLICIT NONE
2115	C
2116	C* - COMMON QQQCOM1
2117	C
2118	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
2119	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
2120	$ AYTMP,CXTMP,CYTMP
2121	INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS,
2122	$ SPOLPTS
2123	INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
2124	INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
2125	INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
2126	$ NPOLMAX,SPOLMAX
2127	LOGICAL PTPOLN, PTPOLS
2128	COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
2129	$YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
2130	$ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
2131	C
2132	C* - COMMON QQQCOM2
2133	C
2134
2135	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
2136	INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
2137	INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
2138	COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
2139	$,LSTIG3,LSTIG4,LSTXOR
2140	C
2141	C* - COMMON QQQCOM3
2142	C
2143
2144	CHARACTER*1 LSTGTYP
2145	COMMON /QQQCOM3/ LSTGTYP
2146	C* -----------------------------------------------------------------------
2147	IF( (XGDPTR.NE.0))THEN
2148	CALL MEMOIRH(XGD, XGDPTR, 0, 0)
2149	XGDPTR = 0
2150	ENDIF
2151	IF( (YGDPTR.NE.0))THEN
2152	CALL MEMOIRH(YGD, YGDPTR, 0, 0)
2153	YGDPTR = 0
2154	ENDIF
2155	IF( (ZPTR.NE.0))THEN
2156	CALL MEMOIRH(ZTMP,ZPTR, 0, 0)
2157	ZPTR = 0
2158	ENDIF
2159	IF( (AXPTR.NE.0))THEN
2160	CALL MEMOIRH(AXTMP,AXPTR, 0, 0)
2161	AXPTR = 0
2162	ENDIF
2163	IF( (AYPTR.NE.0))THEN
2164	CALL MEMOIRH(AYTMP,AYPTR, 0, 0)
2165	AYPTR = 0
2166	ENDIF
2167	IF( (CXPTR.NE.0))THEN
2168	CALL MEMOIRH(CXTMP, CXPTR, 0, 0)
2169	CXPTR = 0
2170	ENDIF
2171	IF( (CYPTR.NE.0))THEN
2172	CALL MEMOIRH(CYTMP, CYPTR, 0, 0)
2173	CYPTR = 0
2174	ENDIF
2175	IF( (NPOLPTR.NE.0))THEN
2176	CALL MEMOIRH(NPOLPTS, NPOLPTR, 0, 0)
2177	NPOLPTR = 0
2178	NPOLMAX = 0
2179	ENDIF
2180	IF( (SPOLPTR.NE.0))THEN
2181	CALL MEMOIRH(SPOLPTS, SPOLPTR, 0, 0)
2182	SPOLPTR = 0
2183	SPOLMAX = 0
2184	ENDIF
2185	LSTLILJ = 0
2186	LSTNI = 0
2187	LSTNJ = 0
2188	LSTGTYP = ' '
2189	LSTIG1 = 0
2190	LSTIG2 = 0
2191	LSTIG3 = 0
2192	LSTIG4 = 0
2193	LSTXOR = 0
2194	RETURN
2195	END
2196	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2197	C
2198	SUBROUTINE RGLL2GD(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
2199	USE mod_kinds_oasis
2200	IMPLICIT NONE
2201	INTEGER (kind=ip_intwp_p) LI,LJ
2202	REAL(kind=ip_realwp_p) SPDO(LI,LJ), PSIO(LI,LJ), XLON(LI,LJ)
2203	CHARACTER*1 GRTYP
2204	INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
2205	EXTERNAL CIGAXG
2206	C
2207	C AUTEUR - Y. CHARTIER - AVRIL 91
2208	C
2209	C OBJET(RGLL2GD)
2210	C - PASSE DE PSIOTESSE ET DIRECTION
2211	C - A VENT DE GRILLE (COMPOSANTES U ET V)
2212	C
2213	C APPEL - CALL RGLL2GD(SPD,PSI,LI,LJ,IYP,XG1,XG2,XG3,XG4)
2214	C
2215	C MODULES - XGAIG
2216	C
2217	C ARGUMENTS
2218	C IN/OUT - SPD - A L'ENTREE CONTIENT LA PSIOTESSE DU VENT ET
2219	C A LA SORTIE LA COMPOSANTE U.
2220	C IN/OUT - PSI - A L'ENTREE CONTIENT LA DIRECTION DU VENT ET
2221	C A LA SORTIE LA COMPOSANTE V.
2222	C IN - LI - PREMIERE DIMENSION DES CHAMPS SPD ET PSI
2223	C IN - LJ - DEUXIEME DIMENSION DES CHAMPS SPD ET PSI
2224	C IN - IGTYP - TYPE DE GRILLE (VOIR OUVRIR)
2225	C IN - XG1 - ** DESCRIPTEUR DE GRILLE (REEL),
2226	C IN - XG2 - IGTYP = 'N', PI, PJ, D60, DGRW
2227	C IN - XG3 - IGTYP = 'L', LAT0, LON0, DLAT, DLON,
2228	C IN - XG4 - IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
2229	C = 1. NORD
2230	C = 2. SUD **
2231	C
2232	CMESSAGES - "ERREUR MAUVAISE GRILLE (RGLL2GD)"
2233	C
2234	C-------------------------------------------------------------
2235	C
2236	C
2237	C
2238	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
2239	C
2240	C RDTODG = 180/PIE, DGTORD = PIE/180
2241
2242	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
2243	DATA PIE /3.14159265358979323846/
2244	DATA RDTODG /57.295779513082/
2245	DATA DGTORD /1.7453292519943E-2/
2246	C
2247
2248	INTEGER (kind=ip_intwp_p) I,J
2249	REAL(kind=ip_realwp_p) PSI,U,V
2250	REAL(kind=ip_realwp_p) XG1, XG2, XG3, XG4
2251	IF( (GRTYP.EQ. 'N'))THEN
2252	CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
2253	DO 23002 I=1,LI
2254	DO 23004 J=1,LJ
2255	PSI =XLON(I,J)+XG4-PSIO(I,J)
2256	U = COS(PSIDGTORD)SPDO(I,J)
2257	V = SIN(PSIDGTORD)SPDO(I,J)
2258	SPDO(I,J) = U
2259	PSIO(I,J) = V
2260	23004 CONTINUE
2261	23002 CONTINUE
2262	RETURN
2263	ENDIF
2264	IF( (GRTYP.EQ. 'S'))THEN
2265	CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
2266	DO 23008 I=1,LI
2267	DO 23010 J=1,LJ
2268	PSI =180.0 - XLON(I,J)+XG4-PSIO(I,J)
2269	U = COS(PSIDGTORD)SPDO(I,J)
2270	V = SIN(PSIDGTORD)SPDO(I,J)
2271	SPDO(I,J) = U
2272	PSIO(I,J) = V
2273	23010 CONTINUE
2274	23008 CONTINUE
2275	RETURN
2276	ENDIF
2277	IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'G'.OR.GRTYP.EQ.
2278	$'L'))THEN
2279	DO 23014 I=1,LI
2280	DO 23016 J=1,LJ
2281	PSI = 270.0 - PSIO(I,J)
2282	U = COS(PSIDGTORD)SPDO(I,J)
2283	V = SIN(PSIDGTORD)SPDO(I,J)
2284	SPDO(I,J) = U
2285	PSIO(I,J) = V
2286	23016 CONTINUE
2287	23014 CONTINUE
2288	RETURN
2289	ENDIF
2290	WRITE(6, 600) GRTYP
2291	600 FORMAT('0',' Error, bad grid (RGLL2GD) - GRTYP = ', A1)
2292	RETURN
2293	END
2294	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2295	C
2296	SUBROUTINE RGOPTC(OP, VAL, FLAG)
2297	USE mod_kinds_oasis
2298	IMPLICIT NONE
2299	CHARACTER () OP, VAL
2300	LOGICAL FLAG
2301	CFLAG = .TRUE. MODE SET
2302	CFLAG = .FALSE. MODE GET
2303
2304	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
2305	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
2306	PARAMETER (VOISIN = 0)
2307	PARAMETER (LINEAIR = 1)
2308	PARAMETER (CUBIQUE = 3)
2309	PARAMETER (OUI = 1)
2310	PARAMETER (MINIMUM = 2)
2311	PARAMETER (MAXIMUM = 3)
2312	PARAMETER (VALEUR = 4)
2313	PARAMETER (ABORT = 13)
2314	LOGICAL FLGXTRAP
2315	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
2316	REAL(kind=ip_realwp_p) VALXTRAP
2317	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
2318	IF( (FLAG))THEN
2319	IF( (OP.EQ. 'EXTRAP'))THEN
2320	IF( (VAL.EQ. 'OUI'))THEN
2321	CODXTRAP = OUI
2322	ELSE
2323	IF( (VAL.EQ. 'ABORT'))THEN
2324	CODXTRAP = ABORT
2325	ELSE
2326	IF( (VAL.EQ. 'MAXIMUM'))THEN
2327	CODXTRAP = MAXIMUM
2328	ELSE
2329	IF( (VAL.EQ. 'MINIMUM'))THEN
2330	CODXTRAP = MINIMUM
2331	ELSE
2332	IF( (VAL.EQ. 'VOISIN'))THEN
2333	CODXTRAP = VOISIN
2334	ELSE
2335	IF( (VAL.EQ. 'VALEUR'))THEN
2336	CODXTRAP = VALEUR
2337	ELSE
2338	WRITE(6,*)
2339	$ ' <RGOPTC>: bad value for VAL'
2340	WRITE(6,*) ' VAL = ', VAL
2341	WRITE(6,*)
2342	$ ' VAL initialized to ''ABORT'''
2343	CODXTRAP = ABORT
2344	ENDIF
2345	ENDIF
2346	ENDIF
2347	ENDIF
2348	ENDIF
2349	ENDIF
2350	ELSE
2351	IF( (OP.EQ. 'INTERP'.OR. OP.EQ. 'INTERP'))THEN
2352	IF( (VAL.EQ. 'VOISIN'))THEN
2353	ORDINT = VOISIN
2354	ELSE
2355	IF( (VAL.EQ. 'LINEAIR'))THEN
2356	ORDINT = LINEAIR
2357	ELSE
2358	IF( (VAL.EQ. 'CUBIQUE'))THEN
2359	ORDINT = CUBIQUE
2360	ELSE
2361	WRITE(6,*)
2362	$ ' <RGOPTC>: bad value for VAL'
2363	WRITE(6,*) ' VAL = ', VAL
2364	WRITE(6,*)
2365	$ ' VAL initialized to ''CUBIQUE'''
2366	ORDINT = CUBIQUE
2367	ENDIF
2368	ENDIF
2369	ENDIF
2370	ELSE
2371	WRITE(6,*) ' <RGOPTC>: bad value for OP'
2372	WRITE(6,*)
2373	$ ' OP should be equal to ''EXTRAP'' or ''INTERP'''
2374	ENDIF
2375	ENDIF
2376	ELSE
2377	IF( (OP.EQ. 'EXTRAP'))THEN
2378	IF( (CODXTRAP.EQ. OUI))THEN
2379	VAL = 'OUI'
2380	ELSE
2381	IF( (CODXTRAP.EQ. ABORT))THEN
2382	VAL = 'ABORT'
2383	ELSE
2384	IF( (CODXTRAP.EQ. MAXIMUM))THEN
2385	VAL = 'MAXIMUM'
2386	ELSE
2387	IF( (CODXTRAP.EQ. MINIMUM))THEN
2388	VAL = 'MINIMUM'
2389	ELSE
2390	IF( (CODXTRAP.EQ. VOISIN))THEN
2391	VAL = 'VOISIN'
2392	ELSE
2393	IF( (CODXTRAP.EQ. VALEUR))THEN
2394	VAL = 'VALEUR'
2395	ELSE
2396	WRITE(6,*)
2397	$ ' <RGOPTC>: bad value for CODXTRAP'
2398	WRITE(6,*) ' CODXTRAP = ',
2399	$ CODXTRAP
2400	VAL = 'SCRAP'
2401	ENDIF
2402	ENDIF
2403	ENDIF
2404	ENDIF
2405	ENDIF
2406	ENDIF
2407	ELSE
2408	IF( (OP.EQ. 'INTERP'))THEN
2409	IF( (ORDINT.EQ. VOISIN))THEN
2410	VAL = 'VOISIN'
2411	ELSE
2412	IF( (ORDINT.EQ. LINEAIR))THEN
2413	VAL = 'LINEAIR'
2414	ELSE
2415	IF( (ORDINT.EQ. CUBIQUE))THEN
2416	VAL = 'CUBIQUE'
2417	ELSE
2418	WRITE(6,*)
2419	$ ' <RGOPTC>: bad value for ORDINT'
2420	WRITE(6,*) ' ORDINT = ', ORDINT
2421	VAL = 'SCRAP'
2422	ENDIF
2423	ENDIF
2424	ENDIF
2425	ELSE
2426	WRITE(6,*) ' <RGOPTC>: bad value for OP'
2427	WRITE(6,*)
2428	$ ' OP should be equal to ''EXTRAP'' OU ''INTERP'''
2429	ENDIF
2430	ENDIF
2431	ENDIF
2432	RETURN
2433	END
2434	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2435	C
2436	BLOCK DATA QQQXTRPBLK
2437	USE mod_kinds_oasis
2438	IMPLICIT NONE
2439	INTEGER (kind=ip_intwp_p) OUI
2440	PARAMETER (OUI = 1)
2441	LOGICAL FLGXTRAP
2442	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
2443	REAL(kind=ip_realwp_p) VALXTRAP
2444	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
2445	DATA CODXTRAP / OUI /
2446	DATA FLGXTRAP / .FALSE. /
2447	DATA ORDINT / 3 /
2448	END
2449	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2450	C
2451	SUBROUTINE RGOPTI(OP, VAL, FLAG)
2452	USE mod_kinds_oasis
2453	IMPLICIT NONE
2454	CHARACTER () OP
2455	INTEGER (kind=ip_intwp_p) VAL
2456	LOGICAL FLAG
2457	CFLAG = .TRUE. MODE SET
2458	CFLAG = .FALSE. MODE GET
2459
2460	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
2461	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
2462	PARAMETER (VOISIN = 0)
2463	PARAMETER (LINEAIR = 1)
2464	PARAMETER (CUBIQUE = 3)
2465	PARAMETER (OUI = 1)
2466	PARAMETER (MINIMUM = 2)
2467	PARAMETER (MAXIMUM = 3)
2468	PARAMETER (VALEUR = 4)
2469	PARAMETER (ABORT = 13)
2470	LOGICAL FLGXTRAP
2471	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
2472	REAL(kind=ip_realwp_p) VALXTRAP
2473	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
2474	EXTERNAL QQQXTRPBLK
2475	IF( (FLAG))THEN
2476	IF( (OP.EQ. 'EXTRAP'))THEN
2477	VALXTRAP = REAL(VAL,ip_realwp_p)
2478	ELSE
2479	IF( (OP.EQ. 'INTERP'))THEN
2480	IF( (VAL.EQ. VOISIN.OR. VAL.EQ. LINEAIR.OR. VAL
2481	$ .EQ. CUBIQUE))THEN
2482	ORDINT = VAL
2483	ELSE
2484	WRITE(6,*) ' <RGOPTI>: bad value for VAL'
2485	WRITE(6,*) ' VAL = ', VAL
2486	WRITE(6,*)
2487	$ ' VAL initialized to ''CUBIQUE'''
2488	ORDINT = CUBIQUE
2489	ENDIF
2490	ELSE
2491	WRITE(6,*) ' <RGOPTI>: bad value for OP'
2492	WRITE(6,*)
2493	$ ' OP should be equal to ''EXTRAP'' or ''INTERP'''
2494	ENDIF
2495	ENDIF
2496	ELSE
2497	IF( (OP.EQ. 'EXTRAP'))THEN
2498	VAL = NINT(VALXTRAP)
2499	ELSE
2500	IF( (OP.EQ. 'INTERP'))THEN
2501	IF( (ORDINT.EQ. VOISIN.OR. ORDINT.EQ. LINEAIR.OR.
2502	$ ORDINT.EQ. CUBIQUE))THEN
2503	VAL = ORDINT
2504	ELSE
2505	WRITE(6,*) ' <RGOPTI>: bad value for ORDINT'
2506	WRITE(6,*) ' ORDINT = ', ORDINT
2507	VAL = -1
2508	ENDIF
2509	ELSE
2510	WRITE(6,*) ' <RGOPTI>: bad value for OP'
2511	WRITE(6,*)
2512	$ ' OP should be equal to ''EXTRAP'' or ''INTERP'''
2513	ENDIF
2514	ENDIF
2515	ENDIF
2516	RETURN
2517	END
2518	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2519	C
2520	SUBROUTINE RGOPTR(OP, VAL, FLAG)
2521	USE mod_kinds_oasis
2522	IMPLICIT NONE
2523	CHARACTER () OP
2524	REAL(kind=ip_realwp_p) VAL
2525	LOGICAL FLAG
2526	C FLAG = .TRUE. MODE SET
2527	C FLAG = .FALSE. MODE GET
2528
2529	INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
2530	INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
2531	PARAMETER (VOISIN = 0)
2532	PARAMETER (LINEAIR = 1)
2533	PARAMETER (CUBIQUE = 3)
2534	PARAMETER (OUI = 1)
2535	PARAMETER (MINIMUM = 2)
2536	PARAMETER (MAXIMUM = 3)
2537	PARAMETER (VALEUR = 4)
2538	PARAMETER (ABORT = 13)
2539	LOGICAL FLGXTRAP
2540	INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
2541	REAL(kind=ip_realwp_p) VALXTRAP
2542	COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
2543	IF( (FLAG))THEN
2544	IF( (OP.EQ. 'EXTRAP'))THEN
2545	VALXTRAP = VAL
2546	ELSE
2547	IF( (OP.EQ. 'INTERP'))THEN
2548	IF( (VAL.EQ. VOISIN.OR. VAL.EQ. LINEAIR.OR. VAL.EQ.
2549	$ CUBIQUE))THEN
2550	ORDINT = NINT(VAL)
2551	ELSE
2552	WRITE(6,*)' <RGOPTR>: bad value for VAL'
2553	WRITE(6,*) ' VAL = ', VAL
2554	WRITE(6,*)
2555	$ ' VAL initialized to ''CUBIQUE'''
2556	ORDINT = CUBIQUE
2557	ENDIF
2558	ELSE
2559	WRITE(6,*) ' <RGOPTR>: bad value for OP'
2560	WRITE(6,*)
2561	$ ' OP should be equal to ''EXTRAP'' or ''INTERP'''
2562	ENDIF
2563	ENDIF
2564	ELSE
2565	IF( (OP.EQ. 'EXTRAP'))THEN
2566	VAL = VALXTRAP
2567	ELSE
2568	IF( (OP.EQ. 'INTERP'))THEN
2569	IF( (ORDINT.EQ. VOISIN.OR. ORDINT.EQ. LINEAIR.OR.
2570	$ ORDINT.EQ. CUBIQUE))THEN
2571	VAL = REAL(ORDINT,ip_realwp_p)
2572	ELSE
2573	WRITE(6,*) ' <RGOPTR>: bad value for ORDINT'
2574	WRITE(6,*) ' ORDINT = ', ORDINT
2575	VAL = REAL(ORDINT,ip_realwp_p)
2576	ENDIF
2577	ELSE
2578	WRITE(6,*) ' <RGOPTR>: bad value for OP'
2579	WRITE(6,*)
2580	$ ' OP should be equal to ''EXTRAP'' OU ''INTERP'''
2581	ENDIF
2582	ENDIF
2583	ENDIF
2584	RETURN
2585	END
2586	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2587	C
2588	SUBROUTINE VLLFXY(DLAT,DLON,X,Y,NI,NJ,D60,DGRW,PI,PJ,NHEM)
2589	C* ----------------------------------------------------------------------
2590	C* S/R VLLFXY(I) - COMPUTES THE GRID CO-ORDINATES OF A POINT
2591	C* ----------------------------------------------------------------------
2592	USE mod_kinds_oasis
2593	IMPLICIT NONE
2594	C* ----------------------------------------------------------------------
2595	C* ARGUMENTS
2596	C OUT - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N).
2597	C - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E).
2598	C IN - X - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
2599	C AS ORIGIN
2600	C - Y - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
2601	C AS ORIGIN
2602	C - D60 - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
2603	C GRID AT 60 DEGREES
2604	C - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
2605	C THE GRID (IN DEGREES)
2606	C - NHEM - 1 FOR NORTHERN HEMISPHERE
2607	C 2 FOR SOUTHERN HEMISPHERE
2608	C
2609	CNOTES - THE COMPANION ROUTINE XYFLL, WHICH COMPUTES THE GRID
2610	C CO-ORDINATES GIVEN THE LATITUDE AND LONGITUDE, IS ALSO
2611	C AVAILABLE.
2612	C
2613	C-----------------------------------------------------------------------
2614	C
2615	C
2616	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
2617	C
2618	C RDTODG = 180/PIE, DGTORD = PIE/180
2619	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
2620	DATA PIE /3.14159265358979323846/
2621	DATA RDTODG /57.295779513082/
2622	DATA DGTORD /1.7453292519943E-2/
2623	C
2624
2625	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
2626	PARAMETER (GLOBAL = 0)
2627	PARAMETER (NORD = 1)
2628	PARAMETER (SUD = 2)
2629	PARAMETER (SUDNORD= 0)
2630	PARAMETER (NORDSUD= 1)
2631	INTEGER (kind=ip_intwp_p) NI,NJ, NHEM
2632	REAL(kind=ip_realwp_p) X(NI,NJ), Y(NI,NJ), DLAT(NI,NJ),
2633	$ DLON(NI,NJ)
2634	REAL(kind=ip_realwp_p) X1, Y1
2635	REAL(kind=ip_realwp_p) D60, DGRW, PI, PJ
2636	REAL(kind=ip_realwp_p) RE,RE2,R2,RLON,RLAT,SINLAT,R
2637	INTEGER (kind=ip_intwp_p) I,J
2638	RE=1.866025*6.371E+6/D60
2639	RE2=RE**2
2640	DO 23000 I=1,NI
2641	DO 23002 J=1,NJ
2642	X1 = X(I,J) - PI
2643	Y1 = Y(I,J) - PJ
2644	IF( (X1.EQ. 0..AND. Y1.EQ. 0.))THEN
2645	DLAT(I,J)=90.
2646	DLON(I,J)=0.
2647
2648	C
2649	C CALCULATE LONGITUDE IN MAP COORDINATES.
2650	C
2651	ENDIF
2652	IF((X1.EQ. 0.))THEN
2653	DLON(I,J)=SIGN(90._ip_realwp_p,Y1)
2654	ENDIF
2655	IF((X1.NE. 0.))THEN
2656	DLON(I,J)=ATAN(Y1/X1)*RDTODG
2657	ENDIF
2658	IF((X1.LT. 0.))THEN
2659	DLON(I,J)=DLON(I,J)+SIGN(180._ip_realwp_p,Y1)
2660	C
2661	C * ADJUST LONGITUDE FOR GRID ORIENTATION.
2662	C
2663
2664	ENDIF
2665	DLON(I,J)=DLON(I,J)-DGRW
2666	IF((DLON(I,J).GT. 180.))THEN
2667	DLON(I,J)=DLON(I,J)-360.
2668	ENDIF
2669	IF((DLON(I,J).LT. -180.))THEN
2670	DLON(I,J)=DLON(I,J)+360.
2671	C
2672	C * CALCULATE LATITUDE.
2673	C
2674
2675	ENDIF
2676	R2=X12+Y12
2677	DLAT(I,J)=(RE2-R2)/(RE2+R2)
2678	DLAT(I,J)= ASIN(DLAT(I,J))*RDTODG
2679	C
2680	C CHANGE SIGNS IF IN SOUTHERN HEMISPHERE.
2681	C
2682
2683	IF((NHEM.EQ. 2))THEN
2684	DLAT(I,J)=-DLAT(I,J)
2685	ENDIF
2686	IF((NHEM.EQ. 2))THEN
2687	DLON(I,J)=-DLON(I,J)
2688	ENDIF
2689	23002 CONTINUE
2690	23000 CONTINUE
2691	RETURN
2692	END
2693	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2694	C
2695	SUBROUTINE VXYFLL(X,Y,DLAT,DLON,NPTS,D60,DGRW,PI,PJ,NHEM)
2696	USE mod_kinds_oasis
2697	IMPLICIT NONE
2698
2699	C AUTHOR - J.D. HENDERSON - FEB 75
2700	C
2701	C OBJECT(VXYFLL)
2702	C - COMPUTES THE GRID CO-ORDINATES MEASURED FROM THE POLE OF A
2703	C POINT, GIVEN THE LATITUDE AND LONGITUDE IN DEGREES.
2704	C
2705	C USAGE - CALL VXYFLL(X,Y,DLAT,DLON,NPTS,D60,DGRW,PI,PJ,NHEM)
2706	C
2707	C ARGUMENTS
2708	C OUT - X - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
2709	C AS ORIGIN
2710	C - Y - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
2711	C AS ORIGIN
2712	C IN - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N)
2713	C - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E)
2714	C - D60 - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
2715	C GRID AT 60 DEGREES
2716	C - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
2717	C THE GRID (IN DEGREES)
2718	C - NHEM - 1 FOR NORTHERN HEMISPHERE
2719	C 2 FOR SOUTHERN HEMISPHERE
2720	C
2721	C NOTES - THE COMPANION ROUTINE LLFXY, WHICH COMPUTES THE LATITUDE
2722	C - AND LONGITUDE GIVEN THE GRID-COORDINATES,
2723	C - IS ALSO AVAILABLE.
2724	C*--------------------------------------------------------------------
2725	C
2726	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
2727	C
2728	C RDTODG = 180/PIE, DGTORD = PIE/180
2729	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
2730	DATA PIE /3.14159265358979323846/
2731	DATA RDTODG /57.295779513082/
2732	DATA DGTORD /1.7453292519943E-2/
2733	C
2734
2735	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
2736	PARAMETER (GLOBAL = 0)
2737	PARAMETER (NORD = 1)
2738	PARAMETER (SUD = 2)
2739	PARAMETER (SUDNORD= 0)
2740	PARAMETER (NORDSUD= 1)
2741	INTEGER (kind=ip_intwp_p) NPTS, NHEM
2742	REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), DLAT(NPTS), DLON(NPTS)
2743	REAL(kind=ip_realwp_p) D60, DGRW, PI, PJ
2744	REAL(kind=ip_realwp_p) RE,RLON,RLAT,SINLAT,R
2745	INTEGER (kind=ip_intwp_p) I
2746	RE=1.866025*6.371E+6/D60
2747	C
2748
2749	IF( (NHEM.EQ.NORD))THEN
2750	DO 23002 I=1,NPTS
2751	RLON=DGTORD*(DLON(I)+DGRW)
2752	RLAT=DGTORD*DLAT(I)
2753	SINLAT=SIN(RLAT)
2754	R=RE*SQRT((1.-SINLAT)/(1.+SINLAT))
2755	X(I)=R*COS(RLON) + PI
2756	Y(I)=R*SIN(RLON) + PJ
2757	23002 CONTINUE
2758	RETURN
2759	ENDIF
2760	IF( (NHEM.EQ.SUD))THEN
2761	DO 23006 I=1,NPTS
2762	RLON = DLON(I)
2763	IF( (RLON.GT. 180.0))THEN
2764	RLON = RLON - 360.0
2765	ENDIF
2766	RLON=DGTORD*(-RLON+DGRW)
2767	RLAT=DGTORD*(-DLAT(I))
2768	SINLAT=SIN(RLAT)
2769	R=RE*SQRT((1.-SINLAT)/(1.+SINLAT))
2770	X(I)=R*COS(RLON)+PI
2771	Y(I)=R*SIN(RLON)+PJ
2772	23006 CONTINUE
2773	ENDIF
2774	RETURN
2775	END
2776	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2777	C
2778	FUNCTION XORCALC(ILAT, ILON, LILJ)
2779	USE mod_kinds_oasis
2780	INTEGER (kind=ip_intwp_p) XORCALC
2781	INTEGER (kind=ip_intwp_p) LILJ
2782	INTEGER (kind=ip_intwp_p) ILAT(LILJ), ILON(LILJ)
2783	XORCALC = 0
2784	DO 23000 I=1,LILJ,17
2785	XORCALC = IEOR(IEOR(XORCALC,ILAT(I)),ILON(I))
2786	23000 CONTINUE
2787	RETURN
2788	END
2789	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2790	C
2791	SUBROUTINE XPNAXEG(NEWAXEX,NEWAXEY,I1,I2,J1,J2,NI,NJ,HEM)
2792	USE mod_kinds_oasis
2793	IMPLICIT NONE
2794	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ,HEM
2795	REAL(kind=ip_realwp_p) NEWAXEX(I1:I2)
2796	REAL(kind=ip_realwp_p) NEWAXEY(J1:J2)
2797	C
2798	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
2799	C
2800	C RDTODG = 180/PIE, DGTORD = PIE/180
2801
2802	REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
2803	DATA PIE /3.14159265358979323846/
2804	DATA RDTODG /57.295779513082/
2805	DATA DGTORD /1.7453292519943E-2/
2806	C
2807	C-- COMMON GAUSSGD
2808	REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
2809	INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
2810	COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
2811	C----------------------
2812
2813	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
2814	PARAMETER (GLOBAL = 0)
2815	PARAMETER (NORD = 1)
2816	PARAMETER (SUD = 2)
2817	PARAMETER (SUDNORD= 0)
2818	PARAMETER (NORDSUD= 1)
2819	INTEGER (kind=ip_intwp_p) I,J
2820	DO 23000 I=1,NI
2821	NEWAXEX(I) = REAL(I,ip_realwp_p)
2822	23000 CONTINUE
2823	NEWAXEX(0) = NEWAXEX(NI) - REAL(NI,ip_realwp_p)
2824	NEWAXEX(NI+1) = NEWAXEX(1) + REAL(NI,ip_realwp_p)
2825	NEWAXEX(NI+2) = NEWAXEX(2) + REAL(NI,ip_realwp_p)
2826	IF( (HEM.EQ. NORD))THEN
2827	DO 23004 J=1,NJ
2828	NEWAXEY(J) = LROOTS(ILROOTS-1+J)
2829	23004 CONTINUE
2830	NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
2831	NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
2832	DO 23006 J=-NJ+1,0
2833	NEWAXEY(J) = -(NEWAXEY(-J+1))
2834	23006 CONTINUE
2835	NEWAXEY(-NJ ) = -180.0 - NEWAXEY(-NJ+1)
2836	NEWAXEY(-NJ-1) = -180.0 - NEWAXEY(-NJ+2)
2837	ENDIF
2838	IF( (HEM.EQ. SUD))THEN
2839	DO 23010 J=1,NJ
2840	NEWAXEY(J) = -(LROOTS(ILROOTS-1+NJ-J+1))
2841	NEWAXEY(NJ+J) = LROOTS(ILROOTS-1+NJ-J+1)
2842	23010 CONTINUE
2843	NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
2844	NEWAXEY(-1) = -180.0 - NEWAXEY(2)
2845	NEWAXEY(NJ+1) = 180.0 - NEWAXEY(2*NJ)
2846	NEWAXEY(NJ+2) = 180.0 - NEWAXEY(2*NJ-1)
2847	ENDIF
2848	IF( (HEM.EQ. GLOBAL))THEN
2849	IF( (0.EQ. MOD(NJ,2)))THEN
2850	DO 23016 J=1,NJ/2
2851	NEWAXEY(J) = -(LROOTS(ILROOTS-1+NJ/2-J+1))
2852	NEWAXEY(NJ/2+J) = LROOTS(ILROOTS-1+J)
2853	23016 CONTINUE
2854	ELSE
2855	DO 23018 J=1,NJ/2
2856	NEWAXEY(J) = -(LROOTS(ILROOTS-1+NJ/2-J+1))
2857	NEWAXEY(NJ/2+J+1) = LROOTS(ILROOTS-1+J)
2858	23018 CONTINUE
2859	NEWAXEY(NJ/2+1) = 0.0
2860	ENDIF
2861	NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
2862	NEWAXEY(-1) = -180.0 - NEWAXEY(2)
2863	NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
2864	NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
2865	ENDIF
2866	RETURN
2867	END
2868	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2869	C
2870	SUBROUTINE XPNAXEZ(NEWAXEX,NEWAXEY,I1,I2,J1,J2,AX,AY,NI,NJ)
2871	USE mod_kinds_oasis
2872	IMPLICIT NONE
2873	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
2874	REAL(kind=ip_realwp_p) NEWAXEX(I1:I2)
2875	REAL(kind=ip_realwp_p) NEWAXEY(J1:J2)
2876	REAL(kind=ip_realwp_p) AX(NI)
2877	REAL(kind=ip_realwp_p) AY(NJ)
2878	INTEGER (kind=ip_intwp_p) I,J
2879	DO 23000 I=1,NI
2880	NEWAXEX(I) = AX(I)
2881	23000 CONTINUE
2882	NEWAXEX(0) = NEWAXEX(NI) - AX(NI)
2883	NEWAXEX(NI+1) = NEWAXEX(1) + AX(NI)
2884	NEWAXEX(NI+2) = NEWAXEX(2) + AX(NI)
2885	DO 23002 J=1,NJ
2886	NEWAXEY(J) = AY(J)
2887	23002 CONTINUE
2888	NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
2889	NEWAXEY(-1) = -180.0 - NEWAXEY(2)
2890	NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
2891	NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
2892	cLT
2893	WRITE(,) NEWAXEY(0 ), NEWAXEY(-1), NEWAXEY(NJ+1),
2894	$ NEWAXEY(NJ+2)
2895	RETURN
2896	END
2897	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2898	C
2899	SUBROUTINE XPNAXEY(NEWAXEX,NEWAXEY,I1,I2,J1,J2,AX,AY,NI,NJ)
2900	USE mod_kinds_oasis
2901	IMPLICIT NONE
2902	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
2903	REAL(kind=ip_realwp_p) NEWAXEX(NI)
2904	REAL(kind=ip_realwp_p) NEWAXEY(NJ)
2905	REAL(kind=ip_realwp_p) AX(NI)
2906	REAL(kind=ip_realwp_p) AY(NJ)
2907	INTEGER (kind=ip_intwp_p) I,J
2908	DO 23000 I=1,NI
2909	NEWAXEX(I) = AX(I)
2910	23000 CONTINUE
2911	DO 23002 J=1,NJ
2912	NEWAXEY(J) = AY(J)
2913	23002 CONTINUE
2914	RETURN
2915	END
2916	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2917	C
2918	SUBROUTINE XPNCOF(I1,I2,J1,J2,NI,NJ,GRTYP,GRREF,IG1,IG2,IG3,
2919	$IG4,SYM,AX,AY)
2920	USE mod_kinds_oasis
2921	IMPLICIT NONE
2922	INTEGER (kind=ip_intwp_p) I1,J1,I2,J2,NI,NJ,IG1,IG2,IG3,IG4
2923	LOGICAL SYM
2924	CHARACTER*1 GRTYP, GRREF
2925	REAL(kind=ip_realwp_p) AX(NI),AY(NJ)
2926	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
2927	PARAMETER (GLOBAL = 0)
2928	PARAMETER (NORD = 1)
2929	PARAMETER (SUD = 2)
2930	PARAMETER (SUDNORD= 0)
2931	PARAMETER (NORDSUD= 1)
2932	IF( (GRTYP.EQ.'L'.OR.GRTYP.EQ.'N'.OR.GRTYP.EQ.'S'.OR.GRTYP.EQ.
2933	$'Y'))THEN
2934	I1 = 1
2935	I2 = NI
2936	J1 = 1
2937	J2 = NJ
2938	RETURN
2939	ENDIF
2940	IF ( (GRTYP .EQ. 'Z')) THEN
2941	I1 = 0
2942	I2 = NI + 2
2943	J1 = -1
2944	J2 = NJ + 2
2945	ENDIF
2946	IF( (GRTYP.EQ.'A'.OR. GRTYP.EQ.'G'))THEN
2947	I1 = 0
2948	I2 = NI + 2
2949	IF( (IG1.EQ.GLOBAL))THEN
2950	J1 = -1
2951	J2 = NJ + 2
2952	ELSE
2953	IF( (IG1.EQ.NORD))THEN
2954	J1 = -NJ - 1
2955	J2 = NJ + 2
2956	ELSE
2957	J1 = -1
2958	J2 = 2 * NJ + 2
2959	ENDIF
2960	ENDIF
2961	RETURN
2962	ENDIF
2963	IF( (GRTYP.EQ.'B'))THEN
2964	I1 = 0
2965	I2 = NI + 1
2966	IF( (IG1.EQ.GLOBAL))THEN
2967	J1 = 0
2968	J2 = NJ + 1
2969	ELSE
2970	IF( (IG1.EQ.NORD))THEN
2971	J1 = -NJ + 1
2972	J2 = NJ + 1
2973	ELSE
2974	J1 = 0
2975	J2 = 2 * NJ
2976	ENDIF
2977	ENDIF
2978	ENDIF
2979	RETURN
2980	END
2981	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2982	C
2983	SUBROUTINE XPNGD(ZOUT,I1,I2,J1,J2,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
2984	$IG4,SYM,VECT)
2985	USE mod_kinds_oasis
2986	IMPLICIT NONE
2987	INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
2988	PARAMETER (GLOBAL = 0)
2989	PARAMETER (NORD = 1)
2990	PARAMETER (SUD = 2)
2991	PARAMETER (SUDNORD= 0)
2992	PARAMETER (NORDSUD= 1)
2993	EXTERNAL PERMUT
2994	INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
2995	INTEGER (kind=ip_intwp_p) IG1, IG2, IG3, IG4
2996	CHARACTER*1 GRTYP
2997	INTEGER (kind=ip_intwp_p) I,J
2998	REAL(kind=ip_realwp_p) ZOUT(I1:I2,J1:J2)
2999	REAL(kind=ip_realwp_p) ZI(NI,NJ)
3000	LOGICAL SYM,VECT
3001	REAL(kind=ip_realwp_p) SIGN
3002	IF( (VECT))THEN
3003	SIGN = -1.0
3004	ELSE
3005	SIGN = 1.0
3006	ENDIF
3007	IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
3008	$ GRTYP.EQ. 'Y'))THEN
3009	DO 23004 J=1,NJ
3010	DO 23006 I=1,NI
3011	ZOUT(I,J) = ZI(I,J)
3012	23006 CONTINUE
3013	23004 CONTINUE
3014	RETURN
3015	ENDIF
3016	IF( (GRTYP.EQ. 'Z') ) THEN
3017	DO J=1,NJ
3018	DO I=1,NI
3019	ZOUT(I,J) = ZI(I,J)
3020	END DO
3021	END DO
3022	DO J=1,NJ
3023	ZOUT(0,J) = ZI(NI,J)
3024	ZOUT(NI+1,J) = ZOUT(1,J)
3025	ZOUT(NI+2,J) = ZOUT(2,J)
3026	END DO
3027	DO I=0,NI/2
3028	ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
3029	ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
3030	ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3031	ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3032	ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
3033	ZOUT(I+NI/2,-1)= SIGN * ZOUT(I, 2)
3034	ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
3035	ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
3036	END DO
3037	DO I=1,2
3038	ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
3039	ZOUT(NI+I,-1)= SIGN * ZOUT(I+NI/2, 2)
3040	ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3041	ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3042	END DO
3043	ENDIF
3044	IF( (GRTYP.EQ.'A'.OR. GRTYP.EQ.'G'))THEN
3045	IF( (IG2.EQ. NORDSUD))THEN
3046	CALL PERMUT(ZI, NI, NJ)
3047	ENDIF
3048	DO 23012 J=1,NJ
3049	DO 23014 I=1,NI
3050	ZOUT(I,J) = ZI(I,J)
3051	23014 CONTINUE
3052	23012 CONTINUE
3053	DO 23016 J=1,NJ
3054	ZOUT(0,J) = ZI(NI,J)
3055	ZOUT(NI+1,J) = ZOUT(1,J)
3056	ZOUT(NI+2,J) = ZOUT(2,J)
3057	23016 CONTINUE
3058	IF( (IG1.EQ. GLOBAL))THEN
3059	DO 23020 I=0,NI/2
3060	ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
3061	ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
3062	ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3063	ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3064	ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
3065	ZOUT(I+NI/2,-1)= SIGN * ZOUT(I, 2)
3066	ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
3067	ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
3068	23020 CONTINUE
3069	DO 23022 I=1,2
3070	ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
3071	ZOUT(NI+I,-1)= SIGN * ZOUT(I+NI/2, 2)
3072	ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3073	ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3074	23022 CONTINUE
3075	ENDIF
3076	IF( (IG1.EQ. NORD))THEN
3077	IF( (SYM))THEN
3078	DO 23028 J=-NJ+1, 0
3079	DO 23030 I=0,NI+2
3080	ZOUT(I,J) = ZOUT(I, -J+1)
3081	23030 CONTINUE
3082	23028 CONTINUE
3083	ELSE
3084	DO 23032 J=-NJ+1, 0
3085	DO 23034 I=0,NI+2
3086	ZOUT(I,J) = -ZOUT(I, -J+1)
3087	23034 CONTINUE
3088	23032 CONTINUE
3089	ENDIF
3090	DO 23036 I=0,NI/2
3091	ZOUT(I,-NJ) = SIGN * ZOUT(I+NI/2, -NJ+1)
3092	ZOUT(I,-NJ-1)= SIGN * ZOUT(I+NI/2,-NJ+2)
3093	ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3094	ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3095	ZOUT(I+NI/2,-NJ) = SIGN * ZOUT(I, -NJ+1)
3096	ZOUT(I+NI/2,-NJ-1)= SIGN * ZOUT(I, -NJ+2)
3097	ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
3098	ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
3099	23036 CONTINUE
3100	DO 23038 I=1,2
3101	ZOUT(NI+I,-NJ) = SIGN * ZOUT(NI/2+I,-NJ+1)
3102	ZOUT(NI+I,-NJ-1)= SIGN * ZOUT(I+NI/2, -NJ+2)
3103	ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
3104	ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
3105
3106	23038 CONTINUE
3107	ENDIF
3108	IF( (IG1.EQ. SUD))THEN
3109	IF( (SYM))THEN
3110	DO 23044 J=1,NJ
3111	DO 23046 I=0,NI+2
3112	ZOUT(I,NJ+J) = ZOUT(I, NJ-J+1)
3113	23046 CONTINUE
3114	23044 CONTINUE
3115	ELSE
3116	DO 23048 J=1,NJ
3117	DO 23050 I=0,NI+2
3118	ZOUT(I,NJ+J) = -ZOUT(I,NJ-J+1)
3119	23050 CONTINUE
3120	23048 CONTINUE
3121	ENDIF
3122	DO 23052 I=0,NI/2
3123	ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
3124	ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
3125	ZOUT(I,2NJ+1)=SIGN ZOUT(I+NI/2, 2*NJ)
3126	ZOUT(I,2NJ+2)=SIGN ZOUT(I+NI/2, 2*NJ-1)
3127	ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
3128	ZOUT(I+NI/2,-1)= SIGN * ZOUT(I,2)
3129	ZOUT(I+NI/2,2NJ+1)=SIGN ZOUT(I, 2*NJ)
3130	ZOUT(I+NI/2,2NJ+2)=SIGN ZOUT(I, 2*NJ-1)
3131	23052 CONTINUE
3132	DO 23054 I=1,2
3133	ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
3134	ZOUT(NI+I,1)= SIGN * ZOUT(I+NI/2, 2)
3135	ZOUT(NI+I,2NJ+1)=SIGN ZOUT(I+NI/2, 2*NJ)
3136	ZOUT(NI+I,2NJ+2)=SIGN ZOUT(I+NI/2, 2*NJ-1)
3137	23054 CONTINUE
3138	ENDIF
3139	ENDIF
3140	IF( (GRTYP.EQ.'B'))THEN
3141	IF( (IG2.EQ. NORDSUD))THEN
3142	CALL PERMUT(ZI, NI, NJ)
3143	ENDIF
3144	DO 23060 J=1,NJ
3145	DO 23062 I=1,NI
3146	ZOUT(I,J) = ZI(I,J)
3147	23062 CONTINUE
3148	23060 CONTINUE
3149	DO 23064 J=1,NJ
3150	ZOUT(0,J) = ZI(NI-1,J)
3151	ZOUT(NI+1,J) = ZOUT(2,J)
3152	23064 CONTINUE
3153	IF( (IG1.EQ. GLOBAL))THEN
3154	DO 23068 I=0,NI/2+1
3155	ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 2)
3156	ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ-1)
3157	ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 2)
3158	ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ-1)
3159	23068 CONTINUE
3160	ZOUT(NI+1,0) = SIGN * ZOUT(NI/2+1,2)
3161	ZOUT(NI+1,NJ+1)=SIGN * ZOUT(NI/2+1, NJ-1)
3162	ENDIF
3163	IF( (IG1.EQ. NORD))THEN
3164	IF( (SYM))THEN
3165	DO 23074 J=-NJ+2, 0
3166	DO 23076 I=0,NI+1
3167	ZOUT(I,J) = ZOUT(I, -J+2)
3168	23076 CONTINUE
3169	23074 CONTINUE
3170	ELSE
3171	DO 23078 J=-NJ+2, 0
3172	DO 23080 I=0,NI+1
3173	ZOUT(I,J) = -ZOUT(I, -J+2)
3174	23080 CONTINUE
3175	23078 CONTINUE
3176	ENDIF
3177	DO 23082 I=0,NI/2+1
3178	ZOUT(I,-NJ+1) = SIGN * ZOUT(I+NI/2, -NJ+3)
3179	ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ-1)
3180	ZOUT(I+NI/2,-NJ+1) = SIGN * ZOUT(I, -NJ+3)
3181	ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ-1)
3182	23082 CONTINUE
3183	ZOUT(NI+1,-NJ+1) = SIGN * ZOUT(NI/2+1,-NJ+3)
3184	ZOUT(NI+1,NJ+1)=SIGN * ZOUT(1+NI/2, NJ-1)
3185	ENDIF
3186	IF( (IG1.EQ. SUD))THEN
3187	IF( (SYM))THEN
3188	DO 23088 J=1,NJ-1
3189	DO 23090 I=0,NI+1
3190	ZOUT(I,NJ+J) = ZOUT(I, NJ-J)
3191	23090 CONTINUE
3192	23088 CONTINUE
3193	ELSE
3194	DO 23092 J=1,NJ-1
3195	DO 23094 I=0,NI+1
3196	ZOUT(I,NJ+J) = -ZOUT(I,NJ-J)
3197	23094 CONTINUE
3198	23092 CONTINUE
3199	ENDIF
3200	DO 23096 I=0,NI/2+1
3201	ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 2)
3202	ZOUT(I,2NJ)=SIGN ZOUT(I+NI/2, 2*NJ-2)
3203	ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 2)
3204	ZOUT(I+NI/2,2NJ)=SIGN ZOUT(I, 2*NJ-2)
3205	23096 CONTINUE
3206	ZOUT(NI+1,0) = SIGN * ZOUT(NI/2+1,2)
3207	ZOUT(NI+1,2NJ)=SIGN ZOUT(1+NI/2, 2*NJ-2)
3208	ENDIF
3209	ENDIF
3210	RETURN
3211	END
3212	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3213	C
3214	SUBROUTINE LLFXY(DLAT,DLON,X,Y,D60,DGRW,NHEM)
3215	USE mod_kinds_oasis
3216	C
3217	C AUTHOR - J.D. HENDERSON - FEB 75
3218	C
3219	C OBJECT(LLFXY)
3220	C - COMPUTES LATITUDE AND LONGITUDE OF A POINT IN A POLAR
3221	C STEREOGRAPHIC GRID FROM CO-ORDINATES IN THE GRID MEASURED
3222	C FROM THE POLE.
3223	C
3224	C USAGE - CALL LLFXY(DLAT,DLON,X,Y,D60,DGRW,NHEM)
3225	C
3226	C ARGUMENTS
3227	C OUT - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N).
3228	C - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E).
3229	C IN - X - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
3230	C AS ORIGIN
3231	C - Y - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
3232	C AS ORIGIN
3233	C - D60 - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
3234	C GRID AT 60 DEGREES
3235	C - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
3236	C THE GRID (IN DEGREES)
3237	C - NHEM - 1 FOR NORTHERN HEMISPHERE
3238	C 2 FOR SOUTHERN HEMISPHERE
3239	C
3240	C NOTES - THE COMPANION ROUTINE XYFLL, WHICH COMPUTES THE GRID
3241	C CO-ORDINATES GIVEN THE LATITUDE AND LONGITUDE, IS ALSO
3242	C AVAILABLE.
3243	C
3244	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
3245	C * RDTODG = 180/PIE, DGTORD = PIE/180
3246	C
3247	DATA PIE /3.14159265358979323846/
3248	DATA RDTODG /57.295779513082/
3249	DATA DGTORD /1.7453292519943E-2/
3250	C* ------------------------------------------------------------------------
3251	RE=1.866025*6.371E+6/D60
3252	RE2=RE**2
3253	C * IF POINT IS AT POLE SET COORD TO (0.,90.).
3254	DLAT=90.
3255	DLON=0.
3256	IF(X.EQ.0. .AND. Y.EQ.0.) GO TO 39
3257	C
3258	C * CALCULATE LONGITUDE IN MAP COORDINATES.
3259	C
3260	IF(X.EQ.0.) DLON=SIGN(90.,Y)
3261	IF(X.NE.0.) DLON=ATAN(Y/X)*RDTODG
3262	IF(X.LT.0.) DLON=DLON+SIGN(180.,Y)
3263	C
3264	C * ADJUST LONGITUDE FOR GRID ORIENTATION.
3265	C
3266	DLON=DLON-DGRW
3267	IF(DLON.GT.+180.) DLON=DLON-360.
3268	IF(DLON.LT.-180.) DLON=DLON+360.
3269	C
3270	C * CALCULATE LATITUDE.
3271	C
3272	R2=X2+Y2
3273	DLAT=(RE2-R2)/(RE2+R2)
3274	DLAT= ASIN(DLAT)*RDTODG
3275	C
3276	C * CHANGE SIGNS IF IN SOUTHERN HEMISPHERE.
3277	C
3278	39 IF(NHEM.EQ.2) DLAT=-DLAT
3279	IF(NHEM.EQ.2) DLON=-DLON
3280	C
3281	RETURN
3282	END
3283	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3284	C
3285	SUBROUTINE XYFLL(X,Y,DLAT,DLON,D60,DGRW,NHEM)
3286	USE mod_kinds_oasis
3287	C
3288	C AUTHOR - J.D. HENDERSON - FEB 75
3289	C
3290	C OBJECT(XYFLL)
3291	C - COMPUTES THE GRID CO-ORDINATES MEASURED FROM THE POLE OF A
3292	C POINT, GIVEN THE LATITUDE AND LONGITUDE IN DEGREES.
3293	C
3294	C USAGE - CALL XYFLL(X,Y,DLAT,DLON,D60,DGRW,NHEM)
3295	C
3296	C ARGUMENTS
3297	C OUT - X - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
3298	C AS ORIGIN
3299	C AS ORIGIN
3300	C - Y - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
3301	C AS ORIGIN
3302	C IN - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N)
3303	C - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E)
3304	C - D60 - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
3305	C GRID AT 60 DEGREES
3306	C - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
3307	C THE GRID (IN DEGREES)
3308	C - NHEM - 1 FOR NORTHERN HEMISPHERE
3309	C 2 FOR SOUTHERN HEMISPHERE
3310	C
3311	C NOTES - THE COMPANION ROUTINE LLFXY, WHICH COMPUTES THE LATITUDE
3312	C AND LONGITUDE GIVEN THE GRID-COORDINATES, IS ALSO AVAILABLE.
3313	C
3314	C---------------------------------------------------------------------------
3315	C
3316	C * 1.866025=(1+SIN60), 6.371E+6=EARTH RADIUS IN METERS.
3317	C
3318	C RDTODG = 180/PIE, DGTORD = PIE/180
3319	DATA PIE /3.14159265358979323846/
3320	DATA RDTODG /57.295779513082/
3321	DATA DGTORD /1.7453292519943E-2/
3322	C
3323	RE=1.866025*6.371E+6/D60
3324	C
3325	GLON=DLON
3326	IF(NHEM.EQ.2) GLON=-DLON
3327	GLAT=DLAT
3328	IF(NHEM.EQ.2) GLAT=-DLAT
3329	C
3330	RLON=DGTORD*(GLON+DGRW)
3331	RLAT=DGTORD*GLAT
3332	SINLAT=SIN(RLAT)
3333	R=RE*SQRT((1.-SINLAT)/(1.+SINLAT))
3334	X=R*COS(RLON)
3335	Y=R*SIN(RLON)
3336	C
3337	RETURN
3338	END
3339	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3340	C
3341	SUBROUTINE ORDLEG(SBX,COA,IR)
3342	USE mod_kinds_oasis
3343	C
3344	C AUTEUR - D. ROBERTSON
3345	C
3346	C OBJET(ORDLEG)
3347	C - THIS ROUTINE IS A SUBSET OF BELOUSOVS ALGORITHM
3348	C USED TO CALCULATE ORDINARY LEGENDRE POLYNOMIALS.
3349	C
3350	C USAGE - CALL ORDLEG(SBX,COA,IR)
3351	C
3352	C ARGUMENTS
3353	C OUT - SBX - LEGENDRE POLYNOMIAL EVALUATED AT COA
3354	C IN - COA - COSINE OF COLATITUDE
3355	C IN - IR - WAVE NUMBER
3356	C
3357	C* -------------------------------------------------------------------------
3358	C
3359	C
3360	C RDTODG = 180/PIE, DGTORD = PIE/180
3361	DATA PIE /3.14159265358979323846/
3362	DATA RDTODG /57.295779513082/
3363	DATA DGTORD /1.7453292519943E-2/
3364	C
3365	PI = PIE
3366	SQR2=SQRT(2.)
3367	IRPP = IR + 1
3368	IRPPM = IRPP - 1
3369	DELTA = ACOS(COA)
3370	SIA = SIN(DELTA)
3371	C
3372	THETA=DELTA
3373	C1=SQR2
3374	C
3375	DO 20 N=1,IRPPM
3376	FN=FLOAT(N)
3377	FN2=2.*FN
3378	FN2SQ=FN2*FN2
3379	C1=C1* SQRT(1.0-1.0/FN2SQ)
3380	20 CONTINUE
3381	C
3382	N=IRPPM
3383	ANG=FN*THETA
3384	S1=0.0
3385	C4=1.0
3386	A=-1.0
3387	B=0.0
3388	N1=N+1
3389	C
3390	DO 27 KK=1,N1,2
3391	K=KK-1
3392	IF (K.EQ.N) C4=0.5*C4
3393	S1=S1+C4* COS(ANG)
3394	A=A+2.0
3395	B=B+1.0
3396	FK=FLOAT(K)
3397	ANG=THETA*(FN-FK-2.0)
3398	C4=(A(FN-B+1.0)/(B(FN2-A)))*C4
3399	27 CONTINUE
3400	C
3401	SBX=S1*C1
3402	C
3403	RETURN
3404	END
3405	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3406	C
3407	LOGICAL FUNCTION VALIDE(NOM, ICHK, MIN, MAX)
3408	C
3409	C AUTEUR - P. SARRAZIN - AVRIL 82
3410	C
3411	C OBJET(VALIDE)
3412	C - VERIFIER SI ICHK => MIN OU =< MAX
3413	C - MESSAGE SI ICHK N'EST PAS DANS LES LIMITES
3414	C
3415	C ARGUMENTS
3416	C IN - NOM - NOM DE LA VARIABLE EMPLOYE PAR LA ROUTINE
3417	C IN - ICHK - VALEUR DE LA VARIABLE POUR VERIFICATION
3418	C IN - MIN - VALEUR MINIMUM DE ICHK
3419	C IN - MAX - VALEUR MAXIMUM DE ICHK
3420	C
3421	C MODULES - SCINT - UVINT
3422	C
3423	C* ----------------------------------------------------------------------
3424	IF( (ICHK.LT.MIN .OR. ICHK.GT.MAX))THEN
3425	WRITE(6,600) NOM,ICHK,MIN,MAX
3426	ENDIF
3427	600 FORMAT("Bad value for",I10,"VALEUR=",I10,"MINIMUM=",
3428	$ I10,"MAXIMUM=",I10)
3429	VALIDE = (ICHK.GE.MIN) .AND. (ICHK.LE.MAX)
3430	RETURN
3431	END
3432	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3433	C
3434	SUBROUTINE PERMUT(Z,NI,NJ)
3435	USE mod_kinds_oasis
3436	REAL(kind=ip_realwp_p) Z(NI,NJ)
3437	C
3438	C AUTEUR - M. VALIN - FEV 82
3439	C
3440	C OBJET(PERMUT)
3441	C - ROTATION D'UNE MATRICE AUTOUR DE LA LIGNE DU MILIEU.
3442	C CETTE ROUTINE EST UTILISEE POUR RE=ARRANGER LES
3443	C DONNEES DANS UN CHAMP. EX: POUR CONTOURER, ON A
3444	C DES DONNEES DANS L'ORDRE SUIVANT, DU SUD AU NORD
3445	C ALORS QUE POUR L'INTERPOLATION ELLES DOIVENT ETRE
3446	C ETALEES DU NORD AU SUD.
3447	C
3448	C APPEL - CALL PERMUT(Z,NI,NJ)
3449	C
3450	C ARGUMENTS
3451	C IN/OUT - Z - CHAMP QUI SUBIT LA ROTATION
3452	C IN - NI - PREMIERE DIMENSION DE Z
3453	C IN - NJ - DEUXIEME DIMENSION DE Z
3454	C* -----------------------------------------------------------------------
3455	C
3456	NCC = NJ/2
3457	DO 23000 J=1,NCC
3458	DO 23002 I=1,NI
3459	T = Z(I,NJ+1-J)
3460	Z(I,NJ+1-J) = Z(I,J)
3461	Z(I,J) = T
3462	23002 CONTINUE
3463	23000 CONTINUE
3464	RETURN
3465	END
3466	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3467	C
3468	SUBROUTINE CIGAXG(CGTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
3469	USE mod_kinds_oasis
3470	CHARACTER*1 CGTYP
3471	C
3472	C AUTEUR - M. VALIN - FEV 82
3473	C
3474	C
3475	C OBJET(CIGAXG)
3476	C - PASSE DES PARAMETRES (ENTIERS) DESCRIPTEURS DE GRILLE
3477	C AUX PARAMETRES REELS.
3478	C
3479	C APPEL - CALL IGAXG(CGTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
3480	C
3481	C ARGUMENTS
3482	C IN - CGTYP - TYPE DE GRILLE (VOIR OUVRIR)
3483	C OUT - XG1 - ** DESCRIPTEUR DE GRILLE (REEL),
3484	C OUT - XG2 - IGTYP = 'N', PI, PJ, D60, DGRW
3485	C OUT - XG3 - IGTYP = 'L', LAT0, LON0, DLAT, DLON,
3486	C OUT - XG4 - IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
3487	C = 1. NORD
3488	C = 2. SUD **
3489	C IN - IG1 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3490	C IN - IG2 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3491	C IN - IG3 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3492	C IN - IG4 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3493	C
3494	CMESSAGES - "ERREUR, MAUVAISE SPECIFICATION DE GRILLE, (TYPE) (IGAXG)"
3495	C
3496	C-------------------------------------------------------------------
3497	C
3498
3499	IF( ((CGTYP.EQ. 'N') .OR. (CGTYP.EQ.'S')))THEN
3500	IF((IG4.LT. 32768))THEN
3501	XG1 = IG2 * 0.1
3502	XG2 = IG1 * 0.1
3503	XG3 = IG4 * 100.
3504	XG4 = IG3 * 0.01
3505	ELSE
3506	JG3 = IG3
3507	JG4 = IG4
3508	JG4 = JG4 - 32768
3509	XG3 = IG1 * 100.
3510	IF((IG3 .GT. 32767))THEN
3511	XG3 = XG3 * 10.
3512	JG3 = JG3 - 32768
3513	ENDIF
3514	XG4 = IG2 * .1
3515	IF((JG4.GT. 16383))THEN
3516	XG4 = 360. - XG4
3517	JG4 = JG4 - 16384
3518	ENDIF
3519	DLAT = 90. -(JG4*180./16383.)
3520	DLON = (JG3*360./32767.)
3521	IHEM = 1
3522	IF(('S'.EQ.CGTYP))THEN
3523	IHEM = 2
3524	ENDIF
3525	CALL XYFLL(XG1,XG2,DLAT,DLON,XG3,XG4,IHEM)
3526	XG1 = 1.0 - XG1
3527	XG2 = 1.0 - XG2
3528	ENDIF
3529	ELSE
3530	IF((CGTYP.EQ. 'C'))THEN
3531	XG1 = IG3 * 0.01 - 90.
3532	XG2 = IG4 * 0.01
3533	XG3 = 180. / IG1
3534	XG4 = 360. / IG2
3535	ELSE
3536	IF( ((CGTYP.EQ. 'A') .OR. (CGTYP.EQ. 'B') .OR. (CGTYP
3537	$ .EQ. 'G')))THEN
3538	XG1 = IG1
3539	XG2 = IG2
3540	XG3 = 0.
3541	XG4 = 0.
3542	ELSE
3543	IF((CGTYP.EQ. 'L'))THEN
3544	XG1 = IG3 * 0.01 - 90.
3545	XG2 = IG4 * 0.01
3546	XG3 = IG1 * 0.01
3547	XG4 = IG2 * 0.01
3548	ELSE
3549	IF((CGTYP.EQ. 'H'))THEN
3550	XG1 = IG3
3551	XG2 = .01*IG4 - 90.
3552	XG3 = 500*IG2
3553	XG4 = IG1*.2
3554	ELSE
3555	WRITE(6,600)
3556	ENDIF
3557	ENDIF
3558	ENDIF
3559	ENDIF
3560	ENDIF
3561	600 FORMAT(' Error bad grid specification , (TYPE)'
3562	$,'(IGAXG)')
3563	RETURN
3564	END
3565	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3566	C
3567	SUBROUTINE CXGAIG(CGTYP,IG1,IG2,IG3,IG4,XG1,XG2,XG3,XG4)
3568	USE mod_kinds_oasis
3569	LOGICAL VALIDE
3570	EXTERNAL VALIDE
3571	CHARACTER * 1 CGTYP
3572	LOGICAL LFLAG
3573	C
3574	C AUTEUR - M. VALIN - FEV 82
3575	C
3576	C OBJET(XGAIG)
3577	C - PASSE DES PARAMETRES (REELS) DESCRIPTEURS DE GRILLE
3578	C AUX PARAMETRES ENTIERS.
3579	C
3580	C APPEL - CALL XGAIG(CGTYP,IG1,IG2,IG3,IG4,XG1,XG2,XG3,XG4)
3581	C
3582	C ARGUMENTS
3583	C IN - CGTYP - TYPE DE GRILLE (VOIR OUVRIR)
3584	C OUT - IG1 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3585	C OUT - IG2 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3586	C OUT - IG3 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3587	C OUT - IG4 - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
3588	C IN - XG1 - ** DESCRIPTEUR DE GRILLE (REEL),
3589	C IN - XG2 - IGTYP = 'N', PI, PJ, D60, DGRW
3590	C IN - XG3 - IGTYP = 'L', LAT0, LON0, DLAT, DLON,
3591	C IN - XG4 - IGTYP = 'A', 'B', 'G', XG1 = 0, GLOBAL
3592	C = 1, NORD
3593	C = 2, SUD **
3594	C
3595	C MESSAGES - "ERREUR DANS LA DESCRIPTION DE LA GRILLE (IG1) (XGAIG)"
3596	C "ERREUR, MAUVAISE SPECIFICATION (LAT0) (XGAIG)"
3597	C "ERREUR, GRILLE INCONNUE (TYPE) (XGAIG)"
3598	C
3599	C* ------------------------------------------------------------------
3600	C
3601	C
3602	IF( (CGTYP.EQ. 'N' .OR. CGTYP.EQ.'S'))THEN
3603	IG1 = NINT(XG2 * 10.)
3604	IG2 = NINT(XG1 * 10.)
3605	IG3 = NINT(XG4 * 100.)
3606	IG4 = NINT(XG3 * 0.01)
3607	23002 IF( (IG3.LT. 0))THEN
3608	IG3 = IG3 + 36000
3609	GOTO 23002
3610	ENDIF
3611	IF((IG1.LT.0.OR.IG2.LT.0.OR.IG1.GT.2047.OR.IG2.GT.2047.OR.
3612	$ IG4.GT.32000))THEN
3613	IG1 = 0
3614	IG2 = 0
3615	IG3 = 0
3616	IG4 = 32768
3617	IF((XG3.GT. 204700))THEN
3618	IG3 = 32768
3619	IG1 = NINT(XG3*.001)
3620	ELSE
3621	IG3 = 0
3622	IG1 = NINT(XG3*.01)
3623	ENDIF
3624	IG2 = NINT(XG4*10)
3625	IF((IG2.LT.0))THEN
3626	IG2 = ABS(IG2)
3627	IG4 = IG4 + 16384
3628	ENDIF
3629	IF((IG2.GT.1800))THEN
3630	IG2 = ABS(IG2 - 3600)
3631	IG4 = IG4 + 16384
3632	ENDIF
3633	IHEM = 1
3634	IF(('S'.EQ.CGTYP))THEN
3635	IHEM = 2
3636	ENDIF
3637	CALL LLFXY(DLAT,DLON,1.-XG1,1.-XG2,XG3,XG4,IHEM)
3638	DLAT = 90. - DLAT
3639	IF((DLON.LT.0))THEN
3640	DLON = DLON + 360.
3641	ENDIF
3642	IG3 = IG3 + NINT(DLON*32767./360.)
3643	IG4 = IG4 + NINT(DLAT*16383./180.)
3644	ENDIF
3645	ELSE
3646	IF( (CGTYP.EQ. 'A' .OR. CGTYP.EQ. 'B' .OR. CGTYP.EQ. 'G')
3647	$ )THEN
3648	IG1 = XG1
3649	IG2 = XG2
3650	IG3 = 0
3651	IG4 = 0
3652	LFLAG=VALIDE("IG1",IG1,0,2)
3653	LFLAG=VALIDE("IG2",IG2,0,1)
3654	ELSE
3655	IF((CGTYP.EQ. 'C'))THEN
3656	IG1 = 180. / XG3 + 0.5
3657	IG2 = 360. / XG4 + 0.5
3658	IG3 = (90. + XG1) * 100. + 0.5
3659	IG4 = XG2 * 100. + 0.5
3660	23020 IF( (IG4.LT. 0))THEN
3661	IG4 = IG4 + 36000
3662	GOTO 23020
3663	ENDIF
3664	IF( (IG3.LT. 0))THEN
3665	WRITE(6,601)
3666	ENDIF
3667	ELSE
3668	IF( (CGTYP.EQ. 'H'))THEN
3669	IG1 = NINT(5.*XG4)
3670	23026 IF( (IG1.LT. 0))THEN
3671	IG1 = IG1 + 1800
3672	GOTO 23026
3673	ENDIF
3674	IG2 = NINT(.002*XG3)
3675	IG3 = NINT(XG1)
3676	IG4 = NINT(100.*(90.+XG2))
3677	ELSE
3678	IF( (CGTYP.EQ. 'L'))THEN
3679	IG1 = XG3 * 100. + 0.5
3680	IG2 = XG4 * 100. + 0.5
3681	IG3 = (90. + XG1) * 100. + 0.5
3682	IG4 = XG2 * 100. + 0.5
3683	23030 IF( (IG4.LT. 0))THEN
3684	IG4 = IG4 + 36000
3685	GOTO 23030
3686	ENDIF
3687	IF( (IG3.LT. 0))THEN
3688	WRITE(6,601)
3689	ENDIF
3690	ELSE
3691	WRITE(6,602)
3692	ENDIF
3693	ENDIF
3694	ENDIF
3695	ENDIF
3696	ENDIF
3697	601 FORMAT(' Error, bad specification (LAT0) (XGAIG)')
3698	602 FORMAT(' Error, unknown grid (TYPE) (XGAIG)')
3699	RETURN
3700	END
3701	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Context Navigation

source: CPL/oasis3/trunk/src/lib/fscint/src/discendo.f @ 1677

Download in other formats: