anchor.diff:        Qc_exact = Qc
         Qr_exact = Qr
   
  -
  -    !-------------------------------------------------------------------------------  
  -    ! Steady state - rotated
  -    !-------------------------------------------------------------------------------      
  -    case(3)
  -      !Field at cell's center
  -      u0    = 20._r8
  -      phi0  = 3._r8*10**4 
  -      w     = omega*erad*u0 + u0*u0*0.5_r8
  -      sigma = w/10._r8
  -      temp0 = phi0*phi0/300._r8
  -      !xsi = 0.01_r8
  -      xsi = 0.000_r8
  -
  -      ! Calculate rotation matrix
  -      lat0 = -45._r8*deg2rad
  -      lon0 = 0._r8!-90._r8*deg2rad
  -      Rmat(1:3,1:3)=0._r8
  -      Rmat(1,1)=dcos(lat0)*dcos(lon0)
  -      Rmat(1,2)=-dsin(lon0)*dcos(lat0)
  -      Rmat(1,3)=-dsin(lat0)
  -      Rmat(2,1)=dsin(lon0)
  -      Rmat(2,2)=dcos(lon0)
  -      Rmat(3,1)=dsin(lat0)*dcos(lon0)
  -      Rmat(3,2)=-dsin(lat0)*dsin(lon0)
  -      Rmat(3,3)=dcos(lat0)
  -        
  -      !Inverse rotation
  -      RmatT=transpose(Rmat)
  -      alpha =-45._r8*deg2rad
  -
  -      fsphere=3
  -
  -      do i=1, mesh%nv
  -        lon = mesh%v(i)%lon
  -        lat = mesh%v(i)%lat
  -        sinn = -dcos(lon)*dcos(lat)*dsin(alpha) + dsin(lat)*dcos(alpha)
  -        coss =  dcos(lat)*dcos(alpha) + dcos(lon)*dsin(lat)*dsin(alpha)
  -
  -        vectmp = mesh%v(i)%p
  -        vectmp = matmul(RmatT,vectmp)
  -        call cart2sph (vectmp(1), vectmp(2), vectmp(3), lon, lat )
  -
  -        sinn = sin(lat)
  -        coss = cos(lat)
  -
  -        h%f(i) = (phi0 -(w+sigma)*sinn**2)*gravi
  -        bt%f(i) = 0._r8 
  -
  -        theta%f(i) = temp0 + sigma*(coss**2)*((w+sigma)*(coss**2) + (phi0-w-sigma)*2._r8 )
  -        theta%f(i) = theta%f(i)/( phi0**2 + ((w+sigma)**2)*sinn**4 -2._r8*phi0*(w+sigma)*sinn**2)
  -        Qv%f(i) = (1._r8-xsi)*qsat(theta%f(i),h%f(i),bt%f(i),1._r8)
  -
  -        !Fluxes
  -        hQv%f(i) = h%f(i)*Qv%f(i)
  -        hQc%f(i) = h%f(i)*Qc%f(i)
  -        hQr%f(i) = h%f(i)*Qr%f(i)
  -        htheta%f(i) = h%f(i)*theta%f(i)
  -      end do
  -
  -      print*,minval(theta%f),maxval(theta%f)
  -      !stop
  -      q0 = 0.02_r8/maxval(qv%f)
  -      print*,maxval(qv%f)
  -      Qv%f = q0*Qv%f
  -      hQv%f = q0*hQv%f
  -
  -      !Velocity
  -      if(useStagHTC)then
  -        do l=1, mesh%ne
  -          lat = mesh%ed(l)%c%lat
  -          lon = mesh%ed(l)%c%lon
  -          utmp=u0*(cos(lat)*cos(alpha) + cos(lon)*sin(lat)*sin(alpha))
  -          vtmp=-u0*(sin(lon)*sin(alpha))
  -          call convert_vec_sph2cart(utmp, vtmp, mesh%ed(l)%c%p, vectmp)
  -          v_ed%p(l)%v=vectmp
  -          u%f(l)=dot_product(vectmp,mesh%ed(l)%tg)
  -        end do
  -
  -      elseif(useStagHC)then
  -        do l=1, mesh%ne
  -          !vectmp=mesh%edhx(l)%c%p
  -          !vectmp = matmul(RmatT,vectmp)
  -          !call cart2sph(vectmp(1),vectmp(2),vectmp(3),lon,lat)
  -          !utmp = u0*dcos(lat)
  -          !vtmp = 0._r8
  -          !call convert_vec_sph2cart(utmp, vtmp, mesh%edhx(l)%c%p, vectmp)
  -          !v_ed%p(l)%v=vectmp
  -          !u%f(l)=dot_product(vectmp,mesh%edhx(l)%nr)
  -       end do
  -      end if
  -
  -      h_exact = h
  -      u_exact = u
  -      theta_exact = theta
  -      Qv_exact = Qv
  -      Qc_exact = Qc
  -      Qr_exact = Qr
  -
       !======================================================================================
  -    ! Flow over Andes mountain
  +    ! Flow over a mountain - from Zerroukat and Allen JCP 2015
       !======================================================================================
  -    case(10)
  -    ! Calculate rotation matrix
  -    lat0 = 0._r8*deg2rad
  -    lon0 = 0._r8!-90._r8*deg2rad
  -    Rmat(1:3,1:3)=0._r8
  -    Rmat(1,1)=dcos(lat0)*dcos(lon0)
  -    Rmat(1,2)=-dsin(lon0)*dcos(lat0)
  -    Rmat(1,3)=-dsin(lat0)
  -    Rmat(2,1)=dsin(lon0)
  -    Rmat(2,2)=dcos(lon0)
  -    Rmat(3,1)=dsin(lat0)*dcos(lon0)
  -    Rmat(3,2)=-dsin(lat0)*dsin(lon0)
  -    Rmat(3,3)=dcos(lat0)
  -    
  -    Rmat(1,1)=1.d0
  -    Rmat(1,2)=0.d0
  -    Rmat(1,3)=0.d0
  -    Rmat(2,1)=0.d0
  -    Rmat(2,2)=dcos(lat0)
  -    Rmat(2,3)=-dsin(lat0)
  -    Rmat(3,1)=0.d0
  -    Rmat(3,2)=dsin(lat0)
  -    Rmat(3,3)=dcos(lat0)
  -    
  -    lat0 = 0._r8*deg2rad
  -    Rmat2(1,1)=dcos(lat0)
  -    Rmat2(1,2)=-dsin(lat0)
  -    Rmat2(1,3)=0.d0
  -    Rmat2(2,1)=dsin(lat0)
  -    Rmat2(2,2)=dcos(lat0)
  -    Rmat2(2,3)=0.d0
  -    Rmat2(3,1)=0.d0
  -    Rmat2(3,2)=0.d0
  -    Rmat2(3,3)=1.d0
  -    
  -    lat0 = 60._r8*deg2rad
  -    Rmat(1,1)=dcos(lat0)
  -    Rmat(1,2)=0.d0
  -    Rmat(1,3)=-dsin(lat0)
  -    Rmat(2,1)=0.d0
  -    Rmat(2,2)=1.d0
  -    Rmat(2,3)=0.d0
  -    Rmat(3,1)=dsin(lat0)
  -    Rmat(3,2)=0.d0
  -    Rmat(3,3)=dcos(lat0)   
  -     
  -    !Inverse rotation
  -    RmatT=matmul(transpose(Rmat),transpose(Rmat2))
  -    !RmatT=transpose(Rmat)
  -    alpha =lat0
  -
  -    !fsphere=3
  +    case(4)
   
       !Parameters
       u0     = 20._r8
  
positive.code:       hQc_old=hQc
        hQr_old=hQr
      end do
    end  subroutine swm_physics_tests
  
  
    subroutine write_evol_file_cswm(time, iniwater, Twater, inimass, Penergy0, Kenergy0, Tenergy0, Availenergy0,&
  
negative.code:     !Variables at Voronoi centers
      do i=1, mesh%nv
          vectmp = mesh%v(i)%p
          lon = mesh%v(i)%lon
          lat = mesh%v(i)%lat
          sinn = sin(lat)
  