!--------------------------------------------------------------------------
   !---------------------------- 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
      !$OMP DO SCHEDULE(STATIC)
      DO ij=1,primal_num
         l=0
         mass_col(ij)=0.
         DO l = 1,llm
            mass_col(ij) = mass_col(ij) + rhodz(l,ij)
         END DO
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = 1, llm
            ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
            ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij)
            ! rhodz(l,ij) = m/g
            rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij)
         END DO
      END DO
      !$OMP END DO
   END IF
   DO iq=1,nqdyn
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = 1, llm
            theta(l,ij,iq) = theta_rhodz(l,ij,iq)/rhodz(l,ij)
         END DO
      END DO
      !$OMP END DO
   END DO
   !---------------------------- theta ----------------------------------
   !--------------------------------------------------------------------------