!
!  mgsolve.f90
!
       subroutine mgsolve(rz,iz,dz,cz,rzp,izp,czp,ierr)
!
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*)
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)
       integer :: ierr
!
       character (len=80) :: cell,cval,zone 
       character (len=4) :: ctype
!
       integer :: i,j,k,p,ncell,ip,length,alength,blength,nzone,maxcell
       integer :: np,npmax
       integer, dimension(0:9) :: nit
       real :: rate  
       real, dimension(0:9) :: sumreso,sumresn,sumresoi
       logical :: lsolve,lstop
       character (len=2), dimension(1:10) :: cnum
       data cnum/'1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ','10'/
!
       real, allocatable, dimension(:,:,:) :: reso,resoi,resn       
       integer, allocatable, dimension(:) :: icell
       integer, allocatable, dimension(:,:,:) :: icvp,jcvp,kcvp       
!           
! Pointers        
       integer, allocatable, dimension(:,:,:) :: aAP,aAE,aAW, &
                                 aAN,aAS,aAT,aAB,aR,aRC,aPhic
       integer, allocatable, dimension(:,:) :: aPhi,aBB                         
       integer :: pzn,pcl
!
       nit=0; sumreso=0.0; sumresn=0.0; sumresoi=0.0
!
!------ Interior nodes 
!
! Determine number of zones and max cells in a zone
       call get_mem('physics','/zone_list',ctype,pzn,nzone,ierr)
       maxcell = 0
       do i =1,nzone
         zone = czp(pzn + i -1) 
         cval = '/'//trim(zone)//'/'//'cell_list' 
         call get_mem('physics',cval,ctype,pcl,ncell,ierr)       
         maxcell = max(ncell,maxcell)
       enddo 
!
! Allocate cell data sizes
       allocate (icell(1:nzone)) 
       allocate (icvp(1:nzone,1:maxcell,0:9))
       allocate (jcvp(1:nzone,1:maxcell,0:9))
       allocate (kcvp(1:nzone,1:maxcell,0:9))       
       icell=0; icvp=0; jcvp=0; kcvp=0
!
! Allocate residual arrays
       allocate (reso(1:nzone,1:maxcell,0:9))
       allocate (resoi(1:nzone,1:maxcell,0:9))
       allocate (resn(1:nzone,1:maxcell,0:9)) 
       reso=0.0; resoi=0.0; resn=0.0
!
! Allocate and initialize pointer space
       allocate(aAP(1:nzone,1:maxcell,0:9)) 
       allocate(aAE(1:nzone,1:maxcell,0:9))
       allocate(aAW(1:nzone,1:maxcell,0:9))
       allocate(aAN(1:nzone,1:maxcell,0:9))
       allocate(aAS(1:nzone,1:maxcell,0:9))
       allocate(aAT(1:nzone,1:maxcell,0:9))
       allocate(aAB(1:nzone,1:maxcell,0:9))
       allocate(aBB(1:nzone,1:maxcell))
       allocate(aR(1:nzone,1:maxcell,0:9))
       allocate(aRC(1:nzone,1:maxcell,0:9))                                                               
       allocate(aPhi(1:nzone,1:maxcell))
       allocate(aPhic(1:nzone,1:maxcell,0:9))
       aAP=0; aAE=0; aAW=0; aAN=0; aAS=0
       aAT=0; aBB=0; aR=0; aRC=0; aPhi=0; aPhic=0
! 
! Loop over zones
       do i =1,nzone
!
! Loop over cells in a zone
         zone = czp(pzn + i -1)        
         cval = '/'//trim(zone)//'/'//'cell_list' 
         call get_mem('physics',cval,ctype,pcl,ncell,ierr)         
         icell (i) = ncell  
!
         do j =1,ncell
           cell =  '/'//trim(zone)//'/'//czp(pcl + j -1)
!
! Cell dimensions      
           call get_int_value('physics',trim(cell)//'/icvs',icvp(i,j,0),izp,ierr)       
           call get_int_value('physics',trim(cell)//'/jcvs',jcvp(i,j,0),izp,ierr)      
           call get_int_value('physics',trim(cell)//'/kcvs',kcvp(i,j,0),izp,ierr)
!
! Dependent variable space           
           call get_mem('coef',trim(cell)//'/energy',ctype,aPhi(i,j),alength,ierr)           
           call get_mem('coef',trim(cell)//'/energyc',ctype,aPhic(i,j,0),alength,ierr)
!
! Coefficient space
           call get_mem('coef',trim(cell)//'/ap',ctype,aAP(i,j,0),alength,ierr)
           call get_mem('coef',trim(cell)//'/ae',ctype,aAE(i,j,0),alength,ierr)
           call get_mem('coef',trim(cell)//'/aw',ctype,aAW(i,j,0),alength,ierr)
           call get_mem('coef',trim(cell)//'/an',ctype,aAN(i,j,0),alength,ierr)
           call get_mem('coef',trim(cell)//'/as',ctype,aAS(i,j,0),alength,ierr)       
           call get_mem('coef',trim(cell)//'/at',ctype,aAT(i,j,0),alength,ierr)              
           call get_mem('coef',trim(cell)//'/ab',ctype,aAB(i,j,0),alength,ierr)
           call get_mem('coef',trim(cell)//'/b',ctype,aBB(i,j),blength,ierr)
           call get_mem('coef',trim(cell)//'/r',ctype,aR(i,j,0),blength,ierr)         
           call get_mem('coef',trim(cell)//'/rc',ctype,aRC(i,j,0),blength,ierr)
!
! Coarse mesh internal coefficients
           call get_int_value('physics',trim(cell)//'/np',np,izp,ierr)
           do p=1,np
             cval = trim(cell)//'/level_'//cnum(p) 
             call get_mem('coef',trim(cval)//'/ap',ctype,aAP(i,j,p),length,ierr) 
             call get_mem('coef',trim(cval)//'/ae',ctype,aAE(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/aw',ctype,aAW(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/an',ctype,aAN(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/as',ctype,aAS(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/at',ctype,aAT(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/ab',ctype,aAB(i,j,p),length,ierr)   
             call get_mem('coef',trim(cval)//'/r',ctype,aR(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/rc',ctype,aRC(i,j,p),length,ierr)
             call get_mem('coef',trim(cval)//'/phic',ctype,aPhic(i,j,p),length,ierr)
!
             call get_int_value('coef',trim(cval)//'/icvs',icvp(i,j,p),iz,ierr)
             call get_int_value('coef',trim(cval)//'/jcvs',jcvp(i,j,p),iz,ierr)
             call get_int_value('coef',trim(cval)//'/kcvs',kcvp(i,j,p),iz,ierr)  
           enddo
         enddo
       enddo
!
!------ Interface control surfaces 
!
! To be developed ....
!
!
! Multigrid solution
       lsolve = .true.; lstop = .false.
       p = 0; ip=0
!     
       do while (lsolve)         
!
! Fine mesh level
!
       if ( p==0 ) then
         if (nit(p)==0) then
!
! Sum over interfaces (not implemented)
!         
           sumresoi(p) = 0.0            
           do i=1,nzone
             do j=1,icell(i)
               call residlp(icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                            rz(aPhi(i,j)),rz(aAP(i,j,0)), &                           
                            rz(aAE(i,j,0)),rz(aAW(i,j,0)),rz(aAN(i,j,0)), &
                            rz(aAS(i,j,0)),rz(aAT(i,j,0)),rz(aAB(i,j,0)), &
                            rz(aBB(i,j)),rz(aR(i,j,0)),resoi(i,j,p)) 
               sumresoi(p) = sumresoi(p) + resoi(i,j,p)
               length = icvp(i,j,p)*jcvp(i,j,p)*kcvp(i,j,p)
               call set_a_0(rz(aPhic(i,j,p)),length)
             enddo             
           enddo  
!           write(*,*) 'reso0i=',sumresoi(p)
         endif 
!
         do i=1,nzone
           do j=1,icell(i)
             call solvegs (rz(aAP(i,j,0)), &                           
                           rz(aAE(i,j,0)),rz(aAW(i,j,0)),rz(aAN(i,j,0)), &
                           rz(aAS(i,j,0)),rz(aAT(i,j,0)),rz(aAB(i,j,0)), &
                           rz(aR(i,j,0)),rz(aPhic(i,j,0)), &
                           icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0))
           enddo             
         enddo 
!
! Sum over interfaces (not implemented)
!
         sumresn(p) = 0.0 
         do i=1,nzone
           do j=1,icell(i)
             call residlp(icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                          rz(aPhic(i,j,0)),rz(aAP(i,j,0)), &                           
                          rz(aAE(i,j,0)),rz(aAW(i,j,0)),rz(aAN(i,j,0)), &
                          rz(aAS(i,j,0)),rz(aAT(i,j,0)),rz(aAB(i,j,0)), &
                          rz(aR(i,j,0)),rz(aRC(i,j,0)), &                    
                          resn(i,j,p))             
             sumresn(p) = sumresn(p) + resn(i,j,p)
           enddo             
         enddo 
!
         if (nit(p) == 0) sumreso(p) = sumresn(p)
!
         rate = sumresn(p)/(sumreso(p)+1E-10)
!         if (nit(p) > 0) write(*,*) 'p=',p,nit(p),'rsn/reso=', &
!                            rate,'resn=',sumresn(p)
!
         if (sumresn(p) <= 0.001*sumresoi(p)) then
           do i=1,nzone
             do j=1,icell(i)                   
               call update_phi(p,icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                       rz(aPhi(i,j)),icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                       rz(aPhic(i,j,0)))
             enddo             
           enddo
           ip =  ip -1
         elseif (nit(p)>0 .and. rate>0.3 .and. np >0 &
                 .and. .not. lstop) then
           ip = ip + 1
           nit(p+1) = 0
           sumreso(p) = sumresn(p)               
         else
           if (nit(p)>200) then
             do i=1,nzone
               do j=1,icell(i)                   
                 call update_phi(p,icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                       rz(aPhi(i,j)),icvp(i,j,0),jcvp(i,j,0),kcvp(i,j,0), &
                       rz(aPhic(i,j,0)))
               enddo             
             enddo           
             ip = ip -1             
           else 
             sumreso(p) = sumresn(p)          
             nit(p) = nit(p) + 1              
           endif
         endif 
!
       elseif ( p>0 ) then
!
         if (nit(p)==0) then
!
! Update interface variables (not implemented)
!
           sumresoi(p) = 0.0  
           do i=1,nzone
             do j=1,icell(i)
               call residlsc(p,icvp(i,j,p-1),jcvp(i,j,p-1),kcvp(i,j,p-1), &
                           rz(aAP(i,j,p-1)),rz(aAE(i,j,p-1)),rz(aAW(i,j,p-1)),rz(aAN(i,j,p-1)), &
                           rz(aAS(i,j,p-1)),rz(aAT(i,j,p-1)),rz(aAB(i,j,p-1)),rz(aR(i,j,p-1)), &
                           rz(aPhic(i,j,p-1)),icvp(i,j,p),jcvp(i,j,p),kcvp(i,j,p),rz(aR(i,j,p)), & 
                           resoi(i,j,p))
               sumresoi(p) = sumresoi(p) + resoi(i,j,p)                                       
               length = icvp(i,j,p)*jcvp(i,j,p)*kcvp(i,j,p)
               call set_a_0(rz(aPhic(i,j,p)),length) 
             enddo
           enddo   
!           write(*,*) 'reso1i=',sumresoi(p)
         endif
!
         do i=1,nzone
           do j=1,icell(i)
             call solvegs(rz(aAP(i,j,p)),rz(aAE(i,j,p)),rz(aAW(i,j,p)),rz(aAN(i,j,p)),  &
                      rz(aAS(i,j,p)),rz(aAT(i,j,p)),rz(aAB(i,j,p)),rz(aR(i,j,p)),rz(aPhic(i,j,p)),   &
                      icvp(i,j,p),jcvp(i,j,p),kcvp(i,j,p))
           enddo
         enddo 
!
! Update interface variables (not implemented)
!
         sumresn(p) = 0.0 
         do i=1,nzone
           do j=1,icell(i)
             call residlp(icvp(i,j,p),jcvp(i,j,p),kcvp(i,j,p),rz(aPhic(i,j,p)),           & 
                      rz(aAP(i,j,p)),rz(aAE(i,j,p)),rz(aAW(i,j,p)),rz(aAN(i,j,p)),rz(aAS(i,j,p)),  &
                      rz(aAT(i,j,p)),rz(aAB(i,j,p)),rz(aR(i,j,p)),rz(aRC(i,j,p)),resn(i,j,p)) 
             sumresn(p) = sumresn(p) + resn(i,j,p)
           enddo 
         enddo

         if (nit(p) == 0) sumreso(p) = sumresn(p)
!
         rate = sumresn(p)/(sumreso(p)+1E-10)
!         if (nit(p) > 0) write(*,*) 'p=',p,nit(p),'rsn/reso=', &
!                               rate,'resn=',sumresn(p)
!
         if (sumresn(p) <= 0.001*sumresoi(p)) then
!
           do i=1,nzone
             do j=1,icell(i)  
               call update_phi(p,icvp(i,j,p-1),jcvp(i,j,p-1),kcvp(i,j,p-1),  &
                         rz(aPhic(i,j,p-1)),icvp(i,j,p),jcvp(i,j,p),kcvp(i,j,p), &
                         rz(aPhic(i,j,p)))             
             enddo
           enddo
           ip = ip - 1            
           nit(p-1) = nit(p-1) + 1 
         elseif (nit(p)>0 .and. rate>0.3) then
           if (np>p) then 
             ip = ip + 1
             nit(p+1) = 0
             sumreso(p) = sumresn(p)               
           else
             lstop = .true. 
             ip = 0
           endif 
         else 
           sumreso(p) = sumresn(p)           
           nit(p) = nit(p) + 1
           if (nit(p).gt.100) stop
         endif
!
       else 
         lsolve = .false. 
         write(*,*) 'p=',p,'multigrid converged'
         do i=1,np
           write(*,*) 'level=',i,'nit=',nit(i),sumresn(i),sumresoi(i)
         enddo
       endif
       p = ip 
!
       end do
!             
       end subroutine mgsolve       
