!--------------------------------------------------------------------------
   !---------------------------- theta ----------------------------------
   IF(caldyn_eta==eta_mass) THEN ! Compute mass
      ! FIXME : here mass_col is computed from rhodz
      ! so that the DOFs are the same whatever caldyn_eta
      ! in DYNAMICO mass_col is prognosed rather than rhodz
      DO ij=ij_omp_begin_ext,ij_omp_end_ext
         mass_col(ij)=0.
      END DO
      DO l = 1,llm
         DO ij=ij_omp_begin_ext,ij_omp_end_ext
            mass_col(ij) = mass_col(ij) + rhodz(ij,l)
         END DO
      END DO
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij=ij_begin_ext, ij_end_ext
            ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
            ! m = mass_dak(l)+(mass_col(ij)*g+ptop)*mass_dbk(l)
            ! rhodz(ij,l) = m/g
            rhodz(ij,l) = mass_dak(l) + mass_col(ij)*mass_dbk(l)
         END DO
      END DO
   END IF
   DO iq=1,nqdyn
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij=ij_begin_ext, ij_end_ext
            theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l)
         END DO
      END DO
   END DO
   !---------------------------- theta ----------------------------------
   !--------------------------------------------------------------------------