anchor.diff:  
         ! Integrate to tabulate h and psi as functions of geographical
         ! latitude
  -      nygg = 4*FLOOR(SQRT(REAL(mesh%nv)))
  +      nygg = 4*FLOOR(SqrT(REAL(mesh%nv)))
         allocate(hgg(nygg+1)) !, psigg(nygg+1))
         dygg = pi/nygg
         hgg(1) = 0.0D0
  
positive.code:     type(scalar_field) :: phi_mass
      real(kind=8):: mass_water_initial
      integer(i4) :: i,j,k, nnb, iter
      real(r8):: mass_hQc, mass_hQv, mass_hQr, modified_mass_hQr, sumareas, summass,summass2
      real(r8):: eps, eps2
  
      phi_mass= phi
  
negative.code:       gamma_v = 1._r8/gamma_v  
  
        !Vapour
        delta_qv%f(i) = gamma_v*(qv%f(i)-qsat(theta%f(i),h%f(i), bt%f(i), q0))
        delta_qv%f(i) = max(0._r8,delta_qv%f(i))/dtime
  
        !Cloud
        delta_qc%f(i) = max(0._r8, -gamma_v*(qv%f(i)-qsat(theta%f(i),h%f(i), bt%f(i), q0)))
        delta_qc%f(i) = min(qc%f(i),delta_qc%f(i))/dtime
  
        !Rain
        delta_qr%f(i) = max(0._r8,gamma_r*(qc%f(i)-q_precip))/dtime
  
        !Computes the sources
        hSqv%f(i) = delta_qc%f(i) - delta_qv%f(i)
        hSqc%f(i) = delta_qv%f(i) - delta_qc%f(i) - delta_qr%f(i)
        hSqr%f(i) = delta_qr%f(i)
        hStheta%f(i) = Lscale*(delta_qv%f(i)-delta_qc%f(i))
  
        hSqv%f(i) = h%f(i)*hSqv%f(i)
        hSqc%f(i) = h%f(i)*hSqc%f(i)
        hSqr%f(i) = h%f(i)*hSqr%f(i)
        hStheta%f(i) =  h%f(i)*hStheta%f(i)
      end do
  
  