Oasis3 4.0.2
discendo.f
Go to the documentation of this file.
00001 C* ----------------------------------------------------------------------
00002 C 
00003 C                            FSCINT
00004 C                            ------
00005 C               *******************************
00006 C               * OASIS INTERPOLATION PACKAGE *
00007 C               * ----- ------------- ------- *
00008 C               *******************************
00009 C 
00010 C* This is the Fast Scalar INTerpolator package initially written by Yves
00011 C* Chartier and colleagues at RPN Canada. This software has been adapted
00012 C* to the OASIS structure and provides several interpolation techniques
00013 C* (see oasis documentation for details). This package has been split in
00014 C* 2 files: discendo.f and semper.f (for people not fluent in latin,
00015 C* "semper discendo" means "always learning" ... ;-)...). discendo.f
00016 C* contains all the interpolation related routines while semper.f has
00017 C* the memory allocation using fortran 90 features.
00018 C
00019 C     History:
00020 C     -------
00021 C       Version   Programmer     Date      Description
00022 C       -------   ----------     ----      -----------
00023 C       1.1       L. Terray      94/06/01  Introduced
00024 C       1.1       P. Braconnot   94/06/15  Extend longitude range for Z grids
00025 C       2.0       L. Terray      96/09/25 
00026 C       2.2       G. Risari      97/12/01  new routine names 
00027 C       2.2       A.P, S.V, L.T  97/12/14  Change memory allocation
00028 C       2.3       L. Terray      99/09/15  Introduction of Y grid and cleaning
00029 C       2.4       E. Maisonnave  00/07/03  REAL NEWAXEY(I1:I2) (line 2806)
00030 C
00031 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00032 C
00033       FUNCTION FINDLON(VAL, TABLEAU, NBELEM, TMP)
00034 C* ----------------------------------------------------------------------
00035 C
00036 C* New function to account for strange AGCM Z-grids (longitude<0 or >360)
00037 C
00038 C* ----------------------------------------------------------------------
00039       USE mod_kinds_oasis
00040       INTEGER (kind=ip_intwp_p) FINDLON
00041       INTEGER (kind=ip_intwp_p) NBELEM
00042       REAL(kind=ip_realwp_p) VAL
00043       REAL(kind=ip_realwp_p) TABLEAU(NBELEM)
00044       REAL(kind=ip_realwp_p) TMP
00045       INTEGER (kind=ip_intwp_p) DEBUT, MILIEU, FIN
00046 C
00047 C* Get longitudes in the right interval
00048 C
00049       IF( VAL .GT. TABLEAU(NBELEM)) THEN
00050           TMP = VAL - 360.
00051       ELSEIF( VAL .LT. TABLEAU(1)) THEN
00052           TMP = VAL + 360.
00053       ELSE
00054           TMP = VAL
00055       ENDIF
00056 C
00057 C* Find grid point coordinates
00058 C
00059       DEBUT = 1
00060       FIN = NBELEM
00061       MILIEU = (DEBUT+FIN)*0.5
00062 23000 IF( (MILIEU.NE. DEBUT))THEN
00063          IF( TMP .LE. TABLEAU(MILIEU) ) THEN
00064             FIN   = MILIEU
00065          ELSE 
00066             DEBUT = MILIEU
00067          ENDIF 
00068          MILIEU = (DEBUT+FIN)*0.5
00069          GOTO 23000
00070       ENDIF 
00071       FINDLON = MILIEU
00072       RETURN
00073       END
00074 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00075 C
00076       FUNCTION CHERCHE(VAL, TABLEAU, NBELEM)
00077       USE mod_kinds_oasis
00078       INTEGER (kind=ip_intwp_p) CHERCHE
00079       INTEGER (kind=ip_intwp_p) NBELEM
00080       REAL(kind=ip_realwp_p) VAL
00081       REAL(kind=ip_realwp_p) TABLEAU(NBELEM)
00082       INTEGER (kind=ip_intwp_p) DEBUT, MILIEU, FIN
00083       DEBUT = 1
00084       FIN = NBELEM
00085       MILIEU = (DEBUT+FIN)*0.5
00086 23000 IF( (MILIEU.NE. DEBUT))THEN
00087          IF( (VAL.LE. TABLEAU(MILIEU)))THEN
00088             FIN   = MILIEU
00089          ELSE 
00090             DEBUT = MILIEU
00091          ENDIF 
00092          MILIEU = (DEBUT+FIN)*0.5
00093          GOTO 23000
00094       ENDIF 
00095       CHERCHE = MILIEU
00096       RETURN
00097       END
00098 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00099 C
00100       SUBROUTINE CHKXTRAP(PX, PY, NPTS, NI, NJ)
00101       USE mod_kinds_oasis
00102       IMPLICIT NONE
00103       INTEGER (kind=ip_intwp_p) NPTS, NI, NJ, OFFL, OFFR
00104       REAL(kind=ip_realwp_p) PX(NPTS), PY(NPTS)
00105       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
00106       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
00107       PARAMETER (VOISIN  =   0)
00108       PARAMETER (LINEAIR =   1)
00109       PARAMETER (CUBIQUE =   3)
00110       PARAMETER (OUI   =   1)
00111       PARAMETER (MINIMUM =   2)
00112       PARAMETER (MAXIMUM =   3)
00113       PARAMETER (VALEUR  =   4)
00114       PARAMETER (ABORT   =  13)
00115       LOGICAL FLGXTRAP
00116       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
00117       REAL(kind=ip_realwp_p) VALXTRAP
00118       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
00119       INTEGER (kind=ip_intwp_p) N, I, J
00120       FLGXTRAP = .FALSE.
00121       IF( (ORDINT.EQ. 3))THEN
00122          OFFR = 3
00123          OFFL = 2
00124       ELSE 
00125          OFFR = 1
00126          OFFL = 0
00127       ENDIF 
00128       DO 23002 N=1, NPTS
00129          I = INT(PX(N))
00130          J = INT(PY(N))
00131          IF( (I.LT. OFFL.OR. J.LT. OFFL.OR. I.GT. (NI-OFFR).OR. J
00132      $   .GT. (NJ-OFFR)))THEN
00133             FLGXTRAP = .TRUE.
00134          ENDIF 
00135 23002 CONTINUE 
00136       IF( (FLGXTRAP.AND. CODXTRAP.EQ. ABORT))THEN
00137          WRITE(6, *)
00138      $
00139 '*****************************************************************     $****'
00140          WRITE(6, *)
00141      $       'target grid contains points outside of source grid'
00142          WRITE(6, *)  'We stop the program'
00143          WRITE(6, *)
00144      $
00145 '*****************************************************************     $****'
00146          STOP
00147       ENDIF 
00148       RETURN
00149       END
00150 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00151 C
00152       SUBROUTINE DGAUSS (N,ROOTS,KASE)
00153       USE mod_kinds_oasis
00154       REAL(kind=ip_realwp_p) ROOTS(*)
00155 C
00156 C AUTEUR- D. ROBERTSON
00157 C
00158 C OBJET(DGAUSS)
00159 C     - CALCULATES THE ZEROES OF THE ORDINARY LEGENDRE
00160 C     POLYNOMIAL OF ORDER N,  I.E. DEFINE GAUSSIAN GRID
00161 C
00162 C ALGORITHME
00163 C     - THE POSITIVE ROOTS ARE APPROXIMATED BY THE BEST
00164 C     ASYMPTOTIC FORMULA AVAILABLE TO THE AUTHOR, FOUND IN
00165 C     ABRAMOWITZ AND STEGUN "HANDBOOK OF MATHEMATICAL FUNCTIONS".
00166 C     CHAPTER 22 FORMULA 22.16.6.
00167 C     NEWTON'S METHOD IS USED TO REFINE THE GUESS TO PRECISION
00168 C     DEFINED BY THE CONSTANT TOL.  SINCE THE ROOTS ARE OF ORDER
00169 C     OF MAGNITUDE UNITY, ABSOLUTE PRECISION IS ADEQUATE, RATHER
00170 C     THAN A RELATIVE TEST.
00171 C     A STANDARD IDENTITY IS USED TO DETERMINE THE DERIVATIVE OF
00172 C     THE POLYNOMIAL IN TERMS OF THE VALUES OF P(N;X),P(N-1;X).
00173 C     (X**2-1.0)*(DP/DX)=N*(X*P(N;,X)-P(N-1;X)).
00174 C     SEE ABRAMOWITZ AND STEGUN FORMULA 22.8.5
00175 C     NOTE THAT IN CONTRAST TO OTHER FORMULAS THIS REQUIRES ONLY
00176 C     2 EVALUATIONS OF A LEGENDRE POLYNOMIAL PER ITERATION.
00177 C     NOTE THAT THE COORDINATE USED IS CONVENTIONALLY REFERRED TO
00178 C     AS MU=COS(THETA), RUNNING FROM +1 TO -1, FOR THETA FROM 0 TO
00179 C     PI. THE NEGATIVE ROOTS ARE  FILLED BY SYMMETRY.
00180 C     FOR KASE=GLOBAL, ALL N ROOTS ARE FOUND, WHILE FOR
00181 C     DASE=NORTH/SOUTH ONLY THE +VE/-VE ROOTS ARE FOUND,
00182 C     (INCLUDING 0 IF N IS ODD)  I.E. N/2+MOD(N,2) ROOTS.
00183 C
00184 C
00185 C APPEL - CALL DGAUSS(N,ROOTS,KASE)
00186 C
00187 C ARGUMENTS
00188 C     IN    - N     - ORDER OF THE POLYNOMIALS
00189 C     OUT   - ROOTS - ARRAY CONTAINING THE ZEROES OF THE
00190 C     ORDINARY LEGENDRE POLYNOMIALS
00191 C     IN    - KASE  - =0, GLOBAL
00192 C     =1, NORTH
00193 C     =2, SOUTH
00194 C 
00195 C MODULES APPELES
00196 C     - ORDLEG
00197 C
00198 C ----------------------------------------------------------------------
00199 C
00200 C     THE ANSWERS ARE RETURNED IN ROOTS.
00201 
00202       REAL(kind=ip_realwp_p) NORMN,NORMNM
00203       INTEGER (kind=ip_intwp_p) GLOBAL,NORTH,SOUTH,NORD,SUD
00204       PARAMETER(GLOBAL=0,NORTH=1,NORD=1,SOUTH=2,SUD=2)
00205       DATA TOL /1.0E-06/
00206 C
00207 C     RDTODG = 180/PI, DGTORD = PI/180
00208 C
00209       DATA PI   /3.14159265358979323846/
00210       DATA RDTODG /57.295779513082/
00211       DATA DGTORD /1.7453292519943E-2/
00212 C
00213 C     ORDLEG RETURNS POLYNOMIALS NORMALIZED TO UNIT INTEGRAL.
00214 C     NORMN,NORMNMN RESTORE THE CONVENTION NORMALIZATION, P(N;1.0)=1.0.
00215 C* ----------------------------------------------------------------------
00216 
00217       NORMN =SQRT(2.0/(2.0*N+1.0))
00218       NORMNM=SQRT(2.0/(2.0*N-1.0))
00219       L = N/2
00220 C
00221 C     CALCULATE ASYMPTOTIC APPROXIMATION
00222 C
00223       DO 23000 I=1,L
00224          IF((KASE.NE.SOUTH))THEN
00225             J = I
00226          ENDIF 
00227          IF((KASE.EQ.SOUTH))THEN
00228             J=I+L+MOD(N,2)
00229          ENDIF 
00230          T = (4*J-1)*PI/FLOAT(4*N+2)
00231          IF((KASE.NE.SOUTH))THEN
00232             IRT = I
00233          ENDIF 
00234          IF((KASE.EQ.SOUTH))THEN
00235             IRT = I + MOD(N,2)
00236          ENDIF 
00237          ROOTS(IRT)=COS(T+1.0/(8.0*FLOAT(N**2)*TAN(T)))
00238 
00239 C
00240 23000 CONTINUE 
00241       DO 23010 I=1,L
00242 
00243 C
00244 C     REPEAT 1 NEWTON ITERATION
00245 C     **  BEGIN
00246 C
00247 6        CALL ORDLEG(G,ROOTS(I),N)
00248          CALL ORDLEG(GM,ROOTS(I),N-1)
00249          PN = NORMN*G
00250          PNM= NORMNM*GM
00251 C
00252 C     **         GUESS(K+1)=GUESS(K)-P/(DP/DX)
00253 C
00254 
00255          RDPDX = (ROOTS(I)**2-1.0)/(N*(ROOTS(I)*PN-PNM))
00256          DELTA = -PN*RDPDX
00257          ROOTS(I) = ROOTS(I)+DELTA
00258          IF((ABS(DELTA).GT.TOL))THEN
00259             GO TO 6
00260 C
00261 C     **  END UNTIL ABS(DELTA).LE.TOL
00262 C
00263 
00264          ENDIF 
00265          ROOTS(N+1-I) = -ROOTS(I)
00266 
00267 C
00268 23010 CONTINUE 
00269       IF((MOD(N,2).EQ.0))THEN
00270          RETURN
00271       ENDIF 
00272       IF((KASE.NE.SOUTH))THEN
00273          IRT = L+1
00274       ENDIF 
00275       IF((KASE.EQ.SOUTH))THEN
00276          IRT = 1
00277       ENDIF 
00278       ROOTS(IRT) = 0.0
00279       RETURN
00280       END
00281 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00282 C
00283       SUBROUTINE GAUSS(NRACP,RACP,PG,SIA,RAD,PGSSIN2,SINM1,SINM2,
00284      $SIN2)
00285 C* -----------------------------------------------------------------------
00286 C     CALCULE LES NRACP RACINES POSITIVES DU POLYNOME DE LEGENDRE DE
00287 C     DEGRE 2*NRACP (ICI-APRES NOTE PN) DEFINI SUR L INTERVALLE DES
00288 C     COLATITUDES ALLANT DE 0 (POLE NORD) A PI (POLE SUD). ON SAIT QUE
00289 C     LES 2*NRACP RACINES SONT ANTI-SYMETRIQUES P/R A L EQUATEUR PI/2,
00290 C     ETANT POSITIVES ENTRE COLAT=0 ET COLAT =PI/2.
00291 C     ON CALCULE ENSUITE LES POIDS DE GAUSS ASSOCIES AUX COLATITUDES
00292 C     GAUSSIENNES (ICI APRES NOTEES CG), AINSI QU UN CERTAIN NOMBRE DE
00293 C     FONCTIONS DE CG DEFINIES PLUS LOIN. ON RAPPELLE ENFIN QUE LA LATI-
00294 C     TUDE LAT=COLAT-PI/2, ET DONC QUE SIN(LAT)=COS(COLAT).
00295 C     NRACP        : NOMBRE DE RACINES POSITIVES DU POLYNOME DE LEGENDRE
00296 C                  : DE DEGRE 2*NRACP.
00297 C     RACP(I)      : RACINES DE PN, =SIN(LG)=COS(CG).
00298 C     PG(I)        : POIDS DE GAUSS CORRESPONDANTS.
00299 C     SIA(I)       : SIN(CG)=COS(LG).
00300 C     RAD(I)       : COLATITUDE CG EN RADIANS.
00301 C     PGSSIN2(I)   : POIDS DE GAUSS / (SIN(CG))**2.
00302 C     SINM1(I)     : (SIN(CG))**-1.
00303 C     SINM2(I)     : (SIN(CG))**-2.
00304 C     VOIR NST 8, CHAP. A, PP.1-7, ET APPENDICE D12, PP. 26-27.
00305 C     VERSION REVISEE PAR MICHEL BELAND, 9 DECEMBRE 1980.
00306 C     *****************************************************************
00307 C
00308 C
00309       USE mod_kinds_oasis
00310 C     -----------------------------------------------------------------
00311 
00312       DIMENSION RACP(1),PG(1),SIA(1),RAD(1),PGSSIN2(1),SINM1(1),
00313      $SINM2(1),SIN2(1)
00314 C     --------------------------------------------------------------
00315 C
00316 C     ON DEMANDE UNE PRECISION DE 1.E-13 POUR LES RACINES DE PN.
00317 C
00318 
00319       XLIM=1.E-6
00320       PI = 3.14159265358979323846
00321       IR = 2*NRACP
00322       FI=FLOAT(IR)
00323       FI1=FI+1.
00324       FN=FLOAT(NRACP)
00325 C
00326 C     ON UTILISE UNE FORMULE ASYMPTOTIQUE POUR OBTENIR LES VALEURS
00327 C     APPROXIMATIVES DES COLATITUDES GAUSSIENNES
00328 C     CG(I) = (PI/2) * (2*I-1)/(2*NRACP).
00329 C     VOIR ABRAMOWITZ AND STEGUN, P. 787, EQU. 22.16.6 .
00330 C
00331 
00332       DO 23000 I=1,NRACP
00333          DOT=FLOAT(I-1)
00334          RACP(I)=-PI*.5*(DOT+.5)/FN + PI*.5
00335          RACP(I) =  SIN(RACP(I))
00336 
00337 C
00338 C     ON CALCULE ENSUITE LES CONSTANTES FACTEURS DE P(N+1) ET P(N-1)
00339 C     DANS L EXPRESSION DE LA PSEUDO-DERIVEE DE PN.
00340 C
00341 23000 CONTINUE 
00342       DN = FI/SQRT(4.*FI*FI-1.)
00343       DN1=FI1/SQRT(4.*FI1*FI1-1.)
00344       A = DN1*FI
00345       B = DN*FI1
00346       IRP = IR + 1
00347       IRM = IR -1
00348 C
00349 C     ON EMPLOIE ENSUITE UNE METHODE DE NEWTON POUR AUGMENTER LA PREC.
00350 C     SI RACTEMP EST UNE SOL. APPROXIMATIVE  DE PN(RACP)=0., ALORS LA
00351 C     SEQUENCE RACTEMP(I+1)=RACTEMP(I)-PN(RACTEMP(I))/DER.PN(RACTEMP(I))
00352 C     CONVERGE VERS RACP DE FACON QUADRATIQUE.
00353 C     VOIR ABRAMOWITZ AND STEGUN, P.18, EQU. 3.9.5.
00354 C     ORDLEG CALCULE LA VALEUR DE PN (RACP) , NORMALISE.
00355 C
00356 
00357       DO 23002 I=1,NRACP
00358 42       CALL ORDLEG(G,RACP(I),IR)
00359          CALL ORDLEG(GM,RACP(I),IRM)
00360          CALL ORDLEG(GP,RACP(I),IRP)
00361          GT = (A*GP-B*GM)/(RACP(I)*RACP(I)-1.)
00362          RACTEMP = RACP(I) - G/GT
00363          GTEMP = RACP(I) - RACTEMP
00364          RACP(I) = RACTEMP
00365          IF(( ABS(GTEMP).GT.XLIM))THEN
00366             GO TO 42
00367          ENDIF 
00368 
00369 C
00370 C     ON CALCULE ENSUITE LES POIDS DE GAUSS SELON L ALGORITHME
00371 C     PG(I) = 2./[(1.-RACP(I)**2)*(DER.PN(RACP(I)))**2].
00372 C     VOIR ABRAMOWITZ AND STEGUN, P.887, EQU. 25.4.29.
00373 C     NOTE: ON DOIT MULTIPLIER LA PRECEDENTE FORMULE PAR UN FACTEUR
00374 C     DE DENORMALISATION, LES PN DONNES PAR ORDLEG ETANT NORMALISES.
00375 C     ON SE SERT D UNE FORMULE DE RECURRENCE POUR LA DERIVEE DE PN.
00376 C
00377 23002 CONTINUE 
00378       DO 23006 I=1,NRACP
00379          A=2.*(1.-RACP(I)**2)
00380          CALL ORDLEG(B,RACP(I),IRM)
00381          B = B*B*FI*FI
00382          PG(I)=A*(FI-.5)/B
00383          RAD(I) =   ACOS(RACP(I))
00384          SIA(I) =  SIN(RAD(I))
00385          C=(SIA(I))**2
00386          SINM1(I) = 1./SIA(I)
00387          SINM2(I) = 1./C
00388          PGSSIN2(I) =PG(I)/C
00389          SIN2(I)=C
00390 23006 CONTINUE 
00391       RETURN
00392       END
00393 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00394 C
00395       SUBROUTINE GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
00396       USE mod_kinds_oasis
00397       IMPLICIT NONE
00398       INTEGER (kind=ip_intwp_p) LI,LJ
00399       REAL(kind=ip_realwp_p) SPDO(LI,LJ), PSIO(LI,LJ), XLON(LI,LJ)
00400       CHARACTER*1 GRTYP
00401       INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
00402       EXTERNAL CIGAXG
00403 C
00404 C AUTEUR   - Y. CHARTIER - AVRIL 91
00405 C
00406 C OBJET(GDW2LLW)
00407 C         - PASSE DE VENT DE GRILLE (COMPOSANTES U ET V)
00408 C         - A VITESSE ET DIRECTION (SPEED, PSI)
00409 C APPEL    - CALL GDW2LLW(SPD,PSI,LI,LJ,IYP,XG1,XG2,XG3,XG4)
00410 C
00411 C MODULES  - XGAIG
00412 C
00413 C ARGUMENTS
00414 C  IN/OUT - SPD   - A L'ENTREE CONTIENT LA PSIOTESSE DU VENT ET
00415 C                   A LA SORTIE LA COMPOSANTE U.
00416 C  IN/OUT - PSI   - A L'ENTREE CONTIENT LA DIRECTION DU VENT ET
00417 C                   A LA SORTIE LA COMPOSANTE V.
00418 C   IN    - LI    - PREMIERE DIMENSION DES CHAMPS SPD ET PSI
00419 C   IN    - LJ    - DEUXIEME DIMENSION DES CHAMPS SPD ET PSI
00420 C   IN    - IGTYP  - TYPE DE GRILLE (VOIR OUVRIR)
00421 C   IN    - XG1   - ** DESCRIPTEUR DE GRILLE (REEL),
00422 C   IN    - XG2   -    IGTYP = 'N', PI, PJ, D60, DGRW
00423 C   IN    - XG3   -    IGTYP = 'L', LAT0, LON0, DLAT, DLON,
00424 C   IN    - XG4   -    IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
00425 C                                                 = 1. NORD
00426 C                                                 = 2. SUD **
00427 C
00428 C MESSAGES - "ERREUR MAUVAISE GRILLE (GDW2LLW)"
00429 C
00430 C-------------------------------------------------------------
00431 C
00432 C
00433 C
00434 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
00435 C
00436 C RDTODG = 180/PIE, DGTORD = PIE/180
00437 
00438       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
00439       DATA PIE   /3.14159265358979323846/
00440       DATA RDTODG /57.295779513082/
00441       DATA DGTORD /1.7453292519943E-2/
00442 C
00443 
00444       INTEGER (kind=ip_intwp_p) I,J
00445       REAL(kind=ip_realwp_p) SPD, PSI
00446       REAL(kind=ip_realwp_p) XG1,XG2,XG3,XG4
00447       IF( (GRTYP.EQ. 'N'))THEN
00448          CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
00449          DO 23002 I=1,LI
00450             DO 23004 J=1,LJ
00451                SPD=SQRT(SPDO(I,J)*SPDO(I,J)+PSIO(I,J)*PSIO(I,J))
00452                IF( (SPD.EQ. 0.0))THEN
00453                   PSI= 0.0
00454                ELSE 
00455                   IF( (SPDO(I,J).EQ. 0.0))THEN
00456                      IF( (PSIO(I,J).GE. 0.0))THEN
00457                         PSI= XLON(I,J)+XG4-90.0
00458                      ELSE 
00459                         PSI= XLON(I,J)+XG4+90.0
00460                      ENDIF 
00461                   ELSE 
00462                      PSI=XLON(I,J)+XG4-RDTODG*ATAN2(PSIO(I,J),SPDO(I
00463      $               ,J))
00464                   ENDIF 
00465                ENDIF 
00466                PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
00467      $             +360.0_ip_realwp_p,360.0_ip_realwp_p)
00468                SPDO(I,J)=SPD
00469                PSIO(I,J)=PSI
00470 23004       CONTINUE 
00471 23002    CONTINUE 
00472          RETURN
00473       ENDIF 
00474       IF( (GRTYP.EQ. 'S'))THEN
00475          CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
00476          DO 23014 I=1,LI
00477             DO 23016 J=1,LJ
00478                SPD=SQRT(SPDO(I,J)*SPDO(I,J)+PSIO(I,J)*PSIO(I,J))
00479                IF( (SPD.EQ. 0.0))THEN
00480                   PSI= 0.0
00481                ELSE 
00482                   IF( (SPDO(I,J).EQ. 0.0))THEN
00483                      IF( (PSIO(I,J).GE. 0.0))THEN
00484                         PSI= 90.0 - XLON(I,J)+XG4
00485                      ELSE 
00486                         PSI= 270.0 - XLON(I,J)+XG4
00487                      ENDIF 
00488                   ELSE 
00489                      PSI= 180.0 - XLON(I,J)+XG4-RDTODG*ATAN2(PSIO(I,
00490      $               J),SPDO(I,J))
00491                   ENDIF 
00492                ENDIF 
00493                PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
00494      $             +360.0_ip_realwp_p,360.0_ip_realwp_p)
00495                SPDO(I,J)=SPD
00496                PSIO(I,J)=PSI
00497 23016       CONTINUE 
00498 23014    CONTINUE 
00499          RETURN
00500       ENDIF 
00501       IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'G'.OR.GRTYP.EQ.
00502      $'L'))THEN
00503          DO 23026 I=1,LI
00504             DO 23028 J=1,LJ
00505                SPD=SQRT(SPDO(I,J)*SPDO(I,J)+PSIO(I,J)*PSIO(I,J))
00506                IF( (SPD.EQ. 0.0))THEN
00507                   PSI= 0.0
00508                ELSE 
00509                   IF( (SPDO(I,J).EQ. 0.0))THEN
00510                      IF( (PSIO(I,J).GE. 0.0))THEN
00511                         PSI= 180.0
00512                      ELSE 
00513                         PSI= 0.0
00514                      ENDIF 
00515                   ELSE 
00516                      PSI=270.0 - RDTODG*ATAN2(PSIO(I,J),SPDO(I,J))
00517                   ENDIF 
00518                ENDIF 
00519                PSI=MOD(MOD(PSI,360.0_ip_realwp_p)
00520      $             +360.0_ip_realwp_p,360.0_ip_realwp_p)
00521                SPDO(I,J)=SPD
00522                PSIO(I,J)=PSI
00523 23028       CONTINUE 
00524 23026    CONTINUE 
00525          RETURN
00526       ENDIF 
00527       WRITE(6, 600) GRTYP
00528 600   FORMAT('0',' Error bad grid (GDW2LLW) - GRTYP = ', A1)
00529       RETURN
00530       END
00531 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00532 C
00533       SUBROUTINE GGGDINT(ZO,PX,PY,NPTS,AY,CY,Z,I1,I2,J1,J2)
00534       USE mod_kinds_oasis
00535       IMPLICIT NONE
00536 C*******
00537 C AUTEUR: Y.CHARTIER, DRPN
00538 C        FEVRIER 1991
00539 C
00540 C OBJET:  INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR
00541 C        D'UNE GRILLE GAUSSIENNE.
00542 C*******
00543       INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
00544       REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
00545       REAL(kind=ip_realwp_p) AY(J1:J2),CY(J1:J2,6)
00546       REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
00547 C
00548 C  NPTS   : NOMBRE DE POINTS A INTERPOLER
00549 C  I1:I2  : DIMENSION DE LA GRILLE SOURCE SELON X
00550 C  J1:J2  : DIMENSION DE LA GRILLE SOURCE SELON Y
00551 C  ZO     : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
00552 C  PX     : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
00553 C           VEUT INTERPOLER
00554 C  PY     : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
00555 C           VEUT INTERPOLER
00556 C  AY     : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES Y.
00557 C  CY     : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON Y.
00558 C  Z      : VALEURS DE LA GRILLE SOURCE.
00559 C
00560 C***************************************************************************
00561 C
00562 C  *   *   *   *
00563 C
00564 C  *   *   *   *
00565 C        #        ==>   PT (X,Y)
00566 C  *  (=)  *   *  ==> = PT (I, J)
00567 C
00568 C  *   *   *   *
00569 C
00570 C
00571 C
00572 C  CY(I,1) = 1.0 / (X2-X1)
00573 C  CY(I,2) = 1.0 / (X3-X1)
00574 C  CY(I,3) = 1.0 / (X3-X2)
00575 C  CY(I,4) = 1.0 / (X4-X1)
00576 C  CY(I,5) = 1.0 / (X4-X2)
00577 C  CY(I,6) = 1.0 / (X4-X3)
00578 C
00579 C  STRUCTURE IDENTIQUE POUR CY(J,1..6)
00580 
00581       INTEGER (kind=ip_intwp_p) I, J, M, N,STRIDE
00582       REAL(kind=ip_realwp_p) X, X1, X2, X3, X4
00583       REAL(kind=ip_realwp_p) B1, B2,  B3,  B4
00584       REAL(kind=ip_realwp_p) B11, B12, B13, B14
00585       REAL(kind=ip_realwp_p) Y,Y1,Y2,Y3,Y4
00586       REAL(kind=ip_realwp_p) Y11, Y12, Y13, Y14
00587       REAL(kind=ip_realwp_p) AY1, AY2, AY3, AY4
00588       REAL(kind=ip_realwp_p) FA, FA2, FA3, FA4
00589       REAL(kind=ip_realwp_p) A1,A2,A3,A4,C1,C2,C3,C4,C5,C6
00590 C  DEFINITION DES FONCTIONS IN-LINE
00591 
00592       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
00593       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
00594       PARAMETER (VOISIN  =   0)
00595       PARAMETER (LINEAIR =   1)
00596       PARAMETER (CUBIQUE =   3)
00597       PARAMETER (OUI   =   1)
00598       PARAMETER (MINIMUM =   2)
00599       PARAMETER (MAXIMUM =   3)
00600       PARAMETER (VALEUR  =   4)
00601       PARAMETER (ABORT   =  13)
00602       LOGICAL FLGXTRAP
00603       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
00604       REAL(kind=ip_realwp_p) VALXTRAP
00605       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
00606       REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
00607       REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
00608       CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)*0.1666666666666 + 0.5*(Z2-Z3)
00609      $)*DX + 0.5*(Z1+Z3)-Z2)*DX + Z3-0.1666666666666*Z4-0.5*Z2-0.
00610      $3333333333333*Z1)*DX+Z2
00611       ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
00612       FA(A1,A2,A3,A4,X,X1,X2,X3)=A1+(X-X1)*(A2+(X-X2)*(A3+A4*(X-X3))
00613      $)
00614       FA2(C1,A1,A2)=C1*(A2-A1)
00615       FA3(C1,C2,C3,A1,A2,A3)=C2*(C3*(A3-A2)-C1*(A2-A1))
00616       FA4(C1,C2,C3,C4,C5,C6,A1,A2,A3,A4)=C4*(C5*(C6*(A4-A3)-C3*(A3-
00617      $A2))   - C2*(C3*(A3-A2)-C1*(A2-A1)))
00618       STRIDE = 1
00619       IF( (ORDINT.EQ. CUBIQUE))THEN
00620          DO 23002 N=1,NPTS
00621             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
00622             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
00623             DX = PX(N) - I
00624             Y1=CUBIC(Z(I-1,J-1),Z(I  ,J-1),Z(I+1,J-1),Z(I+2,J-1),DX)
00625             Y2=CUBIC(Z(I-1,J  ),Z(I  ,J  ),Z(I+1,J  ),Z(I+2,J  ),DX)
00626             Y3=CUBIC(Z(I-1,J+1),Z(I  ,J+1),Z(I+1,J+1),Z(I+2,J+1),DX)
00627             Y4=CUBIC(Z(I-1,J+2),Z(I  ,J+2),Z(I+1,J+2),Z(I+2,J+2),DX)
00628             Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
00629 C           INTERPOLATION FINALE SELON Y
00630 
00631             AY1=AY(J-1)
00632             AY2=AY(J)
00633             AY3=AY(J+1)
00634             AY4=AY(J+2)
00635             Y11 = Y1
00636             Y12 = FA2(CY(J,1),Y1,Y2)
00637             Y13 = FA3(CY(J,1),CY(J,2),CY(J,3),Y1,Y2,Y3)
00638             Y14 = FA4(CY(J,1),CY(J,2),CY(J,3),CY(J,4),CY(J,5),CY(J,6
00639      $      ),Y1,Y2,Y3,Y4)
00640             ZO(N) = FA(Y11,Y12,Y13,Y14,Y,AY1,AY2,AY3)
00641 23002    CONTINUE 
00642       ENDIF 
00643       IF( (ORDINT.EQ. LINEAIR))THEN
00644          DO 23006 N=1,NPTS
00645             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
00646             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
00647             DX = PX(N) - I
00648             Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
00649             DY = (Y -AY(J))/(AY(J+1)-AY(J))
00650             Y2=ZLIN(Z(I,J  ),Z(I+1,J  ),DX)
00651             Y3=ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
00652             ZO(N)=ZLIN(Y2,Y3,DY)
00653 23006    CONTINUE 
00654       ENDIF 
00655       IF( (ORDINT.EQ. VOISIN))THEN
00656          DO 23010 N=1,NPTS
00657             I = MIN(I2,MAX(I1,NINT(PX(N))))
00658             J = MIN(J2,MAX(J1,NINT(PY(N))))
00659             ZO(N) = Z(I,J)
00660 23010    CONTINUE 
00661       ENDIF 
00662       RETURN
00663       END
00664 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00665 C
00666       SUBROUTINE GGLL2GD(X,Y,XLAT,XLON,NPTS,NI,NJ,HEM,LROOTS)
00667 C* -------------------------------------------------------------------
00668 C* S/R GGLL2GD - COMPUTES THE GRID CO-ORDINATES OF A POINT ON
00669 C                A GAUSSIAN GRID
00670 C* -------------------------------------------------------------------
00671       USE mod_kinds_oasis
00672 C* -------------------------------------------------------------------
00673       IMPLICIT NONE
00674 C* --------------------------------------------------------------------
00675       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
00676       PARAMETER (GLOBAL = 0)
00677       PARAMETER (NORD   = 1)
00678       PARAMETER (SUD   = 2)
00679       PARAMETER (SUDNORD= 0)
00680       PARAMETER (NORDSUD= 1)
00681       INTEGER (kind=ip_intwp_p) NPTS, NI, NJ
00682       REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), XLAT(NPTS), XLON(NPTS)
00683       REAL(kind=ip_realwp_p) LROOTS(NJ)
00684       REAL(kind=ip_realwp_p) ABS,DEL,EPSPHI
00685       INTEGER (kind=ip_intwp_p) HEM
00686       INTEGER (kind=ip_intwp_p) I,J
00687       INTEGER (kind=ip_intwp_p) GUESS
00688       INTEGER (kind=ip_intwp_p) NJHEM
00689       REAL(kind=ip_realwp_p) TMPLAT, DELLAT, DELLON, XLAT0, XLON0
00690       IF( (HEM.EQ. GLOBAL))THEN
00691          NJHEM = NJ / 2
00692       ELSE 
00693          NJHEM = NJ
00694       ENDIF 
00695       DELLON = 360.0 / REAL(NI,kind=ip_realwp_p)
00696       XLON0 = 0.0
00697       DELLAT = 90.0 / REAL(NJHEM,kind=ip_realwp_p)
00698       XLAT0 = DELLAT * 0.5
00699 CVDIR NOVECTOR
00700 
00701       DO 23002 I = 1, NPTS
00702          X(I) = (XLON(I) - XLON0)/DELLON + 1.0
00703          TMPLAT = XLAT(I)
00704          IF( (TMPLAT.LT. 0))THEN
00705             TMPLAT = -TMPLAT
00706          ENDIF 
00707          IF( (TMPLAT.GT. LROOTS(NJHEM)))THEN
00708             Y(I)=NJHEM+0.5*(TMPLAT-LROOTS(NJHEM))/(90.0-LROOTS(NJHEM
00709      $      ))
00710          ELSE 
00711             IF( (TMPLAT.LT. LROOTS(1)))THEN
00712                Y(I)=0.5+(0.5*TMPLAT/LROOTS(1))
00713             ELSE 
00714                GUESS=INT((TMPLAT-XLAT0)/DELLAT + 1.0)
00715 23010          IF((LROOTS(GUESS+1).LT.TMPLAT))THEN
00716                   GUESS = GUESS+1
00717                   GOTO 23010
00718                ENDIF 
00719 23012          IF((LROOTS(GUESS).GE.TMPLAT))THEN
00720                   GUESS = GUESS-1
00721                   GOTO 23012
00722                ENDIF 
00723                Y(I)=GUESS+(TMPLAT-LROOTS(GUESS))/(LROOTS(GUESS+1)-
00724      $            LROOTS(GUESS))
00725             ENDIF 
00726          ENDIF 
00727          IF( (HEM.EQ. GLOBAL))THEN
00728             IF( (0.EQ. MOD(NJ, 2)))THEN
00729                IF( (XLAT(I).GE. 0.0))THEN
00730                   Y(I) = Y(I) + NJHEM
00731                ELSE 
00732                   Y(I) = NJHEM - Y(I) + 1
00733                ENDIF 
00734             ELSE 
00735                IF( (XLAT(I).GE. 0.0))THEN
00736                   Y(I) = Y(I) + NJHEM + 1
00737                ELSE 
00738                   Y(I) = NJHEM - Y(I) + 1
00739                ENDIF 
00740             ENDIF 
00741          ENDIF 
00742          IF( (HEM.EQ. NORD))THEN
00743             IF( (XLAT(I).LT. 0.0))THEN
00744                Y(I) = -Y(I) + 1
00745             ENDIF 
00746          ENDIF 
00747          IF( (HEM.EQ. SUD))THEN
00748             IF( (XLAT(I).GE. 0.0))THEN
00749                Y(I) = Y(I) + NJ
00750             ELSE 
00751                Y(I) = NJ - Y(I) + 1.0
00752             ENDIF 
00753          ENDIF 
00754 23002 CONTINUE 
00755       RETURN
00756       END
00757 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00758 C
00759       BLOCK DATA QQQCOMBLK
00760       USE mod_kinds_oasis
00761       IMPLICIT NONE
00762 C* ----------------------------------------------------------------------
00763 C* COMMON QQQCOM1
00764 C
00765       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
00766       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
00767      $     AYTMP,CXTMP,CYTMP
00768       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
00769      $     SPOLPTS
00770       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
00771       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
00772       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
00773      $     NPOLMAX,SPOLMAX
00774       LOGICAL PTPOLN, PTPOLS
00775       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
00776      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
00777      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
00778       DATA XGDPTR, YGDPTR, ZPTR   /   0,   0,   0  /
00779       DATA AXPTR, AYPTR, CXPTR, CYPTR   /   0,   0,   0,   0 /
00780       DATA NPOLPTR,SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX   /   0,
00781      $   0,   0,   0,   0,   0  /
00782       DATA PTPOLN, PTPOLS / .FALSE., .FALSE. /
00783 C
00784 C-- COMMON QQQCOM2, QQQCOM3
00785 C-
00786 
00787       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
00788       INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
00789       INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
00790       COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
00791      $,LSTIG3,LSTIG4,LSTXOR
00792       CHARACTER*1 LSTGTYP
00793       COMMON /QQQCOM3/ LSTGTYP
00794       DATA LSTLILJ, LSTNI, LSTNJ, LSTIG1, LSTIG2, LSTIG3, LSTIG4   /
00795      $   0,   0,   0,   0,   0,   0,   0 /
00796       DATA LSTXOR, LSTGTYP   /   0,   ' ' /
00797       END
00798 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00799 C
00800       SUBROUTINE IGSCINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,GRREF,IG1
00801      $,IG2,IG3,IG4,SYM,AX,AY,VECFLG)
00802       USE mod_kinds_oasis
00803       IMPLICIT NONE
00804       INTEGER (kind=ip_intwp_p) LI,LJ,NI,NJ,IG1,IG2,IG3,IG4,BIDON
00805       REAL(kind=ip_realwp_p) ZO(LI,LJ),XLAT(LI,LJ),XLON(LI,LJ)
00806       REAL(kind=ip_realwp_p) ZI(NI,NJ),AX(NI),AY(NJ)
00807       LOGICAL SYM
00808       CHARACTER*1 GRTYP, GRREF
00809       LOGICAL QQQCMP
00810       EXTERNAL QQQCMP, RGOPTI
00811       EXTERNAL RGFREE, QQQCHK, QQQALOC, QQQKEEP, LL2RGD, LL2IGD
00812       EXTERNAL XPNCOF, XPNGD, NWTNCOF, IRGDINT, RGDINT
00813       EXTERNAL XPNAXEZ,XPNAXEG,GGGDINT,QQQGLAT,QQQXTRAP,POLRINT
00814       EXTERNAL XPNAXEY
00815 C
00816 CDEFINITION DES VARIABLES LOCALES
00817 C
00818 C-- COMMON QQQCOM1
00819 C--
00820 C
00821       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
00822       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
00823      $     AYTMP,CXTMP,CYTMP
00824       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
00825      $     SPOLPTS
00826       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
00827       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
00828       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
00829      $     NPOLMAX,SPOLMAX
00830       LOGICAL PTPOLN, PTPOLS
00831       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
00832      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
00833      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
00834 C----------------------
00835 C-- COMMON QQQCOM2
00836 C-
00837 
00838       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
00839       INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
00840       INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
00841       COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
00842      $,LSTIG3,LSTIG4,LSTXOR
00843 C-
00844 C------------
00845 C* COMMON QQQCOM3
00846 C
00847 
00848       CHARACTER*1 LSTGTYP
00849       COMMON /QQQCOM3/ LSTGTYP
00850 C------------
00851 
00852       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
00853       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
00854       PARAMETER (VOISIN  =   0)
00855       PARAMETER (LINEAIR =   1)
00856       PARAMETER (CUBIQUE =   3)
00857       PARAMETER (OUI   =   1)
00858       PARAMETER (MINIMUM =   2)
00859       PARAMETER (MAXIMUM =   3)
00860       PARAMETER (VALEUR  =   4)
00861       PARAMETER (ABORT   =  13)
00862       LOGICAL FLGXTRAP
00863       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
00864       REAL(kind=ip_realwp_p) VALXTRAP
00865       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
00866       EXTERNAL QQQCOMBLK
00867       LOGICAL OK
00868       INTEGER (kind=ip_intwp_p) XORSUM
00869       LOGICAL VECTEUR
00870 C
00871 C* Test if field is scalar or vector
00872 C* ---------------------------------
00873 C
00874       CHARACTER*6 VECFLG
00875       IF (VECFLG .EQ. 'SCALAR') THEN
00876           VECTEUR = .FALSE.
00877         ELSEIF (VECFLG .EQ. 'VECTOR') THEN
00878           VECTEUR = .TRUE.
00879         ELSE
00880           WRITE(6,*) ' WARNING !!! Wrong value for VECFLG '
00881           WRITE(6,*) ' WE STOP IN igscint                 '
00882           WRITE(6,*) ' '
00883           STOP
00884       ENDIF
00885 C*----------------------
00886 C
00887 CVERIFICATION DU TYPE DE GRILLE
00888 C
00889 
00890       IF( (GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y')
00891      $    .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
00892      $    GRREF.EQ. 'L')) THEN
00893          GOTO 10
00894       ELSE 
00895          GOTO 101
00896       ENDIF 
00897       ENTRY IGVINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,GRREF,IG1,IG2,
00898      $IG3,IG4,SYM,AX,AY)
00899       VECTEUR = .TRUE.
00900       IF( (GRTYP.EQ.'Z' .OR. GRTYP .EQ. 'Y')
00901      $    .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
00902      $    GRREF.EQ. 'L')) THEN
00903          GOTO 10
00904       ELSE 
00905          GOTO 101
00906       ENDIF 
00907 101   CONTINUE
00908       IF( GRTYP .NE. 'Z' .OR. GRTYP .NE. 'Y')THEN
00909          WRITE(6,*) 'Error, bad grid type (GRTYP)(IGSCINT)'
00910          WRITE(6,*) '''GRTYP'' should be equal to ''Z'' or ''Y'''
00911       ELSE 
00912          WRITE(6,*) 
00913      $    'Error, bad grid type (GRREF)(IGSCINT)'
00914          WRITE(6,*)
00915      $    '''GRREF'' should be equal to ''N'', ''S'' or ''L'''
00916       ENDIF 
00917       RETURN
00918       ENTRY RGVINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,IG4
00919      $,SYM)
00920       VECTEUR = .TRUE.
00921       GOTO 3
00922       ENTRY RGSCINT(ZO,LI,LJ,XLAT,XLON,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
00923      $IG4,SYM,VECFLG)
00924 C
00925 C* Test if field is scalar or vector
00926 C* ---------------------------------
00927 C
00928       IF (VECFLG .EQ. 'SCALAR') THEN
00929           VECTEUR = .FALSE.
00930         ELSEIF (VECFLG .EQ. 'VECTOR') THEN
00931           VECTEUR = .TRUE.
00932         ELSE
00933           WRITE(6,*) ' WARNING !!! Wrong value for VECFLG '
00934           WRITE(6,*) ' WE STOP IN igscint                 '
00935           WRITE(6,*) ' '
00936           STOP
00937       ENDIF
00938 C*
00939 3     CONTINUE
00940       IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'L'.OR.GRTYP.EQ.
00941      $'N'.OR.GRTYP.EQ.'S'   .OR.GRTYP.EQ.'G'))THEN
00942          GOTO 10
00943       ENDIF 
00944       WRITE (6,200)
00945 200   FORMAT('0','Error bad grid type (GRTYP)(RGSCINT)')
00946       RETURN
00947 10    CALL RGOPTI('INTERP', BIDON, .FALSE.)
00948       OK = QQQCMP(XLAT,XLON,LI*LJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM
00949      $)
00950       IF( (.NOT.OK))THEN
00951          CALL XPNCOF(I1,I2,J1,J2,NI,NJ,GRTYP,GRREF,IG1,IG1,IG3,IG4,
00952      $   SYM,AX,AY)
00953          CALL RGFREE
00954          CALL QQQALOC(XLAT, XLON, LI, LJ, I1, I2, J1, J2)
00955          CALL QQQKEEP(LI*LJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM)
00956          CALL QQQCHK(XLAT, XLON, LI*LJ)
00957          IF( GRTYP.EQ. 'Z' .OR. GRTYP.EQ. 'Y')THEN
00958          IF( GRTYP.EQ. 'Z') CALL XPNAXEZ(AXTMP(AXPTR),AYTMP(AYPTR),
00959      $           I1,I2,J1,J2,AX,AY,NI,NJ)
00960          IF( GRTYP.EQ. 'Y') CALL XPNAXEY(AXTMP(AXPTR),AYTMP(AYPTR),
00961      $           I1,I2,J1,J2,AX,AY,NI,NJ)
00962             CALL NWTNCOF(CXTMP(CXPTR),CYTMP(CYPTR),AXTMP(AXPTR),
00963      $      AYTMP(AYPTR),I1,I2,J1,J2)
00964             CALL LL2IGD(XGD(XGDPTR),YGD(YGDPTR),XLAT,XLON,LI*LJ,NI,
00965      $      NJ,GRTYP,GRREF,IG1,IG2,IG3,IG4,SYM,AX,AY)
00966          ELSE 
00967             IF( (GRTYP.EQ. 'G'))THEN
00968                CALL QQQGLAT(NJ,IG1)
00969                CALL XPNAXEG(AXTMP(AXPTR),AYTMP(AYPTR),I1,I2,J1,J2,NI
00970      $         ,NJ,IG1)
00971                CALL NWTNCOF(CXTMP(CXPTR),CYTMP(CYPTR),AXTMP(AXPTR),
00972      $         AYTMP(AYPTR),I1,I2,J1,J2)
00973             ENDIF 
00974             CALL LL2RGD(XGD(XGDPTR),YGD(YGDPTR),XLAT,XLON,LI*LJ,NI,
00975      $      NJ,GRTYP,IG1,IG2,IG3,IG4,SYM)
00976          ENDIF 
00977       ENDIF 
00978       CALL XPNGD(ZTMP(ZPTR),I1,I2,J1,J2,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
00979      $IG4,SYM,VECTEUR)
00980       IF( (GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
00981          CALL IRGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,AXTMP(AXPTR),
00982      $   AYTMP(AYPTR),CXTMP(CXPTR),CYTMP(CYPTR),ZTMP(ZPTR),I1,I2,J1,
00983      $   J2)
00984       ENDIF 
00985       IF( (GRTYP.EQ. 'G'))THEN
00986          CALL GGGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,AYTMP(AYPTR),
00987      $   CYTMP(CYPTR),ZTMP(ZPTR),I1,I2,J1,J2)
00988       ENDIF 
00989       IF( (GRTYP.NE. 'G'.AND. GRTYP.NE. 'Z' 
00990      $    .AND. GRTYP .NE. 'Y'))THEN
00991          CALL RGDINT(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,ZTMP(ZPTR),I1,
00992      $   I2,J1,J2)
00993       ENDIF 
00994       IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
00995      $ GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
00996          CALL QQQXTRAP(ZO,XGD(XGDPTR),YGD(YGDPTR),LI*LJ,ZTMP(ZPTR),
00997      $   I1,I2,J1,J2)
00998       ENDIF 
00999       IF( (PTPOLN.OR. PTPOLS))THEN
01000          CALL POLRINT(ZO,XLAT,XLON,LI*LJ,ZTMP(ZPTR),NI,NJ,I1,I2,J1,
01001      $   J2,GRTYP,IG1,IG2,IG3,IG4,VECTEUR)
01002       ENDIF 
01003       RETURN
01004       END
01005 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01006 C
01007       SUBROUTINE IGUVINT(SPDO,PSIO,LI,LJ,XLAT,XLON,UI,VI,NI,NJ,GRTYP
01008      $,GRREF,IG1,IG2,IG3,IG4,SWS,AX,AY)
01009       USE mod_kinds_oasis
01010       IMPLICIT NONE
01011       INTEGER (kind=ip_intwp_p) LI,LJ,NI,NJ,IG1,IG2,IG3,IG4,BIDON
01012       REAL(kind=ip_realwp_p) SPDO(LI,LJ),PSIO(LI,LJ),XLAT(LI,LJ),
01013      $     XLON(LI,LJ)
01014       REAL(kind=ip_realwp_p) UI(NI,NJ),VI(NI,NJ),AX(NI),AY(NJ)
01015       LOGICAL SWS
01016       CHARACTER*1 GRTYP, GRREF
01017       EXTERNAL GDW2LLW,IGSCINT,RGVINT,MODULE,RGOPTI,IGVINT
01018       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
01019       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
01020       PARAMETER (VOISIN  =   0)
01021       PARAMETER (LINEAIR =   1)
01022       PARAMETER (CUBIQUE =   3)
01023       PARAMETER (OUI   =   1)
01024       PARAMETER (MINIMUM =   2)
01025       PARAMETER (MAXIMUM =   3)
01026       PARAMETER (VALEUR  =   4)
01027       PARAMETER (ABORT   =  13)
01028       LOGICAL FLGXTRAP
01029       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
01030       REAL(kind=ip_realwp_p) VALXTRAP
01031       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
01032 C* -------------------------------------------------------------------
01033 C
01034 CVERIFICATION DU TYPE DE GRILLE
01035 C
01036 
01037       CALL RGOPTI('INTERP', BIDON, .FALSE.)
01038       IF( (GRTYP.EQ.'Z' .OR. GRTYP .EQ. 'Y')
01039      $    .AND. (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR.
01040      $ GRREF.EQ. 'L')) THEN
01041          CALL IGVINT(SPDO,LI,LJ,XLAT,XLON,UI,NI,NJ,GRTYP,GRREF,IG1,
01042      $   IG2,IG3,IG4,.TRUE.,AX,AY)
01043          CALL IGVINT(PSIO,LI,LJ,XLAT,XLON,VI,NI,NJ,GRTYP,GRREF,IG1,
01044      $   IG2,IG3,IG4,.FALSE.,AX,AY)
01045          IF( (SWS))THEN
01046             CALL GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRREF,IG1,IG2,IG3,IG4)
01047          ELSE 
01048             CALL MODULE(SPDO,PSIO,LI,LJ)
01049          ENDIF 
01050       ELSE 
01051          IF( (GRTYP.NE. 'Z' .AND. GRTYP .NE. 'Y'))THEN
01052             WRITE(6,*)
01053      $       'Error, Bad grid type (GRTYP)(IGSCINT)'
01054             WRITE(6,*) 
01055      $          '''GRTYP'' should be equal to ''Z'' or ''Y'''
01056          ELSE 
01057             WRITE(6,*)
01058      $ 'Error, Bad reference grid type(GRREF)(IGSCINT)'
01059             WRITE(6,*)
01060      $       '''GRREF'' should be equal to ''N'', ''S'' or ''L'''
01061          ENDIF 
01062       ENDIF 
01063       RETURN
01064       ENTRY RGUVINT(SPDO,PSIO,LI,LJ,XLAT,XLON,UI,VI,NI,NJ,GRTYP,IG1,
01065      $IG2,IG3,IG4,SWS)
01066       CALL RGOPTI('INTERP', BIDON, .FALSE.)
01067       IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'L'.OR.GRTYP.EQ.
01068      $'N'.OR.GRTYP.EQ.'S'   .OR.GRTYP.EQ.'G'))THEN
01069          CALL RGVINT(SPDO,LI,LJ,XLAT,XLON,UI,NI,NJ,GRTYP,IG1,IG2,IG3
01070      $   ,IG4,.TRUE.)
01071          CALL RGVINT(PSIO,LI,LJ,XLAT,XLON,VI,NI,NJ,GRTYP,IG1,IG2,IG3
01072      $   ,IG4,.FALSE.)
01073          IF( (SWS))THEN
01074             CALL GDW2LLW(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
01075          ELSE 
01076             CALL MODULE(SPDO,PSIO,LI,LJ)
01077          ENDIF 
01078       ELSE 
01079          WRITE (6,200)
01080 200      FORMAT('0','Error bad grid type (GRTYP)(RGSCINT)')
01081       ENDIF 
01082       RETURN
01083       END
01084 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01085 C
01086       SUBROUTINE IRGDINT(ZO,PX,PY,NPTS,AX,AY,CX,CY,Z,I1,I2,J1,J2)
01087       USE mod_kinds_oasis
01088       IMPLICIT NONE
01089 C*******
01090 C AUTEUR: Y.CHARTIER, DRPN
01091 C        FEVRIER 1991
01092 C
01093 C OBJET:  INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR
01094 C        D'UNE GRILLE SOURCE IRREGULIERE.
01095 C*******
01096       INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
01097       REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
01098       REAL(kind=ip_realwp_p) FA, FA2, FA3, FA4
01099       REAL(kind=ip_realwp_p) AX(I1:I2),AY(J1:J2),CX(I1:I2,6)
01100       REAL(kind=ip_realwp_p) CY(J1:J2,6)
01101       REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
01102 C
01103 C  NPTS   : NOMBRE DE POINTS A INTERPOLER
01104 C  I1:I2  : DIMENSION DE LA GRILLE SOURCE SELON X
01105 C  J1:J2  : DIMENSION DE LA GRILLE SOURCE SELON Y
01106 C  ZO     : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
01107 C  PX     : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
01108 C           VEUT INTERPOLER
01109 C  PY     : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
01110 C           VEUT INTERPOLER
01111 C  AX     : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES X.
01112 C  AY     : VECTEUR CONTENANT LA POS. DES POINTS SUR L'AXE DES Y.
01113 C  CX     : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON X.
01114 C  CY     : VECTEUR CONTENANT UNE TABLE DE DIFFERENCES SELON Y.
01115 C  Z      : VALEURS DE LA GRILLE SOURCE.
01116 C
01117 C***************************************************************************
01118 C
01119 C  *   *   *   *
01120 C
01121 C  *   *   *   *
01122 C        #        ==>   PT (X,Y)
01123 C  *  (=)  *   *  ==> = PT (I, J)
01124 C
01125 C  *   *   *   *
01126 C
01127 C
01128 C
01129 C  CX(I,1) = 1.0 / (X2-X1)
01130 C  CX(I,2) = 1.0 / (X3-X1)
01131 C  CX(I,3) = 1.0 / (X3-X2)
01132 C  CX(I,4) = 1.0 / (X4-X1)
01133 C  CX(I,5) = 1.0 / (X4-X2)
01134 C  CX(I,6) = 1.0 / (X4-X3)
01135 C
01136 C  STRUCTURE IDENTIQUE POUR CY(J,1..6)
01137 
01138       REAL(kind=ip_realwp_p) A11,A12,A13,A14,A21,A22,A23,A24
01139       REAL(kind=ip_realwp_p) A31,A32,A33,A34,A41,A42,A43,A44
01140       REAL(kind=ip_realwp_p) B1,B2,B3,B4,B11,B12,B13,B14
01141       REAL(kind=ip_realwp_p) X1,X2,X3,X4,Y1,Y2,Y3,Y4
01142       INTEGER (kind=ip_intwp_p) I, J, N
01143       REAL(kind=ip_realwp_p) A1,A2,A3,A4,X,Y,C1,C2,C3,C4,C5,C6
01144 C  DEFINITION DES FONCTIONS IN-LINE
01145 
01146       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
01147       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
01148       PARAMETER (VOISIN  =   0)
01149       PARAMETER (LINEAIR =   1)
01150       PARAMETER (CUBIQUE =   3)
01151       PARAMETER (OUI   =   1)
01152       PARAMETER (MINIMUM =   2)
01153       PARAMETER (MAXIMUM =   3)
01154       PARAMETER (VALEUR  =   4)
01155       PARAMETER (ABORT   =  13)
01156       LOGICAL FLGXTRAP
01157       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
01158       REAL(kind=ip_realwp_p) VALXTRAP
01159       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
01160       REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
01161       REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
01162       CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)*0.1666666666666 + 0.5*(Z2-Z3)
01163      $)*DX + 0.5*(Z1+Z3)-Z2)*DX + Z3-0.1666666666666*Z4-0.5*Z2-0.
01164      $3333333333333*Z1)*DX+Z2
01165       ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
01166       FA(A1,A2,A3,A4,X,X1,X2,X3)=A1+(X-X1)*(A2+(X-X2)*(A3+A4*(X-X3))
01167      $)
01168       FA2(C1,A1,A2)=C1*(A2-A1)
01169       FA3(C1,C2,C3,A1,A2,A3)=C2*(C3*(A3-A2)-C1*(A2-A1))
01170       FA4(C1,C2,C3,C4,C5,C6,A1,A2,A3,A4)=C4*(C5*(C6*(A4-A3)-C3*(A3-
01171      $A2))   - C2*(C3*(A3-A2)-C1*(A2-A1)))
01172       IF( (ORDINT.EQ. CUBIQUE))THEN
01173          DO 23002 N=1,NPTS
01174             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
01175             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
01176             X = AX(I) + (AX(I+1)-AX(I))*(PX(N)-I)
01177             Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
01178             X1=AX(I-1)
01179             X2=AX(I)
01180             X3=AX(I+1)
01181             X4=AX(I+2)
01182             Y1=AY(J-1)
01183             Y2=AY(J)
01184             Y3=AY(J+1)
01185             Y4=AY(J+2)
01186 C        INTERPOLATION 1ERE RANGEE SELON X
01187 
01188             Z1=Z(I-1,J-1)
01189             Z2=Z(I  ,J-1)
01190             Z3=Z(I+1,J-1)
01191             Z4=Z(I+2,J-1)
01192             A11 = Z1
01193             A12 = FA2(CX(I,1),Z1,Z2)
01194             A13 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
01195             A14 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
01196      $      ),Z1,Z2,Z3,Z4)
01197             B1  = FA(A11,A12,A13,A14,X,X1,X2,X3)
01198 C        INTERPOLATION 2EME RANGEE SELON X
01199 
01200             Z1=Z(I-1,J)
01201             Z2=Z(I  ,J)
01202             Z3=Z(I+1,J)
01203             Z4=Z(I+2,J)
01204             A21 = Z1
01205             A22 = FA2(CX(I,1),Z1,Z2)
01206             A23 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
01207             A24 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
01208      $      ),Z1,Z2,Z3,Z4)
01209             B2  = FA(A21,A22,A23,A24,X,X1,X2,X3)
01210 C     INTERPOLATION 3EME RANGEE SELON X
01211 
01212             Z1=Z(I-1,J+1)
01213             Z2=Z(I  ,J+1)
01214             Z3=Z(I+1,J+1)
01215             Z4=Z(I+2,J+1)
01216             A31 = Z1
01217             A32 = FA2(CX(I,1),Z1,Z2)
01218             A33 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
01219             A34 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
01220      $      ),Z1,Z2,Z3,Z4)
01221             B3  = FA(A31,A32,A33,A34,X,X1,X2,X3)
01222 C        INTERPOLATION 4EME RANGEE SELON X
01223 
01224             Z1=Z(I-1,J+2)
01225             Z2=Z(I  ,J+2)
01226             Z3=Z(I+1,J+2)
01227             Z4=Z(I+2,J+2)
01228             A41 = Z1
01229             A42 = FA2(CX(I,1),Z1,Z2)
01230             A43 = FA3(CX(I,1),CX(I,2),CX(I,3),Z1,Z2,Z3)
01231             A44 = FA4(CX(I,1),CX(I,2),CX(I,3),CX(I,4),CX(I,5),CX(I,6
01232      $      ),Z1,Z2,Z3,Z4)
01233             B4  = FA(A41,A42,A43,A44,X,X1,X2,X3)
01234 C        INTERPOLATION FINALE SELON Y
01235 
01236             B11 = B1
01237             B12 = FA2(CY(J,1),B1,B2)
01238             B13 = FA3(CY(J,1),CY(J,2),CY(J,3),B1,B2,B3)
01239             B14 = FA4(CY(J,1),CY(J,2),CY(J,3),CY(J,4),CY(J,5),CY(J,6
01240      $      ),B1,B2,B3,B4)
01241             ZO(N) = FA(B11,B12,B13,B14,Y,Y1,Y2,Y3)
01242 23002    CONTINUE 
01243       ENDIF 
01244       IF( (ORDINT.EQ. LINEAIR))THEN
01245          DO 23006 N=1,NPTS
01246             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
01247             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
01248             X = AX(I) + (AX(I+1)-AX(I))*(PX(N)-I)
01249             Y = AY(J) + (AY(J+1)-AY(J))*(PY(N)-J)
01250             DX = (X - AX(I))/(AX(I+1)-AX(I))
01251             DY = (Y - AY(J))/(AY(J+1)-AY(J))
01252             Y1 = ZLIN(Z(I,J),Z(I+1,J),DX)
01253             Y2 = ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
01254             ZO(N) = ZLIN(Y1,Y2,DY)
01255 23006    CONTINUE 
01256       ENDIF 
01257       IF( (ORDINT.EQ. VOISIN))THEN
01258          DO 23010 N=1,NPTS
01259             I = MIN(I2,MAX(I1,NINT(PX(N))))
01260             J = MIN(J2,MAX(J1,NINT(PY(N))))
01261             ZO(N) = Z(I,J)
01262 23010    CONTINUE 
01263       ENDIF 
01264       RETURN
01265       END
01266 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01267 C
01268       SUBROUTINE LL2IGD(PX,PY,XLAT,XLON,NPTS,NI,NJ,GRTYP,GRREF,IG1,
01269      $IG2,IG3,IG4,SYM,AX,AY)
01270 C* ----------------------------------------------------------------------
01271 C**S/R LL2IGD - CONVERSION DE COORDONNEES LAT-LON A PTS DE GRILLE
01272 C* ----------------------------------------------------------------------
01273       USE mod_kinds_oasis
01274       IMPLICIT NONE
01275       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
01276       PARAMETER (GLOBAL = 0)
01277       PARAMETER (NORD   = 1)
01278       PARAMETER (SUD   = 2)
01279       PARAMETER (SUDNORD= 0)
01280       PARAMETER (NORDSUD= 1)
01281       EXTERNAL CIGAXG,VXYFLL,LLLL2GD,PERMUT,CHKXTRAP
01282       INTEGER (kind=ip_intwp_p) I, NPTS, NI, NJ
01283       REAL(kind=ip_realwp_p) PX(NPTS),PY(NPTS),XLAT(NPTS),XLON(NPTS)
01284       CHARACTER*1 GRTYP, GRREF
01285       INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
01286       LOGICAL SYM
01287       REAL(kind=ip_realwp_p) AX(NI), AY(NJ)
01288       REAL(kind=ip_realwp_p) PI,PJ,DGRW,D60,TMP
01289       REAL(kind=ip_realwp_p) DLAT, DLON, XLAT0, XLON0
01290       INTEGER (kind=ip_intwp_p) INDX, INDY
01291       INTEGER (kind=ip_intwp_p) CHERCHE, FINDLON
01292       EXTERNAL CHERCHE, FINDLON 
01293       IF( (GRREF.EQ. 'N'.OR. GRREF.EQ. 'S'.OR. GRREF.EQ. 'L')
01294      $)THEN
01295          IF( (GRREF.EQ. 'N'))THEN
01296             CALL CIGAXG(GRREF,  PI, PJ, D60, DGRW, IG1, IG2, IG3,
01297      $       IG4)
01298             CALL VXYFLL(PX, PY, XLAT, XLON, NPTS, D60,DGRW,PI,PJ,1)
01299          ELSE 
01300             IF( (GRREF.EQ. 'S'))THEN
01301                CALL CIGAXG(GRREF,  PI, PJ, D60, DGRW, IG1, IG2, IG3,
01302      $          IG4)
01303                CALL VXYFLL(PX, PY, XLAT, XLON, NPTS, D60,DGRW,PI,PJ,
01304      $         2)
01305             ELSE 
01306                CALL CIGAXG(GRREF, XLAT0, XLON0, DLAT, DLON, IG1, IG2
01307      $         , IG3, IG4)
01308                CALL LLLL2GD(PX, PY, XLAT, XLON, NPTS, XLAT0, XLON0,
01309      $          DLAT, DLON)
01310             ENDIF 
01311          ENDIF 
01312          DO 23006 I=1,NPTS
01313             INDX = FINDLON(PX(I), AX, NI, TMP)
01314             INDY = CHERCHE(PY(I), AY, NJ)
01315             IF( (INDX.GE. NI))THEN
01316                INDX = NI - 1
01317             ENDIF 
01318             IF( (INDY.GE. NJ))THEN
01319                INDY = NJ - 1
01320             ENDIF 
01321             PX(I)=REAL(INDX,kind=ip_realwp_p)+
01322      $          (TMP-AX(INDX))/(AX(INDX+1)-AX(INDX))
01323             PY(I) = REAL(INDY,kind=ip_realwp_p) 
01324      $          + (PY(I) - AY(INDY))/(AY(INDY+1)-AY(INDY))
01325 23006    CONTINUE 
01326          CALL CHKXTRAP(PX, PY, NPTS, NI, NJ)
01327       ENDIF 
01328       RETURN
01329       END
01330 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01331 C
01332       SUBROUTINE LL2RGD(PX,PY,XLAT,XLON,NPTS,NI,NJ,GRTYP,IG1,IG2,IG3
01333      $,IG4,SYM)
01334 C* -----------------------------------------------------------------
01335 C* S/R LL2RGD - CONVERSION DE COORDONNEES LAT-LON A PTS DE GRILLE
01336 C* ----------------------------------------------------------------------
01337       USE mod_kinds_oasis
01338       IMPLICIT NONE
01339       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
01340       PARAMETER (GLOBAL = 0)
01341       PARAMETER (NORD   = 1)
01342       PARAMETER (SUD   = 2)
01343       PARAMETER (SUDNORD= 0)
01344       PARAMETER (NORDSUD= 1)
01345 C
01346 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
01347 C
01348 C RDTODG = 180/PIE, DGTORD = PIE/180
01349 
01350       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
01351       DATA PIE   /3.14159265358979323846/
01352       DATA RDTODG /57.295779513082/
01353       DATA DGTORD /1.7453292519943E-2/
01354 C
01355 C-- COMMON GAUSSGD
01356 C      REAL(kind=ip_realwp_p) ROOTS(1),LROOTS(1)
01357       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
01358       INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
01359       COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
01360 C----------------------
01361 
01362       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
01363       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
01364       PARAMETER (VOISIN  =   0)
01365       PARAMETER (LINEAIR =   1)
01366       PARAMETER (CUBIQUE =   3)
01367       PARAMETER (OUI   =   1)
01368       PARAMETER (MINIMUM =   2)
01369       PARAMETER (MAXIMUM =   3)
01370       PARAMETER (VALEUR  =   4)
01371       PARAMETER (ABORT   =  13)
01372       LOGICAL FLGXTRAP
01373       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
01374       REAL(kind=ip_realwp_p) VALXTRAP
01375       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
01376       EXTERNAL CIGAXG,VXYFLL,LLLL2GD,PERMUT
01377       EXTERNAL QQQGLAT,GGLL2GD,CHKXTRAP
01378       INTEGER (kind=ip_intwp_p) NPTS, NI, NJ
01379       REAL(kind=ip_realwp_p) PX(NPTS),PY(NPTS),XLAT(NPTS),XLON(NPTS)
01380       CHARACTER*1 GRTYP
01381       INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
01382       LOGICAL SYM
01383       REAL(kind=ip_realwp_p) PI,PJ,DGRW,D60
01384       REAL(kind=ip_realwp_p) DELLAT, DELLON, XLAT0, XLON0
01385       INTEGER (kind=ip_intwp_p) I,J
01386       IF( (GRTYP.EQ.'N'))THEN
01387          CALL CIGAXG(GRTYP,PI,PJ,D60,DGRW,IG1,IG2,IG3,IG4)
01388          CALL VXYFLL(PX,PY,XLAT,XLON,NPTS,D60,DGRW,PI,PJ,NORD)
01389          CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
01390       ENDIF 
01391       IF( (GRTYP.EQ.'S'))THEN
01392          CALL CIGAXG(GRTYP,PI,PJ,D60,DGRW,IG1,IG2,IG3,IG4)
01393          CALL VXYFLL(PX,PY,XLAT,XLON,NPTS,D60,DGRW,PI,PJ,SUD)
01394          CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
01395       ENDIF 
01396       IF( (GRTYP.EQ.'A'))THEN
01397          DELLON = 360.0 / REAL(NI,ip_realwp_p)
01398          XLON0 = 0.0
01399          IF( (IG1.EQ. GLOBAL))THEN
01400             DELLAT = 180.0 / REAL(NJ,ip_realwp_p)
01401             XLAT0 = -90.0 + DELLAT * 0.5
01402          ENDIF 
01403          IF( (IG1.EQ. NORD))THEN
01404             DELLAT = 90.0 / REAL(NJ,ip_realwp_p)
01405             XLAT0 =  DELLAT * 0.5
01406          ENDIF 
01407          IF( (IG1.EQ. SUD))THEN
01408             DELLAT = 90.0 / REAL(NJ,ip_realwp_p)
01409             XLAT0 = -90.0 + DELLAT * 0.5
01410          ENDIF 
01411          FLGXTRAP = .FALSE.
01412          CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
01413      $    DELLON)
01414       ENDIF 
01415       IF( (GRTYP.EQ.'B'))THEN
01416          DELLON = 360.0 / REAL(NI-1,ip_realwp_p)
01417          XLON0 = 0.0
01418          IF( (IG1.EQ. GLOBAL))THEN
01419             DELLAT = 180.0 / REAL(NJ-1,ip_realwp_p)
01420             XLAT0 = -90.0
01421          ENDIF 
01422          IF( (IG1.EQ. NORD))THEN
01423             DELLAT = 90.0 / REAL(NJ-1,ip_realwp_p)
01424             XLAT0 =  0.0
01425          ENDIF 
01426          IF( (IG1.EQ. SUD))THEN
01427             DELLAT = 90.0 / REAL(NJ-1,ip_realwp_p)
01428             XLAT0 = -90.0
01429          ENDIF 
01430          FLGXTRAP = .FALSE.
01431          CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
01432      $    DELLON)
01433       ENDIF 
01434       IF( (GRTYP.EQ. 'G'))THEN
01435          FLGXTRAP = .FALSE.
01436          CALL GGLL2GD(PX,PY,XLAT,XLON,NPTS,NI,NJ,IG1,LROOTS(ILROOTS)
01437      $   )
01438       ENDIF 
01439       IF( (GRTYP.EQ. 'L'))THEN
01440          CALL CIGAXG(GRTYP,XLAT0,XLON0,DELLAT,DELLON,IG1,IG2,IG3,IG4
01441      $   )
01442          DO 23024 I=1,NPTS
01443             IF( (XLON(I).LT. XLON0))THEN
01444                XLON(I) = XLON(I) + 360.0
01445             ENDIF 
01446             IF( (XLON(I).GT. (XLON0 + NI*DELLON)))THEN
01447                XLON(I) = XLON(I) - 360.0
01448             ENDIF 
01449 23024    CONTINUE 
01450          CALL LLLL2GD(PX,PY,XLAT,XLON,NPTS,XLAT0, XLON0, DELLAT,
01451      $    DELLON)
01452          CALL CHKXTRAP(PX,PY,NPTS,NI,NJ)
01453       ENDIF 
01454       RETURN
01455       END
01456 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01457 C
01458       SUBROUTINE LLLL2GD(X,Y,DLAT,DLON,NPTS,XLAT0,XLON0, DELLAT,
01459      $ DELLON)
01460 C* ----------------------------------------------------------------------
01461 C* S/R LLLL2GD - COMPUTES THE GRID CO-ORDINATES OF A POINT
01462 C* ----------------------------------------------------------------------
01463       USE mod_kinds_oasis
01464       IMPLICIT NONE
01465       INTEGER (kind=ip_intwp_p) NPTS
01466       REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), DLAT(NPTS), DLON(NPTS)
01467       REAL(kind=ip_realwp_p) XLAT0, XLON0, DELLAT, DELLON
01468       INTEGER (kind=ip_intwp_p) I
01469       DO 23000 I=1,NPTS
01470          X(I) = (DLON(I) - XLON0)/DELLON + 1.0
01471          Y(I) = (DLAT(I) - XLAT0)/DELLAT + 1.0
01472 23000 CONTINUE 
01473       RETURN
01474       END
01475 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01476 C
01477       SUBROUTINE MODULE(U,V,LI,LJ)
01478       USE mod_kinds_oasis
01479       INTEGER (kind=ip_intwp_p) LI,LJ
01480       REAL(kind=ip_realwp_p) U(LI,LJ), V(LI,LJ)
01481       INTEGER (kind=ip_intwp_p) I,J
01482       DO 23000 I=1,LI
01483          DO 23002 J=1,LJ
01484             U(I,J)=SQRT(U(I,J)*U(I,J)+V(I,J)*V(I,J))
01485 23002    CONTINUE 
01486 23000 CONTINUE 
01487       RETURN
01488       END
01489 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01490 C
01491       SUBROUTINE NWTNCOF(CX,CY,AX,AY,I1,I2,J1,J2)
01492 C*******
01493 C AUTEUR: Y. CHARTIER, DRPN
01494 C        FEVRIER 1991
01495 C
01496 C OBJET:  CALCUL DE COEFFICIENTS UTI1ISES DANS LA FORME NEWTONIENNE
01497 C        DE L'INTERPOLATION DE LAGRANGE.
01498 C
01499 C===========================================
01500 C
01501 C   -----*-------------*------#------*----------*------->
01502 C        X1            X2     X      X3         X4
01503 C
01504 C===========================================
01505 C     CX(I,1) = 1.0 / (X2-X1)
01506 C     CX(I,2) = 1.0 / (X3-X1)
01507 C     CX(I,3) = 1.0 / (X3-X2)
01508 C     CX(I,4) = 1.0 / (X4-X1)
01509 C     CX(I,5) = 1.0 / (X4-X2)
01510 C     CX(I,6) = 1.0 / (X4-X3)
01511 C
01512 C     STRUCTURE IDENTIQUE POUR CY(J,1..6)
01513 C*******
01514       USE mod_kinds_oasis
01515       IMPLICIT NONE
01516       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
01517       REAL(kind=ip_realwp_p) CX(I1:I2,6),CY(J1:J2,6),AX(I1:I2),AY(J1:J2)
01518       INTEGER (kind=ip_intwp_p) I,J
01519       DO 23000 I=I1+1,I2-2
01520          CX(I,1) = 1. / (AX(I  ) - AX(I-1))
01521          CX(I,2) = 1. / (AX(I+1) - AX(I-1))
01522          CX(I,3) = 1. / (AX(I+1) - AX(I  ))
01523          CX(I,4) = 1. / (AX(I+2) - AX(I-1))
01524          CX(I,5) = 1. / (AX(I+2) - AX(I  ))
01525          CX(I,6) = 1. / (AX(I+2) - AX(I+1))
01526 23000 CONTINUE 
01527       DO 23002 J=J1+1,J2-2
01528          CY(J,1) = 1. / (AY(J  ) - AY(J-1))
01529          CY(J,2) = 1. / (AY(J+1) - AY(J-1))
01530          CY(J,3) = 1. / (AY(J+1) - AY(J  ))
01531          CY(J,4) = 1. / (AY(J+2) - AY(J-1))
01532          CY(J,5) = 1. / (AY(J+2) - AY(J  ))
01533          CY(J,6) = 1. / (AY(J+2) - AY(J+1))
01534 23002 CONTINUE 
01535       RETURN
01536       END
01537 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01538 C
01539       SUBROUTINE POLRINT(ZO,XLAT,XLON,LILJ,ZI,NI,NJ,I1,I2,J1,J2,
01540      $GRTYP,IG1,IG2,IG3,IG4,VECTEUR)
01541       USE mod_kinds_oasis
01542       IMPLICIT NONE
01543       INTEGER (kind=ip_intwp_p) LILJ,NI,NJ,IG1,IG2,IG3,IG4
01544       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
01545       REAL(kind=ip_realwp_p) ZO(LILJ),XLAT(LILJ),XLON(LILJ)
01546       REAL(kind=ip_realwp_p) ZI(I1:I2,J1:J2)
01547       CHARACTER*1 GRTYP
01548       LOGICAL VECTEUR
01549       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
01550       PARAMETER (GLOBAL = 0)
01551       PARAMETER (NORD   = 1)
01552       PARAMETER (SUD   = 2)
01553       PARAMETER (SUDNORD= 0)
01554       PARAMETER (NORDSUD= 1)
01555 C-- COMMON QQQCOM1
01556 C--
01557       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
01558       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
01559      $     AYTMP,CXTMP,CYTMP
01560       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
01561      $     SPOLPTS
01562       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
01563       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
01564       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
01565      $     NPOLMAX,SPOLMAX
01566       LOGICAL PTPOLN, PTPOLS
01567       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
01568      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
01569      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
01570 C----------------------
01571 C
01572 CDEFINITION DES VARIABLES LOCALES
01573 C
01574 
01575       INTEGER (kind=ip_intwp_p) I,J
01576       INTEGER (kind=ip_intwp_p) N1, N2, S1, S2
01577       REAL(kind=ip_realwp_p) VPOLNOR, VPOLSUD
01578       IF( (VECTEUR))THEN
01579          RETURN
01580       ENDIF 
01581       IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
01582      $ GRTYP.EQ. 'Z' .OR. GRTYP .EQ. 'Y'))THEN
01583          RETURN
01584       ENDIF 
01585       IF( (GRTYP.EQ. 'B'))THEN
01586          IF( (IG1.EQ. GLOBAL))THEN
01587             VPOLNOR = ZI(0, NJ)
01588             VPOLSUD = ZI(0, 1)
01589          ENDIF 
01590          IF( (IG1.EQ. NORD))THEN
01591             VPOLNOR = ZI(0, NJ)
01592             VPOLSUD = ZI(0, -NJ+1)
01593          ENDIF 
01594          IF( (IG1.EQ. SUD))THEN
01595             VPOLNOR = ZI(0, 2*NJ-1)
01596             VPOLSUD = ZI(0, 1)
01597          ENDIF 
01598       ENDIF 
01599       IF( (GRTYP.EQ. 'A'.OR. GRTYP.EQ. 'G'))THEN
01600          VPOLNOR = 0.0
01601          VPOLSUD = 0.0
01602          IF( (IG1.EQ. GLOBAL))THEN
01603             N1 = NJ
01604             N2 = NJ - 1
01605             S1 = 1
01606             S2 = 2
01607          ENDIF 
01608          IF( (IG1.EQ. NORD))THEN
01609             N1 = NJ
01610             N2 = NJ - 1
01611             S1 = -NJ+1
01612             S2 = -NJ+2
01613          ENDIF 
01614          IF( (IG1.EQ. SUD))THEN
01615             N1 = 2*NJ
01616             N2 = 2*NJ - 1
01617             S1 = 1
01618             S2 = 2
01619          ENDIF 
01620          IF( (PTPOLN))THEN
01621             DO 23022 I=1,NI
01622                VPOLNOR = VPOLNOR + 9.0 * ZI(I, N1) + ZI(I,N2)
01623 23022       CONTINUE 
01624             VPOLNOR = 0.1 * VPOLNOR / REAL(NI,ip_realwp_p)
01625          ENDIF 
01626          IF( (PTPOLS))THEN
01627             DO 23026 I=1,NI
01628                VPOLSUD = VPOLSUD + 9.0 * ZI(I, S1) + ZI(I,S2)
01629 23026       CONTINUE 
01630             VPOLSUD = 0.1 * VPOLSUD / REAL(NI,ip_realwp_p)
01631          ENDIF 
01632       ENDIF 
01633       DO 23028 I=1,NPOLNUM
01634          ZO(NPOLPTS(NPOLPTR+I-1)) = VPOLNOR
01635 23028 CONTINUE 
01636       DO 23030 I=1,SPOLNUM
01637          ZO(SPOLPTS(SPOLPTR+I-1)) = VPOLSUD
01638 23030 CONTINUE 
01639       RETURN
01640       END
01641 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01642 C
01643       SUBROUTINE QQQALOC(XLAT, XLON, LI, LJ, I1, I2, J1, J2)
01644       USE mod_kinds_oasis
01645       USE memoir
01646       IMPLICIT NONE
01647 C      EXTERNAL MEMOIRH
01648       INTEGER (kind=ip_intwp_p) LI,LJ,I1,I2,J1,J2
01649       REAL(kind=ip_realwp_p) XLAT(LI,LJ), XLON(LI,LJ)
01650 C-- COMMON QQQCOM1
01651 C--
01652       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
01653       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
01654      $     AYTMP,CXTMP,CYTMP
01655       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
01656      $     SPOLPTS
01657       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
01658       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
01659       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
01660      $     NPOLMAX,SPOLMAX
01661       LOGICAL PTPOLN, PTPOLS
01662       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
01663      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
01664      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
01665 C----------------------
01666       CALL MEMOIRH(XGD, XGDPTR, LI*LJ, 0)
01667       CALL MEMOIRH(YGD, YGDPTR, LI*LJ, 0)
01668       CALL MEMOIRH(ZTMP,ZPTR,(I2-I1+1)*(J2-J1+1), 0)
01669       CALL MEMOIRH(AXTMP,AXPTR,(I2-I1+1), 0)
01670       CALL MEMOIRH(AYTMP,AYPTR,(J2-J1+1), 0)
01671       CALL MEMOIRH(CXTMP,CXPTR, (I2-I1+1) * 6, 0)
01672       CALL MEMOIRH(CYTMP,CYPTR, (J2-J1+1) * 6, 0)
01673       CALL MEMOIRH(NPOLPTS, NPOLPTR, 128, 0)
01674       CALL MEMOIRH(SPOLPTS, SPOLPTR, 128, 0)
01675       NPOLMAX = 128
01676       SPOLMAX = 128
01677       RETURN
01678       END
01679 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01680 C
01681       SUBROUTINE QQQCHK(XLAT, XLON, LILJ)
01682       USE mod_kinds_oasis
01683       USE memoir
01684       IMPLICIT NONE
01685 
01686 C-- COMMON QQQCOM1
01687 C--
01688       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
01689       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
01690      $     AYTMP,CXTMP,CYTMP
01691       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
01692      $     SPOLPTS
01693       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
01694       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
01695       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
01696      $     NPOLMAX,SPOLMAX
01697       LOGICAL PTPOLN, PTPOLS
01698       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
01699      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
01700      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
01701 C----------------------
01702 
01703       INTEGER (kind=ip_intwp_p) LILJ
01704       REAL(kind=ip_realwp_p) XLAT(LILJ), XLON(LILJ)
01705       REAL(kind=ip_realwp_p) TMP
01706       LOGICAL BADLAT
01707       INTEGER (kind=ip_intwp_p) NPOLPT2,SPOLPT2
01708       INTEGER (kind=ip_intwp_p) I,J
01709 C  VERIF POUR POINT AU POLE
01710 
01711       REAL(kind=ip_realwp_p) EPSILON
01712       DATA EPSILON /1.0E-6/
01713       BADLAT = .FALSE.
01714       DO 23000 I=1,LILJ
01715          TMP = XLON(I)
01716          XLON(I)=MOD(MOD(XLON(I),360.0_ip_realwp_p)
01717      $       +360.0_ip_realwp_p,360.0_ip_realwp_p)
01718          TMP = XLAT(I)
01719          XLAT(I)=MAX(-90.0_ip_realwp_p,
01720      $       MIN(90.0_ip_realwp_p,XLAT(I)))
01721          IF( (TMP.NE. XLAT(I)))THEN
01722             BADLAT = .TRUE.
01723          ENDIF 
01724 23000 CONTINUE 
01725       IF( (BADLAT))THEN
01726          WRITE(6,*) '<QQQGLAT> Latitudes gotten > 90.0 or < -90.0'
01727          WRITE(6,*)
01728      $    'Put latitudes back into -90.0 and +90.0'
01729          WRITE(6,*) 'LAT. > 90.0 = 90.0, LAT < -90.0 = -90.0'
01730       ENDIF 
01731       PTPOLN = .FALSE.
01732       NPOLNUM = 0
01733       DO 23006 I=1,LILJ
01734          IF( (EPSILON.GT. ABS(90.0 - XLAT(I))))THEN
01735             IF( (NPOLNUM.EQ. NPOLMAX))THEN
01736                CALL MEMOIRH(NPOLPTS, NPOLPT2, NPOLMAX+128, NPOLMAX)
01737                NPOLPTR = NPOLPT2
01738                NPOLMAX = NPOLMAX + 128
01739             ENDIF 
01740             PTPOLN = .TRUE.
01741             NPOLPTS(NPOLPTR+NPOLNUM) = I
01742             NPOLNUM = NPOLNUM + 1
01743          ENDIF 
01744 23006 CONTINUE 
01745       PTPOLS = .FALSE.
01746       SPOLNUM = 0
01747       DO 23016 I=1,LILJ
01748          IF( (EPSILON.GT. ABS(-90.0 - XLAT(I))))THEN
01749             IF( (SPOLNUM.EQ. SPOLMAX))THEN
01750                CALL MEMOIRH(SPOLPTS, SPOLPT2, SPOLMAX+128, SPOLMAX)
01751                SPOLPTR = SPOLPT2
01752                SPOLMAX = SPOLMAX + 128
01753             ENDIF 
01754             PTPOLS = .TRUE.
01755             SPOLPTS(SPOLPTR+SPOLNUM) = I
01756             SPOLNUM = SPOLNUM + 1
01757          ENDIF 
01758 23016 CONTINUE 
01759       RETURN
01760       END
01761 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01762 C
01763       LOGICAL FUNCTION QQQCMP(XLAT, XLON, LILJ, NI, NJ,GRTYP, IG1,
01764      $ IG2, IG3, IG4, XORSUM)
01765       USE mod_kinds_oasis
01766       INTEGER (kind=ip_intwp_p) LILJ
01767       REAL(kind=ip_realwp_p) XLAT(LILJ), XLON(LILJ)
01768       INTEGER (kind=ip_intwp_p) NI, NJ, IG1, IG2, IG3, IG4, XORSUM
01769       CHARACTER*1 GRTYP
01770       INTEGER (kind=ip_intwp_p) XORCALC
01771       CHARACTER*4 GRTMP
01772 C-- COMMON QQQCOM1
01773 C--
01774       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
01775       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
01776      $     AYTMP,CXTMP,CYTMP
01777       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
01778      $     SPOLPTS
01779       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
01780       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
01781       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
01782      $     NPOLMAX,SPOLMAX
01783       LOGICAL PTPOLN, PTPOLS
01784       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
01785      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
01786      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
01787 C----------------------
01788 C-- COMMON QQQCOM2
01789 C-
01790       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
01791       INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
01792       INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
01793       COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
01794      $,LSTIG3,LSTIG4,LSTXOR
01795 C-
01796 C------------
01797 C* COMMON QQQCOM3
01798 C
01799 
01800       CHARACTER*1 LSTGTYP
01801       COMMON /QQQCOM3/ LSTGTYP
01802 C------------
01803 
01804       QQQCMP = .TRUE.
01805       IF( (LILJ.NE.LSTLILJ.OR.NI.NE.LSTNI.OR.NJ.NE.LSTNJ))THEN
01806          QQQCMP = .FALSE.
01807       ENDIF 
01808       GRTMP = GRTYP
01809       IF( (IG1.NE.LSTIG1.OR.IG2.NE.LSTIG2.OR.IG3.NE.LSTIG3.OR.IG4
01810      $.NE.LSTIG4.OR.GRTMP.NE.LSTGTYP))THEN
01811          QQQCMP = .FALSE.
01812       ENDIF 
01813       XORSUM = XORCALC(XLAT, XLON, LILJ)
01814       IF( (XORSUM.NE.LSTXOR))THEN
01815          QQQCMP = .FALSE.
01816       ENDIF 
01817       RETURN
01818       END
01819 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01820 C
01821       BLOCK DATA GAUSSGDBLK
01822       USE mod_kinds_oasis
01823       IMPLICIT NONE
01824 C
01825 C-- COMMON GAUSSGD
01826       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
01827       INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
01828       COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
01829 C----------------------
01830       DATA IROOTS, ILROOTS /0, 0/
01831       END
01832 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01833 C
01834       SUBROUTINE QQQGLAT(NJ,HEM)
01835 C* ----------------------------------------------------------------------
01836 C* S/R QQQGLAT - CALCUL DES LATITUDES D'UNE GRILLE GAUSSIENNE
01837 C*
01838 C  AUTEUR: YVES CHARTIER. MARS 1991.
01839 C* ----------------------------------------------------------------------
01840       USE mod_kinds_oasis
01841       USE memoir
01842       IMPLICIT NONE
01843       INTEGER (kind=ip_intwp_p) NJ, HEM
01844       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
01845       PARAMETER (GLOBAL = 0)
01846       PARAMETER (NORD   = 1)
01847       PARAMETER (SUD   = 2)
01848       PARAMETER (SUDNORD= 0)
01849       PARAMETER (NORDSUD= 1)
01850 C
01851 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
01852 C
01853 C RDTODG = 180/PIE, DGTORD = PIE/180
01854 
01855       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
01856       DATA PIE   /3.14159265358979323846/
01857       DATA RDTODG /57.295779513082/
01858       DATA DGTORD /1.7453292519943E-2/
01859 C
01860 C-- COMMON GAUSSGD
01861       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
01862       INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
01863       COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
01864 C----------------------
01865       EXTERNAL DGAUSS
01866       EXTERNAL GAUSSGDBLK
01867       INTEGER (kind=ip_intwp_p) J,NPOLY
01868       INTEGER (kind=ip_intwp_p) NJ2
01869       IF( (IROOTS.NE. 0))THEN
01870          CALL MEMOIRH(ROOTS,  IROOTS, 0, 0)
01871       ENDIF 
01872       IF( (ILROOTS.NE. 0))THEN
01873          CALL MEMOIRH(LROOTS,ILROOTS, 0, 0)
01874       ENDIF 
01875       IF( (HEM.NE. GLOBAL))THEN
01876          NPOLY = NJ * 2
01877       ELSE 
01878          NPOLY = NJ
01879       ENDIF 
01880       CALL MEMOIRH(ROOTS, IROOTS, NPOLY, 0)
01881       CALL MEMOIRH(LROOTS,ILROOTS,NPOLY, 0)
01882       CALL DGAUSS(NPOLY, ROOTS(IROOTS), HEM)
01883       IF( (HEM.EQ. GLOBAL))THEN
01884           NJ2 = NJ / 2
01885           DO 23008 J=1,NJ2
01886             LROOTS(J-1+ILROOTS)=90.0-RDTODG*ACOS(ROOTS(IROOTS+NJ2-J)
01887      $          )
01888 23008     CONTINUE 
01889       ENDIF 
01890       IF( (HEM.EQ. NORD))THEN
01891          DO 23012 J=1,NJ
01892             LROOTS(J-1+ILROOTS)=90.0-RDTODG*ACOS(ROOTS(IROOTS+NJ-J))
01893 23012    CONTINUE 
01894       ENDIF 
01895       IF( (HEM.EQ. SUD))THEN
01896          DO 23016 J=1,NJ
01897             LROOTS(ILROOTS-1+J)=-(90.0-RDTODG*ACOS(ROOTS(IROOTS-1+J)
01898      $      ))
01899 23016    CONTINUE 
01900       ENDIF 
01901       RETURN
01902       END
01903 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01904 C
01905       SUBROUTINE QQQKEEP(LILJ,NI,NJ,GRTYP,IG1,IG2,IG3,IG4,XORSUM)
01906       USE mod_kinds_oasis
01907       IMPLICIT NONE
01908       INTEGER (kind=ip_intwp_p) LILJ, NI, NJ, IG1, IG2, IG3, IG4, XORSUM
01909       CHARACTER*1 GRTYP
01910 C-- COMMON QQQCOM2
01911 C-
01912 
01913       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
01914       INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
01915       INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
01916       COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
01917      $,LSTIG3,LSTIG4,LSTXOR
01918 C-
01919 C------------
01920 C* COMMON QQQCOM3
01921 C
01922 
01923       CHARACTER*1 LSTGTYP
01924       COMMON /QQQCOM3/ LSTGTYP
01925 C------------
01926 
01927       LSTLILJ = LILJ
01928       LSTNI   = NI
01929       LSTNJ   = NJ
01930       LSTGTYP = GRTYP
01931       LSTIG1  = IG1
01932       LSTIG2  = IG2
01933       LSTIG3  = IG3
01934       LSTIG4  = IG4
01935       LSTXOR  = XORSUM
01936       RETURN
01937       END
01938 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01939 C
01940       SUBROUTINE QQQXTRAP(ZO, PX, PY, NPTS, Z, I1, I2, J1, J2)
01941       USE mod_kinds_oasis
01942       IMPLICIT NONE
01943       INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
01944       REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
01945       REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
01946       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
01947       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
01948       PARAMETER (VOISIN  =   0)
01949       PARAMETER (LINEAIR =   1)
01950       PARAMETER (CUBIQUE =   3)
01951       PARAMETER (OUI   =   1)
01952       PARAMETER (MINIMUM =   2)
01953       PARAMETER (MAXIMUM =   3)
01954       PARAMETER (VALEUR  =   4)
01955       PARAMETER (ABORT   =  13)
01956       LOGICAL FLGXTRAP
01957       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
01958       REAL(kind=ip_realwp_p) VALXTRAP
01959       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
01960       INTEGER (kind=ip_intwp_p) N, I, J, OFFL, OFFR
01961       REAL(kind=ip_realwp_p) RMIN, RMAX, TEMPMIN, TEMPMAX
01962       IF( (.NOT.FLGXTRAP))THEN
01963          RETURN
01964       ENDIF 
01965       IF( (CODXTRAP.EQ. OUI))THEN
01966          RETURN
01967       ENDIF 
01968       IF( (ORDINT.EQ. 3))THEN
01969          OFFR = 2
01970          OFFL = 1
01971       ELSE 
01972          OFFR = 0
01973          OFFL = 0
01974       ENDIF 
01975       RMIN = Z(I1, J1)
01976       RMAX = Z(I1, J1)
01977       DO 23006 J=J1, J2
01978          DO 23008 I=I1, I2
01979             IF( (Z(I,J).LT. RMIN))THEN
01980                RMIN = Z(I,J)
01981             ENDIF 
01982             IF( (Z(I,J).GT. RMAX))THEN
01983                RMAX = Z(I,J)
01984             ENDIF 
01985 23008    CONTINUE 
01986 23006 CONTINUE 
01987       TEMPMIN = RMIN - 0.25*(RMAX - RMIN)
01988       TEMPMAX = RMAX + 0.25*(RMAX - RMIN)
01989       RMIN = TEMPMIN
01990       RMAX = TEMPMAX
01991       DO 23014 N=1, NPTS
01992          I = INT(PX(N))
01993          J = INT(PY(N))
01994          IF( (I.LT.(I1+OFFL).OR. J.LT.(J1+OFFL).OR. I.GT. (I2-OFFR)
01995      $   .OR. J.GT. (J2-OFFR)))THEN
01996             IF( (CODXTRAP.EQ. VOISIN))THEN
01997                I = MIN(I2, MAX(I1, NINT(PX(N))))
01998                J = MIN(J2, MAX(J1, NINT(PY(N))))
01999                ZO(N) = Z(I,J)
02000             ENDIF 
02001             IF( (CODXTRAP.EQ. MINIMUM))THEN
02002                ZO(N) = RMIN
02003             ENDIF 
02004             IF( (CODXTRAP.EQ. MAXIMUM))THEN
02005                ZO(N) = RMAX
02006             ENDIF 
02007             IF( (CODXTRAP.EQ. VALEUR))THEN
02008                ZO(N) = VALXTRAP
02009             ENDIF 
02010          ENDIF 
02011 23014 CONTINUE 
02012       RETURN
02013       END
02014 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02015 C
02016       SUBROUTINE RGDINT(ZO,PX,PY,NPTS,Z,I1,I2,J1,J2)
02017 C*******
02018 C AUTEUR: Y.CHARTIER, DRPN
02019 C        FEVRIER 1991
02020 C
02021 C OBJET:  INTERPOLATION BI-CUBIQUE DE POINTS A PARTIR D'UNE GRILLE
02022 C        SOURCE REGULIERE.
02023 C
02024 C*******
02025       USE mod_kinds_oasis
02026       IMPLICIT NONE
02027       INTEGER (kind=ip_intwp_p) NPTS,I1,I2,J1,J2
02028       REAL(kind=ip_realwp_p) ZO(NPTS),PX(NPTS),PY(NPTS)
02029       REAL(kind=ip_realwp_p) Z(I1:I2,J1:J2)
02030 C
02031 C  NPTS   : NOMBRE DE POINTS A INTERPOLER
02032 C  I1:I2  : DIMENSION DE LA GRILLE SOURCE SELON X
02033 C  J1:J2  : DIMENSION DE LA GRILLE SOURCE SELON Y
02034 C  ZO     : VECTEUR DE SORTIE CONTENANT LES VALEURS INTERPOLEES
02035 C  PX     : VECTEUR CONTENANT LA POSITION X DES POINTS QUE L'ON
02036 C         : VEUT INTERPOLER
02037 C  PY     : VECTEUR CONTENANT LA POSITION Y DES POINTS QUE L'ON
02038 C         : VEUT INTERPOLER
02039 C  Z      : VALEURS DE LA GRILLE SOURCE.
02040 C
02041 C===========================================
02042 C
02043 C     *   *   *   *
02044 C
02045 C     *   *   *   *
02046 C           #        ==>   PT (X,Y)
02047 C     *  (=)  *   *  ==> = PT (IIND, JIND)
02048 C
02049 C     *   *   *   *
02050 C
02051 C===========================================
02052 
02053       REAL(kind=ip_realwp_p) Y1,Y2,Y3,Y4
02054       INTEGER (kind=ip_intwp_p) M,N,I,J,STRIDE
02055       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
02056       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
02057       PARAMETER (VOISIN  =   0)
02058       PARAMETER (LINEAIR =   1)
02059       PARAMETER (CUBIQUE =   3)
02060       PARAMETER (OUI   =   1)
02061       PARAMETER (MINIMUM =   2)
02062       PARAMETER (MAXIMUM =   3)
02063       PARAMETER (VALEUR  =   4)
02064       PARAMETER (ABORT   =  13)
02065       LOGICAL FLGXTRAP
02066       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
02067       REAL(kind=ip_realwp_p) VALXTRAP
02068       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
02069       REAL(kind=ip_realwp_p) CUBIC, DX,DY,Z1,Z2,Z3,Z4
02070       REAL(kind=ip_realwp_p) ZLIN, ZZ1, ZZ2, ZDX
02071       CUBIC(Z1,Z2,Z3,Z4,DX)=((((Z4-Z1)*0.1666666666666 + 0.5*(Z2-Z3)
02072      $)*DX + 0.5*(Z1+Z3)-Z2)*DX + Z3-0.1666666666666*Z4-0.5*Z2-0.
02073      $3333333333333*Z1)*DX+Z2
02074       ZLIN(ZZ1,ZZ2,ZDX) = ZZ1 + (ZZ2 - ZZ1) * ZDX
02075       STRIDE = 1
02076       IF( (ORDINT.EQ. CUBIQUE))THEN
02077          DO 23002 N=1,NPTS
02078             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
02079             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
02080             DX = PX(N) - I
02081             DY = PY(N) - J
02082             Y1=CUBIC(Z(I-1,J-1),Z(I  ,J-1),Z(I+1,J-1),Z(I+2,J-1),DX)
02083             Y2=CUBIC(Z(I-1,J  ),Z(I  ,J  ),Z(I+1,J  ),Z(I+2,J  ),DX)
02084             Y3=CUBIC(Z(I-1,J+1),Z(I  ,J+1),Z(I+1,J+1),Z(I+2,J+1),DX)
02085             Y4=CUBIC(Z(I-1,J+2),Z(I  ,J+2),Z(I+1,J+2),Z(I+2,J+2),DX)
02086             ZO(N)=CUBIC(Y1,Y2,Y3,Y4,DY)
02087 23002    CONTINUE 
02088       ENDIF 
02089       IF( (ORDINT.EQ. LINEAIR))THEN
02090          DO 23006 N=1,NPTS
02091             I = MIN(I2-2,MAX(I1+1,INT(PX(N))))
02092             J = MIN(J2-2,MAX(J1+1,INT(PY(N))))
02093             DX = PX(N) - I
02094             DY = PY(N) - J
02095             Y2=ZLIN(Z(I,J  ),Z(I+1,J  ),DX)
02096             Y3=ZLIN(Z(I,J+1),Z(I+1,J+1),DX)
02097             ZO(N)=ZLIN(Y2,Y3,DY)
02098 23006    CONTINUE 
02099       ENDIF 
02100       IF( (ORDINT.EQ. VOISIN))THEN
02101          DO 23010 N=1,NPTS
02102             I = MIN(I2,MAX(I1,NINT(PX(N))))
02103             J = MIN(J2,MAX(J1,NINT(PY(N))))
02104             ZO(N)=Z(I,J)
02105 23010    CONTINUE 
02106       ENDIF 
02107       RETURN
02108       END
02109 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02110 C
02111       SUBROUTINE RGFREE()
02112       USE mod_kinds_oasis
02113       USE memoir
02114       IMPLICIT NONE
02115 C
02116 C* - COMMON QQQCOM1
02117 C
02118       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: XGD, YGD
02119       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ZTMP,AXTMP,
02120      $     AYTMP,CXTMP,CYTMP
02121       INTEGER (kind=ip_intwp_p), DIMENSION(:), POINTER :: NPOLPTS, 
02122      $     SPOLPTS
02123       INTEGER (kind=ip_intwp_p) XGDPTR,YGDPTR
02124       INTEGER (kind=ip_intwp_p) ZPTR, AXPTR, AYPTR, CXPTR, CYPTR
02125       INTEGER (kind=ip_intwp_p) NPOLPTR, SPOLPTR, NPOLNUM,SPOLNUM,
02126      $     NPOLMAX,SPOLMAX
02127       LOGICAL PTPOLN, PTPOLS
02128       COMMON /QQQCOM1/ XGD,YGD,ZTMP,AXTMP,AYTMP,CXTMP,CYTMP,XGDPTR,
02129      $YGDPTR,ZPTR,AXPTR,AYPTR,CXPTR,CYPTR,NPOLPTS,SPOLPTS,NPOLPTR,
02130      $ SPOLPTR,NPOLNUM,SPOLNUM,NPOLMAX,SPOLMAX,PTPOLN, PTPOLS
02131 C
02132 C* - COMMON QQQCOM2
02133 C
02134 
02135       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2
02136       INTEGER (kind=ip_intwp_p) LSTLILJ, LSTNI, LSTNJ, LSTXOR
02137       INTEGER (kind=ip_intwp_p) LSTIG1, LSTIG2, LSTIG3, LSTIG4
02138       COMMON /QQQCOM2/ I1,I2,J1,J2,LSTLILJ,LSTNI,LSTNJ,LSTIG1,LSTIG2
02139      $,LSTIG3,LSTIG4,LSTXOR
02140 C
02141 C* - COMMON QQQCOM3
02142 C
02143 
02144       CHARACTER*1 LSTGTYP
02145       COMMON /QQQCOM3/ LSTGTYP
02146 C* -----------------------------------------------------------------------
02147       IF( (XGDPTR.NE.0))THEN
02148          CALL MEMOIRH(XGD, XGDPTR, 0, 0)
02149          XGDPTR = 0
02150       ENDIF 
02151       IF( (YGDPTR.NE.0))THEN
02152          CALL MEMOIRH(YGD, YGDPTR, 0, 0)
02153          YGDPTR = 0
02154       ENDIF 
02155       IF( (ZPTR.NE.0))THEN
02156          CALL MEMOIRH(ZTMP,ZPTR, 0, 0)
02157          ZPTR = 0
02158       ENDIF 
02159       IF( (AXPTR.NE.0))THEN
02160          CALL MEMOIRH(AXTMP,AXPTR, 0, 0)
02161          AXPTR = 0
02162       ENDIF 
02163       IF( (AYPTR.NE.0))THEN
02164          CALL MEMOIRH(AYTMP,AYPTR, 0, 0)
02165          AYPTR = 0
02166       ENDIF 
02167       IF( (CXPTR.NE.0))THEN
02168          CALL MEMOIRH(CXTMP, CXPTR, 0, 0)
02169          CXPTR = 0
02170       ENDIF 
02171       IF( (CYPTR.NE.0))THEN
02172          CALL MEMOIRH(CYTMP, CYPTR, 0, 0)
02173          CYPTR = 0
02174       ENDIF 
02175       IF( (NPOLPTR.NE.0))THEN
02176          CALL MEMOIRH(NPOLPTS, NPOLPTR, 0, 0)
02177          NPOLPTR = 0
02178          NPOLMAX = 0
02179       ENDIF 
02180       IF( (SPOLPTR.NE.0))THEN
02181          CALL MEMOIRH(SPOLPTS, SPOLPTR, 0, 0)
02182          SPOLPTR = 0
02183          SPOLMAX = 0
02184       ENDIF 
02185       LSTLILJ = 0
02186       LSTNI   = 0
02187       LSTNJ   = 0
02188       LSTGTYP = ' '
02189       LSTIG1  = 0
02190       LSTIG2  = 0
02191       LSTIG3  = 0
02192       LSTIG4  = 0
02193       LSTXOR  = 0
02194       RETURN
02195       END
02196 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02197 C
02198       SUBROUTINE RGLL2GD(SPDO,PSIO,XLON,LI,LJ,GRTYP,IG1,IG2,IG3,IG4)
02199       USE mod_kinds_oasis
02200       IMPLICIT NONE
02201       INTEGER (kind=ip_intwp_p) LI,LJ
02202       REAL(kind=ip_realwp_p) SPDO(LI,LJ), PSIO(LI,LJ), XLON(LI,LJ)
02203       CHARACTER*1 GRTYP
02204       INTEGER (kind=ip_intwp_p) IG1,IG2,IG3,IG4
02205       EXTERNAL CIGAXG
02206 C
02207 C AUTEUR   - Y. CHARTIER - AVRIL 91
02208 C
02209 C OBJET(RGLL2GD)
02210 C         - PASSE DE PSIOTESSE ET DIRECTION
02211 C         - A VENT DE GRILLE (COMPOSANTES U ET V)
02212 C
02213 C APPEL    - CALL RGLL2GD(SPD,PSI,LI,LJ,IYP,XG1,XG2,XG3,XG4)
02214 C
02215 C MODULES  - XGAIG
02216 C
02217 C ARGUMENTS
02218 C  IN/OUT - SPD   - A L'ENTREE CONTIENT LA PSIOTESSE DU VENT ET
02219 C                   A LA SORTIE LA COMPOSANTE U.
02220 C  IN/OUT - PSI   - A L'ENTREE CONTIENT LA DIRECTION DU VENT ET
02221 C                   A LA SORTIE LA COMPOSANTE V.
02222 C   IN    - LI    - PREMIERE DIMENSION DES CHAMPS SPD ET PSI
02223 C   IN    - LJ    - DEUXIEME DIMENSION DES CHAMPS SPD ET PSI
02224 C   IN    - IGTYP  - TYPE DE GRILLE (VOIR OUVRIR)
02225 C   IN    - XG1   - ** DESCRIPTEUR DE GRILLE (REEL),
02226 C   IN    - XG2   -    IGTYP = 'N', PI, PJ, D60, DGRW
02227 C   IN    - XG3   -    IGTYP = 'L', LAT0, LON0, DLAT, DLON,
02228 C   IN    - XG4   -    IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
02229 C                                                 = 1. NORD
02230 C                                                 = 2. SUD **
02231 C
02232 CMESSAGES - "ERREUR MAUVAISE GRILLE (RGLL2GD)"
02233 C
02234 C-------------------------------------------------------------
02235 C
02236 C
02237 C
02238 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
02239 C
02240 C RDTODG = 180/PIE, DGTORD = PIE/180
02241 
02242       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
02243       DATA PIE   /3.14159265358979323846/
02244       DATA RDTODG /57.295779513082/
02245       DATA DGTORD /1.7453292519943E-2/
02246 C
02247 
02248       INTEGER (kind=ip_intwp_p) I,J
02249       REAL(kind=ip_realwp_p) PSI,U,V
02250       REAL(kind=ip_realwp_p) XG1, XG2, XG3, XG4
02251       IF( (GRTYP.EQ. 'N'))THEN
02252          CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
02253          DO 23002 I=1,LI
02254             DO 23004 J=1,LJ
02255                PSI =XLON(I,J)+XG4-PSIO(I,J)
02256                U = COS(PSI*DGTORD)*SPDO(I,J)
02257                V = SIN(PSI*DGTORD)*SPDO(I,J)
02258                SPDO(I,J) = U
02259                PSIO(I,J) = V
02260 23004       CONTINUE 
02261 23002    CONTINUE 
02262          RETURN
02263       ENDIF 
02264       IF( (GRTYP.EQ. 'S'))THEN
02265          CALL CIGAXG(GRTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
02266          DO 23008 I=1,LI
02267             DO 23010 J=1,LJ
02268                PSI =180.0 - XLON(I,J)+XG4-PSIO(I,J)
02269                U = COS(PSI*DGTORD)*SPDO(I,J)
02270                V = SIN(PSI*DGTORD)*SPDO(I,J)
02271                SPDO(I,J) = U
02272                PSIO(I,J) = V
02273 23010       CONTINUE 
02274 23008    CONTINUE 
02275          RETURN
02276       ENDIF 
02277       IF( (GRTYP.EQ.'A'.OR.GRTYP.EQ.'B'.OR.GRTYP.EQ.'G'.OR.GRTYP.EQ.
02278      $'L'))THEN
02279          DO 23014 I=1,LI
02280             DO 23016 J=1,LJ
02281                PSI = 270.0 - PSIO(I,J)
02282                U = COS(PSI*DGTORD)*SPDO(I,J)
02283                V = SIN(PSI*DGTORD)*SPDO(I,J)
02284                SPDO(I,J) = U
02285                PSIO(I,J) = V
02286 23016       CONTINUE 
02287 23014    CONTINUE 
02288          RETURN
02289       ENDIF 
02290       WRITE(6, 600) GRTYP
02291 600   FORMAT('0',' Error, bad grid (RGLL2GD) - GRTYP = ', A1)
02292       RETURN
02293       END
02294 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02295 C
02296       SUBROUTINE RGOPTC(OP, VAL, FLAG)
02297       USE mod_kinds_oasis
02298       IMPLICIT NONE
02299       CHARACTER *(*) OP, VAL
02300       LOGICAL FLAG
02301 CFLAG = .TRUE.  MODE SET
02302 CFLAG = .FALSE. MODE GET
02303 
02304       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
02305       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
02306       PARAMETER (VOISIN  =   0)
02307       PARAMETER (LINEAIR =   1)
02308       PARAMETER (CUBIQUE =   3)
02309       PARAMETER (OUI   =   1)
02310       PARAMETER (MINIMUM =   2)
02311       PARAMETER (MAXIMUM =   3)
02312       PARAMETER (VALEUR  =   4)
02313       PARAMETER (ABORT   =  13)
02314       LOGICAL FLGXTRAP
02315       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
02316       REAL(kind=ip_realwp_p) VALXTRAP
02317       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
02318       IF( (FLAG))THEN
02319          IF( (OP.EQ.  'EXTRAP'))THEN
02320             IF( (VAL.EQ.  'OUI'))THEN
02321                CODXTRAP = OUI
02322             ELSE 
02323                IF( (VAL.EQ.  'ABORT'))THEN
02324                   CODXTRAP = ABORT
02325                ELSE 
02326                   IF( (VAL.EQ.  'MAXIMUM'))THEN
02327                      CODXTRAP = MAXIMUM
02328                   ELSE 
02329                      IF( (VAL.EQ.  'MINIMUM'))THEN
02330                         CODXTRAP = MINIMUM
02331                      ELSE 
02332                         IF( (VAL.EQ.  'VOISIN'))THEN
02333                            CODXTRAP = VOISIN
02334                         ELSE 
02335                            IF( (VAL.EQ.  'VALEUR'))THEN
02336                               CODXTRAP = VALEUR
02337                            ELSE 
02338                               WRITE(6,*)
02339      $                         ' <RGOPTC>: bad value for VAL'
02340                               WRITE(6,*) '           VAL = ', VAL
02341                               WRITE(6,*)
02342      $                         '           VAL initialized to ''ABORT'''
02343                               CODXTRAP = ABORT
02344                            ENDIF 
02345                         ENDIF 
02346                      ENDIF 
02347                   ENDIF 
02348                ENDIF 
02349             ENDIF 
02350          ELSE 
02351             IF( (OP.EQ.  'INTERP'.OR.  OP.EQ.  'INTERP'))THEN
02352                IF( (VAL.EQ.  'VOISIN'))THEN
02353                   ORDINT = VOISIN
02354                ELSE 
02355                   IF( (VAL.EQ.  'LINEAIR'))THEN
02356                      ORDINT = LINEAIR
02357                   ELSE 
02358                      IF( (VAL.EQ.  'CUBIQUE'))THEN
02359                         ORDINT = CUBIQUE
02360                      ELSE 
02361                         WRITE(6,*)
02362      $                   ' <RGOPTC>: bad value for VAL'
02363                         WRITE(6,*) '           VAL = ', VAL
02364                         WRITE(6,*)
02365      $                   '           VAL initialized to ''CUBIQUE'''
02366                         ORDINT = CUBIQUE
02367                      ENDIF 
02368                   ENDIF 
02369                ENDIF 
02370             ELSE 
02371                WRITE(6,*) ' <RGOPTC>: bad value for OP'
02372                WRITE(6,*)
02373      $ '           OP should be equal to ''EXTRAP'' or ''INTERP'''
02374             ENDIF 
02375          ENDIF 
02376       ELSE 
02377          IF( (OP.EQ.  'EXTRAP'))THEN
02378             IF( (CODXTRAP.EQ.  OUI))THEN
02379                VAL = 'OUI'
02380             ELSE 
02381                IF( (CODXTRAP.EQ.  ABORT))THEN
02382                   VAL = 'ABORT'
02383                ELSE 
02384                   IF( (CODXTRAP.EQ.  MAXIMUM))THEN
02385                      VAL = 'MAXIMUM'
02386                   ELSE 
02387                      IF( (CODXTRAP.EQ.  MINIMUM))THEN
02388                         VAL = 'MINIMUM'
02389                      ELSE 
02390                         IF( (CODXTRAP.EQ.  VOISIN))THEN
02391                            VAL = 'VOISIN'
02392                         ELSE 
02393                            IF( (CODXTRAP.EQ.  VALEUR))THEN
02394                               VAL = 'VALEUR'
02395                            ELSE 
02396                               WRITE(6,*)
02397      $ ' <RGOPTC>: bad value for CODXTRAP'
02398                               WRITE(6,*) '           CODXTRAP = ',
02399      $                         CODXTRAP
02400                               VAL = 'SCRAP'
02401                            ENDIF 
02402                         ENDIF 
02403                      ENDIF 
02404                   ENDIF 
02405                ENDIF 
02406             ENDIF 
02407          ELSE 
02408             IF( (OP.EQ.  'INTERP'))THEN
02409                IF( (ORDINT.EQ.  VOISIN))THEN
02410                   VAL = 'VOISIN'
02411                ELSE 
02412                   IF( (ORDINT.EQ.  LINEAIR))THEN
02413                      VAL = 'LINEAIR'
02414                   ELSE 
02415                      IF( (ORDINT.EQ.  CUBIQUE))THEN
02416                         VAL = 'CUBIQUE'
02417                      ELSE 
02418                         WRITE(6,*)
02419      $                   ' <RGOPTC>: bad value for ORDINT'
02420                         WRITE(6,*) '           ORDINT = ', ORDINT
02421                         VAL = 'SCRAP'
02422                      ENDIF 
02423                   ENDIF 
02424                ENDIF 
02425             ELSE 
02426                WRITE(6,*) ' <RGOPTC>: bad value for OP'
02427                WRITE(6,*)
02428      $ '           OP should be equal to ''EXTRAP'' OU ''INTERP'''
02429             ENDIF 
02430          ENDIF 
02431       ENDIF 
02432       RETURN
02433       END
02434 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02435 C
02436       BLOCK DATA QQQXTRPBLK
02437       USE mod_kinds_oasis
02438       IMPLICIT NONE
02439       INTEGER (kind=ip_intwp_p) OUI
02440       PARAMETER (OUI   =   1)
02441       LOGICAL FLGXTRAP
02442       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
02443       REAL(kind=ip_realwp_p) VALXTRAP
02444       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
02445       DATA CODXTRAP / OUI /
02446       DATA FLGXTRAP / .FALSE. /
02447       DATA ORDINT   / 3 /
02448       END
02449 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02450 C
02451       SUBROUTINE RGOPTI(OP, VAL, FLAG)
02452       USE mod_kinds_oasis
02453       IMPLICIT NONE
02454       CHARACTER *(*) OP
02455       INTEGER (kind=ip_intwp_p) VAL
02456       LOGICAL FLAG
02457 CFLAG = .TRUE.  MODE SET
02458 CFLAG = .FALSE. MODE GET
02459 
02460       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
02461       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
02462       PARAMETER (VOISIN  =   0)
02463       PARAMETER (LINEAIR =   1)
02464       PARAMETER (CUBIQUE =   3)
02465       PARAMETER (OUI   =   1)
02466       PARAMETER (MINIMUM =   2)
02467       PARAMETER (MAXIMUM =   3)
02468       PARAMETER (VALEUR  =   4)
02469       PARAMETER (ABORT   =  13)
02470       LOGICAL FLGXTRAP
02471       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
02472       REAL(kind=ip_realwp_p) VALXTRAP
02473       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
02474       EXTERNAL QQQXTRPBLK
02475       IF( (FLAG))THEN
02476          IF( (OP.EQ.  'EXTRAP'))THEN
02477             VALXTRAP = REAL(VAL,ip_realwp_p)
02478          ELSE 
02479             IF( (OP.EQ.  'INTERP'))THEN
02480                IF( (VAL.EQ.  VOISIN.OR.  VAL.EQ.  LINEAIR.OR.  VAL
02481      $         .EQ. CUBIQUE))THEN
02482                   ORDINT = VAL
02483                ELSE 
02484                   WRITE(6,*) ' <RGOPTI>: bad value for VAL'
02485                   WRITE(6,*) '           VAL = ', VAL
02486                   WRITE(6,*)
02487      $             '           VAL initialized to ''CUBIQUE'''
02488                   ORDINT = CUBIQUE
02489                ENDIF 
02490             ELSE 
02491                WRITE(6,*) ' <RGOPTI>: bad value for OP'
02492                WRITE(6,*)
02493      $ '           OP should be equal to ''EXTRAP'' or ''INTERP'''
02494             ENDIF 
02495          ENDIF 
02496       ELSE 
02497          IF( (OP.EQ.  'EXTRAP'))THEN
02498             VAL = NINT(VALXTRAP)
02499          ELSE 
02500             IF( (OP.EQ.  'INTERP'))THEN
02501                IF( (ORDINT.EQ.  VOISIN.OR.  ORDINT.EQ.  LINEAIR.OR.
02502      $          ORDINT.EQ.  CUBIQUE))THEN
02503                   VAL = ORDINT
02504                ELSE 
02505                   WRITE(6,*) ' <RGOPTI>: bad value for ORDINT'
02506                   WRITE(6,*) '           ORDINT = ', ORDINT
02507                   VAL = -1
02508                ENDIF 
02509             ELSE 
02510                WRITE(6,*) ' <RGOPTI>: bad value for OP'
02511                WRITE(6,*)
02512      $ '           OP should be equal to ''EXTRAP'' or ''INTERP'''
02513             ENDIF 
02514          ENDIF 
02515       ENDIF 
02516       RETURN
02517       END
02518 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02519 C
02520       SUBROUTINE RGOPTR(OP, VAL, FLAG)
02521       USE mod_kinds_oasis
02522       IMPLICIT NONE
02523       CHARACTER *(*) OP
02524       REAL(kind=ip_realwp_p) VAL
02525       LOGICAL FLAG
02526 C  FLAG = .TRUE.  MODE SET
02527 C  FLAG = .FALSE. MODE GET
02528 
02529       INTEGER (kind=ip_intwp_p) VOISIN, LINEAIR, CUBIQUE
02530       INTEGER (kind=ip_intwp_p) OUI, ABORT, VALEUR, MAXIMUM, MINIMUM
02531       PARAMETER (VOISIN  =   0)
02532       PARAMETER (LINEAIR =   1)
02533       PARAMETER (CUBIQUE =   3)
02534       PARAMETER (OUI   =   1)
02535       PARAMETER (MINIMUM =   2)
02536       PARAMETER (MAXIMUM =   3)
02537       PARAMETER (VALEUR  =   4)
02538       PARAMETER (ABORT   =  13)
02539       LOGICAL FLGXTRAP
02540       INTEGER (kind=ip_intwp_p) CODXTRAP, ORDINT
02541       REAL(kind=ip_realwp_p) VALXTRAP
02542       COMMON /QQQXTRP/ ORDINT, FLGXTRAP, CODXTRAP, VALXTRAP
02543       IF( (FLAG))THEN
02544          IF( (OP.EQ. 'EXTRAP'))THEN
02545             VALXTRAP = VAL
02546          ELSE 
02547             IF( (OP.EQ. 'INTERP'))THEN
02548                IF( (VAL.EQ. VOISIN.OR. VAL.EQ. LINEAIR.OR. VAL.EQ.
02549      $          CUBIQUE))THEN
02550                   ORDINT = NINT(VAL)
02551                ELSE 
02552                   WRITE(6,*)' <RGOPTR>: bad value for VAL'
02553                   WRITE(6,*) '           VAL = ', VAL
02554                   WRITE(6,*)
02555      $             '           VAL initialized to ''CUBIQUE'''
02556                   ORDINT = CUBIQUE
02557                ENDIF 
02558             ELSE 
02559                WRITE(6,*) ' <RGOPTR>: bad value for OP'
02560                WRITE(6,*)
02561      $ '           OP should be equal to ''EXTRAP'' or ''INTERP'''
02562             ENDIF 
02563          ENDIF 
02564       ELSE 
02565          IF( (OP.EQ. 'EXTRAP'))THEN
02566             VAL = VALXTRAP
02567          ELSE 
02568             IF( (OP.EQ. 'INTERP'))THEN
02569                IF( (ORDINT.EQ. VOISIN.OR. ORDINT.EQ. LINEAIR.OR.
02570      $          ORDINT.EQ. CUBIQUE))THEN
02571                   VAL = REAL(ORDINT,ip_realwp_p)
02572                ELSE 
02573                   WRITE(6,*) ' <RGOPTR>: bad value for ORDINT'
02574                   WRITE(6,*) '           ORDINT = ', ORDINT
02575                   VAL = REAL(ORDINT,ip_realwp_p)
02576                ENDIF 
02577             ELSE 
02578                WRITE(6,*) ' <RGOPTR>: bad value for OP'
02579                WRITE(6,*)
02580      $ '           OP should be equal to ''EXTRAP'' OU ''INTERP'''
02581             ENDIF 
02582          ENDIF 
02583       ENDIF 
02584       RETURN
02585       END
02586 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02587 C
02588       SUBROUTINE VLLFXY(DLAT,DLON,X,Y,NI,NJ,D60,DGRW,PI,PJ,NHEM)
02589 C* ----------------------------------------------------------------------
02590 C* S/R VLLFXY(I) - COMPUTES THE GRID CO-ORDINATES OF A POINT
02591 C* ----------------------------------------------------------------------
02592       USE mod_kinds_oasis
02593       IMPLICIT NONE
02594 C* ----------------------------------------------------------------------
02595 C* ARGUMENTS
02596 C   OUT   - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N).
02597 C         - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E).
02598 C   IN    - X    - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
02599 C                  AS ORIGIN
02600 C         - Y    - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
02601 C                  AS ORIGIN
02602 C         - D60  - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
02603 C                  GRID AT 60 DEGREES
02604 C         - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
02605 C                  THE GRID (IN DEGREES)
02606 C         - NHEM - 1 FOR NORTHERN HEMISPHERE
02607 C                  2 FOR SOUTHERN HEMISPHERE
02608 C
02609 CNOTES    - THE COMPANION ROUTINE XYFLL, WHICH COMPUTES THE GRID
02610 C           CO-ORDINATES GIVEN THE LATITUDE AND LONGITUDE, IS ALSO
02611 C           AVAILABLE.
02612 C
02613 C-----------------------------------------------------------------------
02614 C
02615 C
02616 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
02617 C
02618 C RDTODG = 180/PIE, DGTORD = PIE/180
02619       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
02620       DATA PIE   /3.14159265358979323846/
02621       DATA RDTODG /57.295779513082/
02622       DATA DGTORD /1.7453292519943E-2/
02623 C
02624 
02625       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
02626       PARAMETER (GLOBAL = 0)
02627       PARAMETER (NORD   = 1)
02628       PARAMETER (SUD   = 2)
02629       PARAMETER (SUDNORD= 0)
02630       PARAMETER (NORDSUD= 1)
02631       INTEGER (kind=ip_intwp_p) NI,NJ, NHEM
02632       REAL(kind=ip_realwp_p) X(NI,NJ), Y(NI,NJ), DLAT(NI,NJ), 
02633      $     DLON(NI,NJ)
02634       REAL(kind=ip_realwp_p) X1, Y1
02635       REAL(kind=ip_realwp_p) D60, DGRW, PI, PJ
02636       REAL(kind=ip_realwp_p) RE,RE2,R2,RLON,RLAT,SINLAT,R
02637       INTEGER (kind=ip_intwp_p) I,J
02638       RE=1.866025*6.371E+6/D60
02639       RE2=RE**2
02640       DO 23000 I=1,NI
02641          DO 23002 J=1,NJ
02642             X1 = X(I,J) - PI
02643             Y1 = Y(I,J) - PJ
02644             IF( (X1.EQ. 0..AND. Y1.EQ. 0.))THEN
02645                DLAT(I,J)=90.
02646                DLON(I,J)=0.
02647 
02648 C
02649 C  CALCULATE LONGITUDE IN MAP COORDINATES.
02650 C
02651             ENDIF 
02652             IF((X1.EQ. 0.))THEN
02653                DLON(I,J)=SIGN(90._ip_realwp_p,Y1)
02654             ENDIF 
02655             IF((X1.NE. 0.))THEN
02656                DLON(I,J)=ATAN(Y1/X1)*RDTODG
02657             ENDIF 
02658             IF((X1.LT.  0.))THEN
02659                DLON(I,J)=DLON(I,J)+SIGN(180._ip_realwp_p,Y1)
02660 C
02661 C  * ADJUST LONGITUDE FOR GRID ORIENTATION.
02662 C
02663 
02664             ENDIF 
02665             DLON(I,J)=DLON(I,J)-DGRW
02666             IF((DLON(I,J).GT.  180.))THEN
02667                DLON(I,J)=DLON(I,J)-360.
02668             ENDIF 
02669             IF((DLON(I,J).LT. -180.))THEN
02670                DLON(I,J)=DLON(I,J)+360.
02671 C
02672 C     * CALCULATE LATITUDE.
02673 C
02674 
02675             ENDIF 
02676             R2=X1**2+Y1**2
02677             DLAT(I,J)=(RE2-R2)/(RE2+R2)
02678             DLAT(I,J)= ASIN(DLAT(I,J))*RDTODG
02679 C
02680 C     CHANGE SIGNS IF IN SOUTHERN HEMISPHERE.
02681 C
02682 
02683             IF((NHEM.EQ. 2))THEN
02684                DLAT(I,J)=-DLAT(I,J)
02685             ENDIF 
02686             IF((NHEM.EQ. 2))THEN
02687                DLON(I,J)=-DLON(I,J)
02688             ENDIF 
02689 23002    CONTINUE 
02690 23000 CONTINUE 
02691       RETURN
02692       END
02693 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02694 C
02695       SUBROUTINE VXYFLL(X,Y,DLAT,DLON,NPTS,D60,DGRW,PI,PJ,NHEM)
02696       USE mod_kinds_oasis
02697       IMPLICIT NONE
02698 
02699 C AUTHOR   - J.D. HENDERSON  -  FEB 75
02700 C
02701 C OBJECT(VXYFLL)
02702 C     - COMPUTES THE GRID CO-ORDINATES MEASURED FROM THE POLE OF A
02703 C       POINT, GIVEN THE LATITUDE AND LONGITUDE IN DEGREES.
02704 C
02705 C USAGE    - CALL VXYFLL(X,Y,DLAT,DLON,NPTS,D60,DGRW,PI,PJ,NHEM)
02706 C
02707 C ARGUMENTS
02708 C   OUT   - X    - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
02709 C                  AS ORIGIN
02710 C         - Y    - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
02711 C                  AS ORIGIN
02712 C   IN    - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N)
02713 C         - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E)
02714 C         - D60  - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
02715 C                  GRID AT 60 DEGREES
02716 C         - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
02717 C                  THE GRID (IN DEGREES)
02718 C         - NHEM - 1 FOR NORTHERN HEMISPHERE
02719 C                  2 FOR SOUTHERN HEMISPHERE
02720 C
02721 C NOTES    - THE COMPANION ROUTINE LLFXY, WHICH COMPUTES THE LATITUDE
02722 C         - AND LONGITUDE GIVEN THE GRID-COORDINATES,
02723 C         - IS ALSO AVAILABLE.
02724 C*--------------------------------------------------------------------
02725 C
02726 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
02727 C
02728 C RDTODG = 180/PIE, DGTORD = PIE/180
02729       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
02730       DATA PIE   /3.14159265358979323846/
02731       DATA RDTODG /57.295779513082/
02732       DATA DGTORD /1.7453292519943E-2/
02733 C
02734 
02735       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
02736       PARAMETER (GLOBAL = 0)
02737       PARAMETER (NORD   = 1)
02738       PARAMETER (SUD   = 2)
02739       PARAMETER (SUDNORD= 0)
02740       PARAMETER (NORDSUD= 1)
02741       INTEGER (kind=ip_intwp_p) NPTS, NHEM
02742       REAL(kind=ip_realwp_p) X(NPTS), Y(NPTS), DLAT(NPTS), DLON(NPTS)
02743       REAL(kind=ip_realwp_p) D60, DGRW, PI, PJ
02744       REAL(kind=ip_realwp_p) RE,RLON,RLAT,SINLAT,R
02745       INTEGER (kind=ip_intwp_p) I
02746       RE=1.866025*6.371E+6/D60
02747 C
02748 
02749       IF( (NHEM.EQ.NORD))THEN
02750          DO 23002 I=1,NPTS
02751             RLON=DGTORD*(DLON(I)+DGRW)
02752             RLAT=DGTORD*DLAT(I)
02753             SINLAT=SIN(RLAT)
02754             R=RE*SQRT((1.-SINLAT)/(1.+SINLAT))
02755             X(I)=R*COS(RLON) + PI
02756             Y(I)=R*SIN(RLON) + PJ
02757 23002    CONTINUE 
02758          RETURN
02759       ENDIF 
02760       IF( (NHEM.EQ.SUD))THEN
02761          DO 23006 I=1,NPTS
02762             RLON = DLON(I)
02763             IF( (RLON.GT. 180.0))THEN
02764                RLON = RLON - 360.0
02765             ENDIF 
02766             RLON=DGTORD*(-RLON+DGRW)
02767             RLAT=DGTORD*(-DLAT(I))
02768             SINLAT=SIN(RLAT)
02769             R=RE*SQRT((1.-SINLAT)/(1.+SINLAT))
02770             X(I)=R*COS(RLON)+PI
02771             Y(I)=R*SIN(RLON)+PJ
02772 23006    CONTINUE 
02773       ENDIF 
02774       RETURN
02775       END
02776 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02777 C
02778       FUNCTION XORCALC(ILAT, ILON, LILJ)
02779       USE mod_kinds_oasis
02780       INTEGER (kind=ip_intwp_p) XORCALC
02781       INTEGER (kind=ip_intwp_p) LILJ
02782       INTEGER (kind=ip_intwp_p) ILAT(LILJ), ILON(LILJ)
02783       XORCALC = 0
02784       DO 23000 I=1,LILJ,17
02785          XORCALC = IEOR(IEOR(XORCALC,ILAT(I)),ILON(I))
02786 23000 CONTINUE 
02787       RETURN
02788       END
02789 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02790 C
02791       SUBROUTINE XPNAXEG(NEWAXEX,NEWAXEY,I1,I2,J1,J2,NI,NJ,HEM)
02792       USE mod_kinds_oasis
02793       IMPLICIT NONE
02794       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ,HEM
02795       REAL(kind=ip_realwp_p) NEWAXEX(I1:I2)
02796       REAL(kind=ip_realwp_p) NEWAXEY(J1:J2)
02797 C
02798 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
02799 C
02800 C RDTODG = 180/PIE, DGTORD = PIE/180
02801 
02802       REAL(kind=ip_realwp_p) PIE,RDTODG,DGTORD
02803       DATA PIE   /3.14159265358979323846/
02804       DATA RDTODG /57.295779513082/
02805       DATA DGTORD /1.7453292519943E-2/
02806 C
02807 C-- COMMON GAUSSGD
02808       REAL(kind=ip_realwp_p), DIMENSION(:), POINTER :: ROOTS,LROOTS
02809       INTEGER (kind=ip_intwp_p) IROOTS, ILROOTS
02810       COMMON /GAUSSGD/ ROOTS,LROOTS,IROOTS,ILROOTS
02811 C----------------------
02812 
02813       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
02814       PARAMETER (GLOBAL = 0)
02815       PARAMETER (NORD   = 1)
02816       PARAMETER (SUD   = 2)
02817       PARAMETER (SUDNORD= 0)
02818       PARAMETER (NORDSUD= 1)
02819       INTEGER (kind=ip_intwp_p) I,J
02820       DO 23000 I=1,NI
02821          NEWAXEX(I) = REAL(I,ip_realwp_p)
02822 23000 CONTINUE 
02823       NEWAXEX(0) = NEWAXEX(NI) - REAL(NI,ip_realwp_p)
02824       NEWAXEX(NI+1) = NEWAXEX(1) + REAL(NI,ip_realwp_p)
02825       NEWAXEX(NI+2) = NEWAXEX(2) + REAL(NI,ip_realwp_p)
02826       IF( (HEM.EQ. NORD))THEN
02827          DO 23004 J=1,NJ
02828             NEWAXEY(J) = LROOTS(ILROOTS-1+J)
02829 23004    CONTINUE 
02830          NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
02831          NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
02832          DO 23006 J=-NJ+1,0
02833             NEWAXEY(J) = -(NEWAXEY(-J+1))
02834 23006    CONTINUE 
02835          NEWAXEY(-NJ  ) = -180.0 - NEWAXEY(-NJ+1)
02836          NEWAXEY(-NJ-1) = -180.0 - NEWAXEY(-NJ+2)
02837       ENDIF 
02838       IF( (HEM.EQ. SUD))THEN
02839          DO 23010 J=1,NJ
02840             NEWAXEY(J) = -(LROOTS(ILROOTS-1+NJ-J+1))
02841             NEWAXEY(NJ+J) = LROOTS(ILROOTS-1+NJ-J+1)
02842 23010    CONTINUE 
02843          NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
02844          NEWAXEY(-1) = -180.0 - NEWAXEY(2)
02845          NEWAXEY(NJ+1) = 180.0 - NEWAXEY(2*NJ)
02846          NEWAXEY(NJ+2) = 180.0 - NEWAXEY(2*NJ-1)
02847       ENDIF 
02848       IF( (HEM.EQ. GLOBAL))THEN
02849          IF( (0.EQ. MOD(NJ,2)))THEN
02850             DO 23016 J=1,NJ/2
02851                NEWAXEY(J)   =  -(LROOTS(ILROOTS-1+NJ/2-J+1))
02852                NEWAXEY(NJ/2+J) = LROOTS(ILROOTS-1+J)
02853 23016       CONTINUE 
02854          ELSE 
02855             DO 23018 J=1,NJ/2
02856                NEWAXEY(J)   =  -(LROOTS(ILROOTS-1+NJ/2-J+1))
02857                NEWAXEY(NJ/2+J+1) = LROOTS(ILROOTS-1+J)
02858 23018       CONTINUE 
02859             NEWAXEY(NJ/2+1) = 0.0
02860          ENDIF 
02861          NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
02862          NEWAXEY(-1) = -180.0 - NEWAXEY(2)
02863          NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
02864          NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
02865       ENDIF 
02866       RETURN
02867       END
02868 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02869 C
02870       SUBROUTINE XPNAXEZ(NEWAXEX,NEWAXEY,I1,I2,J1,J2,AX,AY,NI,NJ)
02871       USE mod_kinds_oasis
02872       IMPLICIT NONE
02873       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
02874       REAL(kind=ip_realwp_p) NEWAXEX(I1:I2)
02875       REAL(kind=ip_realwp_p) NEWAXEY(J1:J2)
02876       REAL(kind=ip_realwp_p) AX(NI)
02877       REAL(kind=ip_realwp_p) AY(NJ)
02878       INTEGER (kind=ip_intwp_p) I,J
02879       DO 23000 I=1,NI
02880          NEWAXEX(I) = AX(I)
02881 23000 CONTINUE 
02882       NEWAXEX(0) = NEWAXEX(NI) - AX(NI)
02883       NEWAXEX(NI+1) = NEWAXEX(1) + AX(NI)
02884       NEWAXEX(NI+2) = NEWAXEX(2) + AX(NI)
02885       DO 23002 J=1,NJ
02886          NEWAXEY(J) = AY(J)
02887 23002 CONTINUE 
02888       NEWAXEY(0 ) = -180.0 - NEWAXEY(1)
02889       NEWAXEY(-1) = -180.0 - NEWAXEY(2)
02890       NEWAXEY(NJ+1) = 180.0 - NEWAXEY(NJ)
02891       NEWAXEY(NJ+2) = 180.0 - NEWAXEY(NJ-1)
02892 cLT
02893          WRITE(*,*) NEWAXEY(0 ), NEWAXEY(-1), NEWAXEY(NJ+1), 
02894      $       NEWAXEY(NJ+2)
02895       RETURN
02896       END
02897 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02898 C
02899       SUBROUTINE XPNAXEY(NEWAXEX,NEWAXEY,I1,I2,J1,J2,AX,AY,NI,NJ)
02900       USE mod_kinds_oasis
02901       IMPLICIT NONE
02902       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
02903       REAL(kind=ip_realwp_p) NEWAXEX(NI)
02904       REAL(kind=ip_realwp_p) NEWAXEY(NJ)
02905       REAL(kind=ip_realwp_p) AX(NI)
02906       REAL(kind=ip_realwp_p) AY(NJ)
02907       INTEGER (kind=ip_intwp_p) I,J
02908       DO 23000 I=1,NI
02909          NEWAXEX(I) = AX(I)
02910 23000 CONTINUE 
02911       DO 23002 J=1,NJ
02912          NEWAXEY(J) = AY(J)
02913 23002 CONTINUE 
02914       RETURN
02915       END
02916 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02917 C
02918       SUBROUTINE XPNCOF(I1,I2,J1,J2,NI,NJ,GRTYP,GRREF,IG1,IG2,IG3,
02919      $IG4,SYM,AX,AY)
02920       USE mod_kinds_oasis
02921       IMPLICIT NONE
02922       INTEGER (kind=ip_intwp_p) I1,J1,I2,J2,NI,NJ,IG1,IG2,IG3,IG4
02923       LOGICAL SYM
02924       CHARACTER*1 GRTYP, GRREF
02925       REAL(kind=ip_realwp_p) AX(NI),AY(NJ)
02926       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
02927       PARAMETER (GLOBAL = 0)
02928       PARAMETER (NORD   = 1)
02929       PARAMETER (SUD   = 2)
02930       PARAMETER (SUDNORD= 0)
02931       PARAMETER (NORDSUD= 1)
02932       IF( (GRTYP.EQ.'L'.OR.GRTYP.EQ.'N'.OR.GRTYP.EQ.'S'.OR.GRTYP.EQ.
02933      $'Y'))THEN
02934          I1 = 1
02935          I2 = NI
02936          J1 = 1
02937          J2 = NJ
02938          RETURN
02939       ENDIF 
02940       IF ( (GRTYP .EQ. 'Z')) THEN
02941          I1 = 0
02942          I2 = NI + 2
02943          J1 = -1
02944          J2 = NJ + 2
02945       ENDIF         
02946       IF( (GRTYP.EQ.'A'.OR. GRTYP.EQ.'G'))THEN
02947          I1 = 0
02948          I2 = NI + 2
02949          IF( (IG1.EQ.GLOBAL))THEN
02950             J1 = -1
02951             J2 = NJ + 2
02952          ELSE 
02953             IF( (IG1.EQ.NORD))THEN
02954                J1 = -NJ - 1
02955                J2 =  NJ + 2
02956             ELSE 
02957                J1 = -1
02958                J2 =  2 * NJ + 2
02959             ENDIF 
02960          ENDIF 
02961          RETURN
02962       ENDIF 
02963       IF( (GRTYP.EQ.'B'))THEN
02964          I1 = 0
02965          I2 = NI + 1
02966          IF( (IG1.EQ.GLOBAL))THEN
02967             J1 = 0
02968             J2 = NJ + 1
02969          ELSE 
02970             IF( (IG1.EQ.NORD))THEN
02971                J1 = -NJ + 1
02972                J2 =  NJ + 1
02973             ELSE 
02974                J1 =  0
02975                J2 =  2 * NJ
02976             ENDIF 
02977          ENDIF 
02978       ENDIF 
02979       RETURN
02980       END
02981 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02982 C
02983       SUBROUTINE XPNGD(ZOUT,I1,I2,J1,J2,ZI,NI,NJ,GRTYP,IG1,IG2,IG3,
02984      $IG4,SYM,VECT)
02985       USE mod_kinds_oasis
02986       IMPLICIT NONE
02987       INTEGER (kind=ip_intwp_p) GLOBAL, NORD, SUD, SUDNORD, NORDSUD
02988       PARAMETER (GLOBAL = 0)
02989       PARAMETER (NORD   = 1)
02990       PARAMETER (SUD   = 2)
02991       PARAMETER (SUDNORD= 0)
02992       PARAMETER (NORDSUD= 1)
02993       EXTERNAL PERMUT
02994       INTEGER (kind=ip_intwp_p) I1,I2,J1,J2,NI,NJ
02995       INTEGER (kind=ip_intwp_p) IG1, IG2, IG3, IG4
02996       CHARACTER*1 GRTYP
02997       INTEGER (kind=ip_intwp_p) I,J
02998       REAL(kind=ip_realwp_p) ZOUT(I1:I2,J1:J2)
02999       REAL(kind=ip_realwp_p) ZI(NI,NJ)
03000       LOGICAL SYM,VECT
03001       REAL(kind=ip_realwp_p) SIGN
03002       IF( (VECT))THEN
03003          SIGN = -1.0
03004       ELSE 
03005          SIGN = 1.0
03006       ENDIF 
03007       IF( (GRTYP.EQ. 'L'.OR. GRTYP.EQ. 'N'.OR. GRTYP.EQ. 'S'.OR.
03008      $ GRTYP.EQ. 'Y'))THEN
03009          DO 23004 J=1,NJ
03010             DO 23006 I=1,NI
03011                ZOUT(I,J) = ZI(I,J)
03012 23006       CONTINUE 
03013 23004    CONTINUE 
03014          RETURN
03015       ENDIF
03016       IF( (GRTYP.EQ. 'Z') ) THEN
03017           DO J=1,NJ
03018             DO I=1,NI
03019               ZOUT(I,J) = ZI(I,J)
03020             END DO 
03021           END DO
03022           DO J=1,NJ
03023             ZOUT(0,J) = ZI(NI,J)
03024             ZOUT(NI+1,J) = ZOUT(1,J)
03025             ZOUT(NI+2,J) = ZOUT(2,J)
03026           END DO 
03027           DO I=0,NI/2
03028             ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
03029             ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
03030             ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03031             ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03032             ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
03033             ZOUT(I+NI/2,-1)= SIGN * ZOUT(I, 2)
03034             ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
03035             ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
03036           END DO 
03037           DO I=1,2
03038             ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
03039             ZOUT(NI+I,-1)= SIGN * ZOUT(I+NI/2, 2)
03040             ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03041             ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03042           END DO 
03043       ENDIF
03044       IF( (GRTYP.EQ.'A'.OR. GRTYP.EQ.'G'))THEN
03045          IF( (IG2.EQ. NORDSUD))THEN
03046             CALL PERMUT(ZI, NI, NJ)
03047          ENDIF 
03048          DO 23012 J=1,NJ
03049             DO 23014 I=1,NI
03050                ZOUT(I,J) = ZI(I,J)
03051 23014       CONTINUE 
03052 23012    CONTINUE 
03053          DO 23016 J=1,NJ
03054             ZOUT(0,J) = ZI(NI,J)
03055             ZOUT(NI+1,J) = ZOUT(1,J)
03056             ZOUT(NI+2,J) = ZOUT(2,J)
03057 23016    CONTINUE 
03058          IF( (IG1.EQ. GLOBAL))THEN
03059             DO 23020 I=0,NI/2
03060                ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
03061                ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
03062                ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03063                ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03064                ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
03065                ZOUT(I+NI/2,-1)= SIGN * ZOUT(I, 2)
03066                ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
03067                ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
03068 23020       CONTINUE 
03069             DO 23022 I=1,2
03070                ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
03071                ZOUT(NI+I,-1)= SIGN * ZOUT(I+NI/2, 2)
03072                ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03073                ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03074 23022       CONTINUE 
03075          ENDIF 
03076          IF( (IG1.EQ. NORD))THEN
03077             IF( (SYM))THEN
03078                DO 23028 J=-NJ+1, 0
03079                   DO 23030 I=0,NI+2
03080                      ZOUT(I,J) = ZOUT(I, -J+1)
03081 23030             CONTINUE 
03082 23028          CONTINUE 
03083             ELSE 
03084                DO 23032 J=-NJ+1, 0
03085                   DO 23034 I=0,NI+2
03086                      ZOUT(I,J) = -ZOUT(I, -J+1)
03087 23034             CONTINUE 
03088 23032          CONTINUE 
03089             ENDIF 
03090             DO 23036 I=0,NI/2
03091                ZOUT(I,-NJ) = SIGN * ZOUT(I+NI/2, -NJ+1)
03092                ZOUT(I,-NJ-1)= SIGN * ZOUT(I+NI/2,-NJ+2)
03093                ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03094                ZOUT(I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03095                ZOUT(I+NI/2,-NJ) = SIGN * ZOUT(I, -NJ+1)
03096                ZOUT(I+NI/2,-NJ-1)= SIGN * ZOUT(I, -NJ+2)
03097                ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ)
03098                ZOUT(I+NI/2,NJ+2)=SIGN * ZOUT(I, NJ-1)
03099 23036       CONTINUE 
03100             DO 23038 I=1,2
03101                ZOUT(NI+I,-NJ) = SIGN * ZOUT(NI/2+I,-NJ+1)
03102                ZOUT(NI+I,-NJ-1)= SIGN * ZOUT(I+NI/2, -NJ+2)
03103                ZOUT(NI+I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ)
03104                ZOUT(NI+I,NJ+2)=SIGN * ZOUT(I+NI/2, NJ-1)
03105 
03106 23038       CONTINUE 
03107          ENDIF 
03108          IF( (IG1.EQ. SUD))THEN
03109             IF( (SYM))THEN
03110                DO 23044 J=1,NJ
03111                   DO 23046 I=0,NI+2
03112                      ZOUT(I,NJ+J) = ZOUT(I, NJ-J+1)
03113 23046             CONTINUE 
03114 23044          CONTINUE 
03115             ELSE 
03116                DO 23048 J=1,NJ
03117                   DO 23050 I=0,NI+2
03118                      ZOUT(I,NJ+J) = -ZOUT(I,NJ-J+1)
03119 23050             CONTINUE 
03120 23048          CONTINUE 
03121             ENDIF 
03122             DO 23052 I=0,NI/2
03123                ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 1)
03124                ZOUT(I,-1)= SIGN * ZOUT(I+NI/2, 2)
03125                ZOUT(I,2*NJ+1)=SIGN * ZOUT(I+NI/2, 2*NJ)
03126                ZOUT(I,2*NJ+2)=SIGN * ZOUT(I+NI/2, 2*NJ-1)
03127                ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 1)
03128                ZOUT(I+NI/2,-1)= SIGN * ZOUT(I,2)
03129                ZOUT(I+NI/2,2*NJ+1)=SIGN * ZOUT(I, 2*NJ)
03130                ZOUT(I+NI/2,2*NJ+2)=SIGN * ZOUT(I, 2*NJ-1)
03131 23052       CONTINUE 
03132             DO 23054 I=1,2
03133                ZOUT(NI+I,0) = SIGN * ZOUT(NI/2+I,1)
03134                ZOUT(NI+I,1)= SIGN * ZOUT(I+NI/2, 2)
03135                ZOUT(NI+I,2*NJ+1)=SIGN * ZOUT(I+NI/2, 2*NJ)
03136                ZOUT(NI+I,2*NJ+2)=SIGN * ZOUT(I+NI/2, 2*NJ-1)
03137 23054       CONTINUE 
03138          ENDIF 
03139       ENDIF 
03140       IF( (GRTYP.EQ.'B'))THEN
03141          IF( (IG2.EQ. NORDSUD))THEN
03142             CALL PERMUT(ZI, NI, NJ)
03143          ENDIF 
03144          DO 23060 J=1,NJ
03145             DO 23062 I=1,NI
03146                ZOUT(I,J) = ZI(I,J)
03147 23062       CONTINUE 
03148 23060    CONTINUE 
03149          DO 23064 J=1,NJ
03150             ZOUT(0,J) = ZI(NI-1,J)
03151             ZOUT(NI+1,J) = ZOUT(2,J)
03152 23064    CONTINUE 
03153          IF( (IG1.EQ. GLOBAL))THEN
03154             DO 23068 I=0,NI/2+1
03155                ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 2)
03156                ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ-1)
03157                ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 2)
03158                ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ-1)
03159 23068       CONTINUE 
03160             ZOUT(NI+1,0) = SIGN * ZOUT(NI/2+1,2)
03161             ZOUT(NI+1,NJ+1)=SIGN * ZOUT(NI/2+1, NJ-1)
03162          ENDIF 
03163          IF( (IG1.EQ. NORD))THEN
03164             IF( (SYM))THEN
03165                DO 23074 J=-NJ+2, 0
03166                   DO 23076 I=0,NI+1
03167                      ZOUT(I,J) = ZOUT(I, -J+2)
03168 23076             CONTINUE 
03169 23074          CONTINUE 
03170             ELSE 
03171                DO 23078 J=-NJ+2, 0
03172                   DO 23080 I=0,NI+1
03173                      ZOUT(I,J) = -ZOUT(I, -J+2)
03174 23080             CONTINUE 
03175 23078          CONTINUE 
03176             ENDIF 
03177             DO 23082 I=0,NI/2+1
03178                ZOUT(I,-NJ+1) = SIGN * ZOUT(I+NI/2, -NJ+3)
03179                ZOUT(I,NJ+1)=SIGN * ZOUT(I+NI/2, NJ-1)
03180                ZOUT(I+NI/2,-NJ+1) = SIGN * ZOUT(I, -NJ+3)
03181                ZOUT(I+NI/2,NJ+1)=SIGN * ZOUT(I, NJ-1)
03182 23082       CONTINUE 
03183             ZOUT(NI+1,-NJ+1) = SIGN * ZOUT(NI/2+1,-NJ+3)
03184             ZOUT(NI+1,NJ+1)=SIGN * ZOUT(1+NI/2, NJ-1)
03185          ENDIF 
03186          IF( (IG1.EQ. SUD))THEN
03187             IF( (SYM))THEN
03188                DO 23088 J=1,NJ-1
03189                   DO 23090 I=0,NI+1
03190                      ZOUT(I,NJ+J) = ZOUT(I, NJ-J)
03191 23090             CONTINUE 
03192 23088          CONTINUE 
03193             ELSE 
03194                DO 23092 J=1,NJ-1
03195                   DO 23094 I=0,NI+1
03196                      ZOUT(I,NJ+J) = -ZOUT(I,NJ-J)
03197 23094             CONTINUE 
03198 23092          CONTINUE 
03199             ENDIF 
03200             DO 23096 I=0,NI/2+1
03201                ZOUT(I,0) = SIGN * ZOUT(I+NI/2, 2)
03202                ZOUT(I,2*NJ)=SIGN * ZOUT(I+NI/2, 2*NJ-2)
03203                ZOUT(I+NI/2,0) = SIGN * ZOUT(I, 2)
03204                ZOUT(I+NI/2,2*NJ)=SIGN * ZOUT(I, 2*NJ-2)
03205 23096       CONTINUE 
03206             ZOUT(NI+1,0) = SIGN * ZOUT(NI/2+1,2)
03207             ZOUT(NI+1,2*NJ)=SIGN * ZOUT(1+NI/2, 2*NJ-2)
03208          ENDIF 
03209       ENDIF 
03210       RETURN
03211       END
03212 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03213 C
03214       SUBROUTINE LLFXY(DLAT,DLON,X,Y,D60,DGRW,NHEM) 
03215       USE mod_kinds_oasis
03216 C
03217 C AUTHOR   - J.D. HENDERSON  -  FEB 75
03218 C 
03219 C OBJECT(LLFXY)
03220 C         - COMPUTES LATITUDE AND LONGITUDE OF A POINT IN A POLAR
03221 C           STEREOGRAPHIC GRID FROM CO-ORDINATES IN THE GRID MEASURED
03222 C           FROM THE POLE.
03223 C
03224 C USAGE    - CALL LLFXY(DLAT,DLON,X,Y,D60,DGRW,NHEM)
03225 C
03226 C ARGUMENTS
03227 C   OUT   - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N).
03228 C         - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E).
03229 C   IN    - X    - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
03230 C                  AS ORIGIN
03231 C         - Y    - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
03232 C                  AS ORIGIN
03233 C         - D60  - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC
03234 C                  GRID AT 60 DEGREES
03235 C         - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO
03236 C                  THE GRID (IN DEGREES)
03237 C         - NHEM - 1 FOR NORTHERN HEMISPHERE
03238 C                  2 FOR SOUTHERN HEMISPHERE 
03239 C
03240 C NOTES    - THE COMPANION ROUTINE XYFLL, WHICH COMPUTES THE GRID
03241 C           CO-ORDINATES GIVEN THE LATITUDE AND LONGITUDE, IS ALSO
03242 C           AVAILABLE.
03243 C
03244 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.
03245 C     * RDTODG = 180/PIE, DGTORD = PIE/180                                     
03246 C
03247       DATA PIE    /3.14159265358979323846/
03248       DATA RDTODG /57.295779513082/
03249       DATA DGTORD /1.7453292519943E-2/
03250 C* ------------------------------------------------------------------------
03251       RE=1.866025*6.371E+6/D60 
03252       RE2=RE**2
03253 C     * IF POINT IS AT POLE SET COORD TO (0.,90.). 
03254       DLAT=90.                 
03255       DLON=0.
03256       IF(X.EQ.0. .AND. Y.EQ.0.) GO TO 39
03257 C                                                 
03258 C     * CALCULATE LONGITUDE IN MAP COORDINATES.
03259 C
03260       IF(X.EQ.0.) DLON=SIGN(90.,Y)
03261       IF(X.NE.0.) DLON=ATAN(Y/X)*RDTODG
03262       IF(X.LT.0.) DLON=DLON+SIGN(180.,Y)
03263 C                  
03264 C     * ADJUST LONGITUDE FOR GRID ORIENTATION. 
03265 C   
03266       DLON=DLON-DGRW 
03267       IF(DLON.GT.+180.) DLON=DLON-360.
03268       IF(DLON.LT.-180.) DLON=DLON+360. 
03269 C              
03270 C     * CALCULATE LATITUDE. 
03271 C                 
03272       R2=X**2+Y**2
03273       DLAT=(RE2-R2)/(RE2+R2) 
03274       DLAT= ASIN(DLAT)*RDTODG 
03275 C   
03276 C     * CHANGE SIGNS IF IN SOUTHERN HEMISPHERE.
03277 C
03278    39 IF(NHEM.EQ.2) DLAT=-DLAT
03279       IF(NHEM.EQ.2) DLON=-DLON 
03280 C 
03281       RETURN
03282       END 
03283 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
03284 C 
03285       SUBROUTINE XYFLL(X,Y,DLAT,DLON,D60,DGRW,NHEM)
03286       USE mod_kinds_oasis
03287 C  
03288 C AUTHOR   - J.D. HENDERSON  -  FEB 75
03289 C    
03290 C OBJECT(XYFLL) 
03291 C         - COMPUTES THE GRID CO-ORDINATES MEASURED FROM THE POLE OF A 
03292 C           POINT, GIVEN THE LATITUDE AND LONGITUDE IN DEGREES.
03293 C                                                                              
03294 C USAGE    - CALL XYFLL(X,Y,DLAT,DLON,D60,DGRW,NHEM) 
03295 C             
03296 C ARGUMENTS   
03297 C   OUT   - X    - X-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE 
03298 C                  AS ORIGIN    
03299 C                  AS ORIGIN
03300 C         - Y    - Y-CO-ORDINATE OF THE POINT AS MEASURED WITH POLE
03301 C                  AS ORIGIN      
03302 C   IN    - DLAT - LATITUDE IN DEGREES (-90 TO +90, POSITIVE N) 
03303 C         - DLON - LONGITUDE IN DEGREES (-180 TO +180, POSITIVE E) 
03304 C         - D60  - GRID LENGTH (IN METRES) OF THE POLAR STEREOGRAPHIC  
03305 C                  GRID AT 60 DEGREES
03306 C         - DGRW - ORIENTATION OF GREENWICH MERIDIAN WITH RESPECT TO 
03307 C                  THE GRID (IN DEGREES)      
03308 C         - NHEM - 1 FOR NORTHERN HEMISPHERE    
03309 C                  2 FOR SOUTHERN HEMISPHERE  
03310 C 
03311 C NOTES    - THE COMPANION ROUTINE LLFXY, WHICH COMPUTES THE LATITUDE   
03312 C           AND LONGITUDE GIVEN THE GRID-COORDINATES, IS ALSO AVAILABLE.
03313 C                                                                
03314 C---------------------------------------------------------------------------
03315 C                                                     
03316 C     * 1.866025=(1+SIN60),   6.371E+6=EARTH RADIUS IN METERS.    
03317 C      
03318 C RDTODG = 180/PIE, DGTORD = PIE/180    
03319       DATA PIE    /3.14159265358979323846/       
03320       DATA RDTODG /57.295779513082/ 
03321       DATA DGTORD /1.7453292519943E-2/  
03322 C 
03323       RE=1.866025*6.371E+6/D60    
03324 C   
03325       GLON=DLON  
03326       IF(NHEM.EQ.2) GLON=-DLON  
03327       GLAT=DLAT     
03328       IF(NHEM.EQ.2) GLAT=-DLAT   
03329 C    
03330       RLON=DGTORD*(GLON+DGRW)      
03331       RLAT=DGTORD*GLAT   
03332       SINLAT=SIN(RLAT)    
03333       R=RE*SQRT((1.-SINLAT)/(1.+SINLAT)) 
03334       X=R*COS(RLON)        
03335       Y=R*SIN(RLON)                                                            
03336 C                                                                              
03337       RETURN                                              
03338       END                                  
03339 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03340 C
03341       SUBROUTINE ORDLEG(SBX,COA,IR) 
03342       USE mod_kinds_oasis
03343 C     
03344 C AUTEUR   - D. ROBERTSON 
03345 C  
03346 C OBJET(ORDLEG)   
03347 C         - THIS ROUTINE IS A SUBSET OF BELOUSOVS ALGORITHM  
03348 C           USED TO CALCULATE ORDINARY LEGENDRE POLYNOMIALS.  
03349 C                                        
03350 C USAGE    - CALL ORDLEG(SBX,COA,IR)
03351 C   
03352 C ARGUMENTS        
03353 C   OUT   - SBX  - LEGENDRE POLYNOMIAL EVALUATED AT COA        
03354 C   IN    - COA - COSINE OF COLATITUDE   
03355 C   IN    - IR  - WAVE NUMBER  
03356 C                                         
03357 C* -------------------------------------------------------------------------
03358 C               
03359 C              
03360 C RDTODG = 180/PIE, DGTORD = PIE/180     
03361       DATA PIE    /3.14159265358979323846/ 
03362       DATA RDTODG /57.295779513082/   
03363       DATA DGTORD /1.7453292519943E-2/    
03364 C      
03365       PI = PIE                        
03366       SQR2=SQRT(2.)                  
03367       IRPP = IR + 1      
03368       IRPPM = IRPP - 1  
03369       DELTA =   ACOS(COA)  
03370       SIA =  SIN(DELTA)   
03371 C                              
03372       THETA=DELTA           
03373       C1=SQR2    
03374 C                                     
03375       DO 20 N=1,IRPPM                              
03376       FN=FLOAT(N)                           
03377       FN2=2.*FN                             
03378       FN2SQ=FN2*FN2                  
03379       C1=C1* SQRT(1.0-1.0/FN2SQ)     
03380    20 CONTINUE                        
03381 C             
03382       N=IRPPM              
03383       ANG=FN*THETA                            
03384       S1=0.0                    
03385       C4=1.0                                          
03386       A=-1.0                               
03387       B=0.0                                       
03388       N1=N+1          
03389 C                        
03390       DO 27 KK=1,N1,2                 
03391       K=KK-1                   
03392       IF (K.EQ.N) C4=0.5*C4           
03393       S1=S1+C4* COS(ANG)         
03394       A=A+2.0        
03395       B=B+1.0         
03396       FK=FLOAT(K)      
03397       ANG=THETA*(FN-FK-2.0)      
03398       C4=(A*(FN-B+1.0)/(B*(FN2-A)))*C4  
03399    27 CONTINUE   
03400 C                               
03401       SBX=S1*C1                                  
03402 C       
03403       RETURN     
03404       END 
03405 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03406 C
03407       LOGICAL FUNCTION VALIDE(NOM, ICHK, MIN, MAX)
03408 C
03409 C AUTEUR   - P. SARRAZIN  - AVRIL 82
03410 C
03411 C OBJET(VALIDE)
03412 C         - VERIFIER SI ICHK => MIN OU =< MAX
03413 C         - MESSAGE SI ICHK N'EST PAS DANS LES LIMITES
03414 C
03415 C ARGUMENTS
03416 C   IN    - NOM   - NOM DE LA VARIABLE EMPLOYE PAR LA ROUTINE
03417 C   IN    - ICHK  - VALEUR DE LA VARIABLE POUR VERIFICATION
03418 C   IN    - MIN   - VALEUR MINIMUM DE ICHK
03419 C   IN    - MAX   - VALEUR MAXIMUM DE ICHK
03420 C
03421 C MODULES - SCINT - UVINT
03422 C
03423 C* ----------------------------------------------------------------------
03424       IF( (ICHK.LT.MIN .OR. ICHK.GT.MAX))THEN
03425          WRITE(6,600) NOM,ICHK,MIN,MAX
03426       ENDIF 
03427 600   FORMAT("Bad value for",I10,"VALEUR=",I10,"MINIMUM=",
03428      $   I10,"MAXIMUM=",I10)
03429       VALIDE = (ICHK.GE.MIN) .AND. (ICHK.LE.MAX)
03430       RETURN
03431       END
03432 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03433 C
03434       SUBROUTINE PERMUT(Z,NI,NJ)
03435       USE mod_kinds_oasis
03436       REAL(kind=ip_realwp_p) Z(NI,NJ)
03437 C
03438 C AUTEUR   - M. VALIN  -  FEV 82
03439 C
03440 C OBJET(PERMUT)
03441 C         - ROTATION D'UNE MATRICE AUTOUR DE LA LIGNE DU MILIEU.
03442 C           CETTE ROUTINE EST UTILISEE POUR RE=ARRANGER LES
03443 C           DONNEES DANS UN CHAMP. EX: POUR CONTOURER, ON A
03444 C           DES DONNEES DANS L'ORDRE SUIVANT, DU SUD AU NORD
03445 C           ALORS QUE POUR L'INTERPOLATION ELLES DOIVENT ETRE
03446 C           ETALEES DU NORD AU SUD.
03447 C
03448 C APPEL    - CALL PERMUT(Z,NI,NJ)
03449 C
03450 C ARGUMENTS
03451 C IN/OUT  - Z  - CHAMP QUI SUBIT LA ROTATION
03452 C   IN    - NI - PREMIERE DIMENSION DE Z
03453 C   IN    - NJ - DEUXIEME DIMENSION DE Z
03454 C* -----------------------------------------------------------------------
03455 C
03456       NCC = NJ/2
03457       DO 23000 J=1,NCC
03458          DO 23002 I=1,NI
03459             T = Z(I,NJ+1-J)
03460             Z(I,NJ+1-J) = Z(I,J)
03461             Z(I,J) = T
03462 23002    CONTINUE 
03463 23000 CONTINUE 
03464       RETURN
03465       END
03466 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
03467 C    
03468       SUBROUTINE CIGAXG(CGTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
03469       USE mod_kinds_oasis
03470       CHARACTER*1 CGTYP
03471 C
03472 C AUTEUR   - M. VALIN  -  FEV 82
03473 C
03474 C
03475 C OBJET(CIGAXG)
03476 C         - PASSE DES PARAMETRES (ENTIERS) DESCRIPTEURS DE GRILLE
03477 C           AUX PARAMETRES REELS.
03478 C
03479 C APPEL    - CALL IGAXG(CGTYP,XG1,XG2,XG3,XG4,IG1,IG2,IG3,IG4)
03480 C
03481 C ARGUMENTS
03482 C   IN    - CGTYP - TYPE DE GRILLE (VOIR OUVRIR)
03483 C   OUT   - XG1   - ** DESCRIPTEUR DE GRILLE (REEL),
03484 C   OUT   - XG2   -    IGTYP = 'N', PI, PJ, D60, DGRW
03485 C   OUT   - XG3   -    IGTYP = 'L', LAT0, LON0, DLAT, DLON,
03486 C   OUT   - XG4   -    IGTYP = 'A', 'B', 'G', XG1 = 0. GLOBAL,
03487 C                                                 = 1. NORD
03488 C                                                 = 2. SUD **
03489 C   IN    - IG1   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03490 C   IN    - IG2   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03491 C   IN    - IG3   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03492 C   IN    - IG4   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03493 C
03494 CMESSAGES - "ERREUR, MAUVAISE SPECIFICATION DE GRILLE, (TYPE) (IGAXG)"
03495 C
03496 C-------------------------------------------------------------------
03497 C
03498 
03499       IF( ((CGTYP.EQ. 'N') .OR. (CGTYP.EQ.'S')))THEN
03500          IF((IG4.LT. 32768))THEN
03501             XG1 = IG2 * 0.1
03502             XG2 = IG1 * 0.1
03503             XG3 = IG4 * 100.
03504             XG4 = IG3 * 0.01
03505          ELSE 
03506             JG3 = IG3
03507             JG4 = IG4
03508             JG4 = JG4 - 32768
03509             XG3 = IG1 * 100.
03510             IF((IG3 .GT. 32767))THEN
03511                XG3 = XG3 * 10.
03512                JG3 = JG3 - 32768
03513             ENDIF 
03514             XG4 = IG2 * .1
03515             IF((JG4.GT. 16383))THEN
03516                XG4 = 360. - XG4
03517                JG4 = JG4 - 16384
03518             ENDIF 
03519             DLAT = 90. -(JG4*180./16383.)
03520             DLON = (JG3*360./32767.)
03521             IHEM = 1
03522             IF(('S'.EQ.CGTYP))THEN
03523                IHEM = 2
03524             ENDIF 
03525             CALL XYFLL(XG1,XG2,DLAT,DLON,XG3,XG4,IHEM)
03526             XG1 = 1.0 - XG1
03527             XG2 = 1.0 - XG2
03528          ENDIF 
03529       ELSE 
03530          IF((CGTYP.EQ. 'C'))THEN
03531             XG1 = IG3 * 0.01 - 90.
03532             XG2 = IG4 * 0.01
03533             XG3 = 180. / IG1
03534             XG4 = 360. / IG2
03535          ELSE 
03536             IF( ((CGTYP.EQ. 'A') .OR. (CGTYP.EQ. 'B') .OR.   (CGTYP
03537      $      .EQ. 'G')))THEN
03538                XG1 = IG1
03539                XG2 = IG2
03540                XG3 = 0.
03541                XG4 = 0.
03542             ELSE 
03543                IF((CGTYP.EQ. 'L'))THEN
03544                   XG1 = IG3 * 0.01 - 90.
03545                   XG2 = IG4 * 0.01
03546                   XG3 = IG1 * 0.01
03547                   XG4 = IG2 * 0.01
03548                ELSE 
03549                   IF((CGTYP.EQ. 'H'))THEN
03550                      XG1 = IG3
03551                      XG2 = .01*IG4 - 90.
03552                      XG3 = 500*IG2
03553                      XG4 = IG1*.2
03554                   ELSE 
03555                      WRITE(6,600)
03556                   ENDIF 
03557                ENDIF 
03558             ENDIF 
03559          ENDIF 
03560       ENDIF 
03561 600   FORMAT(' Error bad grid specification , (TYPE)'
03562      $,'(IGAXG)')
03563       RETURN
03564       END
03565 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03566 C
03567       SUBROUTINE CXGAIG(CGTYP,IG1,IG2,IG3,IG4,XG1,XG2,XG3,XG4)
03568       USE mod_kinds_oasis
03569       LOGICAL VALIDE
03570       EXTERNAL VALIDE
03571       CHARACTER * 1 CGTYP
03572       LOGICAL LFLAG
03573 C
03574 C AUTEUR   - M. VALIN  -  FEV 82
03575 C
03576 C OBJET(XGAIG)
03577 C         - PASSE DES PARAMETRES (REELS) DESCRIPTEURS DE GRILLE
03578 C           AUX PARAMETRES ENTIERS.
03579 C
03580 C APPEL    - CALL XGAIG(CGTYP,IG1,IG2,IG3,IG4,XG1,XG2,XG3,XG4)
03581 C
03582 C ARGUMENTS
03583 C   IN    - CGTYP - TYPE DE GRILLE (VOIR OUVRIR)
03584 C   OUT   - IG1   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03585 C   OUT   - IG2   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03586 C   OUT   - IG3   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03587 C   OUT   - IG4   - DESCRIPTEUR DE GRILLE (ENTIER) VOIR OUVRIR
03588 C   IN    - XG1   - ** DESCRIPTEUR DE GRILLE (REEL),
03589 C   IN    - XG2   -    IGTYP = 'N', PI, PJ, D60, DGRW
03590 C   IN    - XG3   -    IGTYP = 'L', LAT0, LON0, DLAT, DLON,
03591 C   IN    - XG4   -    IGTYP = 'A', 'B', 'G', XG1 = 0, GLOBAL
03592 C                                                 = 1, NORD
03593 C                                                 = 2, SUD **
03594 C
03595 C MESSAGES - "ERREUR DANS LA DESCRIPTION DE LA GRILLE (IG1) (XGAIG)"
03596 C           "ERREUR, MAUVAISE SPECIFICATION (LAT0) (XGAIG)"
03597 C           "ERREUR, GRILLE INCONNUE (TYPE) (XGAIG)"
03598 C
03599 C* ------------------------------------------------------------------
03600 C
03601 C
03602       IF( (CGTYP.EQ. 'N' .OR. CGTYP.EQ.'S'))THEN
03603          IG1 = NINT(XG2 * 10.)
03604          IG2 = NINT(XG1 * 10.)
03605          IG3 = NINT(XG4 * 100.)
03606          IG4 = NINT(XG3 * 0.01)
03607 23002    IF( (IG3.LT. 0))THEN
03608             IG3 = IG3 + 36000
03609             GOTO 23002
03610          ENDIF 
03611          IF((IG1.LT.0.OR.IG2.LT.0.OR.IG1.GT.2047.OR.IG2.GT.2047.OR.
03612      $   IG4.GT.32000))THEN
03613             IG1 = 0
03614             IG2 = 0
03615             IG3 = 0
03616             IG4 = 32768
03617             IF((XG3.GT. 204700))THEN
03618                IG3 = 32768
03619                IG1 = NINT(XG3*.001)
03620             ELSE 
03621                IG3 = 0
03622                IG1 = NINT(XG3*.01)
03623             ENDIF 
03624             IG2 = NINT(XG4*10)
03625             IF((IG2.LT.0))THEN
03626                IG2 = ABS(IG2)
03627                IG4 = IG4 + 16384
03628             ENDIF 
03629             IF((IG2.GT.1800))THEN
03630                IG2 = ABS(IG2 - 3600)
03631                IG4 = IG4 + 16384
03632             ENDIF 
03633             IHEM = 1
03634             IF(('S'.EQ.CGTYP))THEN
03635                IHEM = 2
03636             ENDIF 
03637             CALL LLFXY(DLAT,DLON,1.-XG1,1.-XG2,XG3,XG4,IHEM)
03638             DLAT = 90. - DLAT
03639             IF((DLON.LT.0))THEN
03640                DLON = DLON + 360.
03641             ENDIF 
03642             IG3 = IG3 + NINT(DLON*32767./360.)
03643             IG4 = IG4 + NINT(DLAT*16383./180.)
03644          ENDIF 
03645       ELSE 
03646          IF( (CGTYP.EQ. 'A' .OR. CGTYP.EQ. 'B' .OR.   CGTYP.EQ. 'G')
03647      $   )THEN
03648             IG1 = XG1
03649             IG2 = XG2
03650             IG3 = 0
03651             IG4 = 0
03652             LFLAG=VALIDE("IG1",IG1,0,2)
03653             LFLAG=VALIDE("IG2",IG2,0,1)
03654          ELSE 
03655             IF((CGTYP.EQ. 'C'))THEN
03656                IG1 = 180. / XG3 + 0.5
03657                IG2 = 360. / XG4 + 0.5
03658                IG3 = (90. + XG1) * 100. + 0.5
03659                IG4 = XG2 * 100. + 0.5
03660 23020          IF( (IG4.LT. 0))THEN
03661                   IG4 = IG4 + 36000
03662                   GOTO 23020
03663                ENDIF 
03664                IF( (IG3.LT. 0))THEN
03665                   WRITE(6,601)
03666                ENDIF 
03667             ELSE 
03668                IF( (CGTYP.EQ. 'H'))THEN
03669                   IG1 = NINT(5.*XG4)
03670 23026             IF( (IG1.LT. 0))THEN
03671                      IG1 = IG1 + 1800
03672                      GOTO 23026
03673                   ENDIF 
03674                   IG2 = NINT(.002*XG3)
03675                   IG3 = NINT(XG1)
03676                   IG4 = NINT(100.*(90.+XG2))
03677                ELSE 
03678                   IF( (CGTYP.EQ. 'L'))THEN
03679                      IG1 = XG3 * 100. + 0.5
03680                      IG2 = XG4 * 100. + 0.5
03681                      IG3 = (90. + XG1) * 100. + 0.5
03682                      IG4 = XG2 * 100. + 0.5
03683 23030                IF( (IG4.LT. 0))THEN
03684                         IG4 = IG4 + 36000
03685                         GOTO 23030
03686                      ENDIF 
03687                      IF( (IG3.LT. 0))THEN
03688                         WRITE(6,601)
03689                      ENDIF 
03690                   ELSE 
03691                      WRITE(6,602)
03692                   ENDIF 
03693                ENDIF 
03694             ENDIF 
03695          ENDIF 
03696       ENDIF 
03697 601   FORMAT(' Error, bad specification (LAT0) (XGAIG)')
03698 602   FORMAT(' Error, unknown grid (TYPE) (XGAIG)')
03699       RETURN
03700       END
03701 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 All Data Structures Namespaces Files Functions Variables Defines