      PROGRAM GSOMEGA                                                                                                               
C     PROGRAM GSOMEGA (SSLNSP,       SSVORT,       SSDIV,       GSOMEG,         J2                                                  
C    1                                             INPUT,       OUTPUT, )       J2                                                  
C    2           TAPE1=SSLNSP, TAPE2=SSVORT, TAPE3=SSDIV, TAPE4=GSOMEG,                                                             
C    3                                       TAPE5=INPUT, TAPE6=OUTPUT)                                                             
C     -----------------------------------------------------------------         J2                                                  
C                                                                               J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       J2                                                  
C     FEB 15/94 - F.MAJAESS (REVISE FOR TRUNCATIONS > 99)                                                                           
C     JAN 13/93 - E. CHAN   (DECODE LEVELS IN 8-WORD LABEL)                                                                         
C     JUL 13/92 - E. CHAN   (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII)                                                                   
C     JAN 24/91 - F.MAJAESS (ADJUST FOR PASSING WORKING ARRAYS TO OMEGA2)                                                           
C     NOV 22/88 - M.LAZARE. (MODIFY DIMENDION TO HANDLE T20-L20 OR T32-L15)                                                         
C     DEC 24/86 - M.LAZARE. (ALLOW L20 OPTION WITH T20)                                                                             
C     FEB 19/85 - B.DUGAS.  (RESET DIMENSION TO MAXIMUM T30-L15 EFFECTIVE)                                                          
C     JAN 26/84 - B.DUGAS, R.LAPRISE, J.D.HENDERSON.                                                                                
C                                                                               J2                                                  
CGSOMEGA - CONVERTS Q,D FILES TO VERTICAL MOTION                        3  1 C GJ1                                                  
C                                                                               J3                                                  
CAUTHOR  - J.D.HENDERSON                                                        J3                                                  
C                                                                               J3                                                  
CPURPOSE - COMPUTES A GRID FILE OF VERTICAL MOTION (DP/DT) FROM                 J3                                                  
C          SPECTRAL FILES OF VORTICITY AND DIVERGENCE.                          J3                                                  
C          NOTE - INPUT IS SPECTRAL, OUTPUT IS GAUSSIAN GRIDS,                  J3                                                  
C                 ALL GLOBAL.                                                   J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      SSLNSP = SPECTRAL LN(SF PRES)                                            J3                                                  
C      SSVORT = SPECTRAL VORTICITY  ON SIGMA LEVELS                             J3                                                  
C      SSDIV  = SPECTRAL DIVERGENCE ON SIGMA LEVELS                             J3                                                  
C                                                                               J3                                                  
COUTPUT FILE...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSOMEG = GRIDS OF VERTICAL MOTION OMEGA (DP/DT) ON SIGMA                 J3                                                  
C               HALF LEVELS. UNITS ARE NEWTONS/M**2/SEC.                        J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      ILG   = LENGTH OF A GAUSSIAN GRID ROW                                    J5                                                  
C      ILAT  = NUMBER OF GAUSSIAN LATITUDES.                                    J5                                                  
C      KPACK = OUTPUT PACKING DENSITY (0 DEFAULT TO 4)                          J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARD...                                                       J5                                                  
C                                                                               J5                                                  
C* GSOMEGA   64   52    0                                                       J5                                                  
C---------------------------------------------------------------------------                                                        
C                                                                                                                                   
C     * COMPLEX MULTILEVEL ARRAY FOR SURFACE PRESSURE,VORTICITY AND                                                                 
C     * DIVERGENCE SPECTRAL COEFFICIENTS, IN THAT ORDER IN SP.                                                                      
C     * SP((1+2*ILEV)*LA),T(LA),ES(LA)                                                                                              
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMPLEX SP,T,ES                                                                                                               
      COMMON /SPEC/ SP(431145),T(2145),ES(2145)                                                                                     
                                                                                                                                    
C     * WORKS ARRAYS.                                                                                                               
C     * ALP(LA+LM),DALP(LA+LM),DELALP(LA+LM),EPSI(LA+LM),                                                                           
C     * GGOM((ILG+1)*ILAT*ILEV)                                                                                                     
                                                                                                                                    
      REAL*8 ALP(2210),DALP(2210),DELALP(2210),EPSI(2210)                                                                           
      COMMON/BLANCK/GGOM(1852800)                                                                                                   
                                                                                                                                    
C     * MULTILEVEL WORK ARRAY FOR SLICE VALUES OF PRESSG, P, C,                                                                     
C     * PSDPG, UG, VG AND PSDLG, IN THAT ORDER IN SLICE.                                                                            
C     * SLICE(MAXLG,(4*ILEV+3))                                                                                                     
C     * TG(MAXLG),ESG(MAXLG),OMEGAG(MAXLG*ILEV)                                                                                     
                                                                                                                                    
      COMMON /SLICE/ SLICE(194,403)                                                                                                 
      COMMON /SLICE/ TG(194),ESG(194),OMEGAG(19400)                                                                                 
                                                                                                                                    
C     * WORK FIELDS.                                                                                                                
C     * TRIGS(MAXLG),LSR(2,LM+1),LS(ILEV),...                                                                                       
                                                                                                                                    
      LOGICAL OK,LVECT                                                                                                              
      REAL TRIGS(194)                                                                                                               
      INTEGER LSR(2,66),LS(100),LH(100),LEVP(100),                                                                                  
     1        LEVC(100),IFAX(10)                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
                                                                                                                                    
C     * WRKS((ILG+1)*64)                                                                                                            
                                                                                                                                    
      REAL*8 SL(96),CL(96),WL(96),WOSSL(96),RADL(96)                                                                                
      COMMON/SCM1/ WRKS(12416)                                                                                                      
      COMMON/PARAMS/ WW,TW,A,ASQ,GRAV,RGAS,RGOCP,RGOASQ,CPRES                                                                       
                                                                                                                                    
C     * WORKING ARRAYS NEEDED IN CALLING OMEGA2.                                                                                    
                                                                                                                                    
      REAL CGH(100),UGH(100),VGH(100),CIGH(100),UIGH(100),VIGH(100),                                                                
     1     SDH(100),WA (100),WB (100)                                                                                               
                                                                                                                                    
C     * THE FOLLOWING ARRAYS ARE FOR FUNCTIONS OF ILEV.                                                                             
C     * STAMS(ILEV*ILEV),STAIS(ILEV*ILEV),...                                                                                       
                                                                                                                                    
      REAL STAMS(10000),STAIS(10000)                                                                                                
      REAL S(100),SH(100),SF(100),DS(100),DEL(100)                                                                                  
      REAL SEXPK(100),SHEXPK(100)                                                                                                   
      COMPLEX P2(100)                                                                                                               
                                                                                                                                    
C     * BETA AND AVERT SHOULD BE THE SAME AS IN THE GCM.                                                                            
                                                                                                                                    
      DATA LVECT/.TRUE./                                                                                                            
      DATA BETA,AVERT/ 1.E0, 0.25E0 /                                                                                               
C     * MAXLG/ILG+1 OR ILG+2/, MAXX/ILG+1*ILAT/, MAXSP/LA/, MAXMSP/LA*ILEV/                                                         
      DATA MAXLG/194/, MAXX/18528/, MAXSP/2145/, MAXMSP/214500/                                                                     
      DATA MAXL/100/                                                                                                                
C---------------------------------------------------------------------                                                              
      NFF=6                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,4,5,6)                                                                                                  
      DO 110 N=1,4                                                                                                                  
  110 REWIND N                                                                                                                      
                                                                                                                                    
C     * READ SIZE OF OUTPUT GRID AND PACKING DENSITY (0 DEFAULTS TO 4).                                                             
                                                                                                                                    
      READ(5,5010,END=905) ILG,ILAT,KPACK                                       J4                                                  
      IF(ILG.GT.MAXLG) ILG=MAXLG                                                                                                    
      IF(KPACK.EQ.0) KPACK=4                                                                                                        
      ILG1=ILG+1                                                                                                                    
      ILH=MAXLG/2                                                                                                                   
      ILATH=ILAT/2                                                                                                                  
      WRITE(6,6010) ILG1,ILAT                                                                                                       
                                                                                                                                    
C     * GET SIGMA LEVELS FROM VORT FILE.                                                                                            
                                                                                                                                    
      CALL FILEV(LS,ILEV,IBUF,2)                                                                                                    
      IF((ILEV.LT.1).OR.(ILEV.GT.MAXL)) CALL       XIT('GSOMEGA',-1)                                                                
      CALL LVDCODE(S,LS,ILEV)                                                                                                       
      DO 120 L=1,ILEV                                                                                                               
  120 S(L)=S(L)*0.001E0                                                                                                             
      IP=2                                                                                                                          
      IC=IP+ILEV                                                                                                                    
      IT=IC+ILEV                                                                                                                    
      IU=IT+1                                                                                                                       
      IV=IU+ILEV                                                                                                                    
      IL=IV+ILEV                                                                                                                    
      CALL ISIGL(LS,LH,S,ILEV)                                                                                                      
                                                                                                                                    
C     * CALCULATE CONSTANTS.                                                                                                        
                                                                                                                                    
      LRLMT=IBUF(7)                                                                                                                 
C     IF(LRLMT.LT.100) LRLMT=1000*IBUF(5)+10*IBUF(6)                                                                                
      IF(LRLMT.LT.100) CALL FXLRLMT (LRLMT,IBUF(5),IBUF(6),0)                                                                       
      CALL DIMGT(LSR,LA,LR,LM,KTR,LRLMT)                                                                                            
      CALL PRLRLMT (LA,LR,LM,KTR,LRLMT)                                                                                             
      IF(LA.GT.MAXSP) CALL                         XIT('GSOMEGA',-2)                                                                
      IF(LA*ILEV.GT.MAXMSP) CALL                   XIT('GSOMEGA',-3)                                                                
      ISP=LA+1                                                                                                                      
      ISC=(ILEV+1)*LA+1                                                                                                             
                                                                                                                                    
      CALL SPMCON(S,SEXPK,SH,SHEXPK,SF,DS,DEL,FVORT,PI,ILEV)                                                                        
      CALL EPSCAL(EPSI,LSR,LM)                                                                                                      
      CALL GAUSSG(ILATH,SL,WL,CL,RADL,WOSSL)                                                                                        
      CALL  TRIGL(ILATH,SL,WL,CL,RADL,WOSSL)                                                                                        
      CALL STMCAL(STAMS,STAIS,1.5E0,-0.5E0,S,SH,SF,DEL,RGOCP,ILEV)                                                                  
      CALL FTSETUP(TRIGS,IFAX,ILG)                                                                                                  
C---------------------------------------------------------------------                                                              
C     * GET LN(SF PRES) FOR THE NEXT STEP. CONVERT MB TO N/M**2.                                                                    
                                                                                                                                    
      NSETS=0                                                                                                                       
  150 CALL GETFLD2(1,SP,NC4TO8("SPEC"),-1,NC4TO8("LNSP"),1,IBUF,MAXX,OK)                                                            
      IF(.NOT.OK)THEN                                                                                                               
        IF(NSETS.EQ.0)THEN                                                                                                          
          CALL                                     XIT('GSOMEGA',-4)                                                                
        ELSE                                                                                                                        
          WRITE(6,6030) NSETS                                                                                                       
          CALL                                     XIT('GSOMEGA',0)                                                                 
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF(NSETS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      SP(1)=SP(1)+LOG(100.E0)*SQRT(2.E0)                                                                                            
                                                                                                                                    
C     * GET ILEV LEVELS OF VORTICITY AND DIVERGENCE.                                                                                
                                                                                                                                    
      NST= IBUF(2)                                                                                                                  
      CALL GETSET2 (2,SP(ISP),LEVP,ILEV,IBUF,MAXX,OK)                                                                               
      IF(.NOT.OK.OR.IBUF(2).NE.NST.OR.ILEV.GT.MAXL)                                                                                 
     1  CALL                                       XIT('GSOMEGA',-5)                                                                
      IF(NSETS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      CALL GETSET2 (3,SP(ISC),LEVC,ILEV,IBUF,MAXX,OK)                                                                               
      IF(.NOT.OK.OR.IBUF(2).NE.NST.OR.ILEV.GT.MAXL)                                                                                 
     1  CALL                                       XIT('GSOMEGA',-6)                                                                
      IF(NSETS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      DO 170 L=1,ILEV                                                                                                               
      IF(LEVP(L).NE.LS(L).OR.LEVC(L).NE.LS(L))                                                                                      
     1       CALL                                  XIT('GSOMEGA',-L-20)                                                             
  170 CONTINUE                                                                                                                      
      IF(NSETS.EQ.0) WRITE(6,6025) ILEV,(LH(L),L=1,ILEV)                                                                            
C---------------------------------------------------------------------                                                              
C     * CALCULATE PSI,CHI FROM VORT,DIV.                                                                                            
                                                                                                                                    
      CALL QDTFPC(SP(ISP),SP(ISC),SP(ISP),SP(ISC),                                                                                  
     1            LA,LSR,LM,ILEV,-1)                                                                                                
                                                                                                                                    
C     * LATITUDE LOOP.                                                                                                              
                                                                                                                                    
      DO 390 J=1,ILAT                                                                                                               
                                                                                                                                    
C     * COMPUTE ALP,DALP AND DELALP.                                                                                                
                                                                                                                                    
      CALL ALPST2(ALP,LSR,LM,SL(J),EPSI)                                                                                            
      CALL ALPDY2(DALP,ALP,LSR,LM,EPSI)                                                                                             
      CALL ALPDL2(DELALP,ALP,LSR,LM)                                                                                                
                                                                                                                                    
C     * COMPUTE GRID POINT VALUES FROM SPECTRAL COEFF.                                                                              
C     * (MOISTURE AND TEMPERATURE ARE TURNED OFF).                                                                                  
                                                                                                                                    
      CALL MHEXP (SLICE(1,IP),SLICE(1,IC),SLICE(1,IU),SLICE(1,IV),                                                                  
     1            TG,ESG, SLICE(1,IL),SLICE(1,IT),SLICE, ILH, ILG,                                                                  
     2            SP(ISP),SP(ISC),T,ES,SP,LA,LSR,LM,ILEV, 0,P2, 0,                                                                  
     3            FVORT,ALP,DALP,DELALP,WRKS,IFAX,TRIGS,LVECT)                                                                      
                                                                                                                                    
      DO 320 IK=1,ILG                                                                                                               
  320 SLICE(IK,1)=EXP(SLICE(IK,1))                                                                                                  
                                                                                                                                    
C     * COMPUTE THE VERTICAL MOTION (NEWTONS PER SQ METER PER SEC).                                                                 
                                                                                                                                    
      CALL OMEGA2 (OMEGAG,SLICE(1,IC),SLICE(1,IU),SLICE(1,IV),                                                                      
     1             SLICE(1,IL),SLICE(1,IT),SLICE,                                                                                   
     2             CGH,UGH,VGH,CIGH,UIGH,VIGH,SDH,WA,WB,                                                                            
     3             MAXLG,ILG,ILEV,S,SH,DS,AVERT,BETA,CL(J))                                                                         
                                                                                                                                    
C     * TRANSFER VERTICAL SLICE FOR LATITUDE IHGG IN OMEGAG                                                                         
C     * TO THE MULTILEVEL GAUSSIAN GRID ARRAY GGOM.                                                                                 
                                                                                                                                    
      CALL IGGSL (GGOM,ILG1,ILAT,ILEV, OMEGAG,MAXLG,ILG, J)                                                                         
  390 CONTINUE                                                                                                                      
C---------------------------------------------------------------------                                                              
C     * WRITE ILEV LEVELS OF OMEGA GAUSSIAN GRIDS TO FILE 4 .                                                                       
                                                                                                                                    
      CALL SETLAB(IBUF,NC4TO8("GRID"),NST,NC4TO8("OMEG"),0,                                                                         
     +                                   ILG1,ILAT,0,KPACK)                                                                         
      CALL PUTSET2(4,GGOM,LH,ILEV,IBUF,MAXX)                                                                                        
      IF(NSETS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      NSETS=NSETS+1                                                                                                                 
      GO TO 150                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  905 CALL                                         XIT('GSOMEGA',-7)                                                                
C--------------------------------------------------------------------                                                               
 5010 FORMAT(10X,3I5)                                                           J4                                                  
 6010 FORMAT('0GRID SIZE =',2I6)                                                                                                    
 6025 FORMAT(' OMEGA ON ',I2,' SIGMA LEVELS',/                                                                                      
     1       ' THEY ARE (*1000) : ',20I5/(20X,20I5))                                                                                
 6030 FORMAT('0',I5,'  SETS OF OMEGA COMPUTED')                                                                                     
      END                                                                                                                           
