C ErPVp1.f
C
C VIS-5D function to calculate COMPONENT 1 OF PERTURBATION 3-D Ertel PV 
C distribution (ErPVp1). 

C See below ROMU: For personal choices on variables etc...
C ==============================================================================

      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
C ARGUMENTS:
      implicit none
      INTEGER NVARS
      INTEGER NR, NC, NL(NVARS), LOWLEV(NVARS), MAXNL
      REAL OUTGRID(NR,NC,MAXNL)
      INTEGER OUTNL, OUTLOWLEV
      REAL INGRID(NR,NC,MAXNL,NVARS)
      CHARACTER*8 NAMES(NVARS)
      INTEGER DATE, TIME
      INTEGER PROJECTION
      REAL PROJ_ARGS(*)
      INTEGER VERTICAL
      REAL VERT_ARGS(*)

C LOCAL VARS:
      integer var, ir, ic, il
      integer irh, ierpvp 
      real grav, po, cp, rd, kap

      parameter (grav=9.81,po=100000.,cp=1004.5,rd=2.*cp/7.,kap=rd/cp)

      integer ErPVpval_test1(200,200,0:22),ErPVpval_test2(200,200,0:22),
     $        ErPVpreg_test1(200,200,0:22),ErPVpreg_test2(200,200,0:22),
     $        ErPVpreg_test3(200,200,0:22),ErPVpreg_test4(200,200,0:22),      
     $        rh_test1(200,200,0:22),rh_test2(200,200,0:22),     
     $        ph_test1(0:22),ph_test2(0:22),
     $        left_test1(200,200,0:22),left_test2(200,200,0:22),
     $        left_test3(200,200,0:22),left_test4(200,200,0:22),
     $        left_test5(200,200,0:22),left_test6(200,200,0:22)     
      real p(22),ph(0:22),exnf(22),exnfh(0:22),dexnfh(0:22),
     $     rh(200,200,0:22),qp(200,200,0:22),qp1(200,200,0:22)

C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C Find the RH, ErPVp (PERTURBATION field) variables
      irh = -1 
      ierpvp = -1
      do var=1,nvars
        if (names(var) .eq. 'RH') irh = var         
        if (names(var) .eq. 'ErPVp') ierpvp = var            
      end do

C If RH, ErPVp not found, return error 1
      if ( irh.eq.-1 .or. ierpvp.eq.-1 ) then
        print *, "Couldn't find RH, ErPVp vars!"
        userfunc = 1
        return
      endif
      
C ==============================================================================

C NOTE: In VIS5D (ir:N-->S, ic:W-->E, il:B-->T)

C Fill p/ph(Pa), exnf/exnfh(J kg-1 K-1)
       
      do il=1,maxnl
        p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
        exnf(il)=cp*((p(il)/po)**kap)  
      enddo
      do il=1,maxnl-1
        exnfh(il)=0.5*(exnf(il)+exnf(il+1))
        ph(il)=(exnfh(il)*(po**kap)/cp)**(1./kap)
      enddo 
      exnfh(0)=2.*exnf(1)-exnfh(1)
      ph(0)=(exnfh(0)*(po**kap)/cp)**(1./kap)      
      exnfh(maxnl)=2.*exnf(maxnl)-exnfh(maxnl-1)      
      ph(maxnl)=(exnfh(maxnl)*(po**kap)/cp)**(1./kap)      

C Definition dexnfh
      dexnfh(0)=exnfh(1)-exnfh(0)
      do 120 il=1,maxnl-1
       dexnfh(il)=exnfh(il+1)-exnfh(il) 
120   continue      

C ==============================================================================

C Retrieve rh (%) and perturbation qp (PVU), from which to define component qp1
C for all vertical levels il: 1 --> maxnl-1 (ph levels)
    
      Do 33 ir=1,nr
      Do 33 ic=1,nc
        qp(ir,ic,0)=ingrid(ir,ic,1,ierpvp)
      Do 33 il=1,maxnl-1
        qp(ir,ic,il)=qp(ir,ic,il-1)+dexnfh(il-1)*
     $               (ingrid(ir,ic,il,ierpvp)-
     $               qp(ir,ic,il-1))/(exnf(il)-exnfh(il-1))
        rh(ir,ic,il)=0.5*(ingrid(ir,ic,il,irh)+ingrid(ir,ic,il+1,irh))
33    Continue
      print*,'Retrieved rh and qp'

C ==============================================================================

C CALCULATION OF COMPONENT 1 OF PERTURBATION ERTEL POTENTIAL VORTICITY qp1 
C (PVU) ==========================> ROMU: Define condition/function of qp1 below 

C ------------------------------------------------------------------------------

C ROMU: Parametros SUPERSTORM
c      Call ErPVpval_test(nr,nc,maxnl,qp,
c     $                   0.,100.,ErPVpval_test1)
c      Call rh_test(nr,nc,maxnl,rh,
c     $             70.,100.,rh_test1)
c      Call ph_test(maxnl,ph,
c     $             700.,1000.,ph_test1)
c      Call ph_test(maxnl,ph,
c     $             500.,1000.,ph_test2)
     
C ------------------------------------------------------------------------------

C ROMU: Parametros VALENCIA
      Call ErPVpval_test(nr,nc,maxnl,qp,
     $                   0.,100.,ErPVpval_test1)
      Call ErPVpreg_test(nr,nc,maxnl,qp,
     $                   52,13,0.,100.,ErPVpreg_test1)
      Call ErPVpreg_test(nr,nc,maxnl,qp,
     $                   35,38,0.,100.,ErPVpreg_test2)  
      Call ErPVpreg_test(nr,nc,maxnl,qp,
     $                   29,49,0.,100.,ErPVpreg_test3)          
      Call rh_test(nr,nc,maxnl,rh,
     $             0.,100.,rh_test1)
      Call ph_test(maxnl,ph,
     $             0.,500.,ph_test1)
      Call left_test(nr,nc,maxnl,ph,
     $               81.00,67.81,14.19,1.00,200.,
     $               81.00,63.68,18.40,1.00,400.,
     $               left_test1)   
      Call left_test(nr,nc,maxnl,ph,
     $               81.00,5.61,5.61,81.00,100.,
     $               81.00,12.10,12.10,81.00,400.,
     $               left_test2)   
      Call left_test(nr,nc,maxnl,ph,
     $               45.82,81.00,1.00,36.18,100.,
     $               53.15,81.00,1.00,28.85,300.,
     $               left_test3)   
      Call left_test(nr,nc,maxnl,ph,
     $               81.00,8.00,8.00,81.00,150.,
     $               81.00,6.45,6.45,81.00,300.,
     $               left_test4) 
      Call left_test(nr,nc,maxnl,ph,
     $               81.00,36.24,1.00,67.12,100.,
     $               81.00,36.24,1.00,67.12,1000.,
     $               left_test5)  
      Call left_test(nr,nc,maxnl,ph,
     $               1.00,72.18,38.01,1.00,100.,
     $               1.00,72.18,38.01,1.00,1000.,
     $               left_test6)   

C ------------------------------------------------------------------------------
                            
      DO 20 ir=1,nr
      DO 20 ic=1,nc
      DO 20 il=1,maxnl-1

C DEFAULT/INITIAL qp1 = 0 (no comment !!!)
        qp1(ir,ic,il)=0.

C qp1 = ErPVp (entire)
c        qp1(ir,ic,il)=qp(ir,ic,il)

C qp1 = ErPVp (upper-level)
c        If (ph(il).Le.50000.) 
c     $  qp1(ir,ic,il)=qp(ir,ic,il)
     
C qp1 = ErPVp (lower-level)
c        If (ph(il).Gt.50000.) 
c     $  qp1(ir,ic,il)=qp(ir,ic,il)     

C ------------------------------------------------------------------------------

C qp1 SUPERSTORM ULev
c         qp1(ir,ic,il)=qp(ir,ic,il)*(1-ph_test1(il))*
c     $                 ( 1-ph_test2(il)*rh_test1(ir,ic,il)*
c     $                 ErPVpval_test1(ir,ic,il) )

C qp1 SUPERSTORM LBef
c         qp1(ir,ic,il)=qp(ir,ic,il)*ph_test1(il)*
c     $                 ( 1-rh_test1(ir,ic,il)*
c     $                 ErPVpval_test1(ir,ic,il) )

C qp1 SUPERSTORM DIAB
c         qp1(ir,ic,il)=qp(ir,ic,il)*ph_test2(il)*
c     $                              rh_test1(ir,ic,il)*
c     $                              ErPVpval_test1(ir,ic,il)    
     
C ------------------------------------------------------------------------------

C qp1 VALENCIA SW
         qp1(ir,ic,il)=qp(ir,ic,il)*ErPVpreg_test1(ir,ic,il)*
     $                              ph_test1(il)*
     $                              left_test1(ir,ic,il)*
     $                              left_test2(ir,ic,il)

C qp1 VALENCIA NE
c         qp1(ir,ic,il)=qp(ir,ic,il)*Min(
c     $                              ErPVpreg_test2(ir,ic,il)+
c     $                              ErPVpreg_test3(ir,ic,il),1)*
c     $                              ph_test1(il)*
c     $                              Abs(left_test1(ir,ic,il)-1)*
c     $                              left_test3(ir,ic,il)*
c     $                              left_test4(ir,ic,il)*
c     $                              left_test5(ir,ic,il)*
c     $                              left_test6(ir,ic,il)
                                    
C qp1 VALENCIA SW+NE
c         qp1(ir,ic,il)=qp(ir,ic,il)*ErPVpreg_test1(ir,ic,il)*
c     $                              ph_test1(il)*
c     $                              left_test1(ir,ic,il)*
c     $                              left_test2(ir,ic,il)
c     $                              +
c     $                 qp(ir,ic,il)*Min(
c     $                              ErPVpreg_test2(ir,ic,il)+
c     $                              ErPVpreg_test3(ir,ic,il),1)*
c     $                              ph_test1(il)*
c     $                              Abs(left_test1(ir,ic,il)-1)*
c     $                              left_test3(ir,ic,il)*
c     $                              left_test4(ir,ic,il)*
c     $                              left_test5(ir,ic,il)*
c     $                              left_test6(ir,ic,il)

C qp1 VALENCIA RESTO
c         qp1(ir,ic,il)=qp(ir,ic,il)-
c     $                 qp(ir,ic,il)*ErPVpreg_test1(ir,ic,il)*
c     $                              ph_test1(il)*
c     $                              left_test1(ir,ic,il)*
c     $                              left_test2(ir,ic,il)
c     $                              -
c     $                 qp(ir,ic,il)*Min(
c     $                              ErPVpreg_test2(ir,ic,il)+
c     $                              ErPVpreg_test3(ir,ic,il),1)*
c     $                              ph_test1(il)*
c     $                              Abs(left_test1(ir,ic,il)-1)*
c     $                              left_test3(ir,ic,il)*
c     $                              left_test4(ir,ic,il)*
c     $                              left_test5(ir,ic,il)*
c     $                              left_test6(ir,ic,il)

C ------------------------------------------------------------------------------
     
C BOTTOM/TOP B.C (no comment !!!)
        If (il.Eq.1) qp1(ir,ic,0)=qp1(ir,ic,il)
        If (il.Eq.maxnl-1) qp1(ir,ic,maxnl)=qp1(ir,ic,il)     
  
20    CONTINUE

      print*,'Calculated qp1'
 
C ==============================================================================

C OUTPUT OF COMPONENT 1 OF PERTURBATION 3-D ERTEL PV DISTRIBUTION ErPVp1
C (in PVU)

      print*,
     $'Output of COMPONENT 1 OF PERTURBATION 3-D Ertel PV ',
     $'distribution ErPVp1 (PVU)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
      DO 9999 il=1,maxnl
        outgrid(ir,ic,il)=qp1(ir,ic,il)-(exnfh(il)-exnf(il))*
     $                    (qp1(ir,ic,il)-qp1(ir,ic,il-1))/
     $                    dexnfh(il-1)
9999  CONTINUE    

C ==============================================================================
               
      USERFUNC = 0
      RETURN
      END

C ==============================================================================


C ==============================================================================
C This subroutine determines, for all vertical levels il: 1 --> maxnl-1,
C grid points where ErPVmin <= qp <= ErPVmax (ErPVmin and ErPVmax in PVU). 
C OUTPUT: Yes (test=1) / No (test=0)

      SUBROUTINE ErPVpval_test(nr,nc,maxnl,qp,
     $                        ErPVmin,ErPVmax,test)

        implicit none
        integer ir, ic, il
        INTEGER NR,NC,MAXNL       
        real ErPVmin, ErPVmax
        integer test(200,200,0:22)         
        real qp(200,200,0:22) 

        Do il=1,maxnl-1
          Do ir=1,nr
          Do ic=1,nc
            test(ir,ic,il)=0
            If (ErPVmin.Le.qp(ir,ic,il).And.qp(ir,ic,il).Le.ErPVmax)
     $      test(ir,ic,il)=1             
          Enddo
          Enddo
        Enddo
        
      Return
      END

C ==============================================================================

C ==============================================================================
C This subroutine determines, for all vertical levels il: 1 --> maxnl-1,
C grid points lying inside continuous horizontal regions in which 
C ErPVmin <= qp <= ErPVmax (ErPVmin and ErPVmax in PVU), the regions defined 
C from point (ir0,ic0). 
C OUTPUT: Yes (test=1) / No (test=0)

      SUBROUTINE ErPVpreg_test(nr,nc,maxnl,qp,
     $                      ir0,ic0,ErPVmin,ErPVmax,test)

        implicit none
        integer ir, ic, il, ir0, ic0 
        INTEGER NR,NC,MAXNL       
        real ErPVmin, ErPVmax
        integer test(200,200,0:22)         
        real qp(200,200,0:22) 

        Do il=1,maxnl-1
          Do ir=1,nr
          Do ic=1,nc
            test(ir,ic,il)=0
          Enddo
          Enddo
        Enddo
        If (ir0.Lt.1.Or.ir0.Gt.nr.Or.ic0.Lt.1.Or.ic0.Gt.nc) Then
          Print *,'  '
          Print *,'ERROR: (ir0,ic0) in ErPVpreg_test outside domain'
          Print *,'test = 0'
          Return
        Endif
        Do il=1,maxnl-1
          If (ErPVmin.Le.qp(ir0,ic0,il).And.
     $        qp(ir0,ic0,il).Le.ErPVmax)
     $    test(ir0,ic0,il)=1
        Enddo

        DO 1000 il=1,maxnl-1
                
          Do ic=ic0-1,1,-1
            If ((test(ir0,ic+1,il).Eq.1).And.
     $          (ErPVmin.Le.qp(ir0,ic,il).And.
     $          qp(ir0,ic,il).Le.ErPVmax)) 
     $      test(ir0,ic,il)=1
          Enddo

          Do ic=ic0+1,nc
            If ((test(ir0,ic-1,il).Eq.1).And.
     $          (ErPVmin.Le.qp(ir0,ic,il).And.
     $          qp(ir0,ic,il).Le.ErPVmax)) 
     $      test(ir0,ic,il)=1
          Enddo          
          
          Do ir=1,nr
          Do ic=1,nc
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  

          Do ir=1,nr
          Do ic=nc,1,-1
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  
          
          Do ic=1,nc
          Do ir=1,nr
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  
          
          Do ic=1,nc
          Do ir=nr,1,-1
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo                      

          Do ir=nr,1,-1
          Do ic=1,nc
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  

          Do ir=nr,1,-1
          Do ic=nc,1,-1
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  
          
          Do ic=nc,1,-1
          Do ir=1,nr
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo  
          
          Do ic=nc,1,-1
          Do ir=nr,1,-1
            If ((test(ir,max(ic-1,1),il).Eq.1.Or.
     $           test(max(ir-1,1),max(ic-1,1),il).Eq.1.Or.     
     $           test(max(ir-1,1),ic,il).Eq.1.Or.
     $           test(max(ir-1,1),min(ic+1,nc),il).Eq.1.Or.
     $           test(ir,min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),min(ic+1,nc),il).Eq.1.Or.
     $           test(min(ir+1,nr),ic,il).Eq.1.Or.
     $           test(min(ir+1,nr),max(ic-1,1),il).Eq.1).And.                                     
     $          (ErPVmin.Le.qp(ir,ic,il).And.
     $           qp(ir,ic,il).Le.ErPVmax)) 
     $      test(ir,ic,il)=1 
          Enddo
          Enddo       

          
1000    CONTINUE

      Return
      END

C ==============================================================================

C ==============================================================================
C This subroutine determines, for all vertical levels il: 1 --> maxnl-1,
C grid points where rhmin <= rh <= rhmax (rhmin and rhmax in %). 
C OUTPUT: Yes (test=1) / No (test=0)

      SUBROUTINE rh_test(nr,nc,maxnl,rh,
     $                   rhmin,rhmax,test)

        implicit none
        integer ir, ic, il
        INTEGER NR,NC,MAXNL       
        real rhmin, rhmax
        integer test(200,200,0:22)         
        real rh(200,200,0:22) 

        Do il=1,maxnl-1
          Do ir=1,nr
          Do ic=1,nc
            test(ir,ic,il)=0
            If (rhmin.Le.rh(ir,ic,il).And.rh(ir,ic,il).Le.rhmax)
     $      test(ir,ic,il)=1             
          Enddo
          Enddo
        Enddo
        
      Return
      END

C ==============================================================================
        
        
C ==============================================================================
C This subroutine determines, for all vertical levels il: 1 --> maxnl-1,
C whether phmin <= ph(il) <= phmax (phmin and phmax in hPa). 
C OUTPUT: Yes (test=1) / No (test=0)

      SUBROUTINE ph_test(maxnl,ph,
     $                   phmin,phmax,test)

        implicit none
        integer il
        INTEGER MAXNL       
        real phmin, phmax
        integer test(0:22)         
        real ph(0:22) 

        Do il=1,maxnl-1
          test(il)=0
          If (100*phmin.Le.ph(il).And.ph(il).Le.100*phmax)
     $    test(il)=1 
        Enddo
     
      Return
      END

C ==============================================================================

C ==============================================================================
C This subroutine determines, for all vertical levels il: 1 --> maxnl-1,
C grid points lying to the left of the line (r1,c1)-----(r2,c2) (coordinates
C in grid units)
C OUTPUT: Yes (test=1) / No (test=0)

      SUBROUTINE left_test(nr,nc,maxnl,ph,
     $                     r1min,c1min,r2min,c2min,phmin,
     $                     r1max,c1max,r2max,c2max,phmax,
     $                     test)

        implicit none
        integer ir, ic, il
        INTEGER NR,NC,MAXNL       
        real r1min,c1min,r2min,c2min,phmin,
     $       r1max,c1max,r2max,c2max,phmax,
     $       r1,c1,r2,c2,
     $       rl,cl,rp,cp,cross
        integer test(200,200,0:22)         
        real ph(0:22) 
                
        DO 1000 il=1,maxnl-1
        
          r1=r1min+(ph(il)/100.-phmin)*(r1max-r1min)/(phmax-phmin)
          c1=c1min+(ph(il)/100.-phmin)*(c1max-c1min)/(phmax-phmin)           
          r2=r2min+(ph(il)/100.-phmin)*(r2max-r2min)/(phmax-phmin)
          c2=c2min+(ph(il)/100.-phmin)*(c2max-c2min)/(phmax-phmin)           
          rl=r2-r1
          cl=c2-c1

          Do ir=1,nr
            rp=Real(ir)-r1   
          Do ic=1,nc
            cp=Real(ic)-c1  
            test(ir,ic,il)=0
            cross=rl*cp-rp*cl
            If (cross.Ge.0.) test(ir,ic,il)=1        
          Enddo          
          Enddo

1000    CONTINUE 

      Return
      END

C ==============================================================================
        
