!--------------------------------------------------------------------------
   !---------------------------- compute_hydrostatic_pressure ----------------------------------
   SELECT CASE(caldyn_thermo)
   CASE(thermo_boussinesq)
      ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
      !$OMP DO SCHEDULE(STATIC)
      DO ij=1,primal_num
         pk(llm,ij) = ptop + .5*g* theta_rhodz(llm,ij,1)
         DO l = llm-1,1,-1
            pk(l,ij) = pk(l+1,ij) + (.5*g)*( theta_rhodz(l,ij,1) + theta_rhodz(l+1,ij,1) )
         END DO
         IF(caldyn_eta == eta_lag) THEN
            ps(ij) = pk(1,ij) + .5*g* theta_rhodz(1,ij,1)
         END IF
      END DO
      !$OMP END DO
   CASE(thermo_theta, thermo_entropy, thermo_variable_Cp)
      !$OMP DO SCHEDULE(STATIC)
      DO ij=1,primal_num
         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
         DO l = llm-1,1,-1
            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
         END DO
         IF(caldyn_eta == eta_lag) THEN
            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
         END IF
      END DO
      !$OMP END DO
   CASE(thermo_moist)
      !$OMP DO SCHEDULE(STATIC)
      DO ij=1,primal_num
         pk(llm,ij) = ptop + .5*g* (rhodz(llm,ij)+theta_rhodz(llm,ij,2))
         DO l = llm-1,1,-1
            pk(l,ij) = pk(l+1,ij) + (.5*g)*( (rhodz(l,ij)+theta_rhodz(l,ij,2)) + (rhodz(l+1,ij)+theta_rhodz(l+1,ij,2)) )
         END DO
         IF(caldyn_eta == eta_lag) THEN
            ps(ij) = pk(1,ij) + .5*g* (rhodz(1,ij)+theta_rhodz(1,ij,2))
         END IF
      END DO
      !$OMP END DO
   END SELECT
   !---------------------------- compute_hydrostatic_pressure ----------------------------------
   !--------------------------------------------------------------------------