--- trunk/Sources/phylmd/CV30_routines/cv30_prelim.f	2016/06/02 15:40:30	200
+++ trunk/Sources/phylmd/CV30_routines/cv30_prelim.f	2016/06/06 17:42:15	201
@@ -7,24 +7,32 @@
   SUBROUTINE cv30_prelim(t1, q1, p1, ph1, lv1, cpn1, tv1, gz1, h1, hm1, th1)
 
     USE cv30_param_m, ONLY: nl
-    USE cv_thermo_m, ONLY: cl, clmcpv, cpd, cpv, eps, rrd, rrv
+    USE cv_thermo_m, ONLY: clmcpv, eps
     USE dimphy, ONLY: klev, klon
-    use SUPHEC_M, only: rlvtt
+    use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rd, rv
 
     ! Calculate arrays of geopotential, heat capacity and static energy
 
-    real, intent(in):: t1(klon, klev)
-    real, intent(in):: q1(klon, klev)
-    real, intent(in):: p1(klon, klev), ph1(klon, klev + 1)
+    real, intent(in):: t1(:, :) ! (klon, klev) temperature, in K
+    real, intent(in):: q1(:, :) ! (klon, klev) specific humidity
+    real, intent(in):: p1(:, :) ! (klon, klev) full level pressure, in hPa
+    real, intent(in):: ph1(:, :) ! (klon, klev + 1) half level pressure, in hPa
 
     ! outputs:
-    real lv1(klon, klev), cpn1(klon, klev), tv1(klon, klev)
+
+    real, intent(out):: lv1(:, :) ! (klon, nl)
+    ! specific latent heat of vaporization of water, in J kg-1
+    
+    real, intent(out):: cpn1(:, :) ! (klon, nl)
+    ! specific heat capacity at constant pressure of humid air, in J K-1 kg-1
+
+    real tv1(:, :) ! (klon, klev)
     real gz1(klon, klev), h1(klon, klev), hm1(klon, klev)
-    real th1(klon, klev) ! potential temperature
+    real, intent(out):: th1(:, :) ! (klon, nl) potential temperature, in K
 
     ! Local:
     integer k, i
-    real rdcp
+    real kappa
     real tvx, tvy 
     real cpx(klon, klev)
 
@@ -33,11 +41,11 @@
     do k = 1, nl 
        do i = 1, klon
           lv1(i, k) =  rlvtt - clmcpv * (t1(i, k) - 273.15)
-          cpn1(i, k) = cpd * (1. - q1(i, k)) + cpv * q1(i, k)
-          cpx(i, k) = cpd * (1. - q1(i, k)) + cl * q1(i, k)
-          tv1(i, k) = t1(i, k) * (1. + q1(i, k)/eps - q1(i, k))
-          rdcp = (rrd * (1. - q1(i, k)) + q1(i, k) * rrv)/cpn1(i, k)
-          th1(i, k) = t1(i, k) * (1000./p1(i, k))**rdcp
+          cpn1(i, k) = rcpd * (1. - q1(i, k)) + rcpv * q1(i, k)
+          cpx(i, k) = rcpd * (1. - q1(i, k)) + rcw * q1(i, k)
+          tv1(i, k) = t1(i, k) * (1. + q1(i, k) / eps - q1(i, k))
+          kappa = (rd * (1. - q1(i, k)) + q1(i, k) * rv) / cpn1(i, k)
+          th1(i, k) = t1(i, k) * (1000. / p1(i, k))**kappa
        end do
     end do
 
@@ -49,10 +57,10 @@
 
     do k = 2, nl 
        do i = 1, klon
-          tvx = t1(i, k) * (1. + q1(i, k)/eps - q1(i, k))
-          tvy = t1(i, k - 1) * (1. + q1(i, k - 1)/eps - q1(i, k - 1))
-          gz1(i, k) = gz1(i, k - 1) + 0.5 * rrd * (tvx + tvy) &
-               * (p1(i, k - 1) - p1(i, k))/ph1(i, k)
+          tvx = t1(i, k) * (1. + q1(i, k) / eps - q1(i, k))
+          tvy = t1(i, k - 1) * (1. + q1(i, k - 1) / eps - q1(i, k - 1))
+          gz1(i, k) = gz1(i, k - 1) + 0.5 * rd * (tvx + tvy) &
+               * (p1(i, k - 1) - p1(i, k)) / ph1(i, k)
        end do
     end do