anchor.diff:      !$omp end parallel do
     end subroutine source
   
  -  subroutine swm_physics_tests(meshtmp)
  +  subroutine moist_swm_tests(meshtmp)
       !-----------------------------------------
       !  Main test routine tests routine
       !-----------------------------------------
  
positive.code:           u%f(l)=dot_product(vectmp,mesh%edhx(l)%nr)
          end do
        end if
      !======================================================================================
      ! Flow over Andes mountain
      !======================================================================================
  
      case(4,5,6)
        !Parameters
        u0     = 20._r8
        phi0   = 5960._r8*grav
        h0     = 2000._r8
        phi_ct = (erad*omega*u0+(u0**2)/2._r8)
        w      = omega*erad*u0 + u0*u0*0.5_r8
        sigma  = 0._r8!w/10._r8
        xsi    = 0.001_r8
        mu1    = 0.05_r8
        mu2    = 0.98_r8
        theta_sp = -40._r8/300._r8 
        theta_eq =  30._r8/300._r8 
        theta_np = -20._r8/300._r8 
        
        if(testcase==4)then  
          lon0  = -pi*0.5_r8
          lat0  = pi/6._r8
        else if(testcase==5)then  
          lon0  = -pi*0.5_r8 - 0.35d0*pi
          lat0  = -pi/6._r8
        else if(testcase==6)then  
          lon0  = -pi*0.5_r8 -0.52777d0*pi
          lon1 = lon0 + 2.d0*pi
          lat0  = -pi/6._r8  
        end if
        rmax   = pi/9._r8
  
        !Variables at Voronoi centers
        do i=1, mesh%nv
          lon    = mesh%v(i)%lon
          lat    = mesh%v(i)%lat
          !print*,lon,lat
          coss = cos(lat)
          sinn = sin(lat) 
          if(testcase==6)then  
            !lon = lon-pi! +0.25d0*pi
          end if
          r = dsqrt((lon-lon0)**2+(lat-lat0)**2)
          r1 = dsqrt((lon-lon1)**2+(lat-lat0)**2)
          h%f(i) = (phi0 -w*sinn**2)*gravi
  
          !Topography
          !if(testcase==6)then !tc 6
          !  bt%f(i)=2000._r8*smooth_andean_mountain([lat,lon], alt_table)
          !else !tc 4 or tc5
            if(r<rmax)then
              bt%f(i)=2000._r8*(1._r8-r/rmax)
            else if(r1<rmax)then
              bt%f(i)=2000._r8*(1._r8-r1/rmax)
            else
              bt%f(i)=0.
            endif
          !endif
          ! Correct h to allow orography
          h%f(i)=h%f(i)-bt%f(i)
  
          !Temperature
          !lon = lon+pi + 0.35d0*pi
          if(testcase==6)then  
            lon = lon+1.d0*pi + 0.52777d0*pi !-95graus
          else
             lon = lon+pi + 0.35d0*pi
          end if   
          !lon = lon + 0.2*pi
          !theta%f(i)=F_quad(theta_sp,(1._r8-mu1)*theta_eq,theta_np,lat) + mu1*theta_eq*dcos(lat)*dsin(lon)
          !theta%f(i)=F_quad(theta_np,(1._r8-mu1)*theta_eq,theta_sp,lat) + mu1*theta_eq*dcos(lat)*dsin(lon)
          theta%f(i)= mu1*theta_eq*dcos(lat)*dsin(lon)
          !Vapour
          Qv%f(i) = mu2*qsat(theta%f(i),h%f(i),bt%f(i),1._r8)
          !Qv%f(i) = 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
  
        q0 = 0.02_r8/maxval(qv%f)
        Qv%f = q0*Qv%f
        hQv%f = q0*hQv%f
  
        print*,maxval(bt%f),minval(bt%f)
       ! stop
  
        if(useStagHTC)then
          do l=1, mesh%ne
            utmp=u0*dcos(mesh%ed(l)%c%lat)
            vtmp=0._r8
            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
            utmp=u0*dcos(mesh%edhx(l)%c%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
  
  
  
      !======================================================================================
      ! galewski
      !======================================================================================
  
      case(7,8)
        !Parameters
        !mu1    = 0.05_r8
  
negative.code:     lond = 0._r8*deg2rad
  
      select case(testcase)
      !======================================================================================
      ! Steady state - from Zerroukat and Allen JCP 2015
      !======================================================================================
      case(2)
        !Field at cell's center
        u0    = 20._r8
  