--- trunk/dyn3d/etat0.f	2014/03/26 17:18:58	91
+++ trunk/dyn3d/etat0.f	2014/09/18 19:56:46	113
@@ -13,7 +13,7 @@
 
 contains
 
-  SUBROUTINE etat0
+  SUBROUTINE etat0(phis)
 
     ! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09
 
@@ -43,18 +43,20 @@
     use nr_util, only: pi, assert
     use paramet_m, only: ip1jm, ip1jmp1
     use phyredem_m, only: phyredem
-    use pressure_var, only: pls, p3d
     use q_sat_m, only: q_sat
     use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
     use regr_pr_o3_m, only: regr_pr_o3
-    use serre, only: alphax
     use startdyn, only: start_init_dyn
     USE start_init_orog_m, only: start_init_orog, mask
     use start_init_phys_m, only: start_init_phys
     use start_inter_3d_m, only: start_inter_3d
     use temps, only: itau_phy, annee_ref, day_ref
+    use test_disvert_m, only: test_disvert
 
-    ! Variables local to the procedure:
+    REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1)
+    ! surface geopotential, in m2 s-2
+
+    ! Local:
 
     REAL latfi(klon), lonfi(klon)
     ! (latitude and longitude of a point of the scalar grid identified
@@ -69,19 +71,17 @@
     ! and pressure level "pls(i, j, l)".)
 
     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
-    REAL tsol(klon), qsol(klon), sn(klon)
-    REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf) 
+    REAL sn(klon)
+    REAL qsolsrf(klon, nbsrf), snsrf(klon, nbsrf) 
     REAL albe(klon, nbsrf), evap(klon, nbsrf)
     REAL alblw(klon, nbsrf)
     REAL tsoil(klon, nsoilmx, nbsrf) 
     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
     REAL solsw(klon), sollw(klon), fder(klon)
     !IM "slab" ocean
-    REAL tslab(klon)
     real seaice(klon) ! kg m-2
     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
     REAL rugmer(klon)
-    REAL phis(iim + 1, jjm + 1) ! surface geopotential, in m2 s-2
     real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
     real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, ps
@@ -113,6 +113,15 @@
     real sig1(klon, llm) ! section adiabatic updraft
     real w01(klon, llm) ! vertical velocity within adiabatic updraft
 
+    real pls(iim + 1, jjm + 1, llm)
+    ! (pressure at mid-layer of LMDZ grid, in Pa)
+    ! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
+    ! for layer "l")
+
+    REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
+    ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
+    ! for interface "l")
+
     !---------------------------------
 
     print *, "Call sequence information: etat0"
@@ -123,6 +132,7 @@
 
     pa = 5e4
     CALL disvert
+    call test_disvert
     CALL inigeom
     CALL inifilr
 
@@ -195,16 +205,13 @@
     if (nqmx >= 5) then
        ! Ozone:
        call regr_lat_time_coefoz
-       call regr_pr_o3(q(:, :, :, 5))
+       call regr_pr_o3(p3d, q(:, :, :, 5))
        ! Convert from mole fraction to mass fraction:
        q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
     end if
 
-    tsol = pack(tsol_2d, dyn_phy)
-    qsol = pack(qsol_2d, dyn_phy)
     sn = 0. ! snow
     radsol = 0.
-    tslab = 0. ! IM "slab" ocean
     seaice = 0.
     rugmer = 0.001
     zmea = pack(zmea_2d, dyn_phy)
@@ -288,8 +295,6 @@
     ! Calcul interm\'ediaire :
     CALL massdair(p3d, masse)
 
-    print *, 'ALPHAX = ', alphax
-
     forall (l = 1:llm)
        masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
        masse(:, jjm + 1, l) = &
@@ -309,10 +314,6 @@
     CALL dynredem1("start.nc", vcov, ucov, teta, q, masse, ps, itau=0)
 
     ! Initialisations :
-    tsolsrf(:, is_ter) = tsol
-    tsolsrf(:, is_lic) = tsol
-    tsolsrf(:, is_oce) = tsol
-    tsolsrf(:, is_sic) = tsol
     snsrf(:, is_ter) = sn
     snsrf(:, is_lic) = sn
     snsrf(:, is_oce) = sn
@@ -323,11 +324,8 @@
     albe(:, is_sic) = 0.6
     alblw = albe
     evap = 0.
-    qsolsrf(:, is_ter) = 150.
-    qsolsrf(:, is_lic) = 150.
-    qsolsrf(:, is_oce) = 150.
-    qsolsrf(:, is_sic) = 150.
-    tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
+    qsolsrf = 150.
+    tsoil = spread(spread(pack(tsol_2d, dyn_phy), 2, nsoilmx), 3, nbsrf)
     rain_fall = 0.
     snow_fall = 0.
     solsw = 165.
@@ -335,8 +333,6 @@
     t_ancien = 273.15
     q_ancien = 0.
     agesno = 0.
-    !IM "slab" ocean
-    tslab = tsolsrf(:, is_oce)
     seaice = 0.
 
     frugs(:, is_oce) = rugmer
@@ -351,11 +347,11 @@
     sig1 = 0.
     w01 = 0.
 
-    call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
-         tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
-         evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
-         agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
-         t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
+    call phyredem("startphy.nc", latfi, lonfi, pctsrf, tsoil(:, 1, :), tsoil, &
+         tsoil(:, 1, is_oce), seaice, qsolsrf, pack(qsol_2d, dyn_phy), snsrf, &
+         albe, alblw, evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, &
+         frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
+         q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
     CALL histclo
 
   END SUBROUTINE etat0