      PROGRAM GLOBAVL                                                                                                               
C     PROGRAM GLOBAVL (LL,       AVG,       INPUT,    OUTPUT,           )       C2                                                  
C    1           TAPE1=LL, TAPE2=AVG, TAPE5=INPUT, TAPE6=OUTPUT)                                                                    
C     ----------------------------------------------------------                C2                                                  
C                                                                               C2                                                  
C     FEB 16/10 - F.MAJAESS (ENSURE "ILAT" IS WITHIN SINGLE ARRAY DIMENSION)    C2                                                  
C     SEP 18/06 - F.MAJAESS (IMPROVE SPECIAL VALUE HANDLING)                                                                        
C     NOV 02/04 - S.KHARIN  (PRINT AVERAGES EXCLUDING SPECIAL VALUES.                                                               
C                            PRINT DIMENSIONS OF THE FIELDS.)                                                                       
C     MAR 05/02 - F.MAJAESS (REVISED FOR "CHAR" KIND RECORDS)                                                                       
C     DEC 19/00 - S. KHARIN (USE EXACT EXPRESSION FOR GRID BOX AREAS)                                                               
C                                                                               C2                                                  
CGLOBAVL - CALCULATES SURFACE MEAN OF FIELD ON A LAT/LON GRID           1  1 C  C1                                                  
C                                                                               C3                                                  
CAUTHOR  - G. FLATO AND S. KHARIN (BASED ON GLOBAVG).                           C3                                                  
C                                                                               C3                                                  
CPURPOSE - CALCULATES THE 2-D GLOBAL MEAN OF VALUES IN FILE LL, WHICH           C3                                                  
C          MAY BE A ZONL CROSS-SECTION OR A LAT/LON GRID (MAY BE SHIFTED).      C3                                                  
C          (I.E. LIKE GLOBAVG, BUT FOR NON-GAUSSIAN GRID FIELDS).               C3                                                  
C                                                                               C3                                                  
C          IF OUTPUT FILE, AVG, IS SPECIFIED, GLOBAL MEAN WRITTEN TO THIS FILE. C3                                                  
C          IF NOT, GLOBAL AND HEMISPHERIC AVERAGES ARE WRITTEN TO STANDARD      C3                                                  
C          OUTPUT.                                                              C3                                                  
C          NOTE - NAMES ARE NOT CHECKED.                                        C3                                                  
C                 "LABL" AND/OR "CHAR" KIND RECORDS ARE SIMPLY SKIPPED.         C3                                                  
C                                                                               C3                                                  
CINPUT FILE...                                                                  C3                                                  
C                                                                               C3                                                  
C      LL    = INPUT FILE OF FIELDS.                                            C3                                                  
C              THE FILE CAN CONTAIN "LABL" AND/OR "CHAR" RECORDS.               C3                                                  
C              THE FIELDS THEMSELVES MUST BE TYPE ZONL OR GRID.                 C3                                                  
C                                                                               C3                                                  
COUTPUT FILE...                                                                 C3                                                  
C                                                                               C3                                                  
C       AVG  = SPECIFICATION OF THIS FILE ON THE PROGRAM CALL IS OPTIONAL.      C3                                                  
C              IF AVG IS SPECIFIED THEN                                         C3                                                  
C                THE COMPUTED GLOBAL MEAN IS WRITTEN TO FILE AVG.               C3                                                  
C                THE TYPE OF FILE AVG IS THE SAME AS THAT FOR THE INPUT FILE    C3                                                  
C                LL. (UNCHANGED IBUF STRUCTURE    FOR EACH RECORD).             C3                                                  
C              OTHERWISE, THE COMPUTED VALUES ARE ONLY WRITTEN TO FILE OUTPUT.  C3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               C5                                                  
C       SHFTLT= THE VALUE IN DEGREES BY WHICH THE LAT/LON GRID IS SHIFTED       C5                                                  
C               FROM THE POLES.                                                 C5                                                  
C       IRPT  = 1 IF GRID HAS REPEATED LONGITUDE,                               C5                                                  
C               0 IF NOT.                                                       C5                                                  
C                                                                               C5                                                  
CEXAMPLE OF INPUT CARD...                                                       C5                                                  
C                                                                               C5                                                  
C  FOR A GRID SHIFTED 2.5 DEGREES IN LATITUDE, WITH A REPEATED LONGITUDE.       C5                                                  
C                                                                               C5                                                  
C* GLOBAVL       2.5         1                                                  C5                                                  
C                                                                               C5                                                  
C-----------------------------------------------------------------------                                                            
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/JCLPNTR/STDOUT                                                                                                         
      COMMON/F/F(65341)                                                                                                             
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
      REAL*8 DLAT,BLAT,SINL0,SINL1,SUMW,D,FN,DSH,DNH,DEQ,DGL,SUMWN                                                                  
      REAL*8 WL(181)                                                                                                                
C                                                                                                                                   
      LOGICAL OK,LISTCV,STDOUT                                                                                                      
      DATA MAXX/65341/,SPVAL/1.E+38/,EPS/1.E-6/                                                                                     
      DATA MAXLAT /181/                                                                                                             
      REAL*8 PI                                                                                                                     
      DATA PI/3.14159265358979323846D0/                                                                                             
C-------------------------------------------------------------------                                                                
                                                                                                                                    
C     * SETUP IN "SPVALT" THE TOLERANCE TO CHECK "SPVAL" AGAINST.                                                                   
                                                                                                                                    
      SPVALT=EPS*SPVAL                                                                                                              
                                                                                                                                    
      NF=4                                                                                                                          
      CALL JCLPNT(NF,1,2,5,6)                                                                                                       
      REWIND 1                                                                                                                      
      IF(NF.GT.3)THEN                                                                                                               
        LISTCV=.FALSE.                                                                                                              
        REWIND 2                                                                                                                    
      ELSE                                                                                                                          
        LISTCV=.TRUE.                                                                                                               
      ENDIF                                                                                                                         
C                                                                                                                                   
C     *  READ INPUT CARD AND ECHO INPUT VARIABLES TO STD. OUTPUT                                                                    
C                                                                                                                                   
      READ(5,5000,END=900) SHFTLT,IRPT                                          C4                                                  
      WRITE(6,6010) SHFTLT,IRPT                                                                                                     
C                                                                                                                                   
C     *  READ IN THE FIRST FIELD IN LL                                                                                              
C                                                                                                                                   
      ILAT0=0                                                                                                                       
      NFL = 0                                                                                                                       
  100 CONTINUE                                                                                                                      
      CALL GETFLD2 (1,F,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                   
      IF (.NOT.OK)THEN                                                                                                              
        IF (STDOUT) WRITE(6,6020) NFL                                                                                               
        IF (NFL.EQ.0)THEN                                                                                                           
                                                                                                                                    
C         * ERROR MESSAGE IF NO FIELD WAS READ                                                                                      
          CALL                                     XIT('GLOBAVL',-1)                                                                
                                                                                                                                    
        ELSE                                                                                                                        
                                                                                                                                    
C         * STANDARD MESSAGE WHEN LAST FIELD REACHED                                                                                
          IF (STDOUT) THEN                                                                                                          
            CALL                                   XIT('GLOBAVL',0)                                                                 
          ELSE                                                                                                                      
            STOP                                                                                                                    
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF (IBUF(1).EQ.NC4TO8("ZONL"))     GO TO 120                                                                                  
      IF (IBUF(1).EQ.NC4TO8("GRID"))     GO TO 160                                                                                  
      IF (IBUF(1).EQ.NC4TO8("LABL").OR.                                                                                             
     +    IBUF(1).EQ.NC4TO8("CHAR")    ) GO TO 100                                                                                  
                                                                                                                                    
C     * ERROR MESSAGE IF FIELD(S) OF WRONG TYPE                                                                                     
      CALL                                         XIT('GLOBAVL',-2)                                                                
C                                                                                                                                   
C     *  ZONAL CASE                                                                                                                 
C                                                                                                                                   
  120 ILG = 2                                                                                                                       
      ILG1 = ILG-1                                                                                                                  
      ILAT = IBUF(5)                                                                                                                
      DO 130 J=1,ILAT                                                                                                               
        JJ = ILAT + 1 - J                                                                                                           
        J2 = JJ * 2                                                                                                                 
        J1 = J2 - 1                                                                                                                 
        F(J2) = F(JJ)                                                                                                               
        F(J1) = F(J2)                                                                                                               
  130 CONTINUE                                                                                                                      
      GO TO 170                                                                                                                     
C                                                                                                                                   
C     *  GRID CASE                                                                                                                  
C                                                                                                                                   
  160 ILG = IBUF(5)                                                                                                                 
      ILAT = IBUF(6)                                                                                                                
C                                                                                                                                   
C     * SET NUMBER OF LONGITUDES TO BE USED IN AVERAGE (ILG1) DEPENDING                                                             
C     * ON WHETHER GRID HAS REPEATED LONGITUDE OR NOT.                                                                              
C                                                                                                                                   
      ILG1 = ILG - 1                                                                                                                
      IF(IRPT.EQ.0) ILG1=ILG                                                                                                        
C                                                                                                                                   
  170 IF (ILAT.GT.MAXLAT)  CALL                    XIT('GLOBAVL',-3)                                                                
C                                                                                                                                   
C     *  COMPUTE WEIGHTS FOR AVERAGING (FRACTION OF ZONAL SECTOR AREA                                                               
C     *  ASSOCIATED WITH GRID CELL AT EACH LATITUDE)                                                                                
C     *  DO IT ONLY IF "ILAT" HAS CHANGED.                                                                                          
C                                                                                                                                   
      IF (ILAT.NE.ILAT0) THEN                                                                                                       
         DLAT=(180.E0-(2.E0*SHFTLT))/FLOAT(ILAT-1)                                                                                  
C                                                                                                                                   
C     * BLAT IS GRID BOX BOUNDARY                                                                                                   
C                                                                                                                                   
         BLAT=-90.E0+SHFTLT+DLAT/2.E0                                                                                               
         SINL0=-1.E0                                                                                                                
         SINL1=SIN((PI/180.D0)*BLAT)                                                                                                
         WL(1)=SINL1-SINL0                                                                                                          
         WL(ILAT)=WL(1)                                                                                                             
         SUMW=WL(1)+WL(ILAT)                                                                                                        
         DO J=2,ILAT-1                                                                                                              
            BLAT=BLAT+DLAT                                                                                                          
            SINL0=SINL1                                                                                                             
            SINL1=SIN((PI/180.D0)*BLAT)                                                                                             
            WL(J)=SINL1-SINL0                                                                                                       
            SUMW=SUMW+WL(J)                                                                                                         
         ENDDO                                                                                                                      
                                                                                                                                    
C        * NORMALIZE WEIGHTS                                                                                                        
         SUMWN=1.E0/(SUMW*FLOAT(ILG1))                                                                                              
         DO J=1,ILAT                                                                                                                
            WL(J)=WL(J)*SUMWN                                                                                                       
         ENDDO                                                                                                                      
         ILAT0=ILAT                                                                                                                 
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * NOW COMPUTE AVERAGES                                                                                                        
C                                                                                                                                   
C     * TEST IF ILAT IS ODD (IODD=1) OR EVEN (IODD=0)                                                                               
      IODD=MOD(ILAT,2)                                                                                                              
C                                                                                                                                   
C     * IF EVEN (POINTS ON EITHER SIDE OF EQUATOR)                                                                                  
      IF(IODD.EQ.0) THEN                                                                                                            
C                                                                                                                                   
C       * SOUTHERN HEMISPHERE FIRST                                                                                                 
        D=0.E0                                                                                                                      
        W=0.E0                                                                                                                      
        NSPSH=0                                                                                                                     
        DO J=1,ILAT/2                                                                                                               
          FN = 0.E0                                                                                                                 
          WN = 0.E0                                                                                                                 
          DO I=1,ILG1                                                                                                               
            IJ = (J - 1) * ILG + I                                                                                                  
            IF(ABS(F(IJ)-SPVAL).GT.SPVALT) THEN                                                                                     
              FN = FN + F(IJ) * WL(J)                                                                                               
              WN = WN + WL(J)                                                                                                       
            ELSE                                                                                                                    
              NSPSH = NSPSH + 1                                                                                                     
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
          D = D + FN                                                                                                                
          W = W + WN                                                                                                                
        ENDDO                                                                                                                       
        DSH=2.E0*D                                                                                                                  
        WSH=2.E0*W                                                                                                                  
C                                                                                                                                   
C       * NORTHERN HEMISPHERE NEXT                                                                                                  
        D=0.E0                                                                                                                      
        W=0.E0                                                                                                                      
        NSPNH=0                                                                                                                     
        DO J=ILAT/2+1,ILAT                                                                                                          
          FN = 0.E0                                                                                                                 
          WN = 0.E0                                                                                                                 
          DO I=1,ILG1                                                                                                               
            IJ = (J - 1) * ILG + I                                                                                                  
            IF(ABS(F(IJ)-SPVAL).GT.SPVALT) THEN                                                                                     
              FN = FN + F(IJ) * WL(J)                                                                                               
              WN = WN + WL(J)                                                                                                       
            ELSE                                                                                                                    
              NSPNH = NSPNH + 1                                                                                                     
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
          D = D + FN                                                                                                                
          W = W + WN                                                                                                                
        ENDDO                                                                                                                       
        DNH=2.E0*D                                                                                                                  
        WNH=2.E0*W                                                                                                                  
C                                                                                                                                   
        NSPEQ=0                                                                                                                     
      ELSE                                                                                                                          
C                                                                                                                                   
C       * IF ODD (POINT AT EQUATOR)                                                                                                 
C                                                                                                                                   
C       * SOUTHERN HEMISPHERE FIRST                                                                                                 
        D=0.E0                                                                                                                      
        W=0.E0                                                                                                                      
        NSPSH=0                                                                                                                     
        DO J=1,(ILAT-1)/2                                                                                                           
          FN = 0.E0                                                                                                                 
          WN = 0.E0                                                                                                                 
          DO I=1,ILG1                                                                                                               
            IJ = (J - 1) * ILG + I                                                                                                  
            IF(ABS(F(IJ)-SPVAL).GT.SPVALT) THEN                                                                                     
              FN = FN + F(IJ) * WL(J)                                                                                               
              WN = WN + WL(J)                                                                                                       
            ELSE                                                                                                                    
              NSPSH = NSPSH + 1                                                                                                     
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
          D = D + FN                                                                                                                
          W = W + WN                                                                                                                
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * EQUATORIAL BAND NEXT                                                                                                      
        J = (ILAT-1)/2 + 1                                                                                                          
        FN = 0.E0                                                                                                                   
        WN = 0.E0                                                                                                                   
        NSPEQ=0                                                                                                                     
        DO I=1,ILG1                                                                                                                 
          IJ = (J - 1) * ILG + I                                                                                                    
          IF(ABS(F(IJ)-SPVAL).GT.SPVALT) THEN                                                                                       
            FN = FN + F(IJ) * WL(J)                                                                                                 
            WN = WN + WL(J)                                                                                                         
          ELSE                                                                                                                      
            NSPEQ = NSPEQ + 1                                                                                                       
          ENDIF                                                                                                                     
        ENDDO                                                                                                                       
        DEQ = FN                                                                                                                    
        WEQ = WN                                                                                                                    
        DSH=2.E0*D+DEQ                                                                                                              
        WSH=2.E0*W+WEQ                                                                                                              
C                                                                                                                                   
C       * NORTHERN HEMISPHERE LAST                                                                                                  
        D=0.E0                                                                                                                      
        W=0.E0                                                                                                                      
        NSPNH=0                                                                                                                     
        DO J=(ILAT+1)/2+1,ILAT                                                                                                      
          FN = 0.E0                                                                                                                 
          WN = 0.E0                                                                                                                 
          DO I=1,ILG1                                                                                                               
            IJ = (J - 1) * ILG + I                                                                                                  
            IF(ABS(F(IJ)-SPVAL).GT.SPVALT) THEN                                                                                     
              FN = FN + F(IJ) * WL(J)                                                                                               
              WN = WN + WL(J)                                                                                                       
            ELSE                                                                                                                    
              NSPNH = NSPNH + 1                                                                                                     
            ENDIF                                                                                                                   
          ENDDO                                                                                                                     
          D = D + FN                                                                                                                
          W = W + WN                                                                                                                
        ENDDO                                                                                                                       
        DNH=2.E0*D+DEQ                                                                                                              
        WNH=2.E0*W+WEQ                                                                                                              
      ENDIF                                                                                                                         
                                                                                                                                    
C     * AVERAGE TWO HEMISPHERES TO GET GLOBAL AVERAGE                                                                               
      IF(WSH+WNH.NE.0.E0) THEN                                                                                                      
        DGL=(DSH + DNH)/(WSH+WNH)                                                                                                   
      ELSE                                                                                                                          
        DGL=SPVAL                                                                                                                   
      ENDIF                                                                                                                         
      IF(WSH.NE.0.E0)THEN                                                                                                           
        DSH=DSH/WSH                                                                                                                 
      ELSE                                                                                                                          
        DSH=SPVAL                                                                                                                   
      ENDIF                                                                                                                         
      IF(WNH.NE.0.E0)THEN                                                                                                           
        DNH=DNH/WNH                                                                                                                 
      ELSE                                                                                                                          
        DNH=SPVAL                                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * WRITE GLOBAL AVERAGE AND HEMISPHERIC AVERAGES TO STANDARD OUTPUT                                                            
C                                                                                                                                   
      NSP=NSPNH+NSPSH+NSPEQ                                                                                                         
      IF(NSP.EQ.0) THEN                                                                                                             
        WRITE(6,6030) (IBUF(N),N=1,6),DGL                                                                                           
        WRITE(6,6040) DNH                                                                                                           
        WRITE(6,6050) DSH                                                                                                           
      ELSE                                                                                                                          
        WRITE(6,6031) (IBUF(N),N=1,6),DGL,NSP                                                                                       
        WRITE(6,6041) DNH,NSPNH+NSPEQ/2                                                                                             
        WRITE(6,6051) DSH,NSPSH+(NSPEQ-NSPEQ/2)                                                                                     
      ENDIF                                                                                                                         
                                                                                                                                    
      IF(.NOT.LISTCV) THEN                                                                                                          
C                                                                                                                                   
C       * PUT THE GLOBAL AVERAGE VALUE INTO EVERY ELEMENT OF ARRAY F                                                                
C                                                                                                                                   
        IF (IBUF(1).EQ.NC4TO8("ZONL")) ILG = 1                                                                                      
        DO J=1,ILAT                                                                                                                 
          DO I=1,ILG                                                                                                                
            IJ = (J - 1) * ILG + I                                                                                                  
            F(IJ) = DGL                                                                                                             
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * WRITE F TO OUTPUT FILE, AVG=TAPE2.                                                                                        
C                                                                                                                                   
        CALL PUTFLD2(2,F,IBUF,MAXX)                                                                                                 
                                                                                                                                    
      ENDIF                                                                                                                         
                                                                                                                                    
      NFL = NFL + 1                                                                                                                 
                                                                                                                                    
      GO TO 100                                                                                                                     
C                                                                                                                                   
C     * WRITE ERROR MESSAGE IF END OF FILE ON INPUT CARD                                                                            
C                                                                                                                                   
  900 CALL                                         XIT('GLOBAVL',-4)                                                                
C----------------------------------------------------------------                                                                   
 5000 FORMAT(10X,E10.0,5X,I5)                                                   C4                                                  
 6010 FORMAT('  GLOBAVL GRID PARAMETERS: SHFTLT=',E12.6,', IRPT=',I5)                                                               
 6020 FORMAT('  GLOBAVL READ',I5,' FIELDS')                                                                                         
 6030 FORMAT(' FROM ',A4,1X,I10,1X,A4,3(1X,I5),                                                                                     
     1       ' SURFACE MEAN IS ',E12.6)                                                                                             
 6040 FORMAT('        THE NH MEAN IS ',E12.6)                                                                                       
 6050 FORMAT('        THE SH MEAN IS ',E12.6)                                                                                       
 6031 FORMAT(' FROM ',A4,1X,I10,1X,A4,3(1X,I5),                                                                                     
     1       ' SURFACE MEAN IS ',E12.6,' NSPVAL=',I10)                                                                              
 6041 FORMAT('        THE NH MEAN IS ',E12.6,' NSPVAL=',I10)                                                                        
 6051 FORMAT('        THE SH MEAN IS ',E12.6,' NSPVAL=',I10)                                                                        
      END                                                                                                                           
