      SUBROUTINE PRE_AREAVER(GAPS_LAMBDA_SRCE,LAMBDA_SRCE               &
     &,GAPS_PHI_SRCE,PHI_SRCE,CYCLIC_SRCE,LROW_SRCE,WANT,MASK_SRCE      &
     &,GAPS_LAMBDA_TARG,LAMBDA_TARG,GAPS_PHI_TARG,PHI_TARG              &
     &,CYCLIC_TARG,SPHERICAL                                            &
     &,MAXL,COUNT_TARG,BASE_TARG,INDEX_SRCE,WEIGHT,ICODE,CMESSAGE)      
                                                                       
!     Subroutine PRE_AREAVER                                            
!                                                                       
!     Calculate weights for area-averaging data on the source grid to   
!     data on the target grid.                                          
!                                                                       
!    J.Gregory  <- programmer of some or all of previous code or changes
!                                                                       
!    Model            Modification history from model version 3.0:      
!   version  Date                                                       
!     3.2    13/07/93 Changed CHARACTER*(*) to CHARACTER*(80) for       
!                     portability.  Author Tracey Smith.                
!     3.3   3.9.93    Correct out-of-bounds errors for IXL and IYT      
!     4.1   22.5.96   J.M.Gregory  Add argument SPHERICAL to produce    
!                     correct weights for a spherical surface.          
!     5.3   22.10.01  Change for south -> north ordering. C.F.Durman.   
!     5.5  28.02.03  Add defined RIVERS. C. Bunton                         
!     6.0  12.09.03  Change DEF from A20 to A26. D. Robinson.              
!                                                                       
!     Standard, paper 4 version 4 (14.12.90)                            
!                                                                       
!                                                                       
!     The source grid and target grid are each specified by a lambda set
!     and a phi set of coordinates delimiting the boxes. These sets are 
!     supplied in 1D arrays aa_bb for coordinate aa=LAMBDA,PHI on grid  
!     bb=SRCE,TARG.  The number of gaps is specified by GAPS_aa_bb,     
!     which is equal to the number of lines IF (CY! IC_bb) and aa is    
!     LAMBDA, otherwise one less. (By "gap" we mean the interval        
!     spanning a box in one coordinate only. The total number of boxes, 
!     or grid points, is the product of the numbers of gaps in the two  
!     coordinates.) Whether the axes are cyclic is not known until      
!     run-time, so the dimensions of the arrays LAMBDA_bb are not known 
!     at compile-time, and they are dimensioned assumed-size. There are 
!     no restrictions on the base meridian of lambda, and it does not   
!     have to be the same for source and target grids. The lambda       
!     coordinates should be in increasing order (running from left to   
!     right), the phi increasing (top to bottom). The coordinates must  
!     be given in degrees for cyclic axes, because a range of 360 is    
!     assumed, or IF (SPHERICAL), when trigonometric functions are      
!     used. IF (SPHERICAL), the weights computed are for a spherical    
!     surface, assuming that LAMBDA is longitude and PHI latitude.      
!     Otherwise, LAMBDA and PHI are treated as Cartesian axes on a      
!     plane.                                                            
!                                                                       
!     The array MASK_SRCE is the land/sea mask for the source grid. The 
!     logical value indicating points to be used should be supplied in  
!     WANT. The first dimension of MASK_SRCE should be supplied in      
!     LROW_SRCE. Specifying this separately allows for the possibility  
!     of extra rows and columns in MASK_SRCE which are to be ignored.   
!                                                                       
!     The arrays COUNT_TARG and BASE_TARG should be dimensioned in the  
!     calling program to the number of boxes in the target array.       
!                                                                       
!     The arrays INDEX_SRCE and WEIGHT are returned in the form which   
!     the area-averaging routine DO_AREAVER expects. They are continuous
!     lists comprising consecutive groups of entries. There is a group  
!     for each target point, for which the number of entries is spec-   
!     ified by COUNT_TARG, and the groups appear in the normal order of 
!     grid points. The size required for INDEX_SRCE and WEIGHT depends  
!     on how many source boxes go into each target box, on average, and 
!     is not known at compile-time. The maximum that could be needed is 
!     (GAPS_LAMBDA_SRCE+GAPS_LAMBDA_TARG)*(GAPS_PHI_SRCE+GAPS_PHI_TARG) 
!     and the size to which the arrays are actually dimensioned should  
!     be supplied in MAXL. The size used is returned in MAXL. It is the 
!     responsibility of the calling routine to provide enough space.    

!-----------------------------------------------------------------------
      USE c_pi, ONLY :  &
!      imported scalar parameters
        pi_over_180
!-----------------------------------------------------------------------
!                                                                       
      IMPLICIT NONE                                                     

      INTEGER ::  &
       GAPS_LAMBDA_SRCE       &!IN number of lambda gaps in source grid 
      ,GAPS_PHI_SRCE          &!IN number of phi gaps in source grid    
      ,GAPS_LAMBDA_TARG       &!IN number of lambda gaps in target grid 
     &,GAPS_PHI_TARG          &!IN number of phi gaps in target grid    
      ,LROW_SRCE              &!IN first dimension of MASK_SRCE         
      ,MAXL                   &!INOUT maximum entries in output lists   
     &,COUNT_TARG(GAPS_LAMBDA_TARG,GAPS_PHI_TARG) &
!                              !OUT no. of source boxes in target box   
     &,BASE_TARG(GAPS_LAMBDA_TARG,GAPS_PHI_TARG)  &
!                              !OUT first index in list for target box  
     &,INDEX_SRCE(MAXL)       &!OUT list of source box indices          
      ,ICODE                   !OUT return code                         

      LOGICAL ::  &
     & CYCLIC_SRCE            &!IN source grid is cyclic                
      ,CYCLIC_TARG            &!IN target grid is cyclic                
      ,WANT                   &!IN indicator of wanted points in mask   
      ,MASK_SRCE(LROW_SRCE,*) &!IN land/sea mask for source grid        
      ,SPHERICAL               !IN calculate weights for a sphere       

      REAL ::  &
       LAMBDA_SRCE(*)         &!IN source lambda line coordinates       
      ,PHI_SRCE(*)            &!IN source phi line coordinates          
      ,LAMBDA_TARG(*)         &!IN target lambda line coordinates       
      ,PHI_TARG(*)            &!IN target phi line coordinates          
      ,WEIGHT(MAXL)            !OUT list of weights for source boxes    

      CHARACTER*(80) :: CMESSAGE !OUT error message                       

      INTEGER ::  &
       LINES_LAMBDA_SRCE      &! number of source lambda lines          
      ,LINES_PHI_SRCE         &! number of source phi lines             
      ,LINES_LAMBDA_TARG      &! number of target lambda lines          
      ,LINES_PHI_TARG         &! number of target phi lines             
      ,COUNT_LAMBDA(GAPS_LAMBDA_TARG)  &                                   
!                              ! number of source lambda gaps per target
      ,COUNT_PHI(GAPS_PHI_TARG)  &                                         
!                              ! number of source phi gaps per target   
      ,INDEX_LAMBDA(GAPS_LAMBDA_SRCE+GAPS_LAMBDA_TARG)  &
!                              ! source lambda gap indices              
      ,INDEX_PHI(GAPS_PHI_SRCE+GAPS_PHI_TARG)  &
!                              ! source phi gap indices                 
      ,IX1,IY1,IX2,IY2        &! working SRCE/TARG LAMBDA/PHI indices   
      ,IX1N,IX1W              &! working indices                        
      ,IXL(GAPS_LAMBDA_TARG+1)&! source line on the left of target line 
      ,IX2N                   &! target gap to the right of IX2         
      ,IYT(GAPS_PHI_TARG+1)   &! source line above target line          
      ,IXP,IYP,IP             &! pointers into lists                    
      ,IX,IY,I                 ! loop indices                           

      REAL  ::  &
       BASLAM                 &! minimum lambda for TEMP coordinates    
      ,BTARG                  &! width of target gap                    
      ,BLO,BHI                &! limits of gap                          
      ,TEMP_SRCE(GAPS_LAMBDA_SRCE+1)  &                                  
!                              ! adjusted version of LAMBDA_SRCE        
      ,TEMP_TARG(GAPS_LAMBDA_TARG+1)  &                                    
!                              ! adjusted version of LAMBDA_TARG        
      ,FRAC_LAMBDA(GAPS_LAMBDA_SRCE+GAPS_LAMBDA_TARG)  &
!                              ! fractions of target lambda gaps        
      ,FRAC_PHI(GAPS_PHI_SRCE+GAPS_PHI_TARG)  &
!                              ! fractions of target phi gaps           
      ,SUM                     ! sum of weights                         
                                                                       
!     1  Set up the lambda coordinates to make them easier to handle.   
!                                                                       
!     1.1  Produce in TEMP_SRCE a monotonically increasing set of angles
!     equivalent to LAMBDA_SRCE i.e. equal under modulo 360.            
                                                                       
      IF (CYCLIC_SRCE) THEN                                             
        LINES_LAMBDA_SRCE=GAPS_LAMBDA_SRCE                              
      ELSE                                                              
        LINES_LAMBDA_SRCE=GAPS_LAMBDA_SRCE+1                            
      ENDIF                                                             
      BASLAM=LAMBDA_SRCE(1)                                             
      DO 10 IX1=1,LINES_LAMBDA_SRCE                                     
        IF (LAMBDA_SRCE(IX1).LT.BASLAM) THEN                            
          TEMP_SRCE(IX1)=LAMBDA_SRCE(IX1)+360.                          
        ELSE                                                            
          TEMP_SRCE(IX1)=LAMBDA_SRCE(IX1)                               
        ENDIF                                                           
   10 CONTINUE ! Labelled DO for the sake of fpp                        
                                                                       
!     1.2  Produce in TEMP_TARG a set of angles equivalent to           
!     LAMBDA_TARG i.e. equal under modulo 360, but all in the range     
!     BASLAM to BASLAM+360, where BASLAM=min(TEMP_LAMBDA).              
                                                                       
      IF (CYCLIC_TARG) THEN                                             
        LINES_LAMBDA_TARG=GAPS_LAMBDA_TARG                              
      ELSE                                                              
        LINES_LAMBDA_TARG=GAPS_LAMBDA_TARG+1                            
      ENDIF                                                             
      DO 20 IX2=1,LINES_LAMBDA_TARG                                     
        TEMP_TARG(IX2)=MOD(LAMBDA_TARG(IX2)-BASLAM,360.)                
        IF (TEMP_TARG(IX2).LT.0.) TEMP_TARG(IX2)=TEMP_TARG(IX2)+360.    
        TEMP_TARG(IX2)=TEMP_TARG(IX2)+BASLAM                            
   20 CONTINUE ! Labelled DO for the sake of fpp                        
                                                                       
!     2  For each target lambda line, find the index of the next source 
!     lambda line to the left.                                          
                                                                       
      DO IX2=1,LINES_LAMBDA_TARG                                        
        DO IX1=1,LINES_LAMBDA_SRCE                                      
          IF (TEMP_TARG(IX2).GE.TEMP_SRCE(IX1)) IXL(IX2)=IX1            
        ENDDO                                                           
      ENDDO                                                             
                                                                       
!     3  Find which source lambda gaps cover each target gap and the    
!     fractions they contribute.                                        
!                                                                       
!     At this point IXL(target_line) gives the index of the next source 
!     lambda line to the left of the target lambda line, wrapping round 
!     if the source grid is cyclic. This is also the index of the source
!     gap in which the target line falls. Similarly, the index of the   
!     target line is also that of the target gap of which it is the     
!     left-hand limit. Therefore also IXL(target_gap+1, wrapping round  
!     if reqd.), is the index of the source gap which contains the      
!     right-hand limit of the target gap. For each target gap, we loop  
!     over all source gaps and find the fraction covered by each. Record
!     the fraction and the source index in cumulative lists. If the     
!     source grid is not cyclic, parts of the target gap lying outside  
!     the source grid are neglected.                                    
                                                                       
      IXP=0                                                             
      DO IX2=1,GAPS_LAMBDA_TARG                                         
        IX=0                                                            
        IX2N=MOD(IX2,LINES_LAMBDA_TARG)+1                               
        BTARG=TEMP_TARG(IX2N)-TEMP_TARG(IX2)                            
        IF (BTARG.LT.0.) THEN                                           
          BTARG=BTARG+360.                                              
          IX1N=IXL(IX2N)+LINES_LAMBDA_SRCE                              
        ELSE                                                            
          IX1N=IXL(IX2N)                                                
        ENDIF                                                           
        DO IX1W=IXL(IX2),IX1N                                           
          IX1=MOD(IX1W-1,LINES_LAMBDA_SRCE)+1                           
          IF (CYCLIC_SRCE.OR.IX1.NE.LINES_LAMBDA_SRCE) THEN             
            IF (IX1W.EQ.IXL(IX2)) THEN                                  
              BLO=TEMP_TARG(IX2)                                        
            ELSE                                                        
              BLO=TEMP_SRCE(IX1)                                        
            ENDIF                                                       
            IF (IX1W.EQ.IX1N) THEN                                      
              BHI=TEMP_TARG(IX2N)                                       
            ELSE                                                        
              BHI=TEMP_SRCE(MOD(IX1,LINES_LAMBDA_SRCE)+1)               
            ENDIF                                                       
            IF (BHI.LT.BLO) BHI=BHI+360.                                
            IF (ABS(BHI-BLO).GT.1E-7*ABS(BTARG)) THEN                   
              IX=IX+1                                                   
              INDEX_LAMBDA(IXP+IX)=IX1                                  
              FRAC_LAMBDA(IXP+IX)=(BHI-BLO)/BTARG                       
            ENDIF                                                       
          ENDIF                                                         
        ENDDO                                                           
        COUNT_LAMBDA(IX2)=IX                                            
        IXP=IXP+COUNT_LAMBDA(IX2)                                       
      ENDDO                                                             
                                                                       
!     4  For each target phi line, find the index of the next source phi
!     line above. Comments as for section 2, without wrap-round. Note   
!     that this assumes that the atmosphere gridbox ordering is from    
!     south to north; a north to south atmosphere would require LE      
!     instead of GE in the following IF test.                           
                                                                       
      LINES_PHI_SRCE=GAPS_PHI_SRCE+1                                    
      LINES_PHI_TARG=GAPS_PHI_TARG+1                                    
      DO IY2=1,LINES_PHI_TARG                                           
        IYT(IY2)=0                                                      
        DO IY1=1,LINES_PHI_SRCE                                         
          IF (PHI_TARG(IY2).GE.PHI_SRCE(IY1)) IYT(IY2)=IY1              
        ENDDO                                                           
      ENDDO                                                             
                                                                       
!     5  Find which source phi gaps cover each target gap and the       
!     fractions they contribute. Comments as for section 3, without     
!     wrap-round.                                                       
                                                                       
      IYP=0                                                             
      DO IY2=1,GAPS_PHI_TARG                                            
        IY=0                                                            
        IF (SPHERICAL) THEN                                             
!     Contain angle between +-90. There is no real area outside         
!     these limits on a sphere.                                         
          BTARG=SIN(MAX(MIN(PHI_TARG(IY2+1),90.),-90.)*PI_OVER_180) &
          -SIN(MAX(MIN(PHI_TARG(IY2),90.),-90.)*PI_OVER_180)            
        ELSE                                                            
           BTARG=PHI_TARG(IY2+1)-PHI_TARG(IY2)                          
        ENDIF                                                           
        DO IY1=IYT(IY2),IYT(IY2+1)                                      
          IF (.NOT.(IY1.EQ.0.OR.IY1.EQ.LINES_PHI_SRCE)) THEN            
            IF (IY1.EQ.IYT(IY2)) THEN                                   
              BLO=PHI_TARG(IY2)                                         
            ELSE                                                        
              BLO=PHI_SRCE(IY1)                                         
            ENDIF                                                       
            IF (IY1.EQ.IYT(IY2+1)) THEN                                 
              BHI=PHI_TARG(IY2+1)                                       
            ELSE                                                        
              BHI=PHI_SRCE(IY1+1)                                       
            ENDIF                                                       
            IF (SPHERICAL) THEN                                         
              BLO=MAX(MIN(BLO,90.),-90.)                                
              BHI=MAX(MIN(BHI,90.),-90.)                                
            ENDIF                                                       
            IF (ABS(BHI-BLO).GT.1E-7*ABS(BTARG)) THEN                   
              IY=IY+1                                                   
              INDEX_PHI(IYP+IY)=IY1                                     
!             Both numerator and denominator in the following are -ve.  
              IF (SPHERICAL) THEN                                       
                FRAC_PHI(IYP+IY) = (SIN(BHI*PI_OVER_180)-SIN(BLO*PI_OVER_180))/BTARG      
              ELSE                                                      
                FRAC_PHI(IYP+IY)=(BHI-BLO)/BTARG                        
              ENDIF                                                     
            ENDIF                                                       
          ENDIF                                                         
        ENDDO                                                           
        COUNT_PHI(IY2)=IY                                               
        IYP=IYP+COUNT_PHI(IY2)                                          
      ENDDO                                                             
                                                                       
!     6  For each target box, loop over contributing source boxes and   
!     calculate the weights for each one, ignoring land boxes.          
!                                                                       
!     After the first pass for each target box, go back and normalise   
!     the weights to compensate for land source boxes and any outside   
!     the source grid. Record the source box index and the weight in    
!     cumulative lists.                                                 
                                                                       
      IP=0                                                              
      IYP=0                                                             
      DO IY2=1,GAPS_PHI_TARG                                            
        IXP=0                                                           
        DO IX2=1,GAPS_LAMBDA_TARG                                       
          I=0                                                           
          SUM=0                                                         
          DO IY=IYP+1,IYP+COUNT_PHI(IY2)                                
            DO IX=IXP+1,IXP+COUNT_LAMBDA(IX2)                           
              IF ( MASK_SRCE(INDEX_LAMBDA(IX),INDEX_PHI(IY)) .EQV. WANT) THEN
                I=I+1                                                   
                INDEX_SRCE(IP+I) = INDEX_LAMBDA(IX)  &
                                  + (INDEX_PHI(IY)-1)*GAPS_LAMBDA_SRCE                     
                WEIGHT(IP+I)=FRAC_LAMBDA(IX)*FRAC_PHI(IY)               
                SUM=SUM+WEIGHT(IP+I)                                    
              ENDIF                                                     
            ENDDO                                                       
          ENDDO                                                         
          COUNT_TARG(IX2,IY2)=I                                         
          BASE_TARG(IX2,IY2)=IP                                         
          DO I=1,COUNT_TARG(IX2,IY2)                                    
            WEIGHT(IP+I)=WEIGHT(IP+I)/SUM                               
          ENDDO                                                         
          IP=IP+COUNT_TARG(IX2,IY2)                                     
          IXP=IXP+COUNT_LAMBDA(IX2)                                     
        ENDDO                                                           
        IYP=IYP+COUNT_PHI(IY2)                                          
      ENDDO                                                             
      MAXL=IP                                                           
                                                                       
      ICODE=0                                                           
      CMESSAGE=' '                                                      
      RETURN                                                            

      END SUBROUTINE pre_areaver
