Index: codes/icosagcm/trunk/src/disvert.f90
===================================================================
--- codes/icosagcm/trunk/src/disvert.f90	(revision 206)
+++ codes/icosagcm/trunk/src/disvert.f90	(revision 207)
@@ -23,4 +23,5 @@
   SUBROUTINE init_disvert
   USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert
+  USE disvert_apbp_mod, ONLY: ap_apbp=>ap, bp_apbp=>bp, presnivs_apbp=>presnivs, init_disvert_apbp=>init_disvert
   USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert
   USE disvert_ncarl30_mod, ONLY: ap_ncarl30=>ap, bp_ncarl30=>bp, presnivs_ncarl30=>presnivs, init_disvert_ncarl30=>init_disvert
@@ -64,4 +65,11 @@
         presnivs=>presnivs_std
       
+      CASE('read_apbp')
+    
+        CALL init_disvert_apbp
+        ap=>ap_apbp
+        bp=>bp_apbp
+        presnivs=>presnivs_apbp
+        
       CASE ('ncar')
 
Index: codes/icosagcm/trunk/src/disvert_apbp.f90
===================================================================
--- codes/icosagcm/trunk/src/disvert_apbp.f90	(revision 207)
+++ codes/icosagcm/trunk/src/disvert_apbp.f90	(revision 207)
@@ -0,0 +1,83 @@
+MODULE disvert_apbp_mod
+  USE icosa
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
+!$OMP THREADPRIVATE(ap)
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:)
+!$OMP THREADPRIVATE(bp)
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
+!$OMP THREADPRIVATE(presnivs)
+
+CONTAINS
+
+  SUBROUTINE init_disvert
+!  USE icosa
+  IMPLICIT NONE
+  
+    ALLOCATE(ap(llm+1))
+    ALLOCATE(bp(llm+1))
+    ALLOCATE(presnivs(llm))
+
+    CALL disvert(ap,bp,presnivs)    
+
+  END SUBROUTINE init_disvert  
+
+
+  SUBROUTINE disvert(ap,bp,presnivs)
+!  USE icosa
+  USE mpipara, ONLY: is_mpi_root
+  USE omp_para, ONLY: omp_in_parallel
+  USE transfert_omp_mod, ONLY: bcast_omp
+  IMPLICIT NONE
+  REAL(rstd),INTENT(OUT) :: ap(:)
+  REAL(rstd),INTENT(OUT) :: bp(:)
+  REAL(rstd),INTENT(OUT) :: presnivs(:)
+  
+  INTEGER,PARAMETER :: unit=42
+  CHARACTER(len=255) :: filename
+  INTEGER :: l,ok
+  
+    filename="apbp.txt" !file to read
+    ! but users may want to use some other file name
+    CALL getin('read_apbp_file',filename)
+    
+!$OMP MASTER
+    OPEN(unit,file=filename,status="old",action="read",iostat=ok)
+    IF (ok/=0) THEN
+      WRITE(*,*) "disvert_ap_bp error: input file ",trim(filename)," not found!"
+      CALL ABORT
+    ENDIF
+    ! read in ap() and b() line by line, starting from surface up
+    ! to model top
+    DO l=1,llm+1
+      READ(unit,fmt=*,iostat=ok) ap(l),bp(l)
+      IF (ok/=0) THEN
+        WRITE(*,*) "disvert_ap_bp error: failed reading ap(l) and bp(l) for l=",l
+        CALL ABORT
+      ENDIF
+    ENDDO
+!$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
+      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 ",scaleheight," (km)"
+      DO l=1,llm
+        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scaleheight,       &
+                   ' DZ ~ ',scaleheight*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
+      ENDDO
+    ENDIF
+  
+  END SUBROUTINE disvert
+  
+END MODULE disvert_apbp_mod
Index: codes/icosagcm/trunk/src/earth_const.f90
===================================================================
--- codes/icosagcm/trunk/src/earth_const.f90	(revision 206)
+++ codes/icosagcm/trunk/src/earth_const.f90	(revision 207)
@@ -11,4 +11,5 @@
   REAL(rstd),SAVE :: preff=101325.
   REAL(rstd),SAVE :: pa=50000.
+  REAL(rstd),SAVE :: scaleheight=8. ! atmospheric scale height (km)
   REAL(rstd),SAVE :: scale_factor=1.
 
@@ -29,4 +30,5 @@
     CALL getin("cpp",cpp)  
     CALL getin("preff",preff)  
+    CALL getin("scaleheight",scaleheight)
     
     boussinesq=.FALSE.
