anchor.diff:        call ode_rk4_moist_swm(time, h_old, u_old, htheta_old, hQv_old, hQc_old, hQr_old, &
                            h, u, htheta, hQv, hQc, hQr, dt)
   
  -      call mass_fixer2(hQv)             
  -      call mass_fixer2(hQr)
  -      call mass_fixer2(hQc)
  +      !Apply the monotonic filter for tracers
  +      call monotonic_filter(hQv)             
  +      call monotonic_filter(hQr)
  +      call monotonic_filter(hQc)
   
         !call mass_fixer(hQv,hQr,hQc, iniwater)
         !call scalar_elem_divide(htheta, h, theta)
  
positive.code: 
        !Calculate erngies
        call calc_energies(Penergy, Kenergy, Tenergy, Availenergy)
        print*,iniwater,Twater,(Twater-iniwater)/iniwater
  
        !update fields
        u_old=u
        h_old=h
  
negative.code:         rel_error_Qr = maxval(abs(Qr_error%f))
          rel_error_theta = maxval(abs(theta_error%f))/maxval(abs(theta_exact%f))
          print*, k, ntime
          print*, "Time (dys) :",   k*dt*sec2day, " of ", ntime*dt*sec2day
          print*, "Step = ", k, " of ", ntime
          print '(a33, 3e16.8)','linf errors of (h, u, theta) = ',rel_error_h,rel_error_u,rel_error_theta
          print '(a33, 3e16.8)','linf errors of  (qv, qc, qr) = ',rel_error_Qv,rel_error_Qc,rel_error_Qr
          !print '(a22, 2e16.8)',' height = ',minval(h%f),maxval(h%f)
          !print '(a22, 2e16.8)',' velocity = ',minval(u%f),maxval(u%f)
          !print '(a22, 2e16.8)',' temperature = ',minval(theta%f),maxval(theta%f)
          !print '(a22, 2e16.8)',' vapour = ',minval(qv%f),maxval(qv%f)
          !print '(a22, 2e16.8)',' cloud = ',minval(qc%f),maxval(qc%f)
          !print '(a22, 2e16.8)',' rain = ',minval(qr%f),maxval(qr%f)
          call write_swmp_error_file(time)
        else 
          print*, "Time (dys) :",   k*dt*sec2day, " of ", ntime*dt*sec2day
          print*, "Step = ", k, " of ", ntime
          print*,'                          min               max'
          print '(a22, 2e16.8)',' height = ',minval(h%f),maxval(h%f)
          print '(a22, 2e16.8)',' velocity = ',minval(u%f),maxval(u%f)
          print '(a22, 2e16.8)',' temperature = ',minval(theta%f),maxval(theta%f)
          print '(a22, 2e16.8)',' vapour = ',minval(qv%f),maxval(qv%f)
          print '(a22, 2e16.8)',' cloud = ',minval(qc%f),maxval(qc%f)
          print '(a22, 2e16.8)',' rain = ',minval(qr%f),maxval(qr%f)
        end if
  
        !print*,'CFL = ',cfl
        !Plot fields
        call plotfields_conv(k, time)
  