! $Id: streamfunction.f90 19193 2012-06-30 12:55:46Z wdobler $
!
!  This module solves for a streamfunction for the 2D velocity field.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of special_dummies.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lspecial = .true.
!
! MVAR CONTRIBUTION 1
! MAUX CONTRIBUTION 1 
! COMMUNICATED AUXILIARIES 0
!
!***************************************************************
module Special
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages, only: svn_id, fatal_error
  use Deriv
!
  implicit none
!
  include '../special.h'
!
!  Global to this method
!
  real, dimension(mx,my,mz) :: psi
!
!  Code constants
!
  real :: amplpsi=1d-5   !amplitude of initial perturbation for psi
  real :: tolerance=1d-6,tolerance_coarsest=impossible
  real :: alpha_sor=impossible
!
!  Variables for calculating the Rayleigh number. Default are Europa values.
!
  real ::  gravity_z = 1.3  ! gravity in m/s^2
  real ::  rho0_bq = 917.   ! density in Kg/m^3, Boussinesq value
  real ::  alpha_thermal = 1.65d-4  ! Thermal expansivity in K^-1
  real ::  kappa = 1d-6     ! Thermal diffusivity in m^2/s
  !real ::  cp = 2000.      ! Specific heat in J Kg^-1 K^-1
  real ::  eta_0 = 1d13     ! Viscosity at melting temperature in Pa s
  !real ::  dslab_km = 20.  ! Ice shell thikness in km
  real :: Tbot=270, Tupp=100
!
!  Variables for variable viscosity
!
  real ::  TT_melt = 270.   ! Melting tempertature in K
  real ::  Avisc = 4.00     ! constants for viscosity prescriptions
  real ::  Bvisc = 0.00
  real ::  Cvisc = 0.00
  real ::  deltaT1, Lz1
!
!  Variables for tidal heating
!
  real :: mu_ice = 4d9     ! Rigidity of ice in Pa s
  real :: epsi_0 = 2.1d-5  ! Amplitude of original tidal flexing
  real :: EuropaPeriod = 3.0682204d5 ! Period of Europa, in seconds  
  real :: mu1_ice, OmegaEuropa, kappa1
!
!  These are the needed internal "pencils".
!
  type InternalPencils
     real, dimension(nx,3)   :: uu
     real, dimension(nx) :: ugTT,u2,qtidal,eta,devsigzz
  endtype InternalPencils
!
  type (InternalPencils) :: q
  !$omp threadprivate(q)
!
  character (len=labellen) :: initpsi='nothing'
  character (len=labellen) :: iconv_viscosity='Newtonian'
  character (len=labellen) :: iorder_sor_solver='high_order'
!
  logical :: lprint_residual=.false.,ltidal_heating=.true.
  logical :: ltemperature_advection=.true.,ltemperature_diffusion=.true.
  logical :: lmultigrid=.true.
!  
  integer :: maxit=1000
!
  real :: Ra=impossible
  logical :: lsave_residual=.true.
  real :: kx_TT=2*pi, kz_TT=pi, ampltt=0.01
  logical :: lpoisson_test=.false.
  logical :: ldirect_solver=.true.
!
  integer :: npost=5,npre=5
  integer :: igamma=1,n_vcycles=1
  integer :: nx_coarsest=3
!
  logical :: lsave_residual_grid=.false.
  logical :: ladimensional=.true.
  logical :: lsplit_temperature=.false.
!
  namelist /special_init_pars/ amplpsi,alpha_sor,lprint_residual,&
       tolerance,maxit,gravity_z,rho0_bq,alpha_thermal,kappa,eta_0,&
       iconv_viscosity,Avisc,Bvisc,Cvisc,&
       Tbot,Tupp,Ra,iorder_sor_solver,lsave_residual,&
       kx_TT,kz_TT,ampltt,initpsi,lmultigrid,lpoisson_test,npost,npre,igamma,n_vcycles,&
       ldirect_solver,nx_coarsest,lsave_residual_grid,ladimensional,lsplit_temperature
!
  namelist /special_run_pars/ amplpsi,alpha_sor,Avisc,lprint_residual,&
       tolerance,maxit,gravity_z,rho0_bq,alpha_thermal,kappa,eta_0,&
       iconv_viscosity,Avisc,Bvisc,Cvisc,&
       Tbot,Tupp,Ra,iorder_sor_solver,lsave_residual,&
       ltidal_heating,ltemperature_advection,ltemperature_diffusion,lmultigrid,&
       npost,npre,igamma,n_vcycles,ldirect_solver,nx_coarsest,lsave_residual_grid,&
       lsplit_temperature
!
!  Declare index of new variables in f array. Surface density, midplane
!  temperature, and mass accretion rate.
!
  integer :: ipsi=0, iresidual=0
!
!  Diagnostic variables (needs to be consistent with reset list below).
!
  integer :: idiag_uqxmin=0, idiag_uqxmax=0, idiag_uqxrms=0, idiag_uqxm=0, idiag_uqx2m=0
  integer :: idiag_uqzmin=0, idiag_uqzmax=0, idiag_uqzrms=0, idiag_uqzm=0, idiag_uqz2m=0
  integer :: idiag_uq2m=0, idiag_uqrms=0, idiag_uqmax=0
  integer :: idiag_qtidalmin=0, idiag_qtidalmax=0, idiag_qtidalm=0
  integer :: idiag_dTdz1=0, idiag_dTdz2=0, idiag_dTdz3=0, idiag_dTdz4=0 
  integer :: idiag_nusselt_num=0, idiag_nusselt_den=0
  integer :: idiag_TTmax_cline=0, idiag_TTmin_cline=0
  integer :: idiag_devsigzz1=0,idiag_devsigzz2=0,idiag_devsigzz3=0,idiag_devsigzz4=0
  integer :: idiag_icount=0,idiag_residual=0
!
  logical :: lsolver_highorder, lsolver_loworder
!
  integer :: icount_save
  real :: residual_save
!  
  logical :: lviscosity_const
  logical :: lviscosity_var_newtonian
  logical :: lviscosity_var_blankenbach
!
  real, dimension(mz) :: Tbar = 0.0, etabar=0.0
  real :: gTT_conductive = 0.0, eta_alpha2=0.0
!
  contains
!***********************************************************************
    subroutine register_special()
!
!  Set up indices for variables in special modules.
!
!  6-oct-03/tony: coded
!  01-aug-11/wlad: adapted
!
      use FArrayManager, only: farray_register_pde,farray_register_auxiliary
!
      if (lroot) call svn_id( &
          "$Id: alphadisk.f90 19193 2012-06-30 12:55:46Z wdobler $")
!
!  Register the streamfunction. The streamfunction needs to be a pde so it
!  can have boundaries. 
!
      call farray_register_pde('psi',ipsi)
      call farray_register_auxiliary('residual',iresidual)
!
! Write to read special variables with pc_varcontent
!
      if (lroot) then
        open(4,file=trim(datadir)//'/index_special.pro',status='replace')
        write(4,*) 'ipsi= ',ipsi
        write(4,*) 'iresidual= ',iresidual
        close(4)
      endif
!
      if (lsave_residual) then
         open(9,file='residuals.dat',status='replace')
      endif
!
    endsubroutine register_special
!***********************************************************************
    subroutine initialize_special(f)
!
!  Called after reading parameters, but before the time loop.
!
!  06-oct-03/tony: coded
!  01-aug-11/wlad: adapted
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: dslab,delta_T,delta
!
      if (Lxyz(1)/(nxgrid-1) .ne. Lxyz(3)/(nzgrid-1)) then
        print*,Lxyz(1)/(nxgrid-1),Lxyz(3)/(nzgrid-1)
        call fatal_error("initialize_special","hx ne hz")
      endif
      !if (Lxyz(1)/nxgrid .ne. Lxyz(3)/nzgrid) then
      !   print*,Lxyz(1)/nxgrid,Lxyz(3)/nzgrid
      !  call fatal_error("initialize_special","dx ne dz")
      !endif
!
!  Pre-calculate Rayleigh number in case it is not given
!
      delta_T=Tbot-Tupp
      deltaT1=1./delta_T
      dslab = Lxyz(3)
      Lz1 = 1./Lxyz(3)
!
!  Conductive (linear) Temperature Gradient, used if the T=Tbar(z) + theta(x,y,z,t)
!  split is used. It similarly defines the conductive viscosity eta = etabar(z) * nu(x,y,z,t).
!
      if (lsplit_temperature) then
         do n=n1,n2
            Tbar(n) = Tbot + (n-n1)*(Tupp-Tbot)/(n2-n1)
         enddo
         etabar = exp(-Bvisc*Tbar*deltaT1)
         gTT_conductive = -delta_T*Lz1
         eta_alpha2 = (Bvisc*deltaT1*gTT_conductive)**2
      endif
!
       if (tolerance_coarsest==impossible) tolerance_coarsest=tolerance
!
!  Stuff for the coefficient of SOR. Should be between 0 and 2 for stability.
!     
      if (alpha_sor == impossible) then
        delta=min(dx,dy,dz)/dslab
        alpha_sor= 2./(1+sin(pi*delta))
      endif
!
!  Europa-specific stuff
!
      mu1_ice=1./mu_ice
      OmegaEuropa=1./EuropaPeriod
!
      kappa1=1./kappa
!
      if (ladimensional) then
         if (eta_0 /= 1.0)   call fatal_error("initialize_special",&
              "You are using dimensionless variables, set eta_0=1 and vary Ra")
         if (kappa /= 1.0)   call fatal_error("initialize_special",&
              "You are using dimensionless variables, set kappa=1 and vary Ra")
         if (delta_T /= 1.0) call fatal_error("initialize_special",&
              "You are using dimensionless variables, set Tupp=0 and Tbot=1 and vary Ra")
         if (dslab /= 1.0)   call fatal_error("initialize_special",&
              "You are using dimensionless variables, set Lxyz(3)=1 and vary Ra")
      else
         if (lstart.and.Ra/=impossible) call fatal_error("initialize_special",&
             "You are using dimensioned variables but you set the Rayleigh number. Uncomment it from your start.in")
      endif
!
      if (Ra==impossible) then
        Ra = (gravity_z*alpha_thermal*rho0_bq*delta_T*dslab**3)/(kappa*eta_0)
      endif
!
      if (lroot) print*,'Rayleigh number=',Ra
!
      select case (iorder_sor_solver)
!
      case ('low_order')
         lsolver_loworder=.true.
         lsolver_highorder=.false.
      case ('high_order') 
         lsolver_loworder=.false.
         lsolver_highorder=.true.
      case default  
        write(unit=errormsg,fmt=*) &
             'initialize_special: No such value for iorder_sor_solver: ', &
             trim(iorder_sor_solver)
        call fatal_error('initialize_special',errormsg)
      endselect
!
      select case (iconv_viscosity)
!
      case ('constant')
         lviscosity_const=.true.
         lviscosity_var_newtonian=.false.
         lviscosity_var_blankenbach=.false.
      case ('Newtonian') 
         lviscosity_const=.false.
         lviscosity_var_newtonian=.true.
         lviscosity_var_blankenbach=.false.
      case ('Blankenbach-variable') 
         lviscosity_const=.false.
         lviscosity_var_newtonian=.false.
         lviscosity_var_blankenbach=.true.
      case default
        write(unit=errormsg,fmt=*) &
             'initialize_special: No such value for iconv_viscosity: ', &
             trim(iconv_viscosity)
        call fatal_error('initialize_special',errormsg)
      endselect
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_special
!***********************************************************************
    subroutine init_special(f)
!
!  Initialise special condition; called from start.f90.
!
!  06-oct-2003/tony: coded
!  01-aug-11/wlad: adapted
!
      use Initcond, only: gaunoise
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: amplpsi_
!
!  Give an initial guess for psi
!      
      select case (initpsi)
!
      case ('noise')
         print*,'init_special: gaussian noise initialization'
         call gaunoise(amplpsi,f,ipsi)
      case ('single-mode')
         if (.not.ladimensional) call fatal_error("init_special",&
              "not yet coded for dimensioned variables. Use initpsi='noise' instead")
         print*,'init_special: single-mode initialization for constant viscosity'
         print*,'derived from the solution of the dispersion relation'
         amplpsi_ = -ampltt * Ra*kx_TT/(kz_TT**2 + kx_TT**2)**2
         do n=n1,n2
           do m=m1,m2
             f(l1:l2,m,n,ipsi) = f(l1:l2,m,n,ipsi) + &
                  amplpsi_*sin(kx_TT*x(l1:l2))*sin(kz_TT*z(n))
           enddo
         enddo
         print*,'init_special: done single-mode initialization'
!
      case ('zero')
         f(l1:l2,m1:m2,n1:n2,ipsi)=0.
!         
      case ('nothing')
                   
      case default
!
!  Catch unknown values.
!
          write(unit=errormsg,fmt=*) 'No such value initpsi = ',trim(initpsi)
          call fatal_error('init_special',errormsg)
!
      endselect
      !
      call special_after_boundary(f)
!
    endsubroutine init_special
!***********************************************************************
    subroutine pencil_criteria_special()
!
!  All pencils that this special module depends on are specified here.
!
!  18-07-06/tony: coded
!
      if (ltemperature_advection) lpenc_requested(i_gTT)=.true.
      if (ltemperature_diffusion) lpenc_requested(i_del2TT)=.true.
      if (ltidal_heating) lpenc_requested(i_TT1)=.true.
!
      if (idiag_dTdz1/=0 .or. &
          idiag_dTdz2/=0 .or. &
          idiag_dTdz3/=0 .or. &
          idiag_dTdz4/=0 .or. &
          idiag_nusselt_num/=0) lpenc_diagnos(i_gTT)=.true.
!
      if (idiag_nusselt_den/=0 .or. &
          idiag_TTmax_cline/=0 .or. &
          idiag_TTmin_cline/=0) lpenc_diagnos(i_TT)=.true.
!
    endsubroutine pencil_criteria_special
!***********************************************************************
    subroutine pencil_interdep_special(lpencil_in)
!
!  Interdependency among pencils provided by this module are specified here.
!
!  18-07-06/tony: coded
!
      logical, dimension(npencils), intent(inout) :: lpencil_in
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_special
!***********************************************************************
    subroutine calc_pencils_special(f,p)
!
!  Calculate Special pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  24-nov-04/tony: coded
!  01-aug-11/wlad: adapted
!
      use Deriv, only: der,derij
      use Sub, only: u_dot_grad
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
      real, dimension (nx) :: derxzpsi
!
      intent(inout) :: f
      intent(inout) :: p
!
      if (ltidal_heating .or. &
          idiag_devsigzz1 /=0 .or. &
          idiag_devsigzz2 /=0 .or. &
          idiag_devsigzz3 /=0 .or. &
          idiag_devsigzz4 /=0) then 
         if (lviscosity_const) then 
            q%eta=eta_0
         else if (lviscosity_var_newtonian) then
            q%eta=eta_0*exp(Avisc*(TT_melt*p%TT1 - 1.))
         else if (lviscosity_var_blankenbach) then
            q%eta=eta_0*exp(-Bvisc*p%TT*deltaT1 + Cvisc*(1-z(n))*Lz1)
         else   
            call fatal_error("calc_pencils_special",&
                 "The world is flat, and we never got here.")
         endif
      endif
!
!  Calculate velocities via streamfunction. 
!
      if (ltemperature_advection) then 
         call der( f,ipsi,q%uu(:,1),3)
         q%uu(:,2)=0.
         call der(-f,ipsi,q%uu(:,3),1)
!
!  Calculate temperature advection term. 
!  
         call u_dot_grad(f,iTT,p%gTT,q%uu,q%ugTT)
      endif
!
      if (ltidal_heating) then
        call fatal_error("calc_pencils_special","these are not yet normalized")
        q%qtidal = .5*q%eta*(epsi_0*OmegaEuropa)**2 / (1+(OmegaEuropa*q%eta*mu1_ice)**2)
      endif
!
      if (idiag_uq2m/=0.or.idiag_uqrms/=0.or.idiag_uqmax/=0) &
           q%u2=q%uu(:,1)**2 + q%uu(:,3)**2
!
!  duz/dz for topography calculation. duz/dz = d/dz (-dpsi/dx) = -d2dxdz psi = derij(psi,x,z)
!      
      if (idiag_devsigzz1/=0.or.idiag_devsigzz2/=0.or.idiag_devsigzz3/=0.or.idiag_devsigzz4/=0) then 
         call derij(f,ipsi,derxzpsi,1,3)
         q%devsigzz = q%eta * derxzpsi
      endif
!
      call keep_compiler_quiet(p)
!
    endsubroutine calc_pencils_special
!***********************************************************************
    subroutine special_after_boundary(f)
!
      !use Boundcond, only: update_ghosts
!
      real, dimension (mx,my,mz,mfarray) :: f   
!    
      call solve_for_psi(f)
!
    endsubroutine special_after_boundary
!***********************************************************************
    subroutine solve_for_psi(f)
!
      real, dimension (mx,my,mz,mfarray) :: f   
      real, dimension (mx,mz) :: psi,psi_old,eta
!
      real, dimension (nx,nz) :: rhs,tmp
      real, dimension (nx,nz) :: alpha_factor,beta_factor
!
      real :: residual,aout,dTTdx
      real :: alpha=impossible,beta=impossible
      real :: Rayleigh_factor
!
      integer :: icount,i
!
      if (ladimensional) then
         Rayleigh_factor=Ra
      else
         Rayleigh_factor=alpha_thermal*rho0_bq*gravity_z
      endif
!
!  Define r.h.s.
!
      if (lpoisson_test) then
         psi=0. 
         rhs=1.
      else
!
         call calc_viscosity(f,eta)
!
         do n=n1,n2; do m=m1,m2
            do i=l1,l2
               call f_function(eta,aout,i,n) ; alpha = aout/eta(i,n)
               call g_function(eta,aout,i,n) ; beta  = aout/eta(i,n)
!
               dTTdx=dx_1(i)/60*(+ 45.0*(f(i+1,m,n,iTT)-f(i-1,m,n,iTT)) &
                                 -  9.0*(f(i+2,m,n,iTT)-f(i-2,m,n,iTT)) &
                                 +      (f(i+3,m,n,iTT)-f(i-3,m,n,iTT)))            
!
               if (lsplit_temperature) then 
                 rhs(i-l1+1,n-n1+1) = Rayleigh_factor*dTTdx/(eta(i,n)*etabar(n))
                 alpha_factor(i-l1+1,n-n1+1)=alpha + eta_alpha2
               else
                 rhs(i-l1+1,n-n1+1) = Rayleigh_factor*dTTdx/eta(i,n)
                 alpha_factor(i-l1+1,n-n1+1)=alpha
               endif
               beta_factor(i-l1+1,n-n1+1)=beta
!
            enddo;enddo
         enddo
!         
         psi(l1:l2,n1:n2)=f(l1:l2,mpoint,n1:n2,ipsi)
      
         call update_bounds_psi(psi)
      endif
!
      if (lmultigrid) then
         call multigrid(f,psi,rhs,alpha_factor,beta_factor)
         call update_bounds_psi(psi)
      else
!
!  Initial residual
!
         residual=1e33
!
!  Start counter
!      
         icount=0
!         
         do while (residual > tolerance)
!
!  Calculate psi via SOR
!
            psi_old=psi
!
            call successive_over_relaxation(psi,rhs,alpha_factor,beta_factor)
!
            call update_bounds_psi(psi)
!
            tmp=(psi(l1:l2,n1:n2) - psi_old(l1:l2,n1:n2))**2
            residual = sqrt(sum(tmp)/sum(psi(l1:l2,n1:n2)**2))
!
            if (lsave_residual_grid) then
               do m=1,my
                  f(:,m,:,iresidual) = psi-psi_old
               enddo
            endif
!
! Increase counter. 
!
            icount=icount+1
            if (lprint_residual) &
                 print*, icount,residual,tolerance
            if (lsave_residual) then
               write(9,*) icount,residual
            endif
!
            if (icount >= maxit) then
               call fatal_error("solve_for_psi","reached limit of iterations: did not converge.")
               if (lsave_residual) close(9)
               goto 200
            endif
!
         enddo
!
200      continue
         
         if (lsave_residual) close(9)
!
         icount_save=icount
         residual_save=residual
      endif
!
      do m=m1,m2
         f(:,m,:,ipsi)=psi
      enddo
!
    endsubroutine solve_for_psi
!***********************************************************************
    subroutine successive_over_relaxation(u,r,alp,bet,h)
!
      real, dimension (:,:) :: u
      real, dimension (:,:) :: r
      real, dimension (:,:) :: alp,bet
      real :: cc,ufactor,vterm,u_new,alpha_sor_
      real :: alpha=impossible,beta=impossible
      real, optional :: h
      integer :: ic,nc,k,l
      integer :: nux,nuz,nrx,nrz,nn1,nn2,ll1,ll2,lcc1,lcc2
!
      nux=size(u,1); nuz=size(u,2)
      nrx=size(r,1); nrz=size(r,2)
      ll1=nghost+1; ll2=nux-nghost
      nn1=nghost+1; nn2=nuz-nghost
!
      if (present(h)) then 
         alpha_sor_=2./(1+sin(pi*h))
      else
         alpha_sor_=alpha_sor
      endif
!
      if (lperi(1)) then
         lcc1=ll1
         lcc2=ll2
      else
         lcc1=ll1+1
         lcc2=ll2-1
      endif
!
      do nc=nn1+1,nn2-1
         k=nc-nn1+1
         do ic=lcc1,lcc2
           l=ic-ll1+1
!     
           cc    =   r(l,k)
           alpha = alp(l,k)
           beta  = bet(l,k)
!     
           if (lsolver_highorder) then
              if (present(h)) then
                 call solve_highorder(u,alpha,beta,ufactor,vterm,ic,nc)
              else
                 call solve_highorder(u,alpha,beta,ufactor,vterm,ic,nc,h)
              endif
           else
              if (present(h)) then
                 call solve_loworder(u,alpha,beta,ufactor,vterm,ic,nc)
              else
                 call solve_loworder(u,alpha,beta,ufactor,vterm,ic,nc,h)
              endif
           endif
!     
           u_new=(cc-vterm)/ufactor
           u(ic,nc) = u(ic,nc)*(1-alpha_sor_) +  alpha_sor_*u_new
!     
        enddo
     enddo
!
    endsubroutine successive_over_relaxation
!***********************************************************************
    subroutine solve_highorder(a,alpha,beta,u,v,i,k,h)
!
!   Solve the momentum equation of mantle convection for the streamfunction
!      
!   F2(psi,i,k,h) + G2(psi,i,k,h) + A*F(psi,i,k,h) + B*G(psi,i,k,h) = C
!         
!   d4/dx4 + d4/dz4 + 2*d4dx2dz2 + A*d2/dz2 - A*d2dx2 + B*d2dxdz = C 
!         
      real, dimension (:,:), intent(in) :: a
      real, intent(in) :: alpha,beta
      real, intent(out) :: u,v
      real :: dx1,dz1,fac,fac2
      real :: d4dx4,d4dz4,d2dx2,d2dz2,d2dxdz,d4dx2dz2
!
      integer :: i,k
!
      real, optional :: h
!
      if (present(h)) then
         dx1 = 1./h
         dz1 = 1./h
      else
         dx1=dx_1(i)
         dz1=dz_1(k)
      endif
!      
!  The quantity u collects all terms that multiply psi[i,k]
!  v collects everything else. 
!      
      u = 0.
      v = 0.
!           
! d4dx4
!
      fac=(1.0/6)*dx1**4
      d4dx4=fac*(+ 56.0* a(i,k)               &
                 - 39.0*(a(i+1,k)+a(i-1,k)) &
                 + 12.0*(a(i+2,k)+a(i-2,k)) &
                 -      (a(i+3,k)+a(i-3,k)))
!
! fac2 is the factor by which this term is multiplied in the momentum equation. 
!
      fac2 = 1.
      u = u + fac2*fac*56.
      v = v + fac2*d4dx4 
!
! d4dz4
!
      fac=(1.0/6)*dz1**4
      d4dz4=fac*(+ 56.0* a(i,k)               &
                 - 39.0*(a(i,k+1)+a(i,k-1)) &
                 + 12.0*(a(i,k+2)+a(i,k-2)) &
                 -      (a(i,k+3)+a(i,k-3)))

      fac2 = 1.
      u = u + fac2*fac*56.
      v = v + fac2*d4dz4
!
! d2dx2         
!
      fac=(1./180)*dx1**2
      d2dx2=fac*(-490.0* a(i,k)               &
                 +270.0*(a(i+1,k)+a(i-1,k)) &
                 - 27.0*(a(i+2,k)+a(i-2,k)) &
                 +  2.0*(a(i+3,k)+a(i-3,k)))
          
      fac2 = -alpha
      u = u + fac2*fac*(-490)
      v = v + fac2*d2dx2
!
! d2dz2
!          
      fac=(1./180)*dz1**2
      d2dz2=fac*(-490.0* a(i,k) &
                 +270.0*(a(i,k+1)+a(i,k-1)) &
                 - 27.0*(a(i,k+2)+a(i,k-2)) &
                 +  2.0*(a(i,k+3)+a(i,k-3)))
          
      fac2 = alpha 
      u = u + fac2*fac*(-490)
      v = v + fac2*d2dz2
!                    
! d2dxdz
!          
      fac=(1./60.**2)*dx1*dz1
      d2dxdz=fac*( &
           45.*((45.*(a(i+1,k+1)-a(i-1,k+1))   &
                 -9.*(a(i+2,k+1)-a(i-2,k+1))   &
                    +(a(i+3,k+1)-a(i-3,k+1)))  &
               -(45.*(a(i+1,k-1)-a(i-1,k-1))   &
                 -9.*(a(i+2,k-1)-a(i-2,k-1))   &
                    +(a(i+3,k-1)-a(i-3,k-1)))) &
           -9.*((45.*(a(i+1,k+2)-a(i-1,k+2))   &
                 -9.*(a(i+2,k+2)-a(i-2,k+2))   &
                    +(a(i+3,k+2)-a(i-3,k+2)))  &
               -(45.*(a(i+1,k-2)-a(i-1,k-2))   &
                 -9.*(a(i+2,k-2)-a(i-2,k-2))   &
                    +(a(i+3,k-2)-a(i-3,k-2)))) &
              +((45.*(a(i+1,k+3)-a(i-1,k+3))   &
                 -9.*(a(i+2,k+3)-a(i-2,k+3))   &
                    +(a(i+3,k+3)-a(i-3,k+3)))  &
               -(45.*(a(i+1,k-3)-a(i-1,k-3))   &
                 -9.*(a(i+2,k-3)-a(i-2,k-3))   &
                    +(a(i+3,k-3)-a(i-3,k-3)))) &
           )         
          
      fac2 = beta
      u = u + fac2*fac*0.
      v = v + fac2*d2dxdz
!
! d4dx2dz2       
!
      fac=(1./180.**2)*dx1**2*dz1**2
      d4dx2dz2=fac*( &
           - 490*(  -490.0* a(i,k)                   &
                    +270.0*(a(i+1,k)+a(i-1,k))       &
                    - 27.0*(a(i+2,k)+a(i-2,k))       &
                    +  2.0*(a(i+3,k)+a(i-3,k)))      &
           + 270.*((-490.0* a(i,k+1)                 &
                    +270.0*(a(i+1,k+1)+a(i-1,k+1))   &
                    - 27.0*(a(i+2,k+1)+a(i-2,k+1))   &
                    +  2.0*(a(i+3,k+1)+a(i-3,k+1)))  &
                  +(-490.0* a(i,k-1)                 &
                    +270.0*(a(i+1,k-1)+a(i-1,k-1))   &
                    - 27.0*(a(i+2,k-1)+a(i-2,k-1))   &
                    +  2.0*(a(i+3,k-1)+a(i-3,k-1)))) &
             -27.*((-490.0* a(i,k+2)                 &
                    +270.0*(a(i+1,k+2)+a(i-1,k+2))   &
                    - 27.0*(a(i+2,k+2)+a(i-2,k+2))   &
                    +  2.0*(a(i+3,k+2)+a(i-3,k+2)))  &
                  +(-490.0* a(i,k-2)                 &
                    +270.0*(a(i+1,k-2)+a(i-1,k-2))   &
                    - 27.0*(a(i+2,k-2)+a(i-2,k-2))   &
                    +  2.0*(a(i+3,k-2)+a(i-3,k-2)))) &
              +2.*((-490.0* a(i,k+3)                 &
                    +270.0*(a(i+1,k+3)+a(i-1,k+3))   &
                    - 27.0*(a(i+2,k+3)+a(i-2,k+3))   &
                    +  2.0*(a(i+3,k+3)+a(i-3,k+3)))  &
                  +(-490.0* a(i,k-3)                 &
                    +270.0*(a(i+1,k-3)+a(i-1,k-3))   &
                    - 27.0*(a(i+2,k-3)+a(i-2,k-3))   &
                    +  2.0*(a(i+3,k-3)+a(i-3,k-3)))) &
           )
!
      fac2 = 2.
      u = u + fac2*fac*(490.**2)
      v = v + fac2*d4dx2dz2
!
! Must remove u*a(i,k) from v 
!
      v = v - u*a(i,k)
!
    endsubroutine solve_highorder
!***********************************************************************
    subroutine solve_loworder(a,alpha,beta,u,v,i,k,h)
!
      real, dimension (:,:), intent(in) :: a
      real, intent(in) :: alpha,beta
      real, intent(out) :: u,v
!
      real :: v_1, v_2, v_3, v_4, v_5
      real :: v51,v52,v53,dx1,dz1
      integer :: i,k
!
      real, optional :: h
!
      if (present(h)) then
         dx1 = 1./h
         dz1 = 1./h
      else
         dx1=dx_1(i)
         dz1=dz_1(k)
      endif
!
! a is updated on the fly. So, do it swiping, so that a(i-1), that is updated, gets used for a(i)
!
      u = 20.*(dx1**2 * dz1**2)
!
      v_1 = alpha*(a(i,k+1) + a(i,k-1) - a(i+1,k) - a(i-1,k))*(dx1*dz1)
!
      v_2 = ((a(i,k+2) + a(i,k-2)) + (a(i+2,k) + a(i-2,k)) - 2.*(a(i+1,k+1) + a(i+1,k-1) + a(i-1,k+1) + a(i-1,k-1)))&
              *(dx1**2 * dz1**2)
!
      v_3 = beta*(a(i+1,k+1) - a(i+1,k-1) - a(i-1,k+1) + a(i-1,k-1))*(dx1*dz1)
!      
      v_4 = 4*(a(i+1,k+1) - 2.*a(i+1,k) + a(i+1,k-1) - 2.*a(i,k+1) - 2.*a(i,k-1) + a(i-1,k+1) - 2.*a(i-1,k) + a(i-1,k-1))&
              *(dx1**2 * dz1**2)

      v = v_1 + v_2 + v_3 + v_4

!      u = 2*alpha*(dx1**2 - dz1**2) + 6.*(dx1**4 + dz1**4) + &
!              4.*beta*dx1*dz1 + 8.*dx1**2*dz1**2
!      v_1 = (a(i,k+1) + a(i,k-1))*dz1**2 - (a(i+1,k) + a(i-1,k)) * dx1**2
!      v_2 = ( a(i,k+2) - 4.*a(i,k+1) - 4.*a(i,k-1) + a(i,k-2))   * dz1**4
!      v_3 = ( a(i+2,k) - 4.*a(i+1,k) - 4.*a(i-1,k) + a(i-2,k))   * dx1**4
!      v_4 = (-a(i-1,k) -    a(i,k-1) +    a(i-1,k-1)) * dx1*dz1
!!
!      v51 = a(i+1,k+1) - 2.*a(i,k+1) + a(i-1,k+1)
!      v52 = a(i+1,k  )               + a(i-1,k  )
!      v53 = a(i+1,k-1) - 2.*a(i,k-1) + a(i-1,k-1)
!      v_5 = (v51 - 2.*v52 + v53) * dx1**2*dz1**2
!!
!      v = alpha*v_1 + v_2 + v_3 + 4.*beta*v_4 + 2.*v_5
!
    endsubroutine solve_loworder
!***********************************************************************      
    subroutine calc_viscosity(f,eta)
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,mz), intent(out) :: eta
!
!  Viscosities normalized by eta_0
!
      select case (iconv_viscosity)
!
      case ('constant')
        eta = eta_0
!
      case ('Netwonian') 
        eta = eta_0*exp(Avisc * (TT_melt/f(:,mpoint,:,iTT) - 1.))
!
      case ('Blankenbach-variable')
         do n=1,mz
            eta(:,n) = eta_0*exp(-Bvisc * f(:,mpoint,n,iTT)*deltaT1 + &
                                  Cvisc * (1-z(n)*Lz1)*Lz1 )
         enddo
!
      case default  
        write(unit=errormsg,fmt=*) &
             'calc_viscosity: No such value for iconv_viscosity: ', &
             trim(iconv_viscosity)
        call fatal_error('calc_viscosity',errormsg)
      endselect
!     
    endsubroutine calc_viscosity
!***********************************************************************
    subroutine f_function(a,aout,i,n)
!
!  Calculate F = dz2 - dx2, 6th order derivatives. 
!
      real, dimension(mx,mz), intent(in) :: a
      integer, intent(in) :: i,n
      real :: df2x,df2z,aout,fac
!
      fac=(1./180)*dx_1(i)**2
      df2x=fac*(-490.0*a(i,n) &
                +270.0*(a(i+1,n)+a(i-1,n)) &
                - 27.0*(a(i+2,n)+a(i-2,n)) &
                +  2.0*(a(i+3,n)+a(i-3,n)))

      fac=(1./180)*dz_1(n)**2
      df2z=fac*(-490.0*a(i,n) &
                +270.0*(a(i,n+1)+a(i,n-1)) &
                - 27.0*(a(i,n+2)+a(i,n-2)) &
                +  2.0*(a(i,n+3)+a(i,n-3)))

      aout = df2z-df2x
!
    endsubroutine f_function
!***********************************************************************
    subroutine g_function(a,aout,i,n)
!
!   Calculate G = dzdx, 6th order derivatives. 
!
      real, dimension (mx,mz), intent(in) :: a
      integer, intent(in) :: i,n
      real :: aout,fac
!
      fac=(1./60.**2)*dz_1(n)*dx_1(i)
      aout=fac*( &
              45.*((45.*(a(i+1,n+1)-a(i-1,n+1))  &
                    -9.*(a(i+2,n+1)-a(i-2,n+1))  &
                       +(a(i+3,n+1)-a(i-3,n+1))) &
                  -(45.*(a(i+1,n-1)-a(i-1,n-1))  &
                    -9.*(a(i+2,n-1)-a(i-2,n-1))  &
                       +(a(i+3,n-1)-a(i-3,n-1))))&
              -9.*((45.*(a(i+1,n+2)-a(i-1,n+2))  &
                    -9.*(a(i+2,n+2)-a(i-2,n+2))  &
                       +(a(i+3,n+2)-a(i-3,n+2))) &
                  -(45.*(a(i+1,n-2)-a(i-1,n-2))  &
                    -9.*(a(i+2,n-2)-a(i-2,n-2))  &
                       +(a(i+3,n-2)-a(i-3,n-2))))&
                 +((45.*(a(i+1,n+3)-a(i-1,n+3))  &
                    -9.*(a(i+2,n+3)-a(i-2,n+3))  &
                       +(a(i+3,n+3)-a(i-3,n+3))) &
                  -(45.*(a(i+1,n-3)-a(i-1,n-3))  &
                    -9.*(a(i+2,n-3)-a(i-2,n-3))  &
                       +(a(i+3,n-3)-a(i-3,n-3))))&
                   )
!
    endsubroutine g_function
!***********************************************************************
    subroutine update_bounds_psi(a)
!
      real, dimension(:,:) :: a
      integer :: i,nn1,nn2,ll1,ll2,mxa,mza,ll1i,ll2i
!
      mxa=size(a,1); mza=size(a,2)
      ll1=nghost+1; ll2=mxa-nghost
      nn1=nghost+1; nn2=mza-nghost
!
!  Set boundary of psi - vertical, zero
!
      a(:,nn1)=0.
      a(:,nn2)=0.
!
!  Zero also the second derivative
!
      do i=1,nghost
        a(:,nn1-i) = 2*a(:,nn1) - a(:,nn1+i)
      enddo
      do i=1,nghost
        a(:,nn2+i) = 2*a(:,nn2) - a(:,nn2-i)
      enddo
!
      if (lperi(1)) then 
!
!  Periodic in the lateral 
!
         ll1i=ll1+nghost-1;ll2i=ll2-nghost+1
         a(1   :ll1-1,:) = a(ll2i:ll2,:)
         a(ll2+1:mxa ,:) = a(ll1:ll1i,:)
      else
         a(ll1,:)=0.
         a(ll2,:)=0.
!
         do i=1,nghost
            a(ll1-i,:) = 2*a(ll1,:) - a(ll1+i,:)
         enddo
         do i=1,nghost
            a(ll2+i,:) = 2*a(ll2,:) - a(ll2-i,:)
         enddo
      endif
!
    endsubroutine update_bounds_psi
!***********************************************************************
    subroutine special_calc_energy(f,df,p)
!
      use Diagnostics, only: max_mn_name,sum_mn_name,integrate_mn_name
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      real, dimension (nx) :: dTdz1,dTdz2,dTdz3,dTdz4,nusselt_num,nusselt_den,TTmin_cline,TTmax_cline
      real, dimension (nx) :: devsigzz1,devsigzz2,devsigzz3,devsigzz4
      real, dimension (nx) :: diffus_special,advec_special
      type (pencil_case) :: p
!      
!  Advection
!
      if (ltemperature_advection) then
         df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) - q%ugTT
         if (lsplit_temperature) then
            df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) - q%uu(:,3)*gTT_conductive
         endif
      endif
!
!  Conduction (diffusion)
!
      if (ltemperature_diffusion) & 
           df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + kappa*p%del2TT
!
!  Tidal heating 
!
      if (ltidal_heating) &
           df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + q%qtidal
!
      if (lfirst.and.ldt) then 
        advec_special=sum(abs(q%uu)*dline_1,2)
        maxadvec=maxadvec+advec_special
!
        diffus_special = kappa*dxyz_2
        maxdiffus=max(maxdiffus,diffus_special)
!
        if (headtt.or.ldebug) then
          print*,'special_calc_energy: max(advec_special)  =',maxval(advec_special)
          print*,'special_calc_energy: max(diffus_special) =',maxval(diffus_special)
        endif
      endif
!
      if (ldiagnos) then 
         if (idiag_uqxmin/=0) call max_mn_name(-q%uu(:,1)   ,idiag_uqxmin,lneg=.true.)
         if (idiag_uqxmax/=0) call max_mn_name( q%uu(:,1)   ,idiag_uqxmax)
         if (idiag_uqxm/=0)   call sum_mn_name( q%uu(:,1)   ,idiag_uqxm)
         if (idiag_uqx2m/=0)  call sum_mn_name( q%uu(:,1)**2,idiag_uqx2m)              
         if (idiag_uqxrms/=0) call sum_mn_name( q%uu(:,1)**2,idiag_uqxrms,lsqrt=.true.)
!
         if (idiag_uqzmin/=0) call max_mn_name(-q%uu(:,3)   ,idiag_uqzmin,lneg=.true.)
         if (idiag_uqzmax/=0) call max_mn_name( q%uu(:,3)   ,idiag_uqzmax)
         if (idiag_uqzm/=0)   call sum_mn_name( q%uu(:,3)   ,idiag_uqzm)
         if (idiag_uqz2m/=0)  call sum_mn_name( q%uu(:,3)**2,idiag_uqz2m)
         if (idiag_uqzrms/=0) call sum_mn_name( q%uu(:,3)**2,idiag_uqzrms,lsqrt=.true.)
!
         if (idiag_uq2m/=0)   call sum_mn_name(q%u2,idiag_uq2m)
         if (idiag_uqrms/=0)  call sum_mn_name(q%u2,idiag_uqrms,lsqrt=.true.)
         if (idiag_uqmax/=0)  call max_mn_name(q%u2,idiag_uqmax,lsqrt=.true.)
!
         if (idiag_qtidalmin/=0) call max_mn_name(-q%qtidal,idiag_qtidalmin,lneg=.true.)
         if (idiag_qtidalmax/=0) call max_mn_name( q%qtidal,idiag_qtidalmax)
         if (idiag_qtidalm/=0)   call sum_mn_name( q%qtidal,idiag_qtidalm)
         if (idiag_icount/=0)    call max_mn_name(0*x(l1:l2)+icount_save,idiag_icount)
         if (idiag_residual/=0)  call max_mn_name(0*x(l1:l2)+residual_save,idiag_residual)
!
!  Calculate for benchmark diagnostic the (negative of the) temperature gradient at the
!  four corners of the grid, labeled as below: 
!         
! 1 .______. 2
!   |      |
!   |      |
!   .______. 
! 4          3         
!         
!  I.e., 1 is (x,z)=( 0,Lz), top     left corner, above upwelling
!        2 is (x,z)=(Lx,Lz), top    right corner, above downwelling
!        3 is (x,z)=(Lx, 0), bottom right corner, below downwelling 
!        4 is (x,z)=( 0, 0), bottom  left corner, below upwelling
!
         if (idiag_dTdz1/=0) then
           if (n == n2) then
             dTdz1 = -p%gTT(1,3) - gTT_conductive
           else
             dTdz1 = -impossible
           endif
           call max_mn_name(dTdz1,idiag_dTdz1)
         endif
!
         if (idiag_dTdz2/=0) then
           if (n == n2) then
             dTdz2 = -p%gTT(nx,3) - gTT_conductive
           else
             dTdz2 = -impossible
           endif
           call max_mn_name(dTdz2,idiag_dTdz2)
         endif
!
         if (idiag_dTdz3/=0) then
           if (n == n1) then
             dTdz3 = -p%gTT(nx,3) - gTT_conductive
           else
             dTdz3 = -impossible
           endif
           call max_mn_name(dTdz3,idiag_dTdz3)
         endif
!
         if (idiag_dTdz4/=0) then
           if (n == n1) then
             dTdz4 = -p%gTT(1,3) - gTT_conductive
           else
             dTdz4 = -impossible
           endif
           call max_mn_name(dTdz4,idiag_dTdz4)
         endif
!
         if (idiag_nusselt_num/=0) then
           if (n == n2) then
             nusselt_num=-p%gTT(:,3) - gTT_conductive
           else
             nusselt_num=0.
           endif
           call integrate_mn_name(nusselt_num,idiag_nusselt_num)
         endif
!         
         if (idiag_nusselt_den/=0) then
           if (n == n1) then
             nusselt_den=p%TT + Tbar(n)
           else
             nusselt_den=0.
           endif
           call integrate_mn_name(nusselt_den,idiag_nusselt_den)
         endif
!         
!  temperature extreme at center line
!  min(abs(dTdz)) and z le 0.5
!
!  min(abs(dTdz)) and z ge 0.5
!
!  record position as well

         if (idiag_TTmax_cline/=0) then
            if (z(n) .ge. 0.5*Lxyz(3)) then
               TTmax_cline=p%TT(nx/2) + Tbar(n)
            else
               TTmax_cline=-impossible
            endif
            call max_mn_name(TTmax_cline,idiag_TTmax_cline)
         endif

         if (idiag_TTmin_cline/=0) then
            if (z(n) .lt. 0.5*Lxyz(3)) then
               TTmin_cline=p%TT(nx/2) + Tbar(n)
            else
               TTmin_cline=impossible
            endif
            call max_mn_name(-TTmin_cline,idiag_TTmin_cline,lneg=.true.)
         endif
!         
         if (idiag_devsigzz1/=0) then
            if (n == n2) then
               devsigzz1 = q%devsigzz(1)
            else
               devsigzz1 = -impossible
            endif
            call max_mn_name(devsigzz1,idiag_devsigzz1)
         endif
!         
         if (idiag_devsigzz2/=0) then
            if (n == n2) then
               devsigzz2 = q%devsigzz(nx)
            else
               devsigzz2 = -impossible
            endif
            call max_mn_name(devsigzz2,idiag_devsigzz2)
         endif
!
         if (idiag_devsigzz3/=0) then
            if (n == n1) then
               devsigzz3 = q%devsigzz(nx)
            else
               devsigzz3 = -impossible
            endif
            call max_mn_name(devsigzz3,idiag_devsigzz3)
         endif
!
         if (idiag_devsigzz4/=0) then
            if (n == n1) then
               devsigzz4 = q%devsigzz(1)
            else
               devsigzz4 = -impossible
            endif
            call max_mn_name(devsigzz4,idiag_devsigzz4)
         endif         
      endif
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(p)
!
    endsubroutine special_calc_energy
!***********************************************************************
    subroutine read_special_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=special_init_pars, IOSTAT=iostat)
!
    endsubroutine read_special_init_pars
!***********************************************************************
    subroutine write_special_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_init_pars)
!
    endsubroutine write_special_init_pars
!***********************************************************************
    subroutine read_special_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=special_run_pars, IOSTAT=iostat)
!
    endsubroutine read_special_run_pars
!***********************************************************************
    subroutine write_special_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_run_pars)
!
    endsubroutine write_special_run_pars
!***********************************************************************
    subroutine rprint_special(lreset,lwrite)
!
!  Reads and registers print parameters relevant to special.
!
!  06-oct-03/tony: coded
!
      use Diagnostics, only: parse_name
      use FArrayManager, only: farray_index_append
!
      logical :: lreset
      logical, optional :: lwrite
!
      logical :: lwr
      integer :: iname, inamex
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  Reset everything in case of reset.
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_uqxmin=0; idiag_uqxmax=0; idiag_uqxrms=0; idiag_uqxm=0; idiag_uqx2m=0
        idiag_uqzmin=0; idiag_uqzmax=0; idiag_uqzrms=0; idiag_uqzm=0; idiag_uqz2m=0
        idiag_uq2m=0; idiag_uqrms=0; idiag_uqmax=0; idiag_qtidalmin=0
        idiag_qtidalmax=0; idiag_qtidalm=0
        idiag_dTdz1=0; idiag_dTdz2=0; idiag_dTdz3=0; idiag_dTdz4=0
        idiag_nusselt_num=0; idiag_nusselt_den=0
        idiag_TTmax_cline=0; idiag_TTmin_cline=0
        idiag_devsigzz1=0; idiag_devsigzz2=0; idiag_devsigzz3=0; idiag_devsigzz4=0
        idiag_icount=0; idiag_residual=0
        cformv=''
      endif
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'uqxmin',idiag_uqxmin)
        call parse_name(iname,cname(iname),cform(iname),'uqxmax',idiag_uqxmax)
        call parse_name(iname,cname(iname),cform(iname),'uqxrms',idiag_uqxrms)
        call parse_name(iname,cname(iname),cform(iname),'uqxm',idiag_uqxm)
        call parse_name(iname,cname(iname),cform(iname),'uqx2m',idiag_uqx2m)
!
        call parse_name(iname,cname(iname),cform(iname),'uqzmin',idiag_uqzmin)
        call parse_name(iname,cname(iname),cform(iname),'uqzmax',idiag_uqzmax)
        call parse_name(iname,cname(iname),cform(iname),'uqzrms',idiag_uqzrms)
        call parse_name(iname,cname(iname),cform(iname),'uqzm',idiag_uqzm)
        call parse_name(iname,cname(iname),cform(iname),'uqz2m',idiag_uqz2m)
!
        call parse_name(iname,cname(iname),cform(iname),'uq2m',idiag_uq2m)
        call parse_name(iname,cname(iname),cform(iname),'uqrms',idiag_uqrms)
        call parse_name(iname,cname(iname),cform(iname),'uqmax',idiag_uqmax)
!
        call parse_name(iname,cname(iname),cform(iname),'qtidalmin',idiag_qtidalmin)
        call parse_name(iname,cname(iname),cform(iname),'qtidalmax',idiag_qtidalmax)
        call parse_name(iname,cname(iname),cform(iname),'qtidalm',idiag_qtidalm)
!
        call parse_name(iname,cname(iname),cform(iname),'dTdz1',idiag_dTdz1)
        call parse_name(iname,cname(iname),cform(iname),'dTdz2',idiag_dTdz2)
        call parse_name(iname,cname(iname),cform(iname),'dTdz3',idiag_dTdz3)
        call parse_name(iname,cname(iname),cform(iname),'dTdz4',idiag_dTdz4)
!
        call parse_name(iname,cname(iname),cform(iname),'nusselt_num',idiag_nusselt_num)
        call parse_name(iname,cname(iname),cform(iname),'nusselt_den',idiag_nusselt_den)  
!
        call parse_name(iname,cname(iname),cform(iname),'TTmax_cline',idiag_TTmax_cline)
        call parse_name(iname,cname(iname),cform(iname),'TTmin_cline',idiag_TTmin_cline)
!
        call parse_name(iname,cname(iname),cform(iname),'devsigzz1',idiag_devsigzz1)
        call parse_name(iname,cname(iname),cform(iname),'devsigzz2',idiag_devsigzz2)
        call parse_name(iname,cname(iname),cform(iname),'devsigzz3',idiag_devsigzz3)
        call parse_name(iname,cname(iname),cform(iname),'devsigzz4',idiag_devsigzz4)
!        
        call parse_name(iname,cname(iname),cform(iname),'icount',idiag_icount)
        call parse_name(iname,cname(iname),cform(iname),'residual',idiag_residual)
     enddo
!
!  Check for those quantities for which we want yz-averages.
!
      do inamex=1,nnamex
        !call parse_name(inamex,cnamex(inamex),cformx(inamex),'sigmamx', &
        !    idiag_sigmamx)
      enddo
!
!  check for those quantities for which we want video slices
!
      if (lwrite_slices) then
        where(cnamev=='psi') cformv='DEFINED'
      endif
!
      if (lwr) then
        call farray_index_append('uqxmin',idiag_uqxmin)
        call farray_index_append('uqxmax',idiag_uqxmax)
        call farray_index_append('uqxrms',idiag_uqxrms)
        call farray_index_append('uqxm',idiag_uqxm)
        call farray_index_append('uqx2m',idiag_uqx2m)
!
        call farray_index_append('uqzmin',idiag_uqzmin)
        call farray_index_append('uqzmax',idiag_uqzmax)
        call farray_index_append('uqzrms',idiag_uqzrms)
        call farray_index_append('uqzm',idiag_uqzm)
        call farray_index_append('uqz2m',idiag_uqz2m)
!
        call farray_index_append('uq2m',idiag_uq2m)
        call farray_index_append('uqrms',idiag_uqrms)
        call farray_index_append('uqmax',idiag_uqmax)
!
        call farray_index_append('qtidalmin',idiag_qtidalmin)
        call farray_index_append('qtidalmax',idiag_qtidalmax)
        call farray_index_append('qtidalm',idiag_qtidalm)
!
        call farray_index_append('dTdz1',idiag_dTdz1)
        call farray_index_append('dTdz2',idiag_dTdz2)
        call farray_index_append('dTdz3',idiag_dTdz3)
        call farray_index_append('dTdz4',idiag_dTdz4)
!
        call farray_index_append('nusselt_num',idiag_nusselt_num)
        call farray_index_append('nusselt_den',idiag_nusselt_den)
!
        call farray_index_append('TTmax_cline',idiag_TTmax_cline)
        call farray_index_append('TTmin_cline',idiag_TTmin_cline)
!
        call farray_index_append('devsigzz1',idiag_devsigzz1)
        call farray_index_append('devsigzz2',idiag_devsigzz2)
        call farray_index_append('devsigzz3',idiag_devsigzz3)
        call farray_index_append('devsigzz4',idiag_devsigzz4)
!
        call farray_index_append('icount',idiag_icount)
        call farray_index_append('residual',idiag_residual)
     endif
!
    endsubroutine rprint_special
!***********************************************************************
    subroutine get_slices_special(f,slices)
!
!  Write slices for animation of Special variables.
!
!  26-jun-06/tony: dummy
!
      use Slices_methods, only: assign_slices_scal
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
!  Loop over slices
!
      select case (trim(slices%name))
!
!  Streamfunction
!
        case('psi'); call assign_slices_scal(slices,f,ipsi)

      endselect
!
    endsubroutine get_slices_special
!***********************************************************************
!***********************************************************************
!***********************************************************************
!
!  Stuff for multigrid
!
!***********************************************************************
!***********************************************************************
!***********************************************************************    
    subroutine multigrid(f,psi,rhs,alpha_factor,beta_factor)
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,mz) :: psi
      real, dimension(nx,nz) :: rhs
      real, dimension(nx,nz) :: alpha_factor,beta_factor
!
      call full_multigrid(f,psi,rhs,alpha_factor,beta_factor,n_vcycles)
!
    endsubroutine multigrid
!***********************************************************************
!***********************************************************************
!***********************************************************************
!
!  Stuff for multigrid, adapted from numerical recipes. 
!
!***********************************************************************
!***********************************************************************
!***********************************************************************    
    subroutine full_multigrid(f,psi,r,alp,bet,ncycle)
!
!  Full Multigrid Algorithm for solution of linear elliptic equations.
!  On input u contains the right-hand side ρ in an N × N array,
!  while on output it returns the solution. The dimension N is related
!  to the number of grid levels used in the solution, ng below, by N = 2**ng+1.
!  ncycle is the number of V-cycles to be used at each level.
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,mz), intent(inout) :: psi
      real, dimension(nx,nz), intent(in) :: r
      real, dimension(nx,nz), intent(in) :: alp,bet
      integer, intent(in) :: ncycle
      integer :: j,jcycle,ng,ngrid
      integer :: mvx,mvz,nvx,nvz,nmin,nxc,nzc,mxc,mzc
      integer:: ll1,ll2,ll1_,ll2_,nn1,nn2,nn1_,nn2_
!
      real :: residuals
      integer :: icount
      integer :: nux,nuz
!
!  Define a type so we can have an array of pointers to arrays of grid variables.
!
      type ptr2d
         real, pointer :: rhs(:,:)
         real, pointer :: alpha(:,:)
         real, pointer :: beta(:,:)
      endtype ptr2d
      type(ptr2d), allocatable :: grid(:)
      real, dimension(:,:), pointer :: uj,uj_1
      real, dimension(:,:), allocatable :: uj_old
!
      mvx=size(psi,1); mvz=size(psi,2)
      nvx=size(r  ,1); nvz=size(r  ,2)
      nmin=min(nvx,nvz)

      ngrid=nint(log(nmin-1.0)/log(2.0))
      if (nmin /= 2**ngrid+1) call fatal_error("full_multigrid",&
           "nxgrid-1 or nzgrid-1 must be a power of 2 in full_multigrid")
!
      allocate(grid(ngrid))
!
      nxc=nvx; nzc=nvz 
      ng=ngrid
!
!  Allocate storage for r.h.s. on grid ng, and fill it with the input r.h.s.
!
      allocate(grid(ng)%rhs(nxc,nzc))
      grid(ng)%rhs=r

      allocate(grid(ng)%alpha(nxc,nzc))
      grid(ng)%alpha=alp

      allocate(grid(ng)%beta(nxc,nzc))
      grid(ng)%beta=bet
      
!
!  Similarly allocate storage and fill r.h.s. on all coarse grids by restricting from finer grids.
!
      do
         if ((nxc <= 3).or.(nzc <= 3)) exit
         nxc=nxc/2+1
         nzc=nzc/2+1
         ng=ng-1
!         
         allocate(grid(ng)%rhs(nxc,nzc))
         grid(ng)%rhs=restrict(grid(ng+1)%rhs)
!
         allocate(grid(ng)%alpha(nxc,nzc))
         grid(ng)%alpha=restrict(grid(ng+1)%alpha)
!         
         allocate(grid(ng)%beta(nxc,nzc))
         grid(ng)%beta=restrict(grid(ng+1)%beta)
      enddo
!
!  Initial solution on coarsest grid.
!
      !nn=3
      mxc=nxc+2*nghost     
      mzc=nzc+2*nghost
      allocate(uj(mxc,mzc))
      uj=0.0
      call solve_coarsest(uj,grid(1)%rhs,grid(1)%alpha,grid(1)%beta)
!
!  Nested iteration loop. 
!
      do j=2,ngrid
         ll1_=nghost+1; ll2_=mxc-nghost
         nn1_=nghost+1; nn2_=mzc-nghost
!         
         nxc=2*nxc-1; nzc=2*nzc-1
         mxc=nxc+2*nghost; mzc=nzc+2*nghost
!
         ll1=nghost+1; ll2=mxc-nghost
         nn1=nghost+1; nn2=mzc-nghost
!
         uj_1=>uj
         allocate(uj(mxc,mzc))
         uj=0.0
!
!  Prolongate (interpolate) from grid j-1 to next finer grid j. 
!
         uj(ll1:ll2,nn1:nn2)=prolongate(uj_1(ll1_:ll2_,nn1_:nn2_))
         deallocate(uj_1)
!
!  V-cycle loop. Iterate until convergence. 
!
         residuals=1e33
         icount=0
         do while (residuals > tolerance)
            nux=size(uj,1); nuz=size(uj,2)
            allocate(uj_old(nux,nuz))
            uj_old=uj
            !do jcycle=1,ncycle
            call vcycle(j,uj,grid(j)%rhs,grid(j)%alpha,grid(j)%beta)
            !enddo
            if (j==ngrid.and.lsave_residual_grid) then
              do m=1,my
                f(:,m,:,iresidual) = uj-uj_old
              enddo
            endif
            residuals=sqrt(sum((uj-uj_old)**2)/(nux*nuz))
            icount=icount+1
            if (lprint_residual) print*,icount,residuals,j
            if (lsave_residual) write(9,*) icount,residuals,j
            deallocate(uj_old)
            if (icount >= maxit) goto 100
            !call fatal_error("full multigrid",&
            !     "reached limit of iterations: did not converge.")
         enddo
100      continue 
      enddo
      if (lsave_residual) close(9)
      icount_save=icount
      residual_save=residuals
!      
!  Return solution in psi.
!
      psi=uj
!
!  Deallocate all grids
!
      deallocate(uj)
      do j=1,ngrid
         deallocate(grid(j)%rhs)
         deallocate(grid(j)%alpha)
         deallocate(grid(j)%beta)
!         
      enddo
      deallocate(grid)
!
      contains
!--------------------------------------------------------------------
      recursive subroutine vcycle(jgrid,u,rrhs,alpp,bett)
!
!  Recursive multigrid iteration. On input, jgrid is the current level,
!  u is the current value of the solution, and rhs is the right-hand
!  side. On output u contains the improved solution at the current level.
!  Parameters: npre and npost are the number of relaxation sweeps before
!  and after the coarse-grid correction is computed.
!
!  13-jan-16/wlad: from numerical recipes.
!        
        implicit none
        integer, intent(in) :: jgrid
        real, dimension(:,:), intent(inout) :: u
        real, dimension(:,:), intent(in) :: rrhs
!
        real, dimension(:,:), intent(in) :: alpp,bett
!        
        integer :: jpost,jpre
        integer :: mmvx,mmvz,mmux,mmuz
        integer :: lv1,lv2,lu1,lu2,nv1,nv2,nu1,nu2
!        
        integer :: icycle
        real, dimension((size(rrhs,1)+1)/2,(size(rrhs,2)+1)/2) :: res
        real, dimension((size(rrhs,1)+1)/2,(size(rrhs,2)+1)/2) :: alpp_,bett_
        real, dimension((size(rrhs,1)+1)/2+2*nghost,(size(rrhs,2)+1)/2+2*nghost) :: v
!        
        if (jgrid == 1) then
!
!  Bottom of V: Solve on coarsest grid.
!
           call solve_coarsest(u,rrhs,alpp,bett)
        else
!
!  On downward stroke of the V.
!
!  Pre-smoothing.              
!
           do jpre=1,npre
              call relaxation(u,rrhs,alpp,bett)
           enddo
!
!  Restriction of the residual is the next r.h.s.
!
           res=restrict(resid(u,rrhs,alpp,bett))
!
           alpp_=restrict(alpp)
           bett_=restrict(bett)
!
!  Zero for initial guess in next relaxation.
!
           v=0.0
!
!  Recursive call for the coarse grid correction.
!
           do icycle=1,merge(igamma,1,jgrid/=2)
              call vcycle(jgrid-1,v,res,alpp_,bett_)
           enddo
!
!  On upward stroke of V.
!
           mmvx=size(v,1); mmvz=size(v,2)
           lv1=nghost+1; lv2=mmvx-nghost
           nv1=nghost+1; nv2=mmvz-nghost
!           
           mmux=size(u,1); mmuz=size(u,2)
           lu1=nghost+1; lu2=mmux-nghost
           nu1=nghost+1; nu2=mmuz-nghost
!
           u(lu1:lu2,nu1:nu2)=u(lu1:lu2,nu1:nu2)+prolongate(v(lv1:lv2,nv1:nv2))
!
!  Post-smoothing.
!
           do jpost=1,npost
              call relaxation(u,rrhs,alpp,bett)
           enddo
        endif
!
      endsubroutine vcycle
!--------------------------------------------------------------------
    endsubroutine full_multigrid
!********************************************************************                
    function restrict(uf)
!
!  Half-weighting restriction. If Nc is the coarse-grid dimension, the
!  fine-grid solution is input in the (2Nc − 1) × (2Nc − 1) array uf,
!  the coarse-grid solution is returned in the Nc × Nc array restrict.
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(:,:), intent(in) :: uf
      real, dimension((size(uf,1)+1)/2,(size(uf,2)+1)/2) :: restrict
      integer :: ncx,ncz,nfx,nfz
!
      nfx=size(uf,1); nfz=size(uf,2)
      ncx=(nfx+1)/2; ncz=(nfz+1)/2
!
      !if (nfx == nfz) then 
!
!  Interior points
!      
         restrict(2:ncx-1,2:ncz-1) =   0.5 *uf(3:nfx-2:2,3:nfz-2:2) + &
                                     0.125*(uf(4:nfx-1:2,3:nfz-2:2) + uf(2:nfx-3:2,3:nfz-2:2)+ &
                                            uf(3:nfx-2:2,4:nfz-1:2) + uf(3:nfx-2:2,2:nfz-3:2))
!
!  Boundary points      
!
         restrict(1:ncx,1    ) = uf(1:nfx:2,1      )
         restrict(1:ncx,ncz  ) = uf(1:nfx:2,nfz    )
         restrict(1    ,1:ncz) = uf(1      ,1:nfz:2)
         restrict(ncx  ,1:ncz) = uf(nfx    ,1:nfz:2)
      !else
      !   restrict=1.
      !   print*,'restricted, go on'
      !endif
!      
    endfunction restrict
!********************************************************************
    function prolongate(uc)
!      
!  Coarse-to-fine prolongation by bilinear interpolation. If Nf is the
!  fine-grid dimension and Nc the coarse-grid dimension, then
!  Nf = 2Nc − 1. The coarse-grid solution is input as uc, the fine-grid
!  solution is returned in prolongate.
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(:,:), intent(in) :: uc
      real, dimension(2*size(uc,1)-1,2*size(uc,2)-1) :: prolongate
      integer :: ncx,ncz,nfx,nfz
!
      ncx=size(uc,1); ncz=size(uc,2)
      nfx=2*ncx-1; nfz=2*ncz-1

      !if (ncx == ncz) then 
!
! Do elements that are copies.
!
         prolongate(1:nfx:2,1:nfz:2)=uc(1:ncx,1:ncz)
!
! Do odd-numbered columns, interpolating vertically.
!
         prolongate(2:nfx-1:2,1:nfz:2) = 0.5*(prolongate(3:nfx:2,1:nfz:2) + prolongate(1:nfx-2:2,1:nfz:2))
!
! Do even-numbered columns, interpolating horizontally.
!
         prolongate(1:nfx,2:nfz-1:2) = 0.5*(prolongate(1:nfx,3:nfz:2) + prolongate(1:nfx,1:nfz-2:2))
      !else
      !   prolongate=1.
      !   print*,'prolongated, go on'
      !endif
!
    endfunction prolongate
!********************************************************************
    subroutine solve_coarsest(u,rhs,alp,bet)
!
!  Solution of the model problem on the coarsest grid. 
!  Input in rhs(1:3,1:3) and the solution is returned in u(1:3,1:3). 
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(:,:), intent(inout) :: u
      real, dimension(:,:), intent(in) :: rhs

      real, dimension(:,:), intent(in) :: alp,bet

      real, dimension(size(rhs,1),size(rhs,2)) :: u_old
      real :: h,hx,hz
      real :: alpha=impossible,beta=impossible
!
      real :: res,ufactor,vterm
      integer :: icount,iu
      integer :: nu,nux,nuz,nrx,nrz,lc1,lc2,nc1,nc2
!
      nux=size(u,1); nuz=size(u,2)
      nrx=size(rhs,1); nrz=size(rhs,2)
!
      u = 0.0
      hx = Lxyz(1)/(nrx-1)
      hz = Lxyz(3)/(nrz-1)
      if (hx==hz) then
        h = hx
      else
        print*,"hx=",hx,"hz=",hz
        call fatal_error("solve_coarsest","grid cells are not square")
      endif
!
      if (ldirect_solver) then
         if (nx_coarsest /= 3) call fatal_error("solve_coarsest",&
              "direct solver only for nx=3 in the coarsest grid")
         if (nux/=nuz) call fatal_error("solve_coarsest",&
              "direct solver only for square grids")
         nu=nux
         iu=(nu+1)/2
         if (lpoisson_test) then 
            u(iu,iu) = -h**2 * rhs(2,2)/4.0
         else
            alpha=alp(2,2)
            beta=bet(2,2)
            if (lsolver_highorder) then 
               call solve_highorder(u,alpha,beta,ufactor,vterm,iu,iu,h)
            else
               call solve_loworder(u,alpha,beta,ufactor,vterm,iu,iu,h)
            endif
            u(iu,iu) =  rhs(2,2)/ufactor 
         endif
      else
!
!  No multigrid, do SOR. 
!         
         lc1=nghost+1; lc2=nux-nghost
         nc1=nghost+1; nc2=nuz-nghost      
!      
         res=1e33
         icount=0
!
         do while (res > tolerance_coarsest)
            icount=icount+1
            u_old=u(lc1:lc2,nc1:nc2)
            !call successive_over_relaxation(u,rhs,alp,bet,h)
            
            call relaxation(u,rhs,alp,bet)
            res = sqrt(sum((u(lc1:lc2,nc1:nc2)-u_old)**2)/(nrx*nrz))
            if (ldebug) print*,'solve_coarsest:',icount,res,tolerance_coarsest
         enddo
!
      endif
!
    endsubroutine solve_coarsest
!********************************************************************
    subroutine relaxation(u,rhs,alp,bet)
!
!  Red-black Gauss-Seidel relaxation for model problem. The current
!  value of the solution u is updated, using the right-hand-side
!  function rhs. u and rhs are square arrays of the same odd dimension.
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(:,:), intent(inout) :: u
      real, dimension(:,:), intent(in) :: rhs

      real, dimension(:,:), intent(in) :: alp,bet

      integer :: l,k
      real :: h,h2
!
      real :: cc,ufactor,vterm
      real :: alpha=impossible,beta=impossible
!
      integer :: mvx,mvz,nvx,nvz
      integer :: ic,nc,icc1,icc2
      real :: hx,hz
!
      mvx=size(u  ,1); mvz=size(u  ,2)
      nvx=size(rhs,1); nvz=size(rhs,2)
!
      hx=Lxyz(1)/(nvx-1)
      hz=Lxyz(3)/(nvz-1)
      if (hx==hz) then
        h = hx
      else
        print*,"hx=",hx,"hz=",hz
        call fatal_error("relaxation","grid cells are not square")
      endif

!
! Convection
!
      if (lperi(1)) then
         icc1=1
         icc2=nvx
      else
         icc1=2
         icc2=nvx-1
      endif
!         
      if (.not.lpoisson_test) then 
         call update_bounds_psi(u)
         alpha=0.0
         beta=0.0
         do ic=icc1,icc2
            l=ic+nghost
            do nc=2,nvz-1
              k=nc+nghost
               cc = rhs(ic,nc)
               alpha=alp(ic,nc)
               beta=bet(ic,nc)
               if (lsolver_highorder) then
                  call solve_highorder(u,alpha,beta,ufactor,vterm,l,k,h)
               else
                  call solve_loworder(u,alpha,beta,ufactor,vterm,l,k,h)
               endif
               u(l,k) = (cc-vterm)/ufactor
            enddo
         enddo
         call update_bounds_psi(u)
      else
!
! Poisson
!     
         h2=h**2
!
         do ic=2,nvx-1
            l=ic+nghost
            do nc=2,nvz-1
               k=nc+nghost
               u(l,k) = 0.25*(u(l+1,k) + u(l-1,k) + &
                              u(l,k+1) + u(l,k-1) - h2*rhs(ic,nc))
            enddo
         enddo
      endif
!
    endsubroutine relaxation
!********************************************************************
    function resid(u,rhs,alp,bet)
!
!  Returns minus the residual for the model problem. Input quantities are u and rhs,
!  while the residual is returned in resid. All three quantities are square arrays
!  with the same odd dimension.
!
!  13-jan-16/wlad: from numerical recipes.
!
      implicit none
      real, dimension(:,:), intent(in) :: u
      real, dimension(:,:), intent(in) :: rhs

      real, dimension(:,:), intent(in) :: alp,bet

      real, dimension(size(rhs,1),size(rhs,2)) :: resid
      integer :: l,k
      real :: h,h2i
!
      real :: ufactor,vterm,lhs
      real :: alpha=impossible,beta=impossible
!
      integer :: mvx,mvz,nvx,nvz
      integer :: ic,nc
      real :: hx,hz
!
      mvx=size(u  ,1); mvz=size(u  ,2)
      nvx=size(rhs,1); nvz=size(rhs,2)
!
      hx=Lxyz(1)/(nvx-1)
      hz=Lxyz(3)/(nvz-1)
      if (hx==hz) then
        h = hx
      else
        print*,"hx=",hx,"hz=",hz
        call fatal_error("resid","grid cells are not square")
      endif
!
!  Interior points.
!
      if (lpoisson_test) then 
         h2i=1.0/(h**2)
         do ic=2,nvx-1
            l=ic+nghost
            do nc=2,nvz-1
               k=nc+nghost
               resid(ic,nc) = -h2i*(u(l+1,k  ) + u(l-1,k  ) + &
                                    u(l  ,k+1) + u(l  ,k-1) - &
                                    4.0*u(l,k))+rhs(ic,nc)
            enddo
         enddo
      else
         do ic=2,nvx-1
            l=ic+nghost
            do nc=2,nvz-1
               k=nc+nghost
               alpha=alp(ic,nc)
               beta=bet(ic,nc)
               if (lsolver_highorder) then
                  call solve_highorder(u,alpha,beta,ufactor,vterm,l,k,h)
               else
                  call solve_loworder(u,alpha,beta,ufactor,vterm,l,k,h)
               endif
               lhs = ufactor*u(l,k) + vterm
               resid(ic,nc) = rhs(ic,nc) - lhs
            enddo
         enddo
      endif
!
!  Boundary points.
!      
      resid(1:nvx,1    )=0.0
      resid(1:nvx,nvz  )=0.0
      resid(1    ,1:nvz)=0.0
      resid(nvx  ,1:nvz)=0.0
!
    endfunction resid
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from nospecial.f90 for any Special      **
!**  routines not implemented in this file                         **
!**                                                                **
    include '../special_dummies.inc'
!********************************************************************
endmodule Special
