C I_Sp1.f
C
C VIS-5D function to calculate COMPONENT 1 OF PERTURBATION streamfunction 
C (I_Sp1) by inversion of COMPONENT 1 OF PERTURBATION 3-D Ertel PV 
C distribution. 

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 ii_hm, ii_hp, ii_sm, ii_sp, ierpvp1, itbp1, ittp1,
     $        it, itm, icor, imap
      integer max_cycles,icycle,max_iters,iter,noconverg,
     $        iter_acum_phi,iter_acum_psi
      real acum_o,acum_f       
      real grav, po, cp, rd, kap, ds, ds2, dssqr

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

      real thrs_phi,thrs_psi,over_phi,over_psi,under,
     $     res,resmod
      real Fndim,PI2ndim,PIndim,Lndim,Sndim
      real factor,
     $     phixx,phiyy,phipp,phixp1,phixp2,phixp,phiyp1,phiyp2,phiyp,
     $     psix,psiy,psixx,psiyy,psixy,psixp,psiyp,relvor       
      real phibarpp(200,200,22),
     $     phibarxp1,phibarxp2,phibarxp(200,200,22),
     $     phibaryp1,phibaryp2,phibaryp(200,200,22),
     $     psibarxx(200,200,22),psibaryy(200,200,22),
     $     psibarxy(200,200,22),
     $     psibarxp(200,200,22),psibaryp(200,200,22),
     $     relvorbar(200,200,22)     
      real p(22),ph(0:22),exnf(22),exnfh(0:22),dexnf(22),dexnfh(0:22),
     $     ppl(22),pplm1(22),pplp1(22),
     $     cor(200,200),corx(200,200),cory(200,200),
     $     map(200,200),map2(200,200),
     $     q(200,200,0:22),
     $     phi0(200,200,0:22),phi(200,200,0:22),phibar(200,200,0:22),
     $     for_phi(200,200,22),coef_phi(200,200,22),
     $     psi0(200,200,0:22),psi(200,200,0:22),psibar(200,200,0:22),
     $     for_psi(200,200,22),
     $     thetaB(200,200),thetaT(200,200),
     $     thetaBm,thetaBp 

C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C Find the I_Hm, I_Hp, I_Sm, I_Sp, ErPVp1, TBp1, TTp1, 
C          T, Tm, CORIOLIS  and MAPFACCR variables 
      ii_hm = -1
      ii_hp = -1
      ii_sm = -1
      ii_sp = -1
      ierpvp1 = -1
      itbp1 = -1
      ittp1 = -1
      it = -1
      itm = -1
      icor = -1
      imap = -1
      do var=1,nvars
        if (names(var) .eq. 'I_Hm') ii_hm = var
        if (names(var) .eq. 'I_Hp') ii_hp = var
        if (names(var) .eq. 'I_Sm') ii_sm = var
        if (names(var) .eq. 'I_Sp') ii_sp = var                
        if (names(var) .eq. 'ErPVp1') ierpvp1 = var
        if (names(var) .eq. 'TBp1') itbp1 = var
        if (names(var) .eq. 'TTp1') ittp1 = var 
        if (names(var) .eq. 'T') it = var   
        if (names(var) .eq. 'Tm') itm = var                        
        if (names(var) .eq. 'CORIOLIS') icor = var
        if (names(var) .eq. 'MAPFACCR') imap = var                  
      end do

C If I_Hm, I_Hp, I_Sm, I_Sp, ErPVp1, TBp1, TTp1, 
C    T, Tm, CORIOLIS or MAPFACCR not found, return error 1
      if ( ii_hm.eq.-1 .or. ii_hp.eq.-1 .or.
     $     ii_sm.eq.-1 .or. ii_sp.eq.-1 .or.
     $     ierpvp1.eq.-1 .or. itbp1.eq.-1 .or. ittp1.eq.-1 .or. 
     $     it.eq.-1 .or. itm.eq.-1 .or. 
     $     icor.eq.-1 .or. imap.eq.-1 ) then
        print *, 
     $  "Couldn't find I_Hm, I_Hp, I_Sm, I_Sp, ErPVp1, TBp1, ", 
     $   "TTp1, T, Tm, CORIOLIS or MAPFACCR vars!"
        userfunc = 1
        return
      endif

C Grid-length ds (m)
      if (projection.Eq.2) then
        ds=proj_args(6)*1000.      
      else if (projection.Eq.3) then  
        ds=proj_args(5)*1000.
      else
        print *, "Uncomprensible projection data!"
        userfunc = 1
        return
      endif

C Definition max_cycles,max_iters,thrs_phi,thrs_psi,over_phi,over_psi,under,
C ds2, dssqr
      max_cycles=100
      max_iters=500      
      thrs_phi=1.
      thrs_psi=1.E4
      over_phi=2.-3.1415927*(2./Real(nc-1)**2.+2./Real(nr-1)**2.+
     $         2./Real(maxnl-1)**2.)**.5
      over_psi=2.-3.1415927*(2./Real(nc-1)**2.+2./Real(nr-1)**2.)**.5
      under=0.7
      ds2=2.*ds
      dssqr=ds*ds
      
C Constants for nondimensionalization of the system
      Fndim=1.E-4
      PI2ndim=500./Float(maxnl)
      PIndim=PI2ndim**0.5
      Lndim=ds
      Sndim=Fndim*Fndim*Lndim*Lndim/PI2ndim  
      
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), cor(s-1) and map(adim.)
       
      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)
               
      do 10 ir=1,nr
      do 10 ic=1,nc
        cor(ir,ic)=ingrid(ir,ic,1,icor)
        map(ir,ic)=ingrid(ir,ic,1,imap) 
        map2(ir,ic)=map(ir,ic)*map(ir,ic)           
10    continue   

C Definition of Neumann boundary conditions thetaB and thetaT (K)
      do 100 ir=1,nr
      do 100 ic=1,nc
        thetaB(ir,ic)=ingrid(ir,ic,1,itbp1)*
     $                cp/exnf(1)
        thetaT(ir,ic)=ingrid(ir,ic,1,ittp1)*
     $                cp/exnf(maxnl)             
100   continue      

C Definition dexnf, dexnfh, pplp1, pplm1 and ppl
      dexnfh(0)=exnfh(1)-exnfh(0)
      do 120 il=1,maxnl-1
       dexnf(il)=exnf(il+1)-exnf(il)
       dexnfh(il)=exnfh(il+1)-exnfh(il) 
       pplp1(il)=1./(dexnfh(il)*dexnf(il))           
       pplm1(il)=1./(dexnfh(il-1)*dexnf(il))                     
       ppl(il)=pplp1(il)+pplm1(il)    
120   continue      

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

C Retrieve I_Hm and I_Hp (phi0 and phi temporaly), used to compute phibar
  
      Do 11 ir=1,nr
      Do 11 ic=1,nc
        thetaBm=(ingrid(ir,ic,1,itm)+273.16)*cp/exnf(1)
        phi0(ir,ic,0)=grav*ingrid(ir,ic,1,ii_hm)+
     $                thetaBm*(exnf(1)-exnfh(0))
        thetaBp=(ingrid(ir,ic,1,it)+273.16)*cp/exnf(1)-thetaBm
        phi(ir,ic,0)=grav*ingrid(ir,ic,1,ii_hp)+
     $               thetaBp*(exnf(1)-exnfh(0)) 
        phibar(ir,ic,0)=phi0(ir,ic,0)+0.5*phi(ir,ic,0)            
      Do 11 il=1,maxnl
        phi0(ir,ic,il)=phi0(ir,ic,il-1)+dexnfh(il-1)*
     $                 (grav*ingrid(ir,ic,il,ii_hm)-
     $                 phi0(ir,ic,il-1))/(exnf(il)-exnfh(il-1))
        phi(ir,ic,il)=phi(ir,ic,il-1)+dexnfh(il-1)*
     $                 (grav*ingrid(ir,ic,il,ii_hp)-
     $                 phi(ir,ic,il-1))/(exnf(il)-exnfh(il-1))     
        phibar(ir,ic,il)=phi0(ir,ic,il)+0.5*phi(ir,ic,il)      
11    Continue
      print*,'Calculated PHIBAR'

C SET OF FIRST-GUESS GEOPOTENTIAL phi0 (m2 s-2)    

      Do 30 ir=1,nr
      Do 30 ic=1,nc
      Do 30 il=0,maxnl

C ROMU: Define type of B.C on GEOPOTENTIAL below

C Homogeneous B.C      
        phi0(ir,ic,il)=0.

C B.C = I_Hp
c        phi0(ir,ic,il)=0.5*phi(ir,ic,il)

C SUPERSTORM
c        phi0(ir,ic,il)=0.
c	If ((ph(il)/100.).Le.700.) phi0(ir,ic,il)=phi(ir,ic,il)

30    Continue

      print*,'Set PHI0'
 
C ==============================================================================

C Retrieve I_Sm and I_Sp (psi0 and psi temporaly), used to compute psibar
  
      Do 22 ir=1,nr
      Do 22 ic=1,nc
        thetaBm=(ingrid(ir,ic,1,itm)+273.16)*cp/exnf(1)
        psi0(ir,ic,0)=(grav/Fndim)*ingrid(ir,ic,1,ii_sm)+
     $                (thetaBm/Fndim)*(exnf(1)-exnfh(0))
        thetaBp=(ingrid(ir,ic,1,it)+273.16)*cp/exnf(1)-thetaBm
        psi(ir,ic,0)=(grav/Fndim)*ingrid(ir,ic,1,ii_sp)+
     $               (thetaBp/Fndim)*(exnf(1)-exnfh(0)) 
        psibar(ir,ic,0)=psi0(ir,ic,0)+0.5*psi(ir,ic,0)            
      Do 22 il=1,maxnl
        psi0(ir,ic,il)=psi0(ir,ic,il-1)+dexnfh(il-1)*
     $                 ((grav/Fndim)*ingrid(ir,ic,il,ii_sm)-
     $                 psi0(ir,ic,il-1))/(exnf(il)-exnfh(il-1))
        psi(ir,ic,il)=psi(ir,ic,il-1)+dexnfh(il-1)*
     $                 ((grav/Fndim)*ingrid(ir,ic,il,ii_sp)-
     $                 psi(ir,ic,il-1))/(exnf(il)-exnfh(il-1))     
        psibar(ir,ic,il)=psi0(ir,ic,il)+0.5*psi(ir,ic,il)      
22    Continue
      print*,'Calculated PSIBAR'

C SET OF FIRST GUESS STREAMFUNCTION psi0 (m2 s-1)

      Do 60 ir=1,nr
      Do 60 ic=1,nc
      Do 60 il=0,maxnl

C ROMU: Define type of B.C on STREAMFUNCTION below

C Homogeneous B.C      
        psi0(ir,ic,il)=0.

C B.C = I_Sp
c        psi0(ir,ic,il)=psi(ir,ic,il)

C SUPERSTORM
c        psi0(ir,ic,il)=0.
c	If ((ph(il)/100.).Le.700.) psi0(ir,ic,il)=psi(ir,ic,il)

60    Continue

      print*,'Set PSI0'

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

C Retrieve ErPVp1 (q, in m2 K s-1 kg-1)
      
      Do 33 ir=1,nr
      Do 33 ic=1,nc
        q(ir,ic,0)=1E-6*ingrid(ir,ic,1,ierpvp1)
      Do 33 il=1,maxnl
        q(ir,ic,il)=q(ir,ic,il-1)+dexnfh(il-1)*
     $              (1E-6*ingrid(ir,ic,il,ierpvp1)-
     $              q(ir,ic,il-1))/(exnf(il)-exnfh(il-1))
33    Continue

      print*,'Retrieved q'

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

C NON-DIMENSIONALIZATION OF SYSTEM and psi=psi0, phi=phi0, corx,cory
      ds=ds/Lndim
      ds2=2.*ds
      dssqr=ds*ds
      thrs_phi=thrs_phi/(Sndim*PI2ndim)
      thrs_psi=Fndim*thrs_psi/(Sndim*PI2ndim)
      
      do 910 ir=1,nr
      do 910 ic=1,nc      
        cor(ir,ic)=cor(ir,ic)/Fndim
        thetaB(ir,ic)=thetaB(ir,ic)/(Sndim*PIndim)
        thetaT(ir,ic)=thetaT(ir,ic)/(Sndim*PIndim)
      do 910 il=0,maxnl
        psibar(ir,ic,il)=Fndim*psibar(ir,ic,il)/(Sndim*PI2ndim)
        psi0(ir,ic,il)=Fndim*psi0(ir,ic,il)/(Sndim*PI2ndim)
        psi(ir,ic,il)=psi0(ir,ic,il)
        phibar(ir,ic,il)=phibar(ir,ic,il)/(Sndim*PI2ndim)        
        phi0(ir,ic,il)=phi0(ir,ic,il)/(Sndim*PI2ndim)
        phi(ir,ic,il)=phi0(ir,ic,il)
        factor=grav*kap*exnfh(il)/ph(il)
        q(ir,ic,il)=q(ir,ic,il)/(factor*Sndim*Fndim)                
910   continue     

      do 920 ir=2,nr-1
      do 920 ic=2,nc-1 
        corx(ir,ic)=map(ir,ic)*(cor(ir,ic+1)-cor(ir,ic-1))/ds2
        cory(ir,ic)=map(ir,ic)*(cor(ir-1,ic)-cor(ir+1,ic))/ds2        
920   continue

      exnfh(0)=exnfh(0)/PIndim
      dexnfh(0)=dexnfh(0)/PIndim
      do 930 il=1,maxnl-1
        exnf(il)=exnf(il)/PIndim
        dexnf(il)=dexnf(il)/PIndim     
        exnfh(il)=exnfh(il)/PIndim
        dexnfh(il)=dexnfh(il)/PIndim        
        pplp1(il)=PI2ndim*pplp1(il)
        pplm1(il)=PI2ndim*pplm1(il)
        ppl(il)=PI2ndim*ppl(il)                         
930   continue 
      exnf(maxnl)=exnf(maxnl)/PIndim     
      exnfh(maxnl)=exnfh(maxnl)/PIndim

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

C Terms of the equations depending on phibar, psibar (constant)

      Do 2000 ir=2,nr-1
      Do 2000 ic=2,nc-1
      Do 2000 il=1,maxnl-1

        phibarpp(ir,ic,il)=pplp1(il)*phibar(ir,ic,il+1)+
     $                     pplm1(il)*phibar(ir,ic,il-1)-
     $                     ppl(il)*phibar(ir,ic,il) 
        phibarpp(ir,ic,il)=Max(1.E-4,phibarpp(ir,ic,il))    
  
        phibarxp1=0.5*(phibar(ir,ic+1,il+1)/dexnfh(il)-
     $            phibar(ir,ic+1,il-1)/dexnfh(il-1)+
     $            phibar(ir,ic+1,il)*
     $            (-1./dexnfh(il)+1./dexnfh(il-1)))    
        phibarxp2=0.5*(phibar(ir,ic-1,il+1)/dexnfh(il)-
     $            phibar(ir,ic-1,il-1)/dexnfh(il-1)+
     $            phibar(ir,ic-1,il)*
     $            (-1./dexnfh(il)+1./dexnfh(il-1)))   
        phibarxp(ir,ic,il)=map(ir,ic)*(phibarxp1-
     $                     phibarxp2)/ds2 

        phibaryp1=0.5*(phibar(ir-1,ic,il+1)/dexnfh(il)-
     $            phibar(ir-1,ic,il-1)/dexnfh(il-1)+
     $            phibar(ir-1,ic,il)*
     $            (-1./dexnfh(il)+1./dexnfh(il-1)))   
        phibaryp2=0.5*(phibar(ir+1,ic,il+1)/dexnfh(il)-
     $            phibar(ir+1,ic,il-1)/dexnfh(il-1)+
     $            phibar(ir+1,ic,il)*
     $            (-1./dexnfh(il)+1./dexnfh(il-1)))   
        phibaryp(ir,ic,il)=map(ir,ic)*(phibaryp1-
     $                     phibaryp2)/ds2

        psibarxx(ir,ic,il)=map2(ir,ic)*(psibar(ir,ic+1,il)+
     $                     psibar(ir,ic-1,il)-
     $                     2.*psibar(ir,ic,il))/dssqr          
        psibaryy(ir,ic,il)=map2(ir,ic)*(psibar(ir-1,ic,il)+
     $                     psibar(ir+1,ic,il)-
     $                     2.*psibar(ir,ic,il))/dssqr    
        relvorbar(ir,ic,il)=psibarxx(ir,ic,il)+
     $                      psibaryy(ir,ic,il)        
        relvorbar(ir,ic,il)=Max(relvorbar(ir,ic,il),
     $                      -cor(ir,ic)+1.E-4)   

        psibarxy(ir,ic,il)=map2(ir,ic)*(psibar(ir-1,ic+1,il)-
     $                     psibar(ir-1,ic-1,il)-
     $                     psibar(ir+1,ic+1,il)+
     $                     psibar(ir+1,ic-1,il))/
     $                     (ds2*ds2)       
        psibarxp(ir,ic,il)=map(ir,ic)*(psibar(ir,ic+1,il+1)-
     $                     psibar(ir,ic-1,il+1)-
     $                     psibar(ir,ic+1,il-1)+
     $                     psibar(ir,ic-1,il-1))/
     $                     (ds2*(exnfh(il+1)-exnfh(il-1)))   
        psibaryp(ir,ic,il)=map(ir,ic)*(psibar(ir-1,ic,il+1)-
     $                     psibar(ir+1,ic,il+1)-
     $                     psibar(ir-1,ic,il-1)+
     $                     psibar(ir+1,ic,il-1))/
     $                     (ds2*(exnfh(il+1)-exnfh(il-1)))      

2000  Continue

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

C START CYCLES FOR PV INVERSION (get phi and psi iteratively)

      acum_f=1000000.
       
      DO 1000 icycle=1,max_cycles

       iter_acum_phi=0
       iter_acum_psi=0
       acum_o=acum_f     
       
       print*
       print*,'CYCLE=',icycle

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

C SOLVE GEOPOTENTIAL phi (m2 s-2)      

C ------------------------------------------------------------------------------
      
C Forcing function for_phi (and coef_phi)   
       Do 2100 ir=2,nr-1
       Do 2100 ic=2,nc-1
       Do 2100 il=1,maxnl-1

         psix=map(ir,ic)*(psi(ir,ic+1,il)-psi(ir,ic-1,il))/                   
     $        ds2 
         psiy=map(ir,ic)*(psi(ir-1,ic,il)-psi(ir+1,ic,il))/                   
     $        ds2   
         psixx=map2(ir,ic)*(psi(ir,ic+1,il)+psi(ir,ic-1,il)-
     $         2.*psi(ir,ic,il))/dssqr      
         psiyy=map2(ir,ic)*(psi(ir-1,ic,il)+psi(ir+1,ic,il)-
     $         2.*psi(ir,ic,il))/dssqr  

         relvor=psixx+psiyy      
           
         psixy=map2(ir,ic)*(psi(ir-1,ic+1,il)-psi(ir-1,ic-1,il)-
     $         psi(ir+1,ic+1,il)+psi(ir+1,ic-1,il))/
     $         (ds2*ds2)       
         psixp=map(ir,ic)*(psi(ir,ic+1,il+1)-psi(ir,ic-1,il+1)-
     $         psi(ir,ic+1,il-1)+psi(ir,ic-1,il-1))/
     $         (ds2*(exnfh(il+1)-exnfh(il-1)))      
         psiyp=map(ir,ic)*(psi(ir-1,ic,il+1)-psi(ir+1,ic,il+1)-
     $         psi(ir-1,ic,il-1)+psi(ir+1,ic,il-1))/
     $         (ds2*(exnfh(il+1)-exnfh(il-1)))         

         phixp1=0.5*(phi0(ir,ic+1,il+1)/dexnfh(il)-
     $          phi0(ir,ic+1,il-1)/dexnfh(il-1)+phi0(ir,ic+1,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))         
         phixp2=0.5*(phi0(ir,ic-1,il+1)/dexnfh(il)-
     $          phi0(ir,ic-1,il-1)/dexnfh(il-1)+phi0(ir,ic-1,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))        
         phixp=map(ir,ic)*(phixp1-phixp2)/ds2 
        
         phiyp1=0.5*(phi0(ir-1,ic,il+1)/dexnfh(il)-
     $          phi0(ir-1,ic,il-1)/dexnfh(il-1)+phi0(ir-1,ic,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))       
         phiyp2=0.5*(phi0(ir+1,ic,il+1)/dexnfh(il)-
     $          phi0(ir+1,ic,il-1)/dexnfh(il-1)+phi0(ir+1,ic,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))               
         phiyp=map(ir,ic)*(phiyp1-phiyp2)/ds2

         for_phi(ir,ic,il)=(q(ir,ic,il)+
     $                     (cor(ir,ic)-phibarpp(ir,ic,il))*relvor+
     $                     corx(ir,ic)*psix+cory(ir,ic)*psiy+
     $                     2.*psibarxx(ir,ic,il)*psiyy+
     $                     2.*psibaryy(ir,ic,il)*psixx-
     $                     4.*psibarxy(ir,ic,il)*psixy+
     $                     psixp*phibarxp(ir,ic,il)+
     $                     psiyp*phibaryp(ir,ic,il)+     
     $                     psibarxp(ir,ic,il)*phixp+
     $                     psibaryp(ir,ic,il)*phiyp)
     $                     /map2(ir,ic)

         coef_phi(ir,ic,il)=(cor(ir,ic)+relvorbar(ir,ic,il))/
     $                      map2(ir,ic)
                   
2100   Continue

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

C SOR solution for phi (3-D)
       Do 2250 iter=1,max_iters
         noconverg=0

         Do 2320 il=1,maxnl-1    
         Do 2310 ir=2,nr-1
         Do 2300 ic=2,nc-1 

           If (il.Eq.1) Then
             res=phi(ir,ic+1,il)+phi(ir,ic-1,il)+
     $           phi(ir-1,ic,il)+phi(ir+1,ic,il)-
     $           4.*phi(ir,ic,il)+
     $           coef_phi(ir,ic,il)*dssqr*     
     $           (pplp1(il)*phi(ir,ic,il+1)+
     $           thetaB(ir,ic)/dexnf(il)-
     $           pplp1(il)*phi(ir,ic,il))-          
     $           dssqr*for_phi(ir,ic,il)                   
             resmod=over_phi*res/
     $              (4.+coef_phi(ir,ic,il)*dssqr*pplp1(il))  
             phi(ir,ic,il)=phi(ir,ic,il)+resmod
           Else If (il.Eq.maxnl-1) Then
             res=phi(ir,ic+1,il)+phi(ir,ic-1,il)+
     $           phi(ir-1,ic,il)+phi(ir+1,ic,il)-
     $           4.*phi(ir,ic,il)+
     $           coef_phi(ir,ic,il)*dssqr*     
     $           (pplm1(il)*phi(ir,ic,il-1)-
     $           thetaT(ir,ic)/dexnf(il)-
     $           pplm1(il)*phi(ir,ic,il))-          
     $           dssqr*for_phi(ir,ic,il)
             resmod=over_phi*res/
     $              (4.+coef_phi(ir,ic,il)*dssqr*pplm1(il))  
             phi(ir,ic,il)=phi(ir,ic,il)+resmod       
           Else
             res=phi(ir,ic+1,il)+phi(ir,ic-1,il)+
     $           phi(ir-1,ic,il)+phi(ir+1,ic,il)-
     $           4.*phi(ir,ic,il)+
     $           coef_phi(ir,ic,il)*dssqr*
     $           (pplp1(il)*phi(ir,ic,il+1)+
     $           pplm1(il)*phi(ir,ic,il-1)-
     $           ppl(il)*phi(ir,ic,il))-
     $           dssqr*for_phi(ir,ic,il)
             resmod=over_phi*res/
     $              (4.+coef_phi(ir,ic,il)*dssqr*ppl(il)) 
             phi(ir,ic,il)=phi(ir,ic,il)+resmod               
           Endif
         
           If (Abs(resmod).Gt.thrs_phi) noconverg=noconverg+1              
                              
2300     Continue
2310     Continue
2320     Continue 

         If (noconverg.Eq.0) Goto 2350
      
2250   Continue       

       print*,'Convergence for PHI not reached'
       iter_acum_phi=iter_acum_phi+iter
       Goto 2375
 
2350   Continue 
       print*,'PHI converged at iter=',iter
       iter_acum_phi=iter_acum_phi+iter       

2375   Continue

C ------------------------------------------------------------------------------              
      
C Top and bottom BCs and underrelaxation
       Do 2400 ir=1,nr
       Do 2400 ic=1,nc
         phi(ir,ic,0)=phi(ir,ic,1)+dexnfh(0)*
     $                thetaB(ir,ic)
         phi(ir,ic,maxnl)=phi(ir,ic,maxnl-1)-dexnfh(maxnl-1)*
     $                    thetaT(ir,ic)
       Do 2410 il=0,maxnl
         phi(ir,ic,il)=under*phi(ir,ic,il)+(1.-under)*phi0(ir,ic,il)
2410   Continue       
2400   Continue              

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

C SOLVE STREAMFUNCTION psi (m2 s-1)

C ------------------------------------------------------------------------------
     
C Forcing function for_psi     
       Do 1100 ir=2,nr-1
       Do 1100 ic=2,nc-1
       Do 1100 il=1,maxnl-1

         phixx=map2(ir,ic)*(phi(ir,ic+1,il)+phi(ir,ic-1,il)-
     $         2.*phi(ir,ic,il))/dssqr
         phiyy=map2(ir,ic)*(phi(ir-1,ic,il)+phi(ir+1,ic,il)-
     $         2.*phi(ir,ic,il))/dssqr

         phipp=pplp1(il)*phi(ir,ic,il+1)+
     $         pplm1(il)*phi(ir,ic,il-1)-
     $         ppl(il)*phi(ir,ic,il)

         phixp1=0.5*(phi(ir,ic+1,il+1)/dexnfh(il)-
     $          phi(ir,ic+1,il-1)/dexnfh(il-1)+phi(ir,ic+1,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))    
         phixp2=0.5*(phi(ir,ic-1,il+1)/dexnfh(il)-
     $          phi(ir,ic-1,il-1)/dexnfh(il-1)+phi(ir,ic-1,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))        
         phixp=map(ir,ic)*(phixp1-phixp2)/ds2 
        
         phiyp1=0.5*(phi(ir-1,ic,il+1)/dexnfh(il)-
     $          phi(ir-1,ic,il-1)/dexnfh(il-1)+phi(ir-1,ic,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))    
         phiyp2=0.5*(phi(ir+1,ic,il+1)/dexnfh(il)-
     $          phi(ir+1,ic,il-1)/dexnfh(il-1)+phi(ir+1,ic,il)*
     $          (-1./dexnfh(il)+1./dexnfh(il-1)))               
         phiyp=map(ir,ic)*(phiyp1-phiyp2)/ds2     
          
         psix=map(ir,ic)*(psi0(ir,ic+1,il)-psi0(ir,ic-1,il))/                   
     $        ds2 
         psiy=map(ir,ic)*(psi0(ir-1,ic,il)-psi0(ir+1,ic,il))/                   
     $        ds2 
         psixx=map2(ir,ic)*(psi0(ir,ic+1,il)+psi0(ir,ic-1,il)-
     $         2.*psi0(ir,ic,il))/dssqr   
         psiyy=map2(ir,ic)*(psi0(ir-1,ic,il)+psi0(ir+1,ic,il)-
     $         2.*psi0(ir,ic,il))/dssqr
         psixy=map2(ir,ic)*(psi0(ir-1,ic+1,il)-psi0(ir-1,ic-1,il)-
     $         psi0(ir+1,ic+1,il)+psi0(ir+1,ic-1,il))/
     $         (ds2*ds2)       
         psixp=map(ir,ic)*(psi0(ir,ic+1,il+1)-psi0(ir,ic-1,il+1)-
     $         psi0(ir,ic+1,il-1)+psi0(ir,ic-1,il-1))/
     $         (ds2*(exnfh(il+1)-exnfh(il-1)))      
         psiyp=map(ir,ic)*(psi0(ir-1,ic,il+1)-psi0(ir+1,ic,il+1)-
     $         psi0(ir-1,ic,il-1)+psi0(ir+1,ic,il-1))/
     $         (ds2*(exnfh(il+1)-exnfh(il-1)))          

         for_psi(ir,ic,il)=(q(ir,ic,il)+
     $                     phixx+phiyy-
     $                     (cor(ir,ic)+relvorbar(ir,ic,il))*phipp-
     $                     corx(ir,ic)*psix-cory(ir,ic)*psiy-
     $                     2.*psibarxx(ir,ic,il)*psiyy-
     $                     2.*psibaryy(ir,ic,il)*psixx+
     $                     4.*psibarxy(ir,ic,il)*psixy+
     $                     psixp*phibarxp(ir,ic,il)+
     $                     psiyp*phibaryp(ir,ic,il)+
     $                     psibarxp(ir,ic,il)*phixp+
     $                     psibaryp(ir,ic,il)*phiyp)/
     $                     (map2(ir,ic)*(cor(ir,ic)+
     $                     phibarpp(ir,ic,il))) 
 
1100   Continue

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

C SOR solution for psi (for each level)
       Do 1200 il=1,maxnl-1

         Do 1250 iter=1,max_iters
           noconverg=0

           Do 1300 ir=2,nr-1
           Do 1300 ic=2,nc-1
             res=psi(ir,ic+1,il)+psi(ir,ic-1,il)+
     $           psi(ir-1,ic,il)+psi(ir+1,ic,il)-
     $           4.*psi(ir,ic,il)-
     $           dssqr*for_psi(ir,ic,il)
             resmod=.25*over_psi*res  
             psi(ir,ic,il)=psi(ir,ic,il)+resmod

             If (Abs(resmod).Gt.thrs_psi) noconverg=noconverg+1            

1300       Continue 

           If (noconverg.Eq.0) Goto 1350
      
1250     Continue       

         print*,'Convergence for PSI(',il,') not reached'
         iter_acum_psi=iter_acum_psi+iter       
         Goto 1200
 
1350     Continue 
         print*,'PSI(',il,') converged at iter=',iter
         iter_acum_psi=iter_acum_psi+iter
      
1200   Continue      

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

C When system convergence reached (or minimum), top and bottom BCs for psi and exit
       acum_f=(iter_acum_phi+iter_acum_psi/Real(maxnl-1))/2.
       If (acum_f.Ge.acum_o.And.icycle.Gt.5) Then   
         Do 1375 ir=1,nr
         Do 1375 ic=1,nc
           psi(ir,ic,0)=psi(ir,ic,1)+dexnfh(0)*
     $                  thetaB(ir,ic)
           psi(ir,ic,maxnl)=psi(ir,ic,maxnl-1)-dexnfh(maxnl-1)*
     $                      thetaT(ir,ic)
1375     Continue        
         print* 
         print*,'SYSTEM CONVERGED'
         print*
         print*
         Goto 1023
       Endif

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

C Top and bottom BCs, underrelaxation and new guess psi0=psi, phi0=phi
       Do 1400 ir=1,nr
       Do 1400 ic=1,nc
         psi(ir,ic,0)=psi(ir,ic,1)+dexnfh(0)*
     $                thetaB(ir,ic)
         psi(ir,ic,maxnl)=psi(ir,ic,maxnl-1)-dexnfh(maxnl-1)*
     $                    thetaT(ir,ic)
       Do 1410 il=0,maxnl
         psi(ir,ic,il)=under*psi(ir,ic,il)+(1.-under)*psi0(ir,ic,il)
         psi0(ir,ic,il)=psi(ir,ic,il)
	 phi0(ir,ic,il)=phi(ir,ic,il)
1410   Continue
1400   Continue              

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

C END CYCLE FOR PV INVERSION
          
1000  CONTINUE      

      print*
      print*,'SYSTEM CONVERGENCE NOT REACHED'
      print*
      print*

1023  Continue

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

C OUTPUT OF COMPONENT 1 OF PERTURBATION STREAMFUNCTION I_Sp1 (=[Fndim/grav]*psi, 
C in [grav/Fndim] m2 s-1)

      print*,
     $'Output of COMPONENT 1 OF PERTURBATION streamfunction ',
     $'I_Sp1 ([grav/Fndim] m2 s-1)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
      DO 9999 il=1,maxnl
        outgrid(ir,ic,il)=Sndim*PI2ndim*
     $                    ( psi0(ir,ic,il)-(exnfh(il)-exnf(il))*
     $                    (psi0(ir,ic,il)-psi0(ir,ic,il-1))/
     $                    dexnfh(il-1) )/grav
9999  CONTINUE     

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

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