      PROGRAM INTRAC2                                                                                                               
C     PROGRAM INTRAC2 (ICTL,       GSINIT,       SSTRAC,       INPUT,           I2                                                  
C    1                                                         OUTPUT,  )       I2                                                  
C    2          TAPE99=ICTL, TAPE1=GSINIT, TAPE2=SSTRAC, TAPE5=INPUT,                                                               
C    3                                                   TAPE6=OUTPUT)                                                              
C     ----------------------------------------------------------------          I2                                                  
C                                                                               I2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       I2                                                  
C     MAY 09/96 - F.MAJAESS (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                
C     FEB 15/94 - F.MAJAESS (REVISE FOR TRUNCATIONS > 99)                                                                           
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII)                                                                   
C     JAN 23/91 - F.MAJAESS (RESTRICT "IBUF" TO "ICOM" COMMON BLOCK)                                                                
C     APR 26/90 - J. DE GRANDPRE.  LIKE INITRAC EXCEPT THAT IT INCLUDES THE                                                         
C                                    POSSIBILITY TO USE N TRACERS AND A                                                             
C                                    DIFFERENT TYPE OF INITIAL PATTERN.                                                             
C                                                                               I2                                                  
CINTRAC2 - SET UP A NUMBER ITRAC OF FIELDS TO BE USED AS A FUNCTION             I1                                                  
C          OF CONCENTRATION IN THE G.C.M..                              2  1 C  I1                                                  
C                                                                               I3                                                  
CAUTHOR  - JEAN DE GRANDPRE                                                     I3                                                  
C                                                                               I3                                                  
CPURPOSE - SET UP "N" TRACER FIELDS TO BE USED AS F(TRACER) IN THE GCM.         I3                                                  
C          THESE FIELDS ARE DEFINED ON HALF-LEVELS, LIKE T AND ES.              I3                                                  
C                                                                               I3                                                  
CINPUT FILES...                                                                 I3                                                  
C                                                                               I3                                                  
C      ICTL   = INITIALIZATION CONTROL DATASET (SEE ICNTRLB).                   I3                                                  
C      GSINIT = SIGMA LEVEL GAUSSIAN GRIDS (SEE INITGSC).                       I3                                                  
C                                                                               I3                                                  
COUTPUT FILE...                                                                 I3                                                  
C                                                                               I3                                                  
C      SSTRAC = SIGMA HALF LEVEL TRACER SPECTRAL FIELDS                         I3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               I5                                                  
C      ITRAC = TYPE OF INITIAL PATTERN OF TRACER CONCENTRATION                  I5                                                  
C              -THE FIRST TRACER IS A DISTURBANCE DECAYING AS 1/R**2 AND        I5                                                  
C               UNIFORM IN THE VERTICAL.                                        I5                                                  
C              -THE VERTICAL PROFILE OF THE SECOND TRACER HAS A QUADRATIQUE     I5                                                  
C               FORM LIKE IN THE HORIZONTAL.                                    I5                                                  
C              -ALL OTHERS ARE IDENTICAL TO THE FIRST TRACER.(MAX=99)           I5                                                  
C     ITRVAR = TRACER VARIABLE MAY BE Q, -1./LN(Q), LN(Q)                       I5                                                  
C      AMP   = AMPLITUDE OF DISTURBANCE AT ITS MAXIMUM.                         I5                                                  
C      X0,Y0 = LONGITUDE AND LATITUDE OF PEAK OF DISTURBANCE.                   I5                                                  
C      WID   = IF (ITRAC.EQ.2), WAVELENGTHS IN LONG. OF DOUBLE COS**2 BLOB,     I5                                                  
C              OTHERWISE,       HALF WIDTH OF DISTURBANCE IN DEGREES.           I5                                                  
C      BKG   = IF (ITRAC.EQ.2) WAVELENGTHS IN LAT.  OF DOUBLE COS**2 OF BLOB,   I5                                                  
C              OTHERWISE, UNIFORM BACKGROUND VALUE OF TRACER ADDED TO PATTERN.  I5                                                  
C                                                                               I5                                                  
CEXAMPLE OF INPUT CARD...                                                       I5                                                  
C                                                                               I5                                                  
C*INTRAC2.    1 RLNQ       20.      235.       40.       20.       10.          I5                                                  
C------------------------------------------------------------------------------                                                     
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMPLEX SSTRAC                                                                                                                
      REAL*8 SL,CL,WL,WOSSL,RAD,ALP,EPSI                                                                                            
      COMMON /SPECSIG/SSTRAC(2145)                                                                                                  
      COMMON /GRIDSIG/TRAC(18528),PS(18528)                                                                                         
      COMMON /WRKSPAC/JBUF(8),WRKS(1092),WRKL(728)                                                                                  
      COMMON /ICOM/IBUF(8),IDAT(18528)                                                                                              
      COMMON /WRKSUP/WRK(18528)                                                                                                     
      COMMON /GAUSSGR/SL(96),CL(96),WL(96),WOSSL(96),                                                                               
     1                RAD(96)                                                                                                       
      COMMON /ALPDATA/ALP(2210),EPSI(2210)                                                                                          
C                                                                                                                                   
      REAL RLON(361),PR(100)                                                                                                        
C                                                                                                                                   
      INTEGER LEV(100),LSR(2,66)                                                                                                    
C                                                                                                                                   
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
      DATA PI/3.141592653589793E0/                                                                                                  
      DATA MAXX/18528/, MAXL/100/                                                                                                   
C-------------------------------------------------------------------------                                                          
      NF=5                                                                                                                          
      CALL JCLPNT(NF,99,1,2,5,6)                                                                                                    
C                                                                                                                                   
C     * READ DATA CARD DIRECTIVES.                                                                                                  
C                                                                                                                                   
      READ(5,5000) ITRAC,ITRVAR,AMP,X0,Y0,WID,BKG                               I4                                                  
      WRITE(6,6000) ITRAC,ITRVAR,AMP,X0,Y0,WID,BKG                                                                                  
      NTRAC = ABS(ITRAC)                                                                                                            
C                                                                                                                                   
C     * GET SPECTRAL RESOLUTION INFORMATION FROM CONTROL FILE ICTL.                                                                 
C                                                                                                                                   
      REWIND 99                                                                                                                     
      READ(99,END=902) LABL                                                                                                         
      READ(99,END=903) LABL,I1,I2,I3,I4,LRLMT                                                                                       
      WRITE(6,6001)LRLMT                                                                                                            
C                                                                                                                                   
C     * GET HALF LEVELS AND ANALYSIS GRID SIZE FROM THE GS TEMP FIELD.                                                              
C                                                                                                                                   
      REWIND 1                                                                                                                      
      CALL FIND(1,NC4TO8("GRID"),-1,NC4TO8("TEMP"),-1,OK)                                                                           
      IF(.NOT.OK) CALL                             XIT('INTRAC2',-1)                                                                
      CALL FILEV(LEV,ILEV,IBUF,-1)                                                                                                  
      IF((ILEV.LT.1).OR.(ILEV.GT.MAXL)) CALL       XIT('INTRAC2',-2)                                                                
      ILG=IBUF(5)-1                                                                                                                 
      JLAT=IBUF(6)                                                                                                                  
      WRITE(6,6005)ILEV,(LEV(L),L=1,ILEV)                                                                                           
      WRITE(6,6006)ILG,JLAT                                                                                                         
C                                                                                                                                   
C     * DECODE LEVELS.                                                                                                              
C                                                                                                                                   
      CALL LVDCODE(PR,LEV,ILEV)                                                                                                     
C                                                                                                                                   
C     * GET THE FIELD OF LN(PS).                                                                                                    
C                                                                                                                                   
      REWIND 1                                                                                                                      
      CALL GETFLD2(1,PS,NC4TO8("GRID"),0,NC4TO8("LNSP"),1,IBUF,MAXX,OK)                                                             
      IF(.NOT.OK) CALL                             XIT('INTRAC2',-3)                                                                
      CALL PRTLAB (IBUF)                                                                                                            
C                                                                                                                                   
C     * INITIALIZE CONSTANTS.                                                                                                       
C                                                                                                                                   
      JLATH=JLAT/2                                                                                                                  
      ILG1 =ILG+1                                                                                                                   
      CALL DIMGT (LSR,LA,LR,LM,KTR,LRLMT)                                                                                           
      CALL EPSCAL(EPSI,LSR,LM)                                                                                                      
      CALL GAUSSG(JLATH,SL,WL,CL,RAD,WOSSL)                                                                                         
      CALL TRIGL (JLATH,SL,WL,CL,RAD,WOSSL)                                                                                         
C----------------------------------------------------------------------                                                             
C     * DEFINE LONGITUDES AND OTHER ANGLES IN RADIANS.                                                                              
C                                                                                                                                   
      DO 250 I=1,ILG1                                                                                                               
        RLON(I)=2.E0*PI*FLOAT(I-1)/FLOAT(ILG)                                                                                       
  250 CONTINUE                                                                                                                      
      WID=2.E0*PI*WID/360.E0                                                                                                        
      X0 =2.E0*PI*X0 /360.E0                                                                                                        
      Y0 =2.E0*PI*Y0 /360.E0                                                                                                        
C                                                                                                                                   
C     * DEFINE THE TRACER FIELD ON ONE LEVEL.                                                                                       
C                                                                                                                                   
      IF(ITRAC.NE.99) THEN                                                                                                          
C                                                                                                                                   
C     ------------------------------------------------------------------                                                            
C     * 1. FOR A DISTURBANCE DECAYING AS 1/R**2 :                 ******                                                            
C     ------------------------------------------------------------------                                                            
C                                                                                                                                   
      DO 300 I=1,ILG                                                                                                                
      DO 300 J=1,JLAT                                                                                                               
      ANGL=ACOS(SIN(RAD(J))*SIN(Y0)+COS(RAD(J))*COS(Y0)*COS(RLON(I)-X0))                                                            
      TRAC((J-1)*ILG1+I) = AMP*WID**2/(WID**2+ANGL**2) + BKG                                                                        
  300 CONTINUE                                                                                                                      
C     CALL PRNTROW(TRAC,'TRAC','LIGNE 122',1,1,ILG1*JLAT,                                                                           
C    1             1,ILG1*JLAT,JLAT)                                                                                                
C                                                                                                                                   
      ELSE                                                                                                                          
C                                                                                                                                   
C     ------------------------------------------------------------------                                                            
C     * 2. OR FOR A DOUBLE COSINE**2 BLOB.                        ******                                                            
C     ------------------------------------------------------------------                                                            
C                                                                                                                                   
      DO 310 I=1,ILG1                                                                                                               
      DO 310 J=1,JLAT                                                                                                               
      TRAC((J-1)*ILG1+I)=0.E0                                                                                                       
  310 CONTINUE                                                                                                                      
C                                                                                                                                   
      RX=2.E0*PI*WID/360.E0                                                                                                         
      RY=2.E0*PI*BKG/360.E0                                                                                                         
      DLON=(2.E0*PI)/FLOAT(ILG)                                                                                                     
C                                                                                                                                   
C     * FIND THE INDICES BOUNDING REGION OF BLOB.                                                                                   
C                                                                                                                                   
      IW=MAX(1,  INT((X0-RX)/DLON)+2)                                                                                               
      IE=MIN(ILG,INT((X0+RX)/DLON)+1)                                                                                               
C                                                                                                                                   
      JS=1                                                                                                                          
      DO 320 J=2,JLAT                                                                                                               
      IF(RAD(J).GT.Y0-RY)THEN                                                                                                       
         JS=J                                                                                                                       
         GO TO 322                                                                                                                  
      ENDIF                                                                                                                         
  320 CONTINUE                                                                                                                      
  322 CONTINUE                                                                                                                      
C                                                                                                                                   
      JN=JLAT                                                                                                                       
      JLATM=JLAT-1                                                                                                                  
      DO 330 J=JLATM,1,-1                                                                                                           
      IF(RAD(J).LT.Y0+RY)THEN                                                                                                       
         JN=J                                                                                                                       
         GO TO 332                                                                                                                  
      ENDIF                                                                                                                         
  330 CONTINUE                                                                                                                      
  332 CONTINUE                                                                                                                      
C                                                                                                                                   
      JS=MAX(JS,1)                                                                                                                  
      JN=MIN(JN,JLAT)                                                                                                               
C                                                                                                                                   
C     * INSERT THE BLOB.                                                                                                            
C                                                                                                                                   
      DO 340 I=IW,IE                                                                                                                
      DO 340 J=JS,JN                                                                                                                
      TRAC((J-1)*ILG1+I) = AMP*COS(PI*(RLON(I)-X0)/(2.E0*RX))**2                                                                    
     1                        *COS(PI*( RAD(J)-Y0)/(2.E0*RY))**2                                                                    
  340 CONTINUE                                                                                                                      
C                                                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * ENSURE PERIODICITY.                                                                                                         
C                                                                                                                                   
      DO 350 J=1,JLAT                                                                                                               
      TRAC(J*ILG1)=TRAC((J-1)*ILG1+1)                                                                                               
  350 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * DISTRIBUTION INITIALE DES DIFFERENTS TRACEURS A TOUS LES NIVEAUX.                                                           
C                                                                                                                                   
      DO 600 N=1,NTRAC                                                                                                              
         CALL NOM(NAM,'X',N)                                                                                                        
                                                                                                                                    
         LEN=ILG1*JLAT                                                                                                              
         IF(N.EQ.2) THEN                                                                                                            
C                                                                                                                                   
C           LE PROFIL VERTICAL DU 2IEME TRACEUR EST DE FORME QUADRATIQUE                                                            
C                                                                                                                                   
            DO 400 L=1,ILEV                                                                                                         
               DO 375 I=1,LEN                                                                                                       
                  WRK(I)=TRAC(I)*1.E0/(1.E0+((PR(L)-430.E0)/400.E0)**2)                                                             
  375          CONTINUE                                                                                                             
                                                                                                                                    
               IF(ITRVAR.EQ.NC4TO8("RLNQ"))THEN                                                                                     
                  DO 385 I=1,LEN                                                                                                    
                     WRK(I)=-1.E0/LOG(WRK(I))                                                                                       
  385             CONTINUE                                                                                                          
                                                                                                                                    
               ELSEIF(ITRVAR.EQ.NC4TO8(" LNQ"))THEN                                                                                 
                  DO 395 I=1,LEN                                                                                                    
                     WRK(I)=LOG(WRK(I))                                                                                             
  395             CONTINUE                                                                                                          
                                                                                                                                    
               ENDIF                                                                                                                
C                                                                                                                                   
C              * CONVERSION DANS L'ESPACE SPECTRAL ET SAUVEGARDE DES CHAMPS                                                         
C              * A CHAQUE NIVEAU.                                                                                                   
C                                                                                                                                   
               CALL GGAST2(SSTRAC,LSR,LM,LA,WRK,ILG1,JLAT,1,SL,WL,ALP,                                                              
     1                     EPSI,WRKS,WRKL)                                                                                          
               CALL SETLAB(IBUF,NC4TO8("SPEC"),0,NAM,LEV(L),                                                                        
     +                                         LA,1,LRLMT,1)                                                                        
               CALL PUTFLD2(2,SSTRAC,IBUF,MAXX)                                                                                     
  400       CONTINUE                                                                                                                
                                                                                                                                    
         ELSE                                                                                                                       
C                                                                                                                                   
C           * LE PROFIL DES AUTRES TRACEURS SONT IDENTIQUES ET UNIFORMES                                                            
C           * DANS LA VERTICAL.                                                                                                     
C                                                                                                                                   
            IF(ITRVAR.EQ.NC4TO8("   Q"))THEN                                                                                        
               DO 404 I=1,LEN                                                                                                       
                  WRK(I)=TRAC(I)                                                                                                    
  404          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(ITRVAR.EQ.NC4TO8("RLNQ"))THEN                                                                                    
               DO 405 I=1,LEN                                                                                                       
                  WRK(I)=-1.E0/LOG(TRAC(I))                                                                                         
  405          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(ITRVAR.EQ.NC4TO8(" LNQ"))THEN                                                                                    
               DO 406 I=1,LEN                                                                                                       
                  WRK(I)=LOG(TRAC(I))                                                                                               
  406          CONTINUE                                                                                                             
                                                                                                                                    
            ELSE                                                                                                                    
               CALL                                XIT('INTRAC2',-89)                                                               
            ENDIF                                                                                                                   
C-----------------------------------------------------------------------                                                            
C     * CONVERT TO SPECTRAL.                                                                                                        
C                                                                                                                                   
            CALL GGAST2(SSTRAC,LSR,LM,LA,WRK,ILG1,JLAT,1,SL,WL,ALP,                                                                 
     1                     EPSI,WRKS,WRKL)                                                                                          
C-----------------------------------------------------------------------                                                            
C           * WRITE-OUT WHOLE TRACER FIELD.                                                                                         
C                                                                                                                                   
            DO 500 L=1,ILEV                                                                                                         
               CALL SETLAB(IBUF,NC4TO8("SPEC"),0,NAM,LEV(L),                                                                        
     +                                         LA,1,LRLMT,1)                                                                        
               CALL PUTFLD2(2,SSTRAC,IBUF,MAXX)                                                                                     
               WRITE(6,6010)IBUF                                                                                                    
  500       CONTINUE                                                                                                                
         ENDIF                                                                                                                      
         IF(N.GE.99)  CALL                         XIT('INTRAC2',-90)                                                               
  600 CONTINUE                                                                                                                      
C                                                                                                                                   
      CALL                                         XIT('INTRAC2',0)                                                                 
C                                                                                                                                   
C     * ERROR EXIT SECTION.                                                                                                         
C                                                                                                                                   
  901 CALL                                         XIT('INTRAC2',-4)                                                                
  902 CALL                                         XIT('INTRAC2',-5)                                                                
  903 CALL                                         XIT('INTRAC2',-6)                                                                
  904 CALL                                         XIT('INTRAC2',-7)                                                                
  905 CALL                                         XIT('INTRAC2',-8)                                                                
C-----------------------------------------------------------------------                                                            
 5000 FORMAT(10X,I5,1X,A4,5E10.0)                                               I4                                                  
 6000 FORMAT(' ITRAC,ITRVAR,AMP,X0,Y0,WID,BKG=',I5,1X,A4,5E13.3)                                                                    
 6001 FORMAT(' LRLMT=',I10)                                                                                                         
 6005 FORMAT(I5,' LEVELS ARE...',/,(19X,10I6))                                                                                      
 6006 FORMAT(' ILG,JLAT=',2I5)                                                                                                      
 6010 FORMAT(1X,60X,A4,I10,2X,A4,I10,4I8)                                                                                           
      END                                                                                                                           
