--- trunk/Sources/phylmd/CV3_routines/cv3_feed.f	2015/04/29 15:47:56	134
+++ trunk/Sources/phylmd/CV30_routines/cv30_feed.f	2016/06/06 17:42:15	201
@@ -1,114 +1,78 @@
-module cv3_feed_m
+module cv30_feed_m
 
   implicit none
 
 contains
 
-  SUBROUTINE cv3_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, &
-       iflag, tnk, qnk, gznk, plcl)
+  SUBROUTINE cv30_feed(t1, q1, qs1, p1, ph1, gz1, icb1, iflag1, tnk1, qnk1, &
+       gznk1, plcl1)
 
-    ! Purpose: CONVECTIVE FEED
+    ! Purpose: convective feed
+    ! Assuming origin level of ascending parcels is minorig.
 
-    ! Main differences with cv_feed:
-    ! - ph added in input
-    ! - here, nk(i)=minorig
-    ! - icb defined differently (plcl compared with ph instead of p)
-
-    ! Main differences with convect3:
-    ! - we do not compute dplcldt and dplcldr of CLIFT anymore
-    ! - values iflag different (but tests identical)
-    ! - A, B explicitely defined (!)
-
-    use cv3_param_m
-
-    ! inputs:
-    integer, intent(in):: len, nd
-    real, intent(in):: t(len, nd)
-    real, intent(in):: q(len, nd), qs(len, nd), p(len, nd)
-    real hm(len, nd), gz(len, nd)
-    real, intent(in):: ph(len, nd+1)
+    use cv30_param_m, only: minorig, nl
+    USE dimphy, ONLY: klev, klon
+    use numer_rec_95, only: locate
+
+    real, intent(in):: t1(:, :) ! (klon, klev)
+    real, intent(in):: q1(:, :), qs1(:, :), p1(:, :) ! (klon, klev)
+    real, intent(in):: ph1(:, :) ! (klon, klev+1)
+    real, intent(in):: gz1(:, :) ! (klon, klev)
 
     ! outputs:
-    integer iflag(len)
-    integer, intent(out):: nk(len), icb(len), icbmax
-    real tnk(len), qnk(len), gznk(len), plcl(len)
-
-    ! local variables:
-    integer i, k
-    integer ihmin(len)
-    real work(len)
-    real pnk(len), qsnk(len), rh(len), chi(len)
-    real A, B ! convect3
 
-    !--------------------------------------------------------------------
-
-    plcl=0.0
-
-    ! --- Origin level of ascending parcels for convect3:
-
-    do i=1, len
-       nk(i)=minorig
-    end do
-
-    ! --- Check whether parcel level temperature and specific humidity
-    ! --- are reasonable
-
-    do i=1, len
-       if ((t(i, nk(i)) < 250. .or. q(i, nk(i)) <= 0.) .and. iflag(i) == 0) &
-            iflag(i)=7
-    end do
-
-    ! --- Calculate lifted condensation level of air at parcel origin level
-    ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV., 1980)
-
-    A = 1669.0 ! convect3
-    B = 122.0 ! convect3
-
-    do i=1, len
-       if (iflag(i).ne.7) then
-          tnk(i)=t(i, nk(i))
-          qnk(i)=q(i, nk(i))
-          gznk(i)=gz(i, nk(i))
-          pnk(i)=p(i, nk(i))
-          qsnk(i)=qs(i, nk(i))
-
-          rh(i)=qnk(i)/qsnk(i)
-          chi(i)=tnk(i)/(A-B*rh(i)-tnk(i)) ! convect3
-          plcl(i)=pnk(i)*(rh(i)**chi(i))
-          if ((plcl(i) < 200. .or. plcl(i) >= 2000.) .and. iflag(i) == 0) &
-               iflag(i) = 8
-       endif
-    end do
-
-    ! --- Calculate first level above lcl (=icb)
-
-    do i=1, len
-       icb(i)=nlm
-    end do
+    integer, intent(out):: icb1(:) ! (klon)
+    ! first level above LCL, 2 <= icb1 <= nl - 2
 
-    ! la modification consiste a comparer plcl a ph et non a p:
-    ! icb est defini par : ph(icb) < plcl < ph(icb - 1)
-    do k=3, nl-1 ! modification pour que icb soit supérieur ou égal à 2
-       do i=1, len
-          if(ph(i, k) < plcl(i)) icb(i) = min(icb(i), k)
-       end do
-    end do
+    integer, intent(out):: iflag1(:) ! (klon)
+    real tnk1(:), qnk1(:), gznk1(:) ! (klon)
+    real, intent(out):: plcl1(:) ! (klon)
+
+    ! Local:
+    integer i
+    real rh(klon)
+    real, parameter:: A = 1669., B = 122.
 
-    do i=1, len
-       if((icb(i) == nlm).and.(iflag(i) == 0))iflag(i)=9
-    end do
-
-    do i=1, len
-       icb(i) = icb(i)-1 ! icb sup ou egal a 2
-    end do
-
-    ! Compute icbmax.
+    !--------------------------------------------------------------------
 
-    icbmax=2
-    do i=1, len
-       if (iflag(i) < 7) icbmax=max(icbmax, icb(i)) ! sb Jun7th02
-    end do
+    ! Calculate lifted condensation level of air at parcel origin level
+    ! (within 0.2 % of formula of Bolton, Mon. Wea. Rev., 1980)
+    where (t1(:, minorig) >= 250. .and. q1(:, minorig) > 0.)
+       ! Parcel level temperature and specific humidity are reasonable.
+       tnk1 = t1(:, minorig)
+       qnk1 = q1(:, minorig)
+       gznk1 = gz1(:, minorig)
+
+       rh = qnk1 / qs1(:, minorig)
+       plcl1 = p1(:, minorig) * rh**(tnk1 / (A - B * rh - tnk1))
+       iflag1 = 0
+    elsewhere
+       plcl1 = 0.
+       iflag1 = 7
+    end where
+
+    where (iflag1 == 0 .and. (plcl1 < 200. .or. plcl1 >= 2000.)) iflag1 = 8
+
+    ! Compute icb1:
+    do i = 1, klon
+       if (plcl1(i) <= ph1(i, nl - 2)) then
+          ! Distinguish this case just so that icb1 = nl - 2, not nl =
+          ! 3, for plcl1 exactly == ph1(i, nl - 2). Maybe not useful.
+          icb1(i) = nl - 2
+       else
+          icb1(i) = locate(- ph1(i, 3:nl - 2), - plcl1(i), my_lbound = 3)
+          ! {2 <= icb1(i) <= nl - 3}
+          ! {ph1(i, icb1(i) + 1) < plcl1(i)}
+          ! {plcl1(i) <= ph1(i, icb1(i)) or icb1(i) == 2}
+       end if
+    end do
+
+    where (icb1 == nl - 2 .and. iflag1 == 0) iflag1 = 9
+
+    ! {(2 <= icb1(i) <= nl - 3 and ph1(i, icb1(i) + 1) < plcl1(i) and
+    ! (plcl1(i) <= ph1(i, icb1(i)) or icb1(i) == 2)) or iflag1(i) /=
+    ! 0}
 
-  end SUBROUTINE cv3_feed
+  end SUBROUTINE cv30_feed
 
-end module cv3_feed_m
+end module cv30_feed_m