Index: trunk/ICOSA_LMDZ/bld.cfg
===================================================================
--- trunk/ICOSA_LMDZ/bld.cfg	(revision 2424)
+++ trunk/ICOSA_LMDZ/bld.cfg	(revision 2452)
@@ -75,4 +75,13 @@
 bld::excl_dep        use::write_field_mod
 bld::excl_dep        use::checksum_mod
+## addition for disvert_icosa_lmdz
+bld::excl_dep        use::prec
+bld::excl_dep        use::getin_mod
+bld::excl_dep        use::grid_param
+bld::excl_dep        use::earth_const
+bld::excl_dep        use::disvert_plugin_mod
+bld::excl_dep        use::abort_mod
+bld::excl_dep        use::free_unit_mod
+bld::excl_dep        use::transfert_omp_mod
 ## additions for Venus:
 bld::excl_dep        use::cpdet_phy_mod
Index: trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90
===================================================================
--- trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90	(revision 2452)
+++ trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90	(revision 2452)
@@ -0,0 +1,184 @@
+module disvert_icosa_lmdz_mod
+
+implicit none
+
+contains
+
+  SUBROUTINE disvert_icosa_lmdz(ap,bp,presnivs)
+  ! Routine to compute ap(),bp() and presnivs() in the same way as done
+  ! in LMDZ for planets (see LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F)
+
+  USE prec, ONLY : rstd
+  USE getin_mod, ONLY : getin
+  USE grid_param, ONLY: llm
+  USE earth_const, ONLY: pa,preff,scale_height
+  USE abort_mod, ONLY: dynamico_abort
+  USE mpipara, ONLY: is_mpi_root
+  USE omp_para, ONLY: omp_in_parallel
+  USE transfert_omp_mod, ONLY: bcast_omp
+  USE free_unit_mod, ONLY : free_unit
+  IMPLICIT NONE
+
+  REAL(rstd),INTENT(OUT) :: ap(:)
+  REAL(rstd),INTENT(OUT) :: bp(:)
+  REAL(rstd),INTENT(OUT) :: presnivs(:)
+
+  INTEGER :: unit
+  CHARACTER(len=255) :: filename
+  INTEGER :: l,ok
+  LOGICAL :: hybrid
+  REAL(rstd) :: scaleheight ! atmospheric scale height (km) from file z2sig.def
+  REAL(rstd) :: zsig(llm)
+  REAL(rstd) :: sig(llm+1)
+  REAL(rstd) :: newsig
+  
+    filename="z2sig.def" !file to read
+    
+    ! User may choose between hybrid (default) or sigma coordinates
+    hybrid=.true. ! default: use hybrid coordinates
+    CALL getin('hybrid',hybrid)
+    ! User must also provide "pa" if in hybrid case
+    if (hybrid) then
+      pa=-999 ! absurd initial value to trigger error if not set
+      CALL getin('pa',pa)
+      IF (pa<0) THEN
+        WRITE(*,*) "disvert_icosa_lmdz error: pa=",pa
+        WRITE(*,*) " set a correct value for the pa parameter"
+        CALL dynamico_abort("disvert_icosa_lmdz : wrong value for pa")
+      ENDIF
+    endif
+    
+!$OMP MASTER
+    unit=free_unit()
+    OPEN(unit,file=filename,status="old",action="read",iostat=ok)
+    IF (ok/=0) THEN
+      WRITE(*,*) "disvert_icosa_lmdz error: input file ",trim(filename)," not found!"
+      CALL dynamico_abort("disvert_icosa_lmdz : could not open z2sig.def file" )
+    ENDIF
+    ! read "scaleheight" on the first line
+    READ(unit,fmt=*,iostat=ok) scaleheight
+    IF (ok/=0) THEN
+      WRITE(*,*) "disvert_icosa_lmdz error: failed reading scaleheight"
+      CALL dynamico_abort("disvert_icosa_lmdz : could not read input file")
+    ENDIF
+    ! While at it check that this "scaleheight" matches the "scale_height"
+    ! from the def files
+    IF ((1.-(scaleheight*1000.)/scale_height)>1.e-6) THEN
+      WRITE(*,*) "disvert_icosa_lmdz error: significant mismatch between"
+      WRITE(*,*) " scale_height (m) from def file:",scale_height
+      WRITE(*,*) " and scaleheight (km) from z2sig.def:",scaleheight
+      CALL dynamico_abort("disvert_icosa_lmdz : atmospheric scale height issue")
+    ENDIF
+    ! read zsig() line by line, starting from the surface up to model top
+    DO l=1,llm
+      READ(unit,fmt=*,iostat=ok) zsig(l)
+      IF (ok/=0) THEN
+        WRITE(*,*) "disvert_icosa_lmdz error: failed reading zsig(l) for l=",l
+        CALL dynamico_abort("disvert_icosa_lmdz : could not read input file" )
+      ENDIF
+    ENDDO
+    CLOSE(unit)
+    
+    ! build sig()
+    sig(1)=1.
+    DO l=2,llm
+      sig(l)=0.5*(exp(-zsig(l)/scaleheight)+exp(-zsig(l-1)/scaleheight))
+    ENDDO
+    sig(llm+1)=0 ! here we asume 0 pressure at top of atmosphere
+
+    ! build  ap() and bp() using sig(), pa and preff
+    IF (hybrid) THEN  ! use hybrid coordinates
+      DO l=1,llm
+        CALL sig_hybrid(sig(l),pa,preff,newsig)
+        bp(l)=exp(1.-1./(newsig**2))
+        ap(l)=pa*(newsig-bp(l))
+      ENDDO
+    ELSE ! use sigma coordinates
+      ap(1:llm)=0
+      bp(1:llm)=sig(1:llm)
+    ENDIF
+    ap(llm+1)=0
+    bp(llm+1)=0
+!$OMP END MASTER
+
+    IF (omp_in_parallel()) THEN
+      CALL bcast_omp(ap)
+      CALL bcast_omp(bp)
+    ENDIF
+    
+    ! build presnivs(), approximative mid-layer pressures
+    DO l=1,llm
+      presnivs(l)=0.5*(ap(l)+bp(l)*preff+ap(l+1)+bp(l+1)*preff)
+    ENDDO
+    
+    ! tell the world about it
+    IF (is_mpi_root) THEN
+!$OMP MASTER
+      WRITE(*,*) "Vertical discretization set by disvert_icosa_lmdz with hybrid=",hybrid
+      WRITE(*,*) "ap()=",ap
+      WRITE(*,*) "bp()=",bp
+      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
+      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
+      DO l=1,llm
+        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
+                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
+      ENDDO
+!$OMP END MASTER
+    ENDIF
+  
+    CONTAINS
+  
+    SUBROUTINE sig_hybrid(sig,pa,preff,newsig)
+    ! F. Forget 2002
+    ! Iterative procedure to find "newsig" the sigma value for a layer
+    ! based on the target values provided in z2sig.def (sig)
+    ! but such that:
+    ! (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
+    ! iterative solving of the equation uses the fact that:
+    ! where exp(1-1./x**2) become << x
+    ! x      exp(1-1./x**2)/x
+    ! 1           1
+    ! 0.68       0.5
+    ! 0.5        1.E-1
+    ! 0.391      1.E-2
+    ! 0.333      1.E-3
+    ! 0.295      1.E-4
+    ! 0.269      1.E-5
+    ! 0.248      1.E-6
+    ! => thus on can assume  newsig = sig*preff/pa if sig*preff/pa < 0.25
+    
+      IMPLICIT NONE
+      REAL(rstd),INTENT(in) :: sig
+      REAL(rstd),INTENT(in) :: pa
+      REAL(rstd),INTENT(in) :: preff
+      REAL(rstd),INTENT(out) :: newsig
+      REAL(rstd) :: F,x1,x2
+      INTEGER :: j
+    
+      newsig=sig
+      x1=0
+      x2=1
+      IF (sig.ge.1) THEN ! handle sigma=1 case
+        newsig=sig
+      ELSE IF (sig*preff/pa.ge.0.25) THEN
+        DO J=1,9999 ! Overkill max number of iterations
+          F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
+          IF (F.gt.1) THEN
+            x2=newsig
+            newsig=(x1+newsig)*0.5
+          ELSE
+            x1=newsig
+            newsig=(x2+newsig)*0.5
+          ENDIF
+          ! Convergence test: stop if sig is approximated to within 0.01m
+          IF (abs(10.*log(F)).LT.1.E-5) EXIT
+        ENDDO
+      ELSE ! IF (sig*preff/pa.le.0.25) THEN
+        newsig=sig*preff/pa
+      ENDIF
+    END SUBROUTINE sig_hybrid
+
+  END SUBROUTINE disvert_icosa_lmdz
+
+
+end module disvert_icosa_lmdz_mod
Index: trunk/ICOSA_LMDZ/src/icosa_lmdz.f90
===================================================================
--- trunk/ICOSA_LMDZ/src/icosa_lmdz.f90	(revision 2424)
+++ trunk/ICOSA_LMDZ/src/icosa_lmdz.f90	(revision 2452)
@@ -1,4 +1,9 @@
 PROGRAM icosa_lmdz
-  USE icosa_init_mod
+  USE icosa_init_mod, ONLY: icosa_init
+  USE disvert_plugin_mod, ONLY: disvert_plugin
+  USE disvert_icosa_lmdz_mod, ONLY: disvert_icosa_lmdz
+
+  ! set up plugins
+  disvert_plugin => disvert_icosa_lmdz
 
   CALL icosa_init
