anchor.diff:      allocate(theta_error%f(1:theta%n), stat=ist)
   
       !Vapour
  -    Qv%n=mesh%nv
  -    Qv%name="vapour"
  -    Qv%pos=0
  -    allocate(Qv%f(1:Qv%n), stat=ist)
  -    allocate(hQv_old%f(1:Qv%n), stat=ist)
  -    allocate(hQv%f(1:Qv%n), stat=ist)
  -    allocate(hSQv%f(1:Qv%n), stat=ist)
  -    allocate(delta_Qv%f(1:Qv%n), stat=ist)
  -    allocate(Qv_exact%f(1:Qv%n), stat=ist)
  -    allocate(Qv_error%f(1:Qv%n), stat=ist)
  +    qv%n=mesh%nv
  +    qv%name="vapour"
  +    qv%pos=0
  +    allocate(qv%f(1:qv%n), stat=ist)
  +    allocate(hqv_old%f(1:qv%n), stat=ist)
  +    allocate(hqv%f(1:qv%n), stat=ist)
  +    allocate(hSqv%f(1:qv%n), stat=ist)
  +    allocate(delta_qv%f(1:qv%n), stat=ist)
  +    allocate(qv_exact%f(1:qv%n), stat=ist)
  +    allocate(qv_error%f(1:qv%n), stat=ist)
   
       !Cloud
  -    Qc%n=mesh%nv
  -    Qc%name="cloud"
  -    Qc%pos=0
  -    allocate(Qc%f(1:Qc%n), stat=ist)
  -    allocate(hQc%f(1:Qc%n), stat=ist)
  -    allocate(hQc_old%f(1:Qc%n), stat=ist)
  -    allocate(hSQc%f(1:Qc%n), stat=ist)
  -    allocate(delta_Qc%f(1:Qc%n), stat=ist)
  -    allocate(Qc_exact%f(1:Qc%n), stat=ist)
  -    allocate(Qc_error%f(1:Qc%n), stat=ist)
  +    qc%n=mesh%nv
  +    qc%name="cloud"
  +    qc%pos=0
  +    allocate(qc%f(1:qc%n), stat=ist)
  +    allocate(hqc%f(1:qc%n), stat=ist)
  +    allocate(hqc_old%f(1:qc%n), stat=ist)
  +    allocate(hSqc%f(1:qc%n), stat=ist)
  +    allocate(delta_qc%f(1:qc%n), stat=ist)
  +    allocate(qc_exact%f(1:qc%n), stat=ist)
  +    allocate(qc_error%f(1:qc%n), stat=ist)
   
       !Rain
  -    Qr%n=mesh%nv
  -    Qr%name="rain"
  -    Qr%pos=0
  -    allocate(Qr%f(1:Qr%n), stat=ist)
  -    allocate(hQr%f(1:Qr%n), stat=ist)
  -    allocate(hQr_old%f(1:Qr%n), stat=ist)
  -    allocate(hSQr%f(1:Qr%n), stat=ist)
  -    allocate(delta_Qr%f(1:Qr%n), stat=ist)
  -    allocate(Qr_exact%f(1:Qr%n), stat=ist)
  -    allocate(Qr_error%f(1:Qr%n), stat=ist)
  +    qr%n=mesh%nv
  +    qr%name="rain"
  +    qr%pos=0
  +    allocate(qr%f(1:qr%n), stat=ist)
  +    allocate(hqr%f(1:qr%n), stat=ist)
  +    allocate(hqr_old%f(1:qr%n), stat=ist)
  +    allocate(hSqr%f(1:qr%n), stat=ist)
  +    allocate(delta_qr%f(1:qr%n), stat=ist)
  +    allocate(qr_exact%f(1:qr%n), stat=ist)
  +    allocate(qr_error%f(1:qr%n), stat=ist)
   
       !Tracer
       tracer%n=mesh%nv
  
positive.code: 
          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
  
        q0 = 0.02_r8/maxval(qv%f)
        Qv%f = q0*Qv%f
        hQv%f = q0*hQv%f
  
        !Velocity
        if(useStagHTC)then
  
negative.code:     type(scalar_field), intent(in):: htheta  !General
  
      !Vapour (defined on voronoi centers)
      type(scalar_field), intent(in):: hqv  !General
  
      !Cloud (defined on voronoi centers)
      type(scalar_field), intent(in):: hqc  !General
  
      !Rain (defined on voronoi centers)
      type(scalar_field), intent(in):: hqr  !General
  
      !Tracer (defined on voronoi centers)
      type(scalar_field), intent(in):: htracer  !General
  