!
!  problem_solve.f90
!
       subroutine problem_solve(rz,iz,dz,cz,rzp,izp,czp,ierr)      
!
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*),ierr
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)       
!
! Declarations
       integer :: nzone,ninter,ncoef,ncell,ntask,p,i,j,k,n,maxcloop
       integer :: length,alength,blength,anlength,mflength,ntmp
       integer :: nloop,icvs,jcvs,kcvs,ind1,ind2,icvp,jcvp,kcvp
       integer :: nint,icvp1,jcvp1,kcvp1,icvp2,jcvp2,kcvp2
       integer :: np,npmax,ncs,ncs1,ncs2
       real :: dtime,simtime,time,rcval,tavg
       character (len=80) :: simtype,zone,task,cell,cell1,cell2,cell1p,cell2p
       character (len=80) :: cval,cvalp,ctask,intf,intp,cprint
       character (len=4)  :: ctype   
       character (len=2), dimension(1:10) :: cnum
       data cnum/'1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ','10'/
!
! Pointers       
       integer :: pAP,pAE,pAW,pAN,pAS,pAT,pAB,pPhi,pB,pR,pPhic,pPhi0
       integer :: pRHO0,pTC0,pCP0,pVol,pAnorm,pX,pY,pZ,pMf,pWRK
       integer :: pnl,pnlx,pnm,pnmx,pnn,pnnx
       integer :: pzn,pcl,ptsk,pXC,pYC,pZC
       integer :: aAP,aAE,aAW,aAN,aAS,aAT,aAB,aB
       integer :: aAP1,aAE1,aAW1,aAN1,aAS1,aAT1,aAB1,aB1
       integer :: aAP2,aAE2,aAW2,aAN2,aAS2,aAT2,aAB2,aB2              
       integer :: pAcoef1,pBcoef1,pAcoef2,pBcoef2,pcind,pind1r,pind2r, &
                  pind1,pind12,pind2,pind21,pint,PIan1,PIan2,PImf1,PImf2,pcbc,pbc
!
! Steady or unsteady
       cval = '/simulation/type'
       call get_char_value('physics',cval,simtype,czp,ierr)
!
! If unsteady get transient settings
       simtime = 1E6; dtime=0.0; 
       if ( simtype(1:8) == 'unsteady' ) then
         cval = '/simulation/time_step'
         call get_real_value('physics',cval,dtime,rzp,ierr)
         cval = '/simulation/total_run_time' 
         call get_real_value('physics',cval,simtime,rzp,ierr)
         cval = '/simulation/max_cloops' 
         call get_int_value('physics',cval,maxcloop,izp,ierr)
       endif
!
! Zone list       
       cval = '/zone_list'
       call get_mem('physics',cval,ctype,pzn,nzone,ierr)
!
! Transient level
       do while (time <= simtime) 
!
! Initialize time level
! 
! 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)         
!
         do j =1,ncell
!
           cell =  '/'//trim(zone)//'/'//czp(pcl + j -1)
           call get_mem('coef',trim(cell)//'/energy0',ctype,pPhi0,alength,ierr)
           call get_mem('coef',trim(cell)//'/energy',ctype,pPhi,alength,ierr)
           if (time < 0.9*dtime) then
             cval = trim(cell)//'/energy0/ic'
             call get_real_value('physics',cval,rcval,rzp,ierr)
             call set_a_r(rz(pPhi0),rcval,alength) 
             call set_a_r(rz(pPhi),rcval,alength)
           else
             call set_a_b(rz(pPhi0),rz(pPhi),alength)
           endif
!
           if (time < 0.9*dtime) then           
             call get_mem('coef',trim(cell)//'/props/tc0',ctype,pTC0,alength,ierr)       
             call get_mem('coef',trim(cell)//'/props/cp0',ctype,pCP0,alength,ierr)
             call get_mem('coef',trim(cell)//'/props/rho0',ctype,pRHO0,alength,ierr)
!           
             cval = trim(cell)//'/props/tc0/ic'
             call get_real_value('physics',cval,rcval,rzp,ierr)
             call set_a_r(rz(pTC0),rcval,alength)
!
             cval = trim(cell)//'/props/cp0/ic'
             call get_real_value('physics',cval,rcval,rzp,ierr)
             call set_a_r(rz(pCP0),rcval,alength)           
!           
             cval = trim(cell)//'/props/rho0/ic'
             call get_real_value('physics',cval,rcval,rzp,ierr)
             call set_a_r(rz(pRHO0),rcval,alength)
           endif
!           
         enddo
       enddo 
!
! Inner loop (coefficient advancement)
       nloop = 1  
       do while (nloop <= maxcloop) 
!
! Fluid flow (not implemented)
!
! Heat transfer
!
! Generate internal fine mesh coefficients                 
!
! 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) 
!
         do j =1,ncell
!
           cell = '/'//trim(zone)//'/'//czp(pcl + j -1) 
!
! Cell dimensions      
           call get_int_value('physics',trim(cell)//'/icvs',icvs,izp,ierr)       
           call get_int_value('physics',trim(cell)//'/jcvs',jcvs,izp,ierr)      
           call get_int_value('physics',trim(cell)//'/kcvs',kcvs,izp,ierr)
!
! Coefficient space
           call get_mem('coef',trim(cell)//'/ap',ctype,pAP,alength,ierr)
           call get_mem('coef',trim(cell)//'/ae',ctype,pAE,alength,ierr)
           call get_mem('coef',trim(cell)//'/aw',ctype,pAW,alength,ierr)
           call get_mem('coef',trim(cell)//'/an',ctype,pAN,alength,ierr)
           call get_mem('coef',trim(cell)//'/as',ctype,pAS,alength,ierr)       
           call get_mem('coef',trim(cell)//'/at',ctype,pAT,alength,ierr)              
           call get_mem('coef',trim(cell)//'/ab',ctype,pAB,alength,ierr)
           call get_mem('coef',trim(cell)//'/b',ctype,pB,blength,ierr)
           call get_mem('coef',trim(cell)//'/r',ctype,pR,blength,ierr)         
           call get_mem('coef',trim(cell)//'/rc',ctype,pWRK,blength,ierr)  
!
! Geometry
           call get_mem('cell',trim(cell)//'/mesh/vol',ctype,pVol,alength,ierr)  
           call get_mem('cell',trim(cell)//'/mesh/anorm',ctype,pAnorm,anlength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/x',ctype,pX,anlength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/y',ctype,pY,anlength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/z',ctype,pZ,anlength,ierr)                      
           call get_mem('cell',trim(cell)//'/mesh/xc',ctype,pXC,alength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/yc',ctype,pYC,alength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/zc',ctype,pZC,alength,ierr)
!
! Dependent variable and properties
           call get_mem('coef',trim(cell)//'/energy0',ctype,pPhi0,alength,ierr)
           call get_mem('coef',trim(cell)//'/props/tc0',ctype,pTC0,alength,ierr)           
           call get_mem('coef',trim(cell)//'/props/cp0',ctype,pCP0,alength,ierr)          
           call get_mem('coef',trim(cell)//'/props/rho0',ctype,pRHO0,alength,ierr)
           call get_mem('coef',trim(cell)//'/energy',ctype,pPhi,alength,ierr)           
           call get_mem('coef',trim(cell)//'/energyc',ctype,pPhic,alength,ierr)
           call get_mem('coef',trim(cell)//'/mass/massfl',ctype,pMf,mflength,ierr)  

!
! Fine mesh internal coefficients
           call coefh(cell,rz(pAP),rz(pAE),rz(pAW),rz(pAN),rz(pAS),rz(pAT),rz(pAB), & 
                      rz(pB),rz(pPhi0),rz(pPhi),rz(pR),rz(pVol),rz(pAnorm), &
                      rz(pXC),rz(pYC),rz(pZC), &
                      rz(pMf),rz(pTC0),rz(pCP0),rz(pRHO0),icvs,jcvs,kcvs, &                        
                      rz,iz,dz,cz,rzp,izp,czp,ierr)
         enddo
       enddo
!
! Generate coefficients at interfaces that are boundary conditions.
!
! Interface list       
       call get_mem('physics','/interface_list',ctype,pint,nint,ierr)
!
       do i=1,nint
         intf = '/'//czp(pint + i -1) 

         call get_mem('physics',trim(intf)//'/cell_index',ctype,pcind,ntmp,ierr) 
         cell1=czp(pcind)
         cell2=czp(pcind+1)
         if (trim(cell1)=='boundary') then     
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pind2r,ntmp,ierr) 
           pind1r = pind2r
         elseif (trim(cell2)=='boundary') then
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pind1r,ntmp,ierr)         
           pind2r = pind1r
         else
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pind1r,ntmp,ierr)         
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pind2r,ntmp,ierr)
         endif
         call get_mem('interface',trim(intf)//'/index1',ctype,pind1,ncs1,ierr) 
         call get_mem('interface',trim(intf)//'/index2',ctype,pind2,ncs2,ierr)
! Only structured cells enabled
         ncs1 = ncs1/4; ncs2=ncs2/4
         call get_mem('interface',trim(intf)//'/acoef1',ctype,pAcoef1,ntmp,ierr)          
         call get_mem('interface',trim(intf)//'/bcoef1',ctype,pBcoef1,ntmp,ierr)
         call get_mem('interface',trim(intf)//'/anorm1',ctype,pIan1,ntmp,ierr)
         call get_mem('interface',trim(intf)//'/massfl1',ctype,pImf1,ntmp,ierr)
         call get_mem('interface',trim(intf)//'/acoef2',ctype,pAcoef2,ntmp,ierr)          
         call get_mem('interface',trim(intf)//'/bcoef2',ctype,pBcoef2,ntmp,ierr) 
         call get_mem('interface',trim(intf)//'/anorm2',ctype,pIan2,ntmp,ierr)
         call get_mem('interface',trim(intf)//'/massfl2',ctype,pImf2,ntmp,ierr)
         call get_mem('physics',trim(intf)//'/cbc',ctype,pcbc,ntmp,ierr)
         call get_mem('physics',trim(intf)//'/bc',ctype,pbc,ntmp,ierr)         
!     
         call coefh_intf(cell1,cell2,izp(pind1r),izp(pind2r),iz(pind1),iz(pind2),& 
                         rz(pAcoef1),rz(pBcoef1),rz(pIan1),rz(pImf1),rz(pAcoef2),rz(pBcoef2),&
                         rz(pIan2),rz(pImf2),czp(pcbc),rzp(pbc),ncs1,ncs2,&
                         rz,iz,dz,cz,rzp,izp,czp,ierr)
!          
       enddo 
!
! Coarse mesh internal coefficients
!
! 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) 
!
         do j =1,ncell
!
           cell = '/'//trim(zone)//'/'//czp(pcl + j -1) 
!
! Initialize mg level information arrays
           call get_mem('physics',trim(cell)//'/nl',ctype,pnl,npmax,ierr)
           call get_mem('physics',trim(cell)//'/nlx',ctype,pnlx,npmax,ierr)        
           call get_mem('physics',trim(cell)//'/nm',ctype,pnm,npmax,ierr)
           call get_mem('physics',trim(cell)//'/nmx',ctype,pnmx,npmax,ierr) 
           call get_mem('physics',trim(cell)//'/nn',ctype,pnn,npmax,ierr)
           call get_mem('physics',trim(cell)//'/nnx',ctype,pnnx,npmax,ierr)	   
!
           call get_int_value('physics','/simulation/npcells',np,izp,ierr)
           do p=1,np
             cvalp = trim(cell)//'/level_'//cnum(p) 
             call get_mem('coef',trim(cvalp)//'/ap',ctype,aAP,length,ierr) 
             call get_mem('coef',trim(cvalp)//'/ae',ctype,aAE,length,ierr)
             call get_mem('coef',trim(cvalp)//'/aw',ctype,aAW,length,ierr)
             call get_mem('coef',trim(cvalp)//'/an',ctype,aAN,length,ierr)
             call get_mem('coef',trim(cvalp)//'/as',ctype,aAS,length,ierr)
             call get_mem('coef',trim(cvalp)//'/at',ctype,aAT,length,ierr)
             call get_mem('coef',trim(cvalp)//'/ab',ctype,aAB,length,ierr)
             call get_mem('coef',trim(cvalp)//'/b',ctype,aB,length,ierr)                                                                                          
!
             call get_int_value('coef',trim(cvalp)//'/icvs',icvp,iz,ierr)
             call get_int_value('coef',trim(cvalp)//'/jcvs',jcvp,iz,ierr)
             call get_int_value('coef',trim(cvalp)//'/kcvs',kcvp,iz,ierr)               
!     
             call coefh_mg(cval,p,rz(pAP),rz(pAE),rz(pAW),rz(pAN),rz(pAS),rz(pAT),rz(pAB),rz(pB), &
                rz(aAP),rz(aAE),rz(aAW),rz(aAN),rz(aAS),rz(aAT),rz(aAB),rz(aB), &                                      
                izp(pnl),izp(pnlx),izp(pnm),izp(pnmx),izp(pnn),izp(pnnx), &
                icvs,jcvs,kcvs,icvp,jcvp,kcvp, &
                rz,iz,dz,cz,rzp,izp,czp,ierr)     
!
           enddo
         enddo
       enddo 
!
! Loop over interfaces at fine and mg levels - move all cell boundary link coefs
! to the interface data areas and zero them
!
       call get_int_value('physics','/simulation/npcells',np,izp,ierr)
       do n=1,nint
         intf = '/'//czp(pint + n -1) 
         call get_mem('physics',trim(intf)//'/cell_index',ctype,pcind,ntmp,ierr)
         cell1=czp(pcind)
         cell2=czp(pcind+1)
!
! Work on the fine mesh
         if (trim(cell1)=='boundary') then     
!        
           call get_int_value('physics',trim(cell2)//'/icvs',icvp2,izp,ierr)       
           call get_int_value('physics',trim(cell2)//'/jcvs',jcvp2,izp,ierr)      
           call get_int_value('physics',trim(cell2)//'/kcvs',kcvp2,izp,ierr)
           icvp1=1;jcvp1=1;kcvp1=1
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pind2r,ntmp,ierr) 
           pind1r = pind2r
!
         elseif (trim(cell2)=='boundary') then
!
           call get_int_value('physics',trim(cell1)//'/icvs',icvp1,izp,ierr)  
           call get_int_value('physics',trim(cell1)//'/jcvs',jcvp1,izp,ierr)      
           call get_int_value('physics',trim(cell1)//'/kcvs',kcvp1,izp,ierr)
           icvp2=1;jcvp2=1;kcvp2=1           
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pind1r,ntmp,ierr)         
           pind2r = pind1r
!
         else
!
           call get_int_value('physics',trim(cell1)//'/icvs',icvp1,izp,ierr)  
           call get_int_value('physics',trim(cell1)//'/jcvs',jcvp1,izp,ierr)      
           call get_int_value('physics',trim(cell1)//'/kcvs',kcvp1,izp,ierr)
           call get_int_value('physics',trim(cell2)//'/icvs',icvp2,izp,ierr)       
           call get_int_value('physics',trim(cell2)//'/jcvs',jcvp2,izp,ierr)      
           call get_int_value('physics',trim(cell2)//'/kcvs',kcvp2,izp,ierr)
!           
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pind1r,ntmp,ierr)         
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pind2r,ntmp,ierr)
           call get_mem('interface',trim(intf)//'/index12','int',pind12,ntmp,ierr)
           call get_mem('interface',trim(intf)//'/index21','int',pind21,ntmp,ierr)
!
         endif
         call get_mem('interface',trim(intf)//'/index1',ctype,pind1,ncs1,ierr) 
         call get_mem('interface',trim(intf)//'/index2',ctype,pind2,ncs2,ierr)
! Only structured cells enabled
         ncs1 = ncs1/4; ncs2=ncs2/4
         call get_mem('interface',trim(intf)//'/acoef1',ctype,pAcoef1,ntmp,ierr)          
         call get_mem('interface',trim(intf)//'/bcoef1',ctype,pBcoef1,ntmp,ierr)
         call get_mem('interface',trim(intf)//'/acoef2',ctype,pAcoef2,ntmp,ierr)          
         call get_mem('interface',trim(intf)//'/bcoef2',ctype,pBcoef2,ntmp,ierr) 
!
! Coefficient space
         aAP1=1;aAE1=1;aAW1=1;aAN1=1;aAS1=1;aAT1=1;aAB1=1;aB1=1
         if (trim(cell1)/='boundary') then    
           call get_mem('coef',trim(cell1)//'/ap',ctype,aAP1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/ae',ctype,aAE1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/aw',ctype,aAW1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/an',ctype,aAN1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/as',ctype,aAS1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/at',ctype,aAT1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/ab',ctype,aAB1,alength,ierr)
           call get_mem('coef',trim(cell1)//'/b',ctype,aB1,blength,ierr) 
         endif 
!
         aAP2=1;aAE2=1;aAW2=1;aAN2=1;aAS2=1;aAT2=1;aAB2=1;aB2=1
         if (trim(cell2)/='boundary') then 
           call get_mem('coef',trim(cell2)//'/ap',ctype,aAP2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/ae',ctype,aAE2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/aw',ctype,aAW2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/an',ctype,aAN2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/as',ctype,aAS2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/at',ctype,aAT2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/ab',ctype,aAB2,alength,ierr)
           call get_mem('coef',trim(cell2)//'/b',ctype,aB2,blength,ierr) 
         endif
!
         call coefh_mg_zero(0,izp(pind1r),izp(pind2r),iz(pind1),iz(pind2),& 
                         rz(pAcoef1),rz(pBcoef1),rz(pAcoef2),rz(pBcoef2),&
                         rz(aAP1),rz(aAE1),rz(aAW1),rz(aAN1),rz(aAS1),rz(aAT1),rz(aAB1),rz(aB1),&
                         rz(aAP2),rz(aAE2),rz(aAW2),rz(aAN2),rz(aAS2),rz(aAT2),rz(aAB2),rz(aB2),&
                         icvp1,jcvp1,kcvp1,icvp2,jcvp2,kcvp2,ncs1,ncs2,rz,iz,dz,cz,rzp,izp,czp,ierr)
!
! Now work on the mg coarse mesh levels
         do p=1,np
           intp  = trim(intf)//'/level_'//cnum(p) 
!
           if (trim(cell1) .ne. 'boundary') then 
             cell1p= trim(cell1)//'/level_'//cnum(p)
           else 
             cell1p= cell1
           endif  
           if (trim(cell2) .ne. 'boundary') then
             cell2p= trim(cell2)//'/level_'//cnum(p)
           else
             cell2p= cell2
           endif
! 
           if (trim(cell1)=='boundary') then
!
             call get_int_value('coef',trim(cell2p)//'/icvs',icvp2,iz,ierr)       
             call get_int_value('coef',trim(cell2p)//'/jcvs',jcvp2,iz,ierr)      
             call get_int_value('coef',trim(cell2p)//'/kcvs',kcvp2,iz,ierr) 
             icvp1=1;jcvp1=1;kcvp1=1
             call get_mem('interface',trim(intp)//'/index2_range',ctype,pind2r,ntmp,ierr) 
             pind1r = pind2r
!
           elseif (trim(cell2)=='boundary') then
!           
             call get_int_value('coef',trim(cell1p)//'/icvs',icvp1,iz,ierr)  
             call get_int_value('coef',trim(cell1p)//'/jcvs',jcvp1,iz,ierr)      
             call get_int_value('coef',trim(cell1p)//'/kcvs',kcvp1,iz,ierr)
             icvp2=1;jcvp2=1;kcvp2=1 
             call get_mem('interface',trim(intp)//'/index1_range',ctype,pind1r,ntmp,ierr)         
             pind2r = pind1r
!
           else
!
             call get_int_value('coef',trim(cell1p)//'/icvs',icvp1,iz,ierr)  
             call get_int_value('coef',trim(cell1p)//'/jcvs',jcvp1,iz,ierr)      
             call get_int_value('coef',trim(cell1p)//'/kcvs',kcvp1,iz,ierr)      
             call get_int_value('coef',trim(cell2p)//'/icvs',icvp2,iz,ierr)       
             call get_int_value('coef',trim(cell2p)//'/jcvs',jcvp2,iz,ierr)      
             call get_int_value('coef',trim(cell2p)//'/kcvs',kcvp2,iz,ierr)           
!
             call get_mem('interface',trim(intp)//'/index1_range',ctype,pind1r,ntmp,ierr)         
             call get_mem('interface',trim(intp)//'/index2_range',ctype,pind2r,ntmp,ierr)
             call get_mem('interface',trim(intp)//'/index12','int',pind12,ntmp,ierr)
             call get_mem('interface',trim(intp)//'/index21','int',pind21,ntmp,ierr)             
!
           endif
           call get_mem('interface',trim(intp)//'/index1',ctype,pind1,ncs1,ierr) 
           call get_mem('interface',trim(intp)//'/index2',ctype,pind2,ncs2,ierr)
! Only structured cells enabled
           ncs1 = ncs1/4; ncs2=ncs2/4
           call get_mem('interface',trim(intp)//'/acoef1',ctype,pAcoef1,ntmp,ierr)          
           call get_mem('interface',trim(intp)//'/bcoef1',ctype,pBcoef1,ntmp,ierr)
           call get_mem('interface',trim(intp)//'/acoef2',ctype,pAcoef2,ntmp,ierr)          
           call get_mem('interface',trim(intp)//'/bcoef2',ctype,pBcoef2,ntmp,ierr)          
!
! Coefficient space
           aAP1=1;aAE1=1;aAW1=1;aAN1=1;aAS1=1;aAT1=1;aAB1=1;aB1=1
           if (trim(cell1)/='boundary') then 
             call get_mem('coef',trim(cell1p)//'/ap',ctype,aAP1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/ae',ctype,aAE1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/aw',ctype,aAW1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/an',ctype,aAN1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/as',ctype,aAS1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/at',ctype,aAT1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/ab',ctype,aAB1,alength,ierr)
             call get_mem('coef',trim(cell1p)//'/b',ctype,aB1,blength,ierr) 
           endif
!
           aAP2=1;aAE2=1;aAW2=1;aAN2=1;aAS2=1;aAT2=1;aAB2=1;aB2=1
           if (trim(cell2)/='boundary') then
             call get_mem('coef',trim(cell2p)//'/ap',ctype,aAP2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/ae',ctype,aAE2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/aw',ctype,aAW2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/an',ctype,aAN2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/as',ctype,aAS2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/at',ctype,aAT2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/ab',ctype,aAB2,alength,ierr)
             call get_mem('coef',trim(cell2p)//'/b',ctype,aB2,blength,ierr)
           endif
!
           call coefh_mg_zero(p,iz(pind1r),iz(pind2r),iz(pind1),iz(pind2),& 
                         rz(pAcoef1),rz(pBcoef1),rz(pAcoef2),rz(pBcoef2),&
                         rz(aAP1),rz(aAE1),rz(aAW1),rz(aAN1),rz(aAS1),rz(aAT1),rz(aAB1),rz(aB1),&
                         rz(aAP2),rz(aAE2),rz(aAW2),rz(aAN2),rz(aAS2),rz(aAT2),rz(aAB2),rz(aB2),&
                         icvp1,jcvp1,kcvp1,icvp2,jcvp2,kcvp2,ncs1,ncs2,rz,iz,dz,cz,rzp,izp,czp,ierr)     
!          
         enddo
       enddo
!       
! Solve mg solution across all zones/cells/interfaces
!       call mgsolve(rz,iz,dz,cz,rzp,izp,czp,ierr)         
       call mgsolve_v(rz,iz,dz,cz,rzp,izp,czp,ierr) 
!
! Check cofficient level conergence
       nloop = nloop + 1
!
! Move to new time level
       enddo
!       
       time = time + dtime
!
! End of timestep postprocessing
       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)         
!
         do j =1,ncell
!
           cell =  '/'//trim(zone)//'/'//czp(pcl + j -1)
           call get_mem('coef',trim(cell)//'/energy',ctype,pPhi,alength,ierr) 
           call get_mem('coef',trim(cell)//'/props/cp0',ctype,pCP0,alength,ierr)
           call get_mem('coef',trim(cell)//'/rc',ctype,pWRK,alength,ierr)                      
!
           call set_a_b(rz(pWRK),rz(pPhi),alength)
           call set_adivb(rz(pWRK),rz(pCP0),alength)
           call cell_avg(rz(pWRK),tavg,alength)
!
           write(*,*)  trim(cell)
           write(*,*) 'time=',time,'temp_avg=',tavg
         enddo
       enddo
!
       enddo
!
! Final postprocessing
!
       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)         
!
         do j =1,ncell
!
           cprint = czp(pcl + j -1)
           cell =  '/'//trim(zone)//'/'//czp(pcl + j -1)

           call get_mem('coef',trim(cell)//'/energy',ctype,pPhi,alength,ierr) 
           call get_mem('coef',trim(cell)//'/props/cp0',ctype,pCP0,alength,ierr)
           call get_mem('coef',trim(cell)//'/rc',ctype,pWRK,alength,ierr)                      
!
           call set_a_b(rz(pWRK),rz(pPhi),alength)
           call set_adivb(rz(pWRK),rz(pCP0),alength)
! 
           call get_int_value('physics',trim(cell)//'/icvs',icvs,izp,ierr)       
           call get_int_value('physics',trim(cell)//'/jcvs',jcvs,izp,ierr)      
           call get_int_value('physics',trim(cell)//'/kcvs',kcvs,izp,ierr)
           call get_mem('cell',trim(cell)//'/mesh/xc',ctype,pXC,alength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/yc',ctype,pYC,alength,ierr)
           call get_mem('cell',trim(cell)//'/mesh/zc',ctype,pZC,alength,ierr)
!
           call write_plt3d(cprint,rz(pXC),rz(pYC),rz(pZC),rz(pWRK),icvs,jcvs,kcvs)
!
         enddo
       enddo      
       end subroutine problem_solve         
