anchor.diff:      call scalar_elem_product(u, hqc_ed, uhqc) !Flux uhqc at edges
   
       !Calculate divergence / cloud eq RHS  
  -    !call div_hx(uhqc, div_uhqc, mesh)
  -    call divhx(u, hqc, hqc_ed, uhqc, div_uhqc, mesh)
  +    call div_hx(uhqc, div_uhqc, mesh)
  +    !call divhx(u, hqc, hqc_ed, uhqc, div_uhqc, mesh)
   
       !Cloud eq. RHS
       cloudeq = -div_uhqc%f
  
positive.code:         print '(a22, 3e16.8)',' rain = ',minval(qr%f),maxval(qr%f), Train
        end if
  
        print '(a22, 2e16.8)',' rain = ',minval(htracer%f),maxval(htracer%f)
        if (rel_error_h>1000.d0) stop
        !print*,'CFL = ',cfl
        !Plot fields
  
negative.code:       Train=sumf_areas(qr)
        Tcloud=sumf_areas(qc)
        Tvapour=sumf_areas(qv)
  
        if(testcase==2)then
          h_error%f = h_exact%f - h%f
          u_error%f = u_exact%f - u%f
  