anchor.diff:      do i=1, mesh%nv
          Penergy=Penergy+mesh%hx(i)%areag*grav*h%f(i)*(h%f(i)*0.5+bt%f(i))
       end do
  -	!rint*, "Potential:",Penergy
  -
  -	UnavailPenergy=0.
  -	!Calculate mean height
  -	mean_height=sumf_areas(hbt)
  +    !rint*, "Potential:",Penergy
  +
  +    UnavailPenergy=0.
  +    !Calculate mean height
  +    mean_height=sumf_areas(hbt)
       !print*, mean_height, maxval(h%f(1:mesh%nv)), mean_height_bt, maxval(bt%f(1:mesh%nv))
       do i=1, mesh%nv
  -    	mean_height_bt=mean_height-bt%f(i)
  -
  -    	!Unavailable energy
  +       mean_height_bt=mean_height-bt%f(i)
  +
  +       !Unavailable energy
          UnavailPenergy=UnavailPenergy+mesh%hx(i)%areag*grav*mean_height_bt*(mean_height_bt*0.5+bt%f(i))
       end do
  -	!print*, "Unavailable:", UnavailPenergy
  +    !print*, "Unavailable:", UnavailPenergy
   
       Kenergy=0.
  -    !if(useReconmtdPerhx)then
  -    ! do i=1, mesh%nv
  -    !   Kenergy=Kenergy+Kin_energy%f(i)*mesh%hx(i)%areag
  -    ! end do
  -    !else
  -
  -    do l=1, mesh%ne
  -       Kenergy=Kenergy+0.5*mesh%ed(l)%leng*mesh%edhx(l)%leng*h_ed%f(l)*u%f(l)**2
  -    end do
  -    !end if
  +    if(useReconmtdPerhx)then
  +       do i=1, mesh%nv
  +          Kenergy=Kenergy+Kin_energy%f(i)*mesh%hx(i)%areag
  +       end do
  +    else       
  +       do l=1, mesh%ne
  +          Kenergy=Kenergy+0.5*mesh%ed(l)%leng*mesh%edhx(l)%leng*h_ed%f(l)*u%f(l)**2
  +       end do
  +    end if
   
       Tenergy=Penergy+Kenergy
  -	Availenergy=Tenergy-UnavailPenergy
  -	!print*, "Available:", Availenergy
  +    Availenergy=Tenergy-UnavailPenergy
  +    !print*, "Available:", Availenergy
   
       return
     end subroutine calc_energies
  
positive.code:           stop
         end if
  
       ! For tc33, h is contant, u is rotation and the f=0
         if(testcase==33)then
            fcte=0.
            fsphere=2
  		  !u0=200
            h_ct=((u0**2)/2._r8)*gravi
            eta_ct=(2.*u0/erad)
         else
         	h_ct=(erad*omega*u0+(u0**2)/2._r8)*gravi
         	eta_ct=(2.*u0/erad+2.*omega)
         end if
  
         grad_ct=u0*eta_ct
  	   print*, "Test case parameters"
  
  	   print*, "u0", u0
  	   print*, "h_ct", h_ct
  
  
         do i=1, mesh%nv
            !print*, i
            h%f(i)=hollgw
            !if(testcase==32)then !Only put bottom topograpy in tc32 - geostrophic balanced flow
               bt%f(i)=h0-h_ct*dsin(mesh%v(i)%lat)**2
            !elseif(testcase==33)then
            !   bt%f(i)=h0-((u0**2)/2._r8)*gravi*dsin(mesh%v(i)%lat)**2
            !end if
  
negative.code:     !File for errors
      filename=trim(datadir)//trim(swmname)//"_mass_energy_"//trim(mesh%name)//".txt"
      buffer="        n        mvdist  tcase       k       nt        dt(s)       "//&
           "          cfl                 time(dys)             mass               Penergy    "//&
           "          Kenergy            Tenergy            Availenergy           RMSdiv"
  
      inquire(file=filename, opened=iopen)
      if(iopen)then
  