Index: unk/LMDZ.GENERIC/libf/phystd/N2broadprof_fn.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/N2broadprof_fn.F90	(revision 304)
+++ 	(revision )
@@ -1,110 +1,0 @@
-      subroutine N2broadprof_fn(tsurf,temp,popsk,play,plev,qvap)
-      implicit none
-
-!     creates a temperature profile for an N2-CO2 atmosphere
-!     a fine resolution grid is used for accuracy
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "callkeys.h"
-
-      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
-      REAL plev(nlayermx+1) ! Intermediate pressure levels (pa)
-      REAL psurf,tsurf      
-!      REAL tcond(nlayermx)
-      REAL qvap(nlayermx)
-      REAL temp(nlayermx)   ! temperature at the middle of the layers
-
-      real, parameter :: Rstar   = 8.314
-      real, parameter :: ptriple = 518000.0
-      real, parameter :: latcond = 5.9e5     ! as in condense_co2
-      real, parameter :: mn      = 0.02801   ! molar mass of background gas (N2)
-      real, parameter :: mv      = 0.04401   ! molar mass of background gas (CO2)
-
-      real popsk(nlayermx)  ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx)
-                            ! but kastprof requires 1D, so its ok I guess
-
-      real gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p
-
-      integer l,iter
-
-      real sigma_loc,q_loc,vmr_loc,ppvar_loc,T_loc,p_loc
-      real Tcond
-      real Told,qold,pold
-
-      integer, parameter :: Niter = 2000
-
-      if(ngridmx.gt.1)then
-         print*,'DO NOT RUN KASTPROF IN 3D!!!'
-         call abort
-      end if
-
-      l=1
-
-      vmr_loc = 0.05 ! value at ground
-      q_loc   = vmr_loc*mv/mn
-      T_loc   = Tsurf
-
-      dln_p = log(pceil/plev(1))/Niter ! please check
-      do iter=1,Niter
-
-         sigma_loc = (pceil/plev(1))**(real(iter)/real(Niter))
-         p_loc     = plev(1)*sigma_loc
-
-         ppvar_loc = vmr_loc*p_loc
-         if(ppvar_loc.lt.ptriple)then
-            Tcond = (-3167.8)/(log(.01*ppvar_loc)-23.23)         ! Fanale's formula
-         else
-            Tcond = 684.2-92.3*log(ppvar_loc)+4.32*log(ppvar_loc)**2 ! liquid-vapour (CRC handbook 2003 data)
-         endif
-
-            if (T_loc.le.Tstrat)then
-               T_loc     = Tstrat            ! isothermal stratosphere
-            elseif(T_loc.le.Tcond)then
-
-               gamma     = cpp - q_loc*latcond/T_loc
-               dqdlnT    = (mv*latcond/(mn*T_loc) - gamma) & 
-                            / (latcond/T_loc + Rstar/(mn*q_loc))
-               dlnqdlnT  = (1./q_loc)*dqdlnT
-               dlnpdlnT  = mv*latcond/(Rstar*T_loc) - dlnqdlnT/(1. + q_loc*mn/mv)
-               dpdT      = (p_loc/T_loc)*dlnpdlnT
-               q_loc     = q_loc  + dqdlnT*dln_p/dlnpdlnT
-               vmr_loc   = q_loc*mn/mv ! needed for Tcond calc.
-               T_loc     = T_loc*(1 + dln_p/dlnpdlnT)
-               print*,T_loc
-            else
-               ! nothing doing; so dry convection
-               T_loc     = tsurf*(p_loc/plev(1))**rcp
-            endif
-
-            ! update true gridpoints when we reach them
-            if((p_loc.le.play(l)) .and. (l.le.nlayermx))then
-               temp(l)=T_loc
-               qvap(l)=q_loc
-               ! do an avg you lazy bastard
-!               temp(l)=T_loc*(p_loc-play(l))/(p_loc-pold)+Told*(play(l)-pold)/(p_loc-pold)
-!               qvap(l)=q_loc*(p_loc-play(l))/(p_loc-pold)+qold*(play(l)-pold)/(p_loc-pold)
-
-!               f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
-
-!               print*,'l=',l
-!               print*,'temp(l)=',temp(l)
-!               print*,'qvap(l)=',qvap(l)
-
-
-               qvap(l)=0.05*mv/mn ! suface amount everywhere to test
-
-               l=l+1
-            endif
-
-            Told=T_loc
-            qold=q_loc
-            pold=p_loc
-
-      enddo ! vertical loop
-
-      print*,'---fin---'
-
-      return
-    end subroutine N2broadprof_fn
Index: /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 305)
@@ -181,9 +181,41 @@
 !       2. Opacity calculation
             if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then
-               do ig=1, ngrid ! to stop the rad tran bug
+               do ig=1,ngrid ! to stop the rad tran bug
                   do l=1,nlayer
                      aerosol(ig,l,iaer) =1.e-9
                   end do
                end do
+
+               ! put cloud at cloudlvl
+               if(kastprof.and.(cloudlvl.ne.0.0))then
+                  ig=1
+                  do l=1,nlayer
+                     if(int(cloudlvl).eq.l)then
+                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
+                        print*,'Inserting cloud at level ',l
+                        !aerosol(ig,l,iaer)=10.0
+
+                        rho_ice=920.0
+
+                        ! the Kasting approximation
+                        aerosol(ig,l,iaer) =                      &
+                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
+                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
+                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
+                          ( 4.0e-4 + 1.E-9 ) *         &
+                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+
+
+                        open(115,file='clouds.out',form='formatted')
+                        write(115,*) l,aerosol(ig,l,iaer)
+                        close(115)
+
+                        return
+                     endif
+                  end do
+
+                  call abort
+               endif
+
             else
                do ig=1, ngrid
@@ -232,6 +264,6 @@
 !       2. Opacity calculation
 
-            do l=1,nlayer
-               do ig=1,ngrid-1 ! to stop the rad tran bug
+            do ig=1,ngrid
+               do l=1,nlayer
                   zp=700./pplay(ig,l)
                   aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
Index: /trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90	(revision 305)
@@ -0,0 +1,80 @@
+!-------------------------------------------------------------------------
+      subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x_in,y_in,f)
+!     Necessary for interpolation of continuum data
+
+      implicit none
+
+      integer nX,nY,i,j,a,b
+
+      real*8 x_in,y_in,x,y,x1,x2,y1,y2
+      real*8 f,f11,f12,f21,f22,fA,fB
+      real*8 x_arr(nX)
+      real*8 y_arr(nY)
+      real*8 f2d_arr(nX,nY)
+      
+      integer strlen
+      character*100 label
+      label='subroutine bilinear'
+
+      x=x_in
+      y=y_in
+
+!     1st check we're within the wavenumber range
+      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
+         f=0.0D+0
+         return
+      else
+         
+!     in the x (wavenumber) direction 1st
+         i=1
+ 10      if (i.lt.(nX+1)) then
+            if (x_arr(i).gt.x) then
+               x1=x_arr(i-1)
+               x2=x_arr(i)
+               a=i-1
+               i=9999
+            endif
+            i=i+1
+            goto 10
+         endif
+      endif
+      
+      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
+         write(*,*) 'Warning from bilinear.for:'
+         write(*,*) 'Outside continuum temperature range!'
+         if(y.lt.y_arr(1))then
+            y=y_arr(1)+0.01
+         endif
+         if(y.gt.y_arr(nY))then
+            y=y_arr(nY)-0.01
+         endif
+      else
+
+!     in the y (temperature) direction 2nd
+         j=1
+ 20      if (j.lt.(nY+1)) then
+            if (y_arr(j).gt.y) then
+               y1=y_arr(j-1)
+               y2=y_arr(j)
+               b=j-1
+               j=9999
+            endif
+            j=j+1
+            goto 20
+         endif
+      endif
+      
+      f11=f2d_arr(a,b)
+      f21=f2d_arr(a+1,b)
+      f12=f2d_arr(a,b+1)
+      f22=f2d_arr(a+1,b+1)
+      
+!     1st in x-direction
+      fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
+      fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
+      
+!     then in y-direction
+      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
+      
+      return
+    end subroutine bilinear
Index: /trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 305)
@@ -46,11 +46,9 @@
 
       ! tau0/p0=tau/p (Hansen 1974)
-      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0)  *  4*delta^2/(g*mugaz*lambda^4)
-      ! Then calculate tau0 = 1.16e-20 * 4*delta^2/(g*mugaz*lambda^4 [in um])
-      ! where delta=n-1
+      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
+      ! where delta=n-1 and N0 is an amagat
 
       typeknown=.false.
       do igas=1,ngasmx
-         !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then
          if(igas.eq.vgas)then
             print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// &
@@ -60,5 +58,5 @@
             'as its mixing ratio is less than 0.05.' 
             ! ignore variable gas in Rayleigh calculation
-            ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation
+            ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
             tauconsti(igas) = 0.0
          else
@@ -67,9 +65,11 @@
             elseif(gnom(igas).eq.'N2_')then
                tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
+            elseif(gnom(igas).eq.'H2O')then
+               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
             elseif(gnom(igas).eq.'H2_')then
-               tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
+               tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
+               !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
                ! uses n=1.000132 from Optics, Fourth Edition.
-               ! around four times more opaque than CO2
-               ! around 5.5 times more opaque than air
+               ! was out by a factor 2!
             elseif(gnom(igas).eq.'He_')then
                print*,'Helium not ready yet!'
@@ -104,8 +104,7 @@
             tauvar=0.0
             do igas=1,ngasmx
-               !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then
                if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then
                   ! ignore variable gas in Rayleigh calculation
-                  ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation
+                  ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
                   tauvari(igas)   = 0.0
                else
@@ -114,6 +113,8 @@
                   elseif(gnom(igas).eq.'N2_')then
                      tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
+                  elseif(gnom(igas).eq.'H2O')then
+                     tauvari(igas) = 1.0/wl**4 ! to be improved...
                   elseif(gnom(igas).eq.'H2_')then
-                     tauvari(igas) = 1.0/wl**4 ! no wvl dependence of ref. index here - improve?
+                     tauvari(igas) = 1.0/wl**4 
                   else
                      print*,'Error in calc_rayleigh: Gas species not recognised!'
@@ -138,8 +139,8 @@
       end do
 
-      print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-5), &
-           ' um, tau_R = ',TAURAY(L_NSPECTV-5)*1013.25
-      print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6), &
-           ' um, tau_R = ',TAURAY(L_NSPECTV-6)*1013.25
+      print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6),' um'
+      print*,'tau_R = ',TAURAY(L_NSPECTV-6)*1013.25
+      print*,'sig_R = ',TAURAY(L_NSPECTV-6)*g*mugaz*1.67e-27*100, &
+             'cm^2 molecule^-1'
 
     end subroutine calc_rayleigh
Index: /trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90	(revision 305)
@@ -0,0 +1,79 @@
+subroutine calcenergy_kcm(Tsurf,T,Play,Plev,Qsurf,Q,muvar,Eatmtot)
+
+
+use params_h, only : Rc
+use watercommon_h, only : mH2O
+implicit none
+
+!     ----------------------------------------------------------------
+!     Purpose: Calculate total energy of the (steam) atmosphere
+!     Authour: R. Wordsworth (2011)
+!     ----------------------------------------------------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+!#include "callkeys.h"
+
+
+
+  ! inputs
+  real Tsurf,T(1:nlayermx)
+  real Play(1:nlayermx),Plev(1:nlayermx+1)
+  real Qsurf,Q(1:nlayermx)
+  real muvar(ngridmx,nlayermx+1)
+
+  ! internal
+  integer il
+  real m_n,m_v,cp_n,cp_v,cp_a,Eatmlay
+  real s_c,rho_v,L
+  double precision p_v, s_v, nul
+  real VMR(1:nlayermx)
+
+  ! output
+  real Eatmtot
+
+  ! functions
+  double precision cp_neutral
+
+  m_n = mugaz/1000.
+  m_v = mH2O/1000.
+
+
+  do il=1,nlayermx
+     VMR(il)=Q(il)*(muvar(1,il+1)/mH2O)
+  end do
+
+
+  Eatmtot = 0.0
+  do il=1,nlayermx
+
+
+     cp_n = cp_neutral(dble(T(il)))
+     cp_v = (32.24+1.923d-3*T(il) + 1.055d-5*T(il)**2-3.511d-9*T(il)**3)/m_v
+     ! / by m_v makes it per kg not per mol
+
+     call psat_H2O(dble(T(il)),p_v)
+     p_v   = p_v*1d6
+     rho_v = real(p_v)*m_v/(real(Rc)*T(il))
+
+     call therm(dble(T(il)),dble(rho_v)*1d-3,nul,nul,nul,nul,nul,nul,nul,&
+                nul,nul,nul,s_v,nul)
+
+     s_c = 2.06 * log(T(il)/273.15)
+     s_v = s_v * 1d3
+     s_c = s_c * 1d3
+     L   = (real(s_v)-s_c)*T(il)
+
+     cp_a    = (1-VMR(il))*cp_n + VMR(il)*cp_v 
+     !cp_a    = (1-Q(il))*cp_n + Q(il)*cp_v 
+
+     Eatmlay = ( cp_a*T(il) + L*VMR(il) ) * (Plev(il)-Plev(il+1))/g
+
+     Eatmtot = Eatmtot + Eatmlay 
+  enddo
+
+  print*,'WARNING: E-calc currently assumes H2O saturation!'
+
+end subroutine calcenergy_kcm
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 305)
@@ -1,9 +1,10 @@
-      subroutine callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,    &
+      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
           albedo,emis,mu0,pplev,pplay,pt,                      & 
-          tsurf,fract,dist_star,aerosol,cpp3D,                 &
+          tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
-          reffrad,tau_col,ptime,pday,cloudfrac,totcloudfrac,   &
+          reffrad,tau_col,cloudfrac,totcloudfrac,   &
           clearsky,firstcall,lastcall)
+
 
       use radinc_h
@@ -117,5 +118,5 @@
       real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
 
-      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component
+      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
       REAL pq(ngridmx,nlayer,nq)
       REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
@@ -159,8 +160,7 @@
       character(len=2)  :: tempOLR
       character(len=30) :: filenomOLR
-      real ptime, pday
+      !real ptime, pday
       logical OLRz
       real OLR_nu(ngrid,L_NSPECTI)
-
       real GSR_nu(ngrid,L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
@@ -177,7 +177,4 @@
       real pqtest(ngridmx,nlayer,nq)
 
-!     for Dave Crisp LBL data comparison
-      logical crisp
-
 !     are particle radii fixed?
       logical radfixed
@@ -187,5 +184,9 @@
       external CBRT
 
-      crisp=.false.
+!     included by RW for runway greenhouse 1D study
+      real muvar(ngridmx,nlayermx+1)
+      real vtmp(nlayermx)
+      REAL*8 muvarrad(L_LEVELS)
+
       radfixed=.false.
 
@@ -245,4 +246,5 @@
          Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
 
+
          if((igcm_h2o_vap.eq.0) .and. varactive)then
             print*,'varactive in callcorrk but no h2o_vap tracer.'
@@ -265,4 +267,8 @@
       enddo
 
+
+      if(kastprof)then
+         radfixed=.true.
+      endif
 
       if(radfixed)then
@@ -276,5 +282,6 @@
             do l=1,nlayer
                do ig=1,ngrid
-                  reffrad(ig,l,2) = 2.e-5 ! H2O ice
+                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
+                  reffrad(ig,l,2) = 5.e-6 ! H2O ice
                enddo
             enddo
@@ -297,8 +304,11 @@
 
          maxrad=0.0
+         !averad=0.0
          minrad=1.0
          do l=1,nlayer
+
+            !masse = (pplev(ig,l) - pplev(ig,l+1))/g
+
             do ig=1,ngrid
-
                if(tracer)then
                   reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
@@ -308,4 +318,5 @@
                reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
 
+               !averad = averad + reffrad(ig,l,1)*area(ig)
                maxrad = max(reffrad(ig,l,1),maxrad)
                minrad = min(reffrad(ig,l,1),minrad)
@@ -345,4 +356,5 @@
       endif
 
+
 !     how much light we get
       fluxtoplanet=0
@@ -360,4 +372,5 @@
            reffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
+
 
 !-----------------------------------------------------------------------
@@ -512,8 +525,6 @@
       endif
 
-      if(ngridmx.eq.1) then       ! fixed zenith angle in 1D
+      if(ngridmx.eq.1) then       ! fixed zenith angle 'szangle' in 1D
          acosz = cos(pi*szangle/180.0)
-         !acosz=mu0(ig)
-         !print*,'oned SZ angle OVERRIDE in callcorrk'
          print*,'acosz=',acosz,', szangle=',szangle
       else
@@ -521,17 +532,15 @@
       endif
 
-!-----------------------------------------------------------------------
-!     Water vapour (to be generalised / ignored for other planets)
+
+!-----------------------------------------------------------------------
+!     Water vapour (to be generalised for other gases eventually)
       
-      if(varactive) then
+      if(varactive)then
 
          i_var=igcm_h2o_vap
-
          do l=1,nlayer
-
             qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
             qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
             ! Average approximation as for temperature...
-
          end do
          qvar(1)=qvar(2)
@@ -575,8 +584,39 @@
       end if
 
-      ! Now convert from kg/kg to mol/mol
+      ! IMPORTANT: Now convert from kg/kg to mol/mol
       do k=1,L_LEVELS
          qvar(k) = qvar(k)/epsi
       end do
+
+      if(kastprof)then
+
+         DO l=1,nlayer
+            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
+            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
+                                muvar(ig,max(nlayer+1-l,1)))/2
+         END DO
+      
+         muvarrad(1) = muvarrad(2)
+         muvarrad(2*nlayermx+1)=muvar(ig,1)
+
+         print*,'Recalculating qvar with VARIABLE epsi for kastprof'
+         i_var=igcm_h2o_vap
+         do l=1,nlayer
+            vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
+         end do
+
+         do l=1,nlayer
+            qvar(2*l)   = vtmp(nlayer+1-l)
+            qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
+         end do
+         qvar(1)=qvar(2)
+         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
+
+         !do k=1,L_LEVELS
+         !   qvar(k) = 1.0
+         !end do
+         !print*,'ASSUMING qH2O=1 EVERYWHERE IN CALLCORRK!!'
+
+      endif
 
       ! Keep values inside limits for which we have radiative transfer coefficients
@@ -607,17 +647,4 @@
       tlevrad(2*nlayermx+1)=tsurf(ig)
       
-
-      if(crisp)then
-         open(111,file='/u/rwlmd/CrispLBL/p.dat')
-         open(112,file='/u/rwlmd/CrispLBL/T.dat')
-         do k=1,L_LEVELS
-            read(111,*) plevrad(k)
-            read(112,*) tlevrad(k)
-            plevrad(k)=plevrad(k)*1000.0
-         enddo
-         close(111)
-         close(112)
-      endif
-
       tmid(1) = tlevrad(2)
       tmid(2) = tlevrad(2)
@@ -651,11 +678,12 @@
             print*,'Minimum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
-            !print*,'WARNING, OVERRIDING FOR TEST'
-            call abort
+            print*,'WARNING, OVERRIDING FOR TEST'
+            !call abort
          elseif(tlevrad(k).gt.tgasmax)then
             print*,'Maximum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
-            !print*,'WARNING, OVERRIDING FOR TEST'
-            call abort
+            print*,'WARNING, OVERRIDING FOR TEST'
+            !print*, 'T=',pt
+            !call abort
          endif
       enddo
@@ -688,5 +716,5 @@
             call optcv(dtauv,tauv,taucumv,plevrad,                 &
                  qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
-                 tmid,pmid,taugsurf,qvar)
+                 tmid,pmid,taugsurf,qvar,muvarrad)
 
 
@@ -704,5 +732,5 @@
          end if
 
-	!print*,'taugsurf=',taugsurf(:,L_NGAUSS-1)
+
 
 
@@ -712,5 +740,5 @@
          call optci(plevrad,tlevrad,dtaui,taucumi,                  &
               qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
-              taugsurfi,qvar)
+              taugsurfi,qvar,muvarrad)
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
@@ -729,4 +757,5 @@
          fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
          fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
+
 
          if(fluxtop_dn(ig).lt.0.0)then
@@ -740,18 +769,4 @@
          endif
 
-         if(crisp)then
-            open(111,file='/u/rwlmd/CrispLBL/GCMfluxdn.dat')
-            open(112,file='/u/rwlmd/CrispLBL/GCMfluxup.dat')
-            do k=1,L_NLAYRAD
-               write(111,*) fluxdni(k)
-               write(112,*) fluxupi(k)
-            enddo
-            close(111)
-            close(112)
-            print*,'fluxdni',fluxdni
-            print*,'fluxupi',fluxupi
-            call abort
-         endif
-
 !     Spectral output, for exoplanet observational comparison
          if(specOLR)then
@@ -809,21 +824,13 @@
       if(specOLR)then
         if(ngrid.ne.1)then
-          call writediagspec(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
-!          call writediagspec(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
+          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
+          call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W m^-2",3,GSR_nu)
         endif
       endif
 
-      if(ngrid.eq.1)then
-         open(113,file='surf_vals_long.out',ACCESS='append')
-         write(113,*) tsurf(1),fluxtop_dn(1),         &
-              real(-nfluxtopv),real(nfluxtopi) 
-         close(113)
-      endif
-
-      !if(kastprof)then ! for kasthop only
-
-        if(ngrid.eq.1)then
-
-           print*,'Saving tsurf,psurf in surf_vals.out...'
+      if(lastcall.and.(ngrid.eq.1))then
+
+           print*,'Saving scalar quantities in surf_vals.out...'
+           print*,'psurf = ', pplev(1,1),' Pa'
            open(116,file='surf_vals.out')
            write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
@@ -843,22 +850,4 @@
                enddo
                close(127)
-           endif
-
-!     mixing ratio: CO2 ice particles
-           if(igcm_co2_ice.ne.0)then
-              open(122,file='qCO2ice.out')
-              do l=1,L_NLAYRAD
-                 write(122,*) pq(1,l,igcm_co2_ice)
-              enddo
-              close(122)
-           endif
-
-!     mixing ratio: H2O vapour
-           if(igcm_h2o_vap.ne.0)then
-              open(123,file='qH2Ovap.out')
-              do l=1,L_NLAYRAD
-                 write(123,*) pq(1,l,igcm_h2o_vap)
-              enddo
-              close(123)
            endif
 
@@ -878,7 +867,6 @@
               close(119)
            endif
-        endif
-
-      !endif
+
+      endif
 
     end subroutine callcorrk
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 305)
@@ -16,5 +16,5 @@
      &   , kastprof                                                     &
      &   , noradsurf                                                    &
-     &   , Tsurfk                                                       &
+     &   , graybody                                                     &
      &   , Tstrat                                                       &
      &   , newtonian                                                    &
@@ -45,4 +45,5 @@
      &   , n2mixratio                                                   &
      &   , co2supsat                                                    &
+     &   , cloudlvl                                                     &
      &   , pceil                                                    
 
@@ -73,4 +74,5 @@
       logical CLFvarying
       logical noradsurf
+      logical graybody
 
       integer iddist
@@ -97,4 +99,4 @@
       real maxicethick
       real Tsaldiff
-      real Tsurfk
       real tau_relax
+      real cloudlvl
Index: /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 305)
@@ -0,0 +1,535 @@
+      subroutine condense_cloud(ngrid,nlayer,nq,ptimestep, &
+          pcapcal,pplay,pplev,ptsrf,pt, &
+          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
+          piceco2,psolaralb,pemisurf, &
+          pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
+          pdqc,reffrad,cpp3D)
+
+      use radinc_h, only : naerkind
+
+      implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Condense and/or sublime CO2 ice on the ground and in the
+!     atmosphere, and sediment the ice.
+!  
+!     Inputs
+!     ------
+!     ngrid                 Number of vertical columns
+!     nlayer                Number of layers
+!     pplay(ngrid,nlayer)   Pressure layers
+!     pplev(ngrid,nlayer+1) Pressure levels 
+!     pt(ngrid,nlayer)      Temperature (in K)
+!     ptsrf(ngrid)          Surface temperature
+!     
+!     pdt(ngrid,nlayermx)   Time derivative before condensation/sublimation of pt
+!     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf
+!     pqsurf(ngrid,nq)      Sedimentation flux at the surface (kg.m-2.s-1)
+!     
+!     Outputs
+!     -------
+!     pdpsrf(ngrid)         \  Contribution of condensation/sublimation
+!     pdtc(ngrid,nlayermx)  /  to the time derivatives of Ps, pt, and ptsrf
+!     pdtsrfc(ngrid)       /
+!     
+!     Both
+!     ----
+!     piceco2(ngrid)        CO2 ice at the surface (kg/m2)
+!     psolaralb(ngrid)      Albedo at the surface
+!     pemisurf(ngrid)       Emissivity of the surface
+!
+!     Authors
+!     ------- 
+!     Francois Forget (1996)
+!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
+!     
+!==================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "surfdat.h"
+#include "comgeomfi.h"
+#include "comvert.h"
+#include "callkeys.h"
+#include "tracer.h"
+#include "gases.h"
+
+!-----------------------------------------------------------------------
+!     Arguments
+
+      INTEGER ngrid, nlayer, nq
+
+      REAL ptimestep 
+      REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
+      REAL pcapcal(ngrid)
+      REAL pt(ngrid,nlayer)
+      REAL ptsrf(ngrid)
+      REAL pphi(ngrid,nlayer)
+      REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
+      REAL pdtsrfc(ngrid),pdpsrf(ngrid)
+!      REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
+      REAL piceco2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
+
+      REAL pu(ngrid,nlayer) ,  pv(ngrid,nlayer)
+      REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
+      REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
+      REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
+      REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq)
+
+      REAL reffrad(ngrid,nlayer,naerkind)
+
+!     Allow variations in cp with location
+      REAL cpp3D(ngridmx,nlayermx)   ! specific heat capacity at const. pressure
+
+
+
+!-----------------------------------------------------------------------
+!     Local variables
+
+      INTEGER l,ig,icap,ilay,it,iq
+
+      REAL*8 zt(ngridmx,nlayermx)
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL zcpi
+      REAL ztcond (ngridmx,nlayermx)
+      REAL ztcondsol(ngridmx) 
+      REAL zdiceco2(ngridmx)
+      REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
+      REAL zfallice(ngridmx), Mfallice(ngridmx) 
+      REAL zmflux(nlayermx+1)
+      REAL zu(nlayermx),zv(nlayermx)
+      REAL ztsrf(ngridmx)
+      REAL ztc(nlayermx), ztm(nlayermx+1) 
+      REAL zum(nlayermx+1) , zvm(nlayermx+1)
+      LOGICAL condsub(ngridmx)
+      REAL subptimestep
+      Integer Ntime
+      real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx)
+      real wq(ngridmx,nlayermx+1)
+      real vstokes,reff
+
+!     Special diagnostic variables
+      real tconda1(ngridmx,nlayermx)
+      real tconda2(ngridmx,nlayermx)
+      real zdtsig (ngridmx,nlayermx)
+      real zdt (ngridmx,nlayermx)
+
+!-----------------------------------------------------------------------
+!     Saved local variables
+
+      REAL emisref(ngridmx)
+      REAL latcond
+      REAL ccond
+      REAL cpice
+      SAVE emisref, cpice
+      SAVE latcond,ccond
+
+      LOGICAL firstcall
+      SAVE firstcall
+      REAL SSUM
+      EXTERNAL SSUM
+
+      DATA latcond /5.9e5/
+      DATA cpice /1000./
+      DATA firstcall/.true./
+
+      REAL CBRT
+      EXTERNAL CBRT
+
+      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
+      CHARACTER(LEN=20) :: tracername ! to temporarily store text
+
+      integer igas
+      integer,save :: igasco2=0
+      character(len=3) :: gasname
+
+      real reffradmin, reffradmax
+
+      real ppco2
+
+!-----------------------------------------------------------------------
+!     Initializations
+
+      call zerophys(ngrid*nlayer*nq,pdqc)
+      call zerophys(ngrid*nlayer*nq,pdtc)
+      call zerophys(ngridmx*nlayermx*nqmx,zq)
+      call zerophys(ngridmx*nlayermx,zt)
+
+      !reffradmin=1.e-7
+      !reffradmax=5.e-4
+      !reffradmin=0.5e-7
+      !reffradmax=0.1e-3 ! FF data
+      reffradmin=0.1e-5
+      reffradmax=0.1e-3 ! JB data
+!     improve this in the future...
+
+!     Initialisation
+      IF (firstcall) THEN
+
+         ! find CO2 ice tracer
+         do iq=1,nqmx
+            tracername=noms(iq)
+            if (tracername.eq."co2_ice") then
+               i_co2ice=iq
+            endif
+         enddo
+         
+         ! find CO2 gas
+         do igas=1,ngasmx
+            gasname=gnom(igas)
+!            gasname=noms(igas) ! was a bug
+            if (gasname.eq."CO2") then
+               igasco2=igas
+            endif
+         enddo
+
+        write(*,*) "condense_co2cloud: i_co2ice=",i_co2ice       
+
+        if((i_co2ice.lt.1))then
+           print*,'In condens_co2cloud but no CO2 ice tracer, exiting.'
+           stop
+        endif
+
+         ccond=cpp/(g*latcond)
+         print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond
+
+!          Prepare special treatment if gas is not pure CO2
+             !if (addn2) then
+             !   m_co2   = 44.01E-3 ! CO2 molecular mass (kg/mol)   
+             !   m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)  
+!               Compute A and B coefficient use to compute
+!               mean molecular mass Mair defined by
+!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
+!               1/Mair = A*q(ico2) + B
+             !   A = (1/m_co2 - 1/m_noco2)
+             !   B = 1/m_noco2
+             !endif
+
+!          Minimum CO2 mixing ratio below which mixing occurs with layer above: 
+           !qco2min =0.75  
+
+           firstcall=.false.
+      ENDIF
+      zcpi=1./cpp
+
+!-----------------------------------------------------------------------
+!     Calculation of CO2 condensation / sublimation 
+!
+!     Variables used:
+!     piceco2(ngrid)       amount of co2 ice on the ground  (kg/m2)
+!     zcondicea(ngrid,l)   condensation rate in layer l     (kg/m2/s)
+!     zcondices(ngrid)     condensation rate on the ground  (kg/m2/s)
+!     zfallice(ngrid)      flux of ice falling on surface   (kg/m2/s)
+!     pdtc(ngrid,nlayermx) dT/dt due to phase changes       (K/s)
+     
+
+!     Tendencies initially set to 0 (except pdtc)
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zcondicea(ig,l) = 0.
+            pduc(ig,l)  = 0
+            pdvc(ig,l)  = 0
+            pdqc(ig,l,i_co2ice)  = 0
+         END DO
+      END DO
+      
+      DO ig=1,ngrid
+         Mfallice(ig) = 0.
+         zfallice(ig) = 0.
+         zcondices(ig) = 0.
+         pdtsrfc(ig) = 0.
+         pdpsrf(ig) = 0.
+         condsub(ig) = .false.
+         zdiceco2(ig) = 0.
+      ENDDO
+
+!-----------------------------------------------------------------------
+!     Atmospheric condensation
+
+
+!     Compute CO2 Volume mixing ratio
+!     -------------------------------
+!      if (addn2) then
+!         DO l=1,nlayer
+!            DO ig=1,ngrid
+!              qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
+!!             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
+!              mmean=1/(A*qco2 +B)
+!              vmr_co2(ig,l) = qco2*mmean/m_co2 
+!            ENDDO
+!         ENDDO
+!      else
+!         DO l=1,nlayer
+!            DO ig=1,ngrid
+!              vmr_co2(ig,l)=0.5
+!            ENDDO
+!         ENDDO
+!      end if
+
+!     Forecast the atmospheric frost temperature 'ztcond'
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            ppco2=gfrac(igasco2)*pplay(ig,l)
+            call get_tcond_co2(ppco2,ztcond(ig,l))
+         ENDDO
+      ENDDO
+      
+!     Initialize zq and zt at the beginning of the sub-timestep loop
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zt(ig,l)=pt(ig,l)
+            zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
+            IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
+               print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
+               if(l.eq.1)then
+                  print*,'Perhaps the atmosphere is collapsing on surface...?'
+               endif
+            END IF
+         ENDDO
+      ENDDO
+
+!     Calculate the mass of each atmospheric layer (kg.m-2)
+      do  ilay=1,nlayer
+         do ig=1, ngrid
+            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
+         end do
+      end do
+
+!     -----------------------------------------------
+!     START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
+!     -----------------------------------------------
+      Ntime =  20               ! number of sub-timestep 
+      subptimestep = ptimestep/float(Ntime)           
+
+      DO it=1,Ntime
+
+!     Add the tendencies from other physical processes at each subtimstep
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zt(ig,l)   = zt(ig,l)   + pdt(ig,l)   * subptimestep
+               zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
+            END DO
+         END DO
+
+
+!     Gravitational sedimentation
+            
+!     sedimentation computed from radius computed from q
+!     assuming that the ice is splitted in Nmix particle
+         do  ilay=1,nlayer
+            do ig=1, ngrid
+
+               reff = CBRT( 3*zq(ig,ilay,i_co2ice)/( 4*Nmix_co2*pi*rho_co2 ))
+
+               ! there should be a more elegant way of doing this...
+               if(reff.lt.1.e-16) reff=1.e-16
+
+               ! update reffrad for radiative transfer
+               if(reff.lt.reffradmin)then
+                  reffrad(ig,ilay,1)=reffradmin
+                  !print*,'reff below optical limit'
+               elseif(reff.gt.reffradmax)then
+                  reffrad(ig,ilay,1)=reffradmax
+                  !print*,'reff above optical limit'
+               else
+                  reffrad(ig,ilay,1)=reff
+               endif
+
+               call stokes                      &
+                   (pplev(ig,ilay),pt(ig,ilay), &
+                    reff,vstokes,rho_co2)
+
+               !w(ig,ilay,i_co2ice) = 0.0
+               w(ig,ilay,i_co2ice) = vstokes *  subptimestep * &
+                   pplev(ig,ilay)/(r*pt(ig,ilay))
+
+            end do
+         end do
+
+!     Computing q after sedimentation
+
+         call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
+
+
+!     Progressively accumulating the flux to the ground
+!     Mfallice is the total amount of ice fallen to the ground
+         do ig=1,ngrid
+            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_co2ice)
+         end do
+
+
+!     Condensation / sublimation in the atmosphere
+!     --------------------------------------------
+!     (calculation of zcondicea, zfallice and pdtc)
+!     (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
+!     of CO2 into tracer i_co2ice)
+
+         DO l=nlayer , 1, -1
+            DO ig=1,ngrid
+               pdtc(ig,l)=0.
+
+               ! tweak ccond if the gas is non-ideal
+               if(nonideal)then    
+                  ccond=cpp3D(ig,l)/(g*latcond)
+               endif
+
+
+               IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN
+                  condsub(ig)=.true.
+                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
+                  pdqc(ig,l,i_co2ice) = pdtc(ig,l)*ccond*g
+
+!     Case when the ice from above sublimes entirely
+                  IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
+                      .AND. (zq(ig,l,i_co2ice).gt.0)) THEN
+
+                     pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
+                     pdtc(ig,l)   =-zq(ig,l,i_co2ice)/(ccond*g*subptimestep)
+
+                  END IF
+
+!     Temperature and q after condensation
+                  zt(ig,l)   = zt(ig,l)   + pdtc(ig,l)   * subptimestep
+                  zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
+               END IF
+
+            ENDDO
+         ENDDO
+      ENDDO                     ! end of subtimestep loop
+
+!     Computing global tendencies after the subtimestep
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            pdtc(ig,l) = &
+              (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
+            pdqc(ig,l,i_co2ice) = &
+              (zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
+         END DO
+      END DO
+      DO ig=1,ngrid
+         zfallice(ig) = Mfallice(ig)/ptimestep
+      END DO
+
+
+!-----------------------------------------------------------------------
+!     Condensation/sublimation on the ground
+!     (calculation of zcondices and pdtsrfc)
+
+!     forecast of ground temperature ztsrf and frost temperature ztcondsol
+      DO ig=1,ngrid
+         ppco2=gfrac(igasco2)*pplay(ig,1)
+         call get_tcond_co2(ppco2,ztcondsol(ig))
+         
+         ztsrf(ig) = ptsrf(ig)
+
+         if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
+            print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
+            print*,'T_surf = ',ztsrf,'K'
+            print*,'T_cond = ',ztcondsol,'K'
+            open(116,file='surf_vals.out')
+            write(116,*) 0.0, pplev(1,1), 0.0, 0.0
+            close(116)
+            call abort
+         endif
+
+         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
+
+      ENDDO
+     
+      DO ig=1,ngrid
+         IF(ig.GT.ngrid/2+1) THEN
+            icap=2
+         ELSE
+            icap=1
+         ENDIF
+
+!     Loop over where we have condensation / sublimation
+         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &        ! ground condensation
+                    (zfallice(ig).NE.0.) .OR.              &        ! falling snow
+                    ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &        ! ground sublimation
+                    ((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
+            condsub(ig) = .true.
+
+!     Condensation or partial sublimation of CO2 ice
+            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
+                /(latcond*ptimestep)         
+            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
+
+!     If the entire CO_2 ice layer sublimes
+!     (including what has just condensed in the atmosphere)
+            IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
+                -zcondices(ig))THEN
+               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
+               pdtsrfc(ig)=(latcond/pcapcal(ig))*       &
+                   (zcondices(ig))
+            END IF
+
+!     Changing CO2 ice amount and pressure
+
+            zdiceco2(ig) = zcondices(ig) + zfallice(ig)
+            piceco2(ig)  = piceco2(ig) + zdiceco2(ig)*ptimestep
+            pdpsrf(ig)   = -zdiceco2(ig)*g
+
+            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
+               PRINT*,'STOP in condens'
+               PRINT*,'condensing more than total mass'
+               PRINT*,'Grid point ',ig
+               PRINT*,'Ps = ',pplev(ig,1)
+               PRINT*,'d Ps = ',pdpsrf(ig)
+               STOP
+            ENDIF
+         END IF
+      ENDDO
+
+!     Surface albedo and emissivity of the ground below the snow (emisref)
+!     --------------------------------------------------------------------
+      do ig =1,ngrid
+         IF(ig.GT.ngrid/2+1) THEN
+            icap=2
+         ELSE
+            icap=1
+         ENDIF
+
+         if(.not.piceco2(ig).ge.0.) THEN
+            if(piceco2(ig).le.-1.e-8) print*,   &
+                'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
+            piceco2(ig)=0.
+         endif
+         if (piceco2(ig).gt.0) then
+            psolaralb(ig) = albedice(icap)
+            emisref(ig)   = emisice(icap)
+         else
+            psolaralb(ig) = albedodat(ig)
+            emisref(ig)   = emissiv
+            pemisurf(ig)  = emissiv
+         end if
+      end do
+
+      return
+    end subroutine condense_cloud
+
+!-------------------------------------------------------------------------
+    subroutine get_tcond_co2(p,tcond)
+!   Calculates the condensation temperature for CO2
+
+      implicit none
+
+#include "callkeys.h"
+
+      real p, peff, tcond
+      real, parameter :: ptriple=518000.0
+
+      peff=p!/co2supsat
+
+      if(peff.lt.ptriple)then
+         tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
+      else
+         tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 
+         ! liquid-vapour transition (based on CRC handbook 2003 data)
+      endif
+      return
+
+    end subroutine get_tcond_co2
Index: /trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90	(revision 305)
@@ -0,0 +1,33 @@
+double precision function cp_neutral(T)
+
+  implicit none
+
+#include "gases.h"
+
+  ! inputs
+  double precision T
+
+
+  ! Cp_n : cf CO2 dans abe&matsui (1988)
+  !cp_neutral = (22.26+5.981d-2*T-3.501d-5*T**2+7.469d-9*T**3)/m_n
+
+  if(gnom(1).eq.'N2_')then
+     cp_neutral = 1040.0
+  elseif(gnom(1).eq.'H2_')then
+     cp_neutral = 14310.0
+  else
+     print*,'Gas not recognised in cp_neutral!'
+     call abort
+  endif
+
+  
+end function cp_neutral
+
+
+
+
+
+
+
+
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/gases.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/gases.h	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/gases.h	(revision 305)
@@ -8,5 +8,5 @@
       real gfrac(ngasmx)
 
-      COMMON/gases/gnom,gfrac,vgas 
+      COMMON/gases/gnom,vgas,gfrac
 
 !-----------------------------------------------------------------------
Index: /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 305)
@@ -0,0 +1,121 @@
+subroutine gradients_kcm(profil_flag,rho_v,rho_n,T,dTdp,dPvdp,dPndp)
+
+  use params_h
+  implicit none
+
+#include "gases.h"
+
+  ! inputs
+  integer profil_flag ! 0 = dry, 1 = moist, 2 = isothermal
+  double precision rho_v,rho_n,T
+
+  ! outputs
+  double precision dTdp,dPndp,dPvdp
+  double precision a_v
+
+  ! internal
+  double precision cp_n,cp_v
+
+  double precision press, rho_plus, rho_minus, dVdT, rho_c
+  double precision dlnr,dlna,dpsat,dsv
+  double precision s_minus,s_plus
+  double precision s_v,s_c,L
+  double precision psat_plus,psat_minus,Pn
+  double precision nul
+
+  ! functions
+  double precision cp_neutral
+
+  cp_n = cp_neutral(T)
+
+
+  select case(profil_flag)
+  case(2) ! isothermal
+
+     dTdp  = 0.
+     a_v   = rho_v/rho_n ! constant here
+     dPndp = 1/(1d0+m_n/m_v*a_v)
+     dPvdp = 1 - dPndp
+
+  case(1) ! moist
+
+     Pn = rho_n*T*rmn
+
+     if(gnom(2).eq.'H2O')then
+
+        call psat_H2O(T-2d-1,psat_minus)
+        call psat_H2O(T+2d-1,psat_plus)
+        call psat_H2O(T,press)
+
+        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
+        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
+
+        call therm(T-2d-1,rho_minus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
+                   nul,nul,press,s_minus,nul)
+        call therm(T+2d-1,rho_plus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
+                   nul,nul,press,s_plus,nul)
+        s_c = 2.06 * log(T/273.15)
+
+        s_plus  = s_plus  * 1d3
+        s_minus = s_minus * 1d3
+        s_c     = s_c     * 1d3 ! convert to SI
+
+        if(T.lt.280.0)then 
+           dpsat = press*1d6 * ( 1730.63*log(10.) / (T-39.714)**2 )
+        else
+           call tdpsdt(T,dpsat)
+           dpsat = dpsat * 1d6 / T
+        endif
+
+     elseif(gnom(2).eq.'NH3')then
+
+        call psat_NH3(T-2d-1,psat_minus)
+        call psat_NH3(T+2d-1,psat_plus)
+        call psat_NH3(T,press)
+
+        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
+        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
+
+        call latheat_NH3(T-2d-1,nul,s_minus)
+        call latheat_NH3(T+2d-1,nul,s_plus)
+        call latheat_NH3(T,s_c,nul)
+
+        dpsat = press*1d6 * (-2*1.5609d-4*T + 0.1236)
+
+     endif
+
+     dsv  = (s_plus-s_minus)/4d-1  ! dsv*T = ds / d ln[T]
+     s_v  = (s_plus+s_minus)/2d0
+     dlnr = T/rho_v * (rho_plus-rho_minus)/4d-1  ! d ln[rho_v] / d ln[T]
+
+     if(rho_n/rho_v.lt.1e-5)then
+        dlna = -T*dsv/(s_v-s_c)
+     else
+        a_v  = rho_v/rho_n
+        dlna = (rmn*dlnr - cp_n + rmn - a_v*T*dsv)/(a_v*(s_v-s_c)+rmn)  
+        ! d ln[alpha_v] / d ln[T]
+        ! note cp_n + rmn = cv_n, which is what's required
+     endif
+
+     dTdp  = 1d0 / (dpsat + rho_n*rmn*(1d0 + dlnr - dlna)) ! c.f. Marcq S2.2.2
+     dPvdp = dTdp * dpsat
+     dPndp = 1d0 - dPvdp ! from p = p_v + p_n
+
+  case(0) ! dry
+
+     if(gnom(2).eq.'H2O')then
+        cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v
+     elseif(gnom(2).eq.'NH3')then
+        cp_v = 2.058d3
+     elseif(gnom(2).eq.'CH4')then
+        cp_v = 2.226d3
+     endif
+
+     dTdp  = 1/(rho_n*cp_n+rho_v*cp_v)
+     dPndp = 1/(1d0+m_n/m_v*rho_v/rho_n)
+     dPvdp = 1/(1d0+m_v/m_n*rho_n/rho_v)
+      
+  end select
+
+
+end subroutine gradients_kcm
Index: /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 305)
@@ -53,5 +53,5 @@
       parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm 
       real oceantime
-      parameter (oceantime=1*3600)
+      parameter (oceantime=10*24*3600)
 
       logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
@@ -101,5 +101,5 @@
 
       oceanbulkavg=.false.
-      activerunoff=.false.
+      activerunoff=.true.
       oceanalbvary=.false.
 
@@ -311,5 +311,5 @@
          END DO
 
-         print*,'Mean ocean temperature = ',tsea
+         print*,'Mean ocean temperature               = ',tsea,' K'
 
       endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 305)
@@ -13,5 +13,5 @@
 !
 !   Initialisation for the physical parametrisations of the LMD 
-!   martian atmospheric general circulation modele.
+!   Generic Model.
 !
 !   author: Frederic Hourdin 15 / 10 /93
@@ -56,4 +56,5 @@
 
 
+
       REAL prad,pg,pr,pcpp,pdaysec,ptimestep
  
@@ -265,8 +266,8 @@
          write(*,*)" kastprof = ",kastprof
 
-         write(*,*)"Surface temperature in kastprof mode?"
-         Tsurfk=273.15
-         call getin("Tsurfk",Tsurfk)
-         write(*,*)" Tsurfk = ",Tsurfk
+         write(*,*)"Uniform absorption coefficient in IR?"
+         graybody=.false.
+         call getin("graybody",graybody)
+         write(*,*)" graybody = ",graybody
 
 ! Test of incompatibility:
@@ -344,8 +345,4 @@
          write(*,*)" Fat1AU = ",Fat1AU
 
-         write(*,*)"Set temperature to 1 K above CO2 condensation?"
-         nearco2cond=.false.
-         call getin("nearco2cond",nearco2cond)
-         write(*,*)" nearco2cond = ",nearco2cond
 
 !     1D solar zenith angle
@@ -394,5 +391,11 @@
          endif
 
+         write(*,*)"Cloud pressure level (with kastprof only):"
+         cloudlvl=0. ! default value
+         call getin("cloudlvl",cloudlvl)
+         write(*,*)" cloudlvl = ",cloudlvl
+
          write(*,*)"Is the variable gas species radiatively active?"
+         Tstrat=167.0
          varactive=.false.
          call getin("varactive",varactive)
Index: unk/LMDZ.GENERIC/libf/phystd/iniwrite_spec.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_spec.F	(revision 304)
+++ 	(revision )
@@ -1,391 +1,0 @@
-      SUBROUTINE iniwrite_spec(nid,idayref,phis)
-
-      use radinc_h, only: L_NSPECTI
-      use radcommon_h, only: WAVEI
-
-      implicit none
-
-c=======================================================================
-c
-c   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   'Initialize' the diagfi_spec.nc file: write down dimensions as well
-c   as time-independent fields (e.g: geopotential, mesh area, ...)
-c
-c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-!include "comconst.h"
-#include "comcstfi.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "temps.h"
-#include "ener.h"
-#include "logic.h"
-#include "netcdf.inc"
-#include "description.h"
-#include "serre.h"
-#include"dimphys.h"
-#include"comsoil.h"
-
-c   Arguments:
-c   ----------
-
-      integer nid        ! NetCDF file ID
-      INTEGER*4 idayref  ! date (initial date for this run)
-      REAL phis(ip1jmp1) ! surface geopotential
-
-c   Local:
-c   ------
-      INTEGER length,l
-      parameter (length = 100)
-      REAL tab_cntrl(length) ! run parameters are stored in this array
-      INTEGER ierr
-
-      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
-      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
- !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
-      integer :: idim_bandsIR ! "bandsIR" dimension ID #
-      integer, dimension(2) :: id  
-
-c-----------------------------------------------------------------------
-
-      DO l=1,length
-         tab_cntrl(l)=0.
-      ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
-      tab_cntrl(4)  = FLOAT(idayref)
-      tab_cntrl(5)  = rad
-      tab_cntrl(6)  = omeg
-      tab_cntrl(7)  = g
-      tab_cntrl(8)  = mugaz
-      tab_cntrl(9)  = rcp
-      tab_cntrl(10) = daysec
-      tab_cntrl(11) = dtphys
-      tab_cntrl(12) = etot0
-      tab_cntrl(13) = ptot0
-      tab_cntrl(14) = ztot0
-      tab_cntrl(15) = stot0
-      tab_cntrl(16) = ang0
-c
-c    ..........    P.Le Van  ( ajout le 8/04/96 )    .........
-c         .....        parametres  pour le zoom          ......   
-      tab_cntrl(17)  = clon
-      tab_cntrl(18)  = clat
-      tab_cntrl(19)  = grossismx
-      tab_cntrl(20)  = grossismy
-c
-c     .....   ajout  le 6/05/97 et le 15/10/97  .......
-c
-      IF ( fxyhypb )   THEN
-        tab_cntrl(21) = 1.
-        tab_cntrl(22) = dzoomx
-        tab_cntrl(23) = dzoomy
-      ELSE
-        tab_cntrl(21) = 0.
-        tab_cntrl(22) = dzoomx
-        tab_cntrl(23) = dzoomy
-        tab_cntrl(24) = 0.
-        IF( ysinus )  tab_cntrl(24) = 1.
-      ENDIF
-
-c    .........................................................
-
-! Define dimensions
-    
-      ierr = NF_REDEF (nid)
-
-      ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-!      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
-!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
-!      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
-!      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
-      ierr = NF_DEF_DIM (nid, "bandsIR",L_NSPECTI,idim_bandsIR)
-
-      ierr = NF_ENDDEF(nid)
-
-c  Contol parameters for this run
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, 
-     .       idim_index,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, 
-     .       idim_index,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
-     .                       "Control parameters")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
-#endif
-
-c --------------------------
-c  longitudes and latitudes
-
-
-c$$$      ierr = NF_REDEF (nid)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
-c$$$#else
-c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
-c$$$#endif
-c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
-c$$$     .                       "Longitudes at u nodes")
-c$$$      ierr = NF_ENDDEF(nid)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
-c$$$#endif
-
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr =NF_DEF_VAR(nid, "latitude", NF_DOUBLE, 1, idim_rlatu,nvarid)
-#else
-      ierr =NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim_rlatu,nvarid)
-#endif
-      ierr =NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north")
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
-     .      "North latitude")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr =NF_DEF_VAR(nid,"longitude", NF_DOUBLE, 1, idim_rlonv,nvarid)
-#else
-      ierr = NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, idim_rlonv,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
-     .      "East longitude")
-      ierr = NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
-#endif
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, 
-     .       idim_llm,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, 
-     .       idim_llm,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",10,"pseudo-alt")
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
-
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
-#endif
-
-
-c$$$c --------------------------
-c$$$      ierr = NF_REDEF (nid)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
-c$$$#else
-c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
-c$$$#endif
-c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
-c$$$     .                       "Latitudes at v nodes")
-c$$$      ierr = NF_ENDDEF(nid)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
-c$$$#endif
-
-
-c$$$c --------------------------
-c$$$c  Vertical levels
-c$$$      call def_var(nid,"aps","hybrid pressure at midlayers ","Pa",
-c$$$     .            1,idim_llm,nvarid,ierr)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
-c$$$#endif
-c$$$
-c$$$      call def_var(nid,"bps","hybrid sigma at midlayers"," ",
-c$$$     .            1,idim_llm,nvarid,ierr)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
-c$$$#endif
-c$$$
-c$$$      call def_var(nid,"ap","hybrid pressure at interlayers","Pa",
-c$$$     .            1,idim_llmp1,nvarid,ierr)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
-c$$$#endif
-c$$$
-c$$$      call def_var(nid,"bp","hybrid sigma at interlayers"," ",
-c$$$     .            1,idim_llmp1,nvarid,ierr)
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
-c$$$#else
-c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
-c$$$#endif
-
-c$$$!-------------------------------
-c$$$! (soil) depth variable mlayer() (known from comsoil.h)
-c$$$!-------------------------------
-c$$$      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
-c$$$      ! define variable
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_DOUBLE,1,idim_nsoilmx,nvarid)
-c$$$#else
-c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_FLOAT,1,idim_nsoilmx,nvarid)
-c$$$#endif
-c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20,
-c$$$     .                        "Soil mid-layer depth")
-c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"m")
-c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"positive",4,"down")
-c$$$      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
-c$$$      ! write variable
-c$$$#ifdef NC_DOUBLE
-c$$$      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
-c$$$#else
-c$$$      ierr=NF_PUT_VAR_REAL (nid,nvarid,mlayer)
-c$$$#endif
-
-!-------------------------------
-! Number of bands in the IR
-!-------------------------------
-
-      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
-      ! define variable
-#ifdef NC_DOUBLE
-      ierr=NF_DEF_VAR(nid,"bandsIR",NF_DOUBLE,1,idim_bandsIR,nvarid)
-#else
-      ierr=NF_DEF_VAR(nid,"bandsIR",NF_FLOAT,1,idim_bandsIR,nvarid)
-#endif
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
-     .                        "Band limits in the infrared")
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
-      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
-      ! write variable
-#ifdef NC_DOUBLE
-      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEI))
-#else
-      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEI))
-#endif
-
-c
-c --------------------------
-c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
-
-      id(1)=idim_rlonu
-      id(2)=idim_rlatu
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
-#endif
-c
-      id(1)=idim_rlonv
-      id(2)=idim_rlatv
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
-#endif
-c
-      id(1)=idim_rlonv
-      id(2)=idim_rlatu
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "aire", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "aire", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
-     .                       "Mesh area")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
-#endif
-
-c  Surface geopotential
-      id(1)=idim_rlonv
-      id(2)=idim_rlatu
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "phisinit", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "phisinit", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
-     .                       "Geopotential at the surface")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
-#endif,
-
-
-!      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
-!      write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',rad,omeg,g,mugaz,rcp
-!      write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
-
-      END
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 305)
@@ -0,0 +1,391 @@
+      SUBROUTINE iniwrite_specIR(nid,idayref,phis)
+
+      use radinc_h, only: L_NSPECTI
+      use radcommon_h, only: WAVEI
+
+      implicit none
+
+c=======================================================================
+c
+c   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
+c   -------
+c
+c   Objet:
+c   ------
+c
+c   'Initialize' the diagfi_spec.nc file: write down dimensions as well
+c   as time-independent fields (e.g: geopotential, mesh area, ...)
+c
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+#include "dimensions.h"
+#include "paramet.h"
+!include "comconst.h"
+#include "comcstfi.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "temps.h"
+#include "ener.h"
+#include "logic.h"
+#include "netcdf.inc"
+#include "description.h"
+#include "serre.h"
+#include"dimphys.h"
+#include"comsoil.h"
+
+c   Arguments:
+c   ----------
+
+      integer nid        ! NetCDF file ID
+      INTEGER*4 idayref  ! date (initial date for this run)
+      REAL phis(ip1jmp1) ! surface geopotential
+
+c   Local:
+c   ------
+      INTEGER length,l
+      parameter (length = 100)
+      REAL tab_cntrl(length) ! run parameters are stored in this array
+      INTEGER ierr
+
+      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
+      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
+ !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
+      integer :: idim_bandsIR ! "bandsIR" dimension ID #
+      integer, dimension(2) :: id  
+
+c-----------------------------------------------------------------------
+
+      DO l=1,length
+         tab_cntrl(l)=0.
+      ENDDO
+      tab_cntrl(1)  = FLOAT(iim)
+      tab_cntrl(2)  = FLOAT(jjm)
+      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(4)  = FLOAT(idayref)
+      tab_cntrl(5)  = rad
+      tab_cntrl(6)  = omeg
+      tab_cntrl(7)  = g
+      tab_cntrl(8)  = mugaz
+      tab_cntrl(9)  = rcp
+      tab_cntrl(10) = daysec
+      tab_cntrl(11) = dtphys
+      tab_cntrl(12) = etot0
+      tab_cntrl(13) = ptot0
+      tab_cntrl(14) = ztot0
+      tab_cntrl(15) = stot0
+      tab_cntrl(16) = ang0
+c
+c    ..........    P.Le Van  ( ajout le 8/04/96 )    .........
+c         .....        parametres  pour le zoom          ......   
+      tab_cntrl(17)  = clon
+      tab_cntrl(18)  = clat
+      tab_cntrl(19)  = grossismx
+      tab_cntrl(20)  = grossismy
+c
+c     .....   ajout  le 6/05/97 et le 15/10/97  .......
+c
+      IF ( fxyhypb )   THEN
+        tab_cntrl(21) = 1.
+        tab_cntrl(22) = dzoomx
+        tab_cntrl(23) = dzoomy
+      ELSE
+        tab_cntrl(21) = 0.
+        tab_cntrl(22) = dzoomx
+        tab_cntrl(23) = dzoomy
+        tab_cntrl(24) = 0.
+        IF( ysinus )  tab_cntrl(24) = 1.
+      ENDIF
+
+c    .........................................................
+
+! Define dimensions
+    
+      ierr = NF_REDEF (nid)
+
+      ierr = NF_DEF_DIM (nid, "index", length, idim_index)
+!      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+!      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+!      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
+      ierr = NF_DEF_DIM (nid, "bandsIR",L_NSPECTI,idim_bandsIR)
+
+      ierr = NF_ENDDEF(nid)
+
+c  Contol parameters for this run
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, 
+     .       idim_index,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, 
+     .       idim_index,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
+     .                       "Control parameters")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+
+c --------------------------
+c  longitudes and latitudes
+
+
+c$$$      ierr = NF_REDEF (nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
+c$$$#else
+c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
+c$$$#endif
+c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+c$$$     .                       "Longitudes at u nodes")
+c$$$      ierr = NF_ENDDEF(nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
+c$$$#endif
+
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr =NF_DEF_VAR(nid, "latitude", NF_DOUBLE, 1, idim_rlatu,nvarid)
+#else
+      ierr =NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim_rlatu,nvarid)
+#endif
+      ierr =NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "North latitude")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+#endif
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr =NF_DEF_VAR(nid,"longitude", NF_DOUBLE, 1, idim_rlonv,nvarid)
+#else
+      ierr = NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, idim_rlonv,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "East longitude")
+      ierr = NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+#endif
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, 
+     .       idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, 
+     .       idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",10,"pseudo-alt")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif
+
+
+c$$$c --------------------------
+c$$$      ierr = NF_REDEF (nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
+c$$$#else
+c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
+c$$$#endif
+c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+c$$$     .                       "Latitudes at v nodes")
+c$$$      ierr = NF_ENDDEF(nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
+c$$$#endif
+
+
+c$$$c --------------------------
+c$$$c  Vertical levels
+c$$$      call def_var(nid,"aps","hybrid pressure at midlayers ","Pa",
+c$$$     .            1,idim_llm,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"bps","hybrid sigma at midlayers"," ",
+c$$$     .            1,idim_llm,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"ap","hybrid pressure at interlayers","Pa",
+c$$$     .            1,idim_llmp1,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"bp","hybrid sigma at interlayers"," ",
+c$$$     .            1,idim_llmp1,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
+c$$$#endif
+
+c$$$!-------------------------------
+c$$$! (soil) depth variable mlayer() (known from comsoil.h)
+c$$$!-------------------------------
+c$$$      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+c$$$      ! define variable
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_DOUBLE,1,idim_nsoilmx,nvarid)
+c$$$#else
+c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_FLOAT,1,idim_nsoilmx,nvarid)
+c$$$#endif
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20,
+c$$$     .                        "Soil mid-layer depth")
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"m")
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"positive",4,"down")
+c$$$      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+c$$$      ! write variable
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
+c$$$#else
+c$$$      ierr=NF_PUT_VAR_REAL (nid,nvarid,mlayer)
+c$$$#endif
+
+!-------------------------------
+! Number of bands in the IR
+!-------------------------------
+
+      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+      ! define variable
+#ifdef NC_DOUBLE
+      ierr=NF_DEF_VAR(nid,"bandsIR",NF_DOUBLE,1,idim_bandsIR,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"bandsIR",NF_FLOAT,1,idim_bandsIR,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
+     .                        "Band limits in the infrared")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
+      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+      ! write variable
+#ifdef NC_DOUBLE
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEI))
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEI))
+#endif
+
+c
+c --------------------------
+c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
+
+      id(1)=idim_rlonu
+      id(2)=idim_rlatu
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+     .             "Conversion coefficients cov <--> natural")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+#endif
+c
+      id(1)=idim_rlonv
+      id(2)=idim_rlatv
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+     .             "Conversion coefficients cov <--> natural")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+#endif
+c
+      id(1)=idim_rlonv
+      id(2)=idim_rlatu
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "aire", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "aire", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
+     .                       "Mesh area")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+#endif
+
+c  Surface geopotential
+      id(1)=idim_rlonv
+      id(2)=idim_rlatu
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "phisinit", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "phisinit", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
+     .                       "Geopotential at the surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
+#endif,
+
+
+!      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+!      write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',rad,omeg,g,mugaz,rcp
+!      write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
+
+      END
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 305)
@@ -0,0 +1,391 @@
+      SUBROUTINE iniwrite_specVI(nid,idayref,phis)
+
+      use radinc_h, only: L_NSPECTV
+      use radcommon_h, only: WAVEV
+
+      implicit none
+
+c=======================================================================
+c
+c   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
+c   -------
+c
+c   Objet:
+c   ------
+c
+c   'Initialize' the diagfi_spec.nc file: write down dimensions as well
+c   as time-independent fields (e.g: geopotential, mesh area, ...)
+c
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+#include "dimensions.h"
+#include "paramet.h"
+!include "comconst.h"
+#include "comcstfi.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "temps.h"
+#include "ener.h"
+#include "logic.h"
+#include "netcdf.inc"
+#include "description.h"
+#include "serre.h"
+#include"dimphys.h"
+#include"comsoil.h"
+
+c   Arguments:
+c   ----------
+
+      integer nid        ! NetCDF file ID
+      INTEGER*4 idayref  ! date (initial date for this run)
+      REAL phis(ip1jmp1) ! surface geopotential
+
+c   Local:
+c   ------
+      INTEGER length,l
+      parameter (length = 100)
+      REAL tab_cntrl(length) ! run parameters are stored in this array
+      INTEGER ierr
+
+      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
+      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
+ !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
+      integer :: idim_bandsVI ! "bandsIR" dimension ID #
+      integer, dimension(2) :: id  
+
+c-----------------------------------------------------------------------
+
+      DO l=1,length
+         tab_cntrl(l)=0.
+      ENDDO
+      tab_cntrl(1)  = FLOAT(iim)
+      tab_cntrl(2)  = FLOAT(jjm)
+      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(4)  = FLOAT(idayref)
+      tab_cntrl(5)  = rad
+      tab_cntrl(6)  = omeg
+      tab_cntrl(7)  = g
+      tab_cntrl(8)  = mugaz
+      tab_cntrl(9)  = rcp
+      tab_cntrl(10) = daysec
+      tab_cntrl(11) = dtphys
+      tab_cntrl(12) = etot0
+      tab_cntrl(13) = ptot0
+      tab_cntrl(14) = ztot0
+      tab_cntrl(15) = stot0
+      tab_cntrl(16) = ang0
+c
+c    ..........    P.Le Van  ( ajout le 8/04/96 )    .........
+c         .....        parametres  pour le zoom          ......   
+      tab_cntrl(17)  = clon
+      tab_cntrl(18)  = clat
+      tab_cntrl(19)  = grossismx
+      tab_cntrl(20)  = grossismy
+c
+c     .....   ajout  le 6/05/97 et le 15/10/97  .......
+c
+      IF ( fxyhypb )   THEN
+        tab_cntrl(21) = 1.
+        tab_cntrl(22) = dzoomx
+        tab_cntrl(23) = dzoomy
+      ELSE
+        tab_cntrl(21) = 0.
+        tab_cntrl(22) = dzoomx
+        tab_cntrl(23) = dzoomy
+        tab_cntrl(24) = 0.
+        IF( ysinus )  tab_cntrl(24) = 1.
+      ENDIF
+
+c    .........................................................
+
+! Define dimensions
+    
+      ierr = NF_REDEF (nid)
+
+      ierr = NF_DEF_DIM (nid, "index", length, idim_index)
+!      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+!      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+!      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
+      ierr = NF_DEF_DIM (nid, "bandsVI",L_NSPECTV,idim_bandsVI)
+
+      ierr = NF_ENDDEF(nid)
+
+c  Contol parameters for this run
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, 
+     .       idim_index,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, 
+     .       idim_index,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
+     .                       "Control parameters")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+
+c --------------------------
+c  longitudes and latitudes
+
+
+c$$$      ierr = NF_REDEF (nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
+c$$$#else
+c$$$      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
+c$$$#endif
+c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+c$$$     .                       "Longitudes at u nodes")
+c$$$      ierr = NF_ENDDEF(nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
+c$$$#endif
+
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr =NF_DEF_VAR(nid, "latitude", NF_DOUBLE, 1, idim_rlatu,nvarid)
+#else
+      ierr =NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim_rlatu,nvarid)
+#endif
+      ierr =NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "North latitude")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+#endif
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr =NF_DEF_VAR(nid,"longitude", NF_DOUBLE, 1, idim_rlonv,nvarid)
+#else
+      ierr = NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, idim_rlonv,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "East longitude")
+      ierr = NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+#endif
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, 
+     .       idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, 
+     .       idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",10,"pseudo-alt")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif
+
+
+c$$$c --------------------------
+c$$$      ierr = NF_REDEF (nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
+c$$$#else
+c$$$      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
+c$$$#endif
+c$$$      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+c$$$     .                       "Latitudes at v nodes")
+c$$$      ierr = NF_ENDDEF(nid)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
+c$$$#endif
+
+
+c$$$c --------------------------
+c$$$c  Vertical levels
+c$$$      call def_var(nid,"aps","hybrid pressure at midlayers ","Pa",
+c$$$     .            1,idim_llm,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"bps","hybrid sigma at midlayers"," ",
+c$$$     .            1,idim_llm,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"ap","hybrid pressure at interlayers","Pa",
+c$$$     .            1,idim_llmp1,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
+c$$$#endif
+c$$$
+c$$$      call def_var(nid,"bp","hybrid sigma at interlayers"," ",
+c$$$     .            1,idim_llmp1,nvarid,ierr)
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
+c$$$#else
+c$$$      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
+c$$$#endif
+
+c$$$!-------------------------------
+c$$$! (soil) depth variable mlayer() (known from comsoil.h)
+c$$$!-------------------------------
+c$$$      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+c$$$      ! define variable
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_DOUBLE,1,idim_nsoilmx,nvarid)
+c$$$#else
+c$$$      ierr=NF_DEF_VAR(nid,"soildepth",NF_FLOAT,1,idim_nsoilmx,nvarid)
+c$$$#endif
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20,
+c$$$     .                        "Soil mid-layer depth")
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"m")
+c$$$      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"positive",4,"down")
+c$$$      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+c$$$      ! write variable
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
+c$$$#else
+c$$$      ierr=NF_PUT_VAR_REAL (nid,nvarid,mlayer)
+c$$$#endif
+
+!-------------------------------
+! Number of bands in the IR
+!-------------------------------
+
+      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+      ! define variable
+#ifdef NC_DOUBLE
+      ierr=NF_DEF_VAR(nid,"bandsVI",NF_DOUBLE,1,idim_bandsVI,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"bandsVI",NF_FLOAT,1,idim_bandsVI,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
+     .                        "Band limits in the visible")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
+      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+      ! write variable
+#ifdef NC_DOUBLE
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEV))
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEV))
+#endif
+
+c
+c --------------------------
+c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
+
+      id(1)=idim_rlonu
+      id(2)=idim_rlatu
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+     .             "Conversion coefficients cov <--> natural")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+#endif
+c
+      id(1)=idim_rlonv
+      id(2)=idim_rlatv
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+     .             "Conversion coefficients cov <--> natural")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+#endif
+c
+      id(1)=idim_rlonv
+      id(2)=idim_rlatu
+c
+c --------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "aire", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "aire", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
+     .                       "Mesh area")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+#endif
+
+c  Surface geopotential
+      id(1)=idim_rlonv
+      id(2)=idim_rlatu
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "phisinit", NF_DOUBLE, 2, id,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "phisinit", NF_FLOAT, 2, id,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
+     .                       "Geopotential at the surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
+#endif,
+
+
+!      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+!      write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',rad,omeg,g,mugaz,rcp
+!      write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
+
+      END
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 305)
@@ -27,5 +27,4 @@
       integer nS,nT
       parameter(nS=1649)
-!      parameter(nT=7)
       parameter(nT=14)
 
@@ -34,5 +33,4 @@
       double precision temp_arr(nT)
       double precision abs_arr(nS,nT)
-      !double precision data_tmp(nT+1)
       double precision data_tmp(nT/2+1)
 
@@ -52,5 +50,5 @@
 
          ! cold
-         dt_file=datafile(1:LEN_TRIM(datafile))//'/H2H2_CIA_LT.dat'
+         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2H2_CIA_LT.dat'
          open(33,file=dt_file,form='formatted',status='old',iostat=ios)
          if (ios.ne.0) then        ! file not found
@@ -69,5 +67,5 @@
 
          ! hot
-         dt_file=datafile(1:LEN_TRIM(datafile))//'/H2H2_CIA_HT.dat'
+         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2H2_CIA_HT.dat'
          open(34,file=dt_file,form='formatted',status='old',iostat=ios)
          if (ios.ne.0) then        ! file not found
@@ -106,5 +104,5 @@
          print*,'   pressure ',pres,' Pa'
 
-         call bilinear(wn_arr,temp_arr,nS,nT,abs_arr,wn,temp,abcoef)
+         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
 
          print*,'the absorption is ',abcoef,' cm^-1 amg^-2'
@@ -117,5 +115,5 @@
       else
   
-         call bilinear(wn_arr,temp_arr,nS,nT,abs_arr,wn,temp,abcoef)
+         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
          abcoef=abcoef*100.0*amagat**2 ! convert to m^-1
          !print*,'The absorption is ',abcoef,' m^-1'
@@ -130,9 +128,6 @@
 
 
-
-! put as separate file eventually
-
 !-------------------------------------------------------------------------
-      subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x_in,y_in,f)
+      subroutine bilinearH2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
 !     Necessary for interpolation of continuum data
 
@@ -140,4 +135,6 @@
 
       integer nX,nY,i,j,a,b
+      parameter(nX=1649)
+      parameter(nY=14)
 
       real*8 x_in,y_in,x,y,x1,x2,y1,y2
@@ -176,5 +173,5 @@
       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
          write(*,*) 'Warning from bilinear.for:'
-         write(*,*) 'Extrapolating outside H2-H2 CIA temperature range!'
+         write(*,*) 'Outside continuum temperature range!'
          if(y.lt.y_arr(1))then
             y=y_arr(1)+0.01
@@ -212,3 +209,3 @@
       
       return
-    end subroutine bilinear
+    end subroutine bilinearH2H2
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90	(revision 305)
@@ -0,0 +1,249 @@
+     subroutine interpolateH2Ocont(wn,temp,presS,presF,abcoef,firstcall)
+      implicit none
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Calculates the H2O continuum opacity, using a lookup table from
+!     Clough (2005)
+!
+!     Authors
+!     -------
+!     R. Wordsworth (2011)
+!     
+!==================================================================
+
+#include "datafile.h"
+
+      ! input
+      double precision wn                 ! wavenumber             (cm^-1)
+      double precision temp               ! temperature            (Kelvin)
+      double precision presS              ! self-pressure          (Pascals)
+      double precision presF              ! foreign (air) pressure (Pascals)
+
+      ! output
+      double precision abcoef             ! absorption coefficient (m^-1)
+
+      integer nS,nT
+      parameter(nS=1001)
+      parameter(nT=11)
+
+      double precision kB
+      parameter(kB=1.3806488e-23)
+
+      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
+      double precision wn_arr(nS)
+      double precision temp_arr(nT)
+      double precision abs_arrS(nS,nT)
+      double precision abs_arrF(nS,nT)
+      double precision data_tmp(nT)
+
+      integer k
+      logical firstcall
+
+      save wn_arr, temp_arr, abs_arrS, abs_arrF
+
+      character*100 dt_file
+      integer strlen,ios
+
+      amagatS=(273.15/temp)*(presS/101325.0)
+      amagatF=(273.15/temp)*(presF/101325.0)
+
+      if(firstcall)then
+
+!     1.1 Open the ASCII files
+
+         ! nu array
+         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2O_CONT_NU.dat'
+         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+            write(*,*) 'Error from interpolateH2O_cont.for'
+            write(*,*) 'Data file could not be found:'
+            write(*,*) dt_file
+            call abort
+         else
+            do k=1,nS
+               read(33,*) wn_arr(k)
+            enddo
+         endif
+         close(33)
+
+         ! self broadening
+         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2O_CONT_SELF.dat'
+         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+            write(*,*) 'Error from interpolateH2O_cont.for'
+            write(*,*) 'Data file could not be found:'
+            write(*,*) dt_file
+            call abort
+         else
+            do k=1,nS
+               read(34,*) data_tmp
+               abs_arrS(k,1:nT)=data_tmp(1:nT)
+            end do
+         endif
+         close(34)
+
+         ! foreign (N2+O2+Ar) broadening
+         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2O_CONT_FOREIGN.dat'
+         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+            write(*,*) 'Error from interpolateH2O_cont.for'
+            write(*,*) 'Data file could not be found:'
+            write(*,*) dt_file
+            call abort
+         else
+            do k=1,nS
+               read(35,*) data_tmp
+               abs_arrF(k,1:nT)=data_tmp(1:nT)
+            end do
+         endif
+         close(35)
+
+         temp_arr(1)  = 200.
+         temp_arr(2)  = 250.
+         temp_arr(3)  = 300.
+         temp_arr(4)  = 350.
+         temp_arr(5)  = 400.
+         temp_arr(6)  = 450.
+         temp_arr(7)  = 500.
+         temp_arr(8)  = 550.
+         temp_arr(9)  = 600.
+         temp_arr(10) = 650.
+         temp_arr(11) = 700.
+
+         print*,'At wavenumber ',wn,' cm^-1'
+         print*,'   temperature ',temp,' K'
+         print*,'   H2O pressure ',presS,' Pa'
+         print*,'   air pressure ',presF,' Pa'
+
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
+         print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
+
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
+         print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
+
+         print*,'We have ',amagatS,' amagats of H2O vapour'
+         print*,'and ',amagatF,' amagats of air'
+
+         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
+
+         Nmolec = (presS+presF)/(kB*temp)
+
+         print*,'Total number of molecules per m^3 is',Nmolec
+
+         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
+         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
+
+         print*,'So the total absorption is ',abcoef,' m^-1'
+         print*,'And optical depth / km : ',1000.0*abcoef
+
+      else
+  
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
+
+         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
+
+         !abcoef = (presS/100000.0)*3.0e-24
+         !print*,'Matsui TEST'
+
+         Nmolec = (presS+presF)/(kB*temp)  ! sure this is correct??
+         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
+         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
+
+         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
+         ! however our bands are normally thin, so this is no big deal.
+
+
+      endif
+
+      return
+    end subroutine interpolateH2Ocont
+
+
+
+
+!-------------------------------------------------------------------------
+      subroutine bilinearH2Ocont(x_arr,y_arr,f2d_arr,x_in,y_in,f)
+!     Necessary for interpolation of continuum data
+
+      implicit none
+
+      integer nX,nY,i,j,a,b
+      parameter(nX=1001)
+      parameter(nY=11)
+
+      real*8 x_in,y_in,x,y,x1,x2,y1,y2
+      real*8 f,f11,f12,f21,f22,fA,fB
+      real*8 x_arr(nX)
+      real*8 y_arr(nY)
+      real*8 f2d_arr(nX,nY)
+      
+      integer strlen
+      character*100 label
+      label='subroutine bilinear'
+
+      x=x_in
+      y=y_in
+
+!     1st check we're within the wavenumber range
+      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
+         f=0.0D+0
+         return
+      else
+         
+!     in the x (wavenumber) direction 1st
+         i=1
+ 10      if (i.lt.(nX+1)) then
+            if (x_arr(i).gt.x) then
+               x1=x_arr(i-1)
+               x2=x_arr(i)
+               a=i-1
+               i=9999
+            endif
+            i=i+1
+            goto 10
+         endif
+      endif
+      
+      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
+         write(*,*) 'Warning from bilinear.for:'
+         write(*,*) 'Outside continuum temperature range!'
+         if(y.lt.y_arr(1))then
+            y=y_arr(1)+0.01
+         endif
+         if(y.gt.y_arr(nY))then
+            y=y_arr(nY)-0.01
+         endif
+      else
+
+!     in the y (temperature) direction 2nd
+         j=1
+ 20      if (j.lt.(nY+1)) then
+            if (y_arr(j).gt.y) then
+               y1=y_arr(j-1)
+               y2=y_arr(j)
+               b=j-1
+               j=9999
+            endif
+            j=j+1
+            goto 20
+         endif
+      endif
+      
+      f11=f2d_arr(a,b)
+      f21=f2d_arr(a+1,b)
+      f12=f2d_arr(a,b+1)
+      f22=f2d_arr(a+1,b+1)
+      
+!     1st in x-direction
+      fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
+      fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
+      
+!     then in y-direction
+      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
+
+      return
+    end subroutine bilinearH2Ocont
Index: unk/LMDZ.GENERIC/libf/phystd/kastprof_fn.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/kastprof_fn.F90	(revision 304)
+++ 	(revision )
@@ -1,142 +1,0 @@
-      subroutine kastprof_fn(tsurf,temp,popsk,play,plev,qvap)
-      implicit none
-
-!     creates a kasting temperature profile
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "callkeys.h"
-
-      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
-      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
-      REAL psurf,tsurf      
-      REAL tcond(nlayermx)
-      REAL qvap(nlayermx)
-      REAL temp(nlayermx)   ! temperature at the middle of the layers
-
-      real, parameter :: Rstar   = 8.314
-      real, parameter :: ptriple = 518000.0
-      real, parameter :: latcond = 5.9e5     ! as in condense_co2
-      real, parameter :: mn      = 0.02801   ! molar mass of background gas (N2)
-      real, parameter :: mv      = 0.04401   ! molar mass of background gas (CO2)
-
-      real popsk(nlayermx)  ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx)
-                            ! but kastprof requires 1D, so its ok I guess
-
-      real ppco2
-      real vmr_local,q_local,gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p
-
-      integer l,iter
-
-      logical co2mixedcond
-
-      co2mixedcond=.true.
-
-      if(ngridmx.gt.1)then
-         print*,'DO NOT RUN KASTPROF IN 3D!!!'
-         call abort
-      end if
-
-      ! 1st, we make pot. temperature constant everywhere
-      do l=1,nlayermx
-            temp(l)=TsurfK*popsk(l)
-      enddo
-
-!     then do all the other stuff
-      vmr_local = 1.-n2mixratio
-      q_local   = vmr_local*mv/mn
-
-      qvap(1)   = q_local
-
-      ppco2=vmr_local*play(1)
-
-      if(play(1).lt.ptriple)then !
-         tcond(1) = (-3167.8)/(log(.01*ppco2)-23.23)         ! Fanale's formula
-      else
-         tcond(1) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data)
-      endif
-
-      do l=1,nlayermx-1
-
-         ppco2=vmr_local*play(l) ! this shld be implicit i.e. vmr_local(l+1) ??
-
-         if(play(l).le.ptriple)then ! ERROR
-            tcond(l) = (-3167.8)/(log(.01*ppco2)-23.23)         ! Fanale's formula
-         else
-            tcond(l) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data)
-         endif
-
-         if(co2mixedcond)then
-            if (temp(l) .le. Tstrat)then
-               temp(l)  = Tstrat            ! isothermal stratosphere
-               temp(l+1) = Tstrat
-               qvap(l+1) = qvap(l)
-               !print*,'strato'
-
-            elseif(temp(l) .le. tcond(l))then
-
-               dln_p     = log(play(l+1))-log(play(l))
-               gamma     = cpp - q_local*latcond/temp(l)
-               dqdlnT    = (mv*latcond/(mn*temp(l)) - gamma) & 
-                            / (latcond/temp(l) + Rstar/(mn*q_local))
-               dlnqdlnT  = (1./q_local)*dqdlnT
-               dlnpdlnT  = mv*latcond/(Rstar*temp(l)) - dlnqdlnT/(1. + q_local*mn/mv)
-               dpdT      = (play(l)/temp(l))*dlnpdlnT
-               q_local   = q_local + dqdlnT*dln_p/dlnpdlnT
-
-               vmr_local = q_local*mn/mv
-
-
-               if(1.eq.1)then
-               print*,'l         =',l
-               print*,'gamma     =',gamma
-               print*,'gamma bit =',q_local*latcond/temp(l)
-               print*,'cpp       =',cpp
-               print*,'latcond   =',latcond
-               print*,'q_local   =',q_local
-               print*,'T         =',temp(l)
-               print*,'p         =',play(l)
-               print*,'dqdlnT    =',dqdlnT
-               print*,'dlnqdlnT  =',dlnqdlnT
-               print*,'dlnpdlnT  =',dlnpdlnT
-               print*,'dpdT      =',dpdT
-               call abort
-               endif
-
-               qvap(l+1) = q_local
-               temp(l+1) = temp(l) + temp(l)*dln_p/(dlnpdlnT)
-               !print*,'CO2co'
-            else
-               !print*,'tropo'
-
-
-               ! CO2 is not condensing; therefore dry convection
-               qvap(l+1) = qvap(l)
-            endif
-
-         else
-
-            if (temp(l) .le. Tstrat)then
-               temp(l) = Tstrat            
-               ! isothermal stratosphere
-            elseif(temp(l) .le. tcond(l))then
-               temp(l)=tcond(l)
-               ! pure CO2 condensation
-            else
-               ! CO2 is not condensing; therefore dry convection 
-            endif
-
-         endif
-
-      enddo ! vertical loop
-
-      temp(nlayermx) = Tstrat
-      qvap(nlayermx) = qvap(nlayermx-1)
-
-      print*,'---fin---'
-
-      call abort
-
-      return
-    end subroutine kastprof_fn
Index: /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 305)
@@ -0,0 +1,321 @@
+program kcm1d
+
+  use radinc_h,      only: NAERKIND
+  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
+  use watercommon_h, only: mH2O
+  use ioipsl_getincom 
+
+  implicit none
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Run the universal model radiative transfer once in a 1D column.
+!     Useful for climate evolution studies etc.
+!     
+!     It can be compiled with a command like (e.g. for 25 layers):
+!                                  "makegcm -p std -d 25 kcm1d"
+!     It requires the files "callphys.def", "gases.def"
+!     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
+!
+!     Authors
+!     -------
+!     R. Wordsworth
+!
+!==================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "advtrac.h"
+
+! --------------------------------------------------------------
+!  Declarations
+! --------------------------------------------------------------
+
+  integer nlayer,nlevel,nq
+  integer ilay,ilev,iq,iw,iter
+  real play(nlayermx)     ! Pressure at the middle of the layers [Pa]
+  real zlay(nlayermx)     ! Altitude at middle of the layers [km]
+  real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
+  real temp(nlayermx)     ! temperature at the middle of the layers [K]
+  real q(nlayermx,nqmx)   ! tracer mixing ratio [kg/kg]
+  real vmr(nlayermx,nqmx) ! tracer mixing ratio [mol/mol]
+  real qsurf(nqmx)        ! tracer surface budget [kg/kg] ????
+  real psurf,psurf_n,tsurf      
+
+  real emis, albedo
+
+  real cpp3D(nlayermx)   ! specific heat capacity at const. pressure
+  real rcp3D(nlayermx)   ! R / specific heat capacity
+  real muvar(nlayermx+1)
+
+  real dtsw(nlayermx) ! heating rate (K/s) due to SW
+  real dtlw(nlayermx) ! heating rate (K/s) due to LW
+  real fluxsurf_lw   ! incident LW flux to surf (W/m2)
+  real fluxtop_lw    ! outgoing LW flux to space (W/m2)
+  real fluxsurf_sw   ! incident SW flux to surf (W/m2)
+  real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
+  real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
+      
+! not used
+  real reffrad(nlayermx,naerkind)
+  real nueffrad(nlayermx,naerkind)
+  real cloudfrac(nlayermx)
+  real totcloudfrac
+  real tau_col
+
+  real dTstrat
+  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
+  real OLR_nu(L_NSPECTI)
+  real GSR_nu(L_NSPECTV)
+  real Eatmtot
+
+  integer ierr
+  logical firstcall,lastcall
+
+! --------------
+! Initialisation
+! --------------
+
+  pi=2.E+0*asin(1.E+0)
+
+  reffrad(:,:)  = 0.0
+  nueffrad(:,:) = 0.0
+  cloudfrac(:)  = 0.0
+  totcloudfrac  = 0.0
+
+!  Load parameters from "run.def" 
+! -------------------------------
+
+! check if 'run1d.def' file is around (otherwise reading parameters
+! from callphys.def via getin() routine won't work.)
+  open(90,file='run.def',status='old',form='formatted',&
+       iostat=ierr)
+  if (ierr.ne.0) then
+     write(*,*) 'Cannot find required file "run.def"'
+     write(*,*) '  (which should contain some input parameters'
+     write(*,*) '   along with the following line:'
+     write(*,*) '   INCLUDEDEF=callphys.def'
+     write(*,*) '   )'
+     write(*,*) ' ... might as well stop here ...'
+     stop
+  else
+     close(90)
+  endif
+
+  nlayer=nlayermx
+  nlevel=nlayer+1
+
+  psurf_n=100000. ! default value for psurf
+  print*,'Dry surface pressure (Pa)?'
+  call getin("psurf",psurf_n)
+  write(*,*) " psurf = ",psurf_n
+      
+  tsurf=300.0 ! default value for tsurf
+  print*,'Surface temperature (K)?'
+  call getin("tref",tsurf)
+  write(*,*) " tsurf = ",tsurf
+
+  g=10.0 ! default value for g
+  print*,'Gravity ?'
+  call getin("g",g)
+  write(*,*) " g = ",g
+  
+  periastr = 1.0
+  apoastr  = 1.0
+  print*,'Periastron (AU)?'
+  call getin("periastr",periastr)
+  write(*,*) "periastron = ",periastr
+  dist_star = periastr 
+  fract     = 0.5
+  print*,'Apoastron (AU)?'
+  call getin("apoastr",apoastr)
+  write(*,*) "apoastron = ",apoastr
+ 
+  periastr = periastr*149.598 ! AU to Mkm
+  apoastr  = apoastr*149.598  ! AU to Mk
+
+  albedo=0.2 ! default value for albedo
+  print*,'Albedo of bare ground?'
+  call getin("albedo",albedo)
+  write(*,*) " albedo = ",albedo
+
+  emis=1.0 ! default value for emissivity
+  PRINT *,'Emissivity of bare ground ?'
+  call getin("emis",emis)
+  write(*,*) " emis = ",emis
+
+  pceil=100.0 ! Pascals
+  PRINT *,'Ceiling pressure (Pa) ?'
+  call getin("pceil",pceil)
+  write(*,*) " pceil = ", pceil
+
+  mugaz=0.
+  cpp=0.
+
+  call su_gases
+  call calc_cpp_mugaz
+
+  call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
+
+! Tracer initialisation
+! ---------------------
+  if (tracer) then
+     ! load tracer names from file 'traceur.def'
+     open(90,file='traceur.def',status='old',form='formatted',&
+          iostat=ierr)
+     if (ierr.eq.0) then
+        write(*,*) "kcm1d: Reading file traceur.def"
+        ! read number of tracers:
+        read(90,*,iostat=ierr) nq
+        if (ierr.ne.0) then
+           write(*,*) "kcm1d: error reading number of tracers"
+           write(*,*) "   (first line of traceur.def) "
+           stop
+        else
+           ! check that the number of tracers is indeed nqmx
+           if (nq.ne.nqmx) then
+              write(*,*) "kcm1d: error, wrong number of tracers:"
+              write(*,*) "nq=",nq," whereas nqmx=",nqmx
+              stop
+           endif
+        endif
+        
+        do iq=1,nq
+           ! minimal version, just read in the tracer names, 1 per line
+           read(90,*,iostat=ierr) tnom(iq)
+           if (ierr.ne.0) then
+              write(*,*) 'rcm1d: error reading tracer names...'
+              stop
+           endif
+        enddo !of do iq=1,nq
+     endif
+  endif
+  
+  call initracer()
+
+  do iq=1,nqmx
+     do ilay=1,nlayer
+        q(ilay,iq) = 0.
+     enddo
+  enddo
+  
+  do iq=1,nqmx
+     qsurf(iq) = 0.
+  enddo
+
+  firstcall = .true.
+  lastcall  = .false.
+
+  iter    = 1
+  Tstrat  = 200.0
+  dTstrat = 1000.0
+
+! ---------
+! Run model
+! ---------
+  do
+     psurf = psurf_n
+
+!    Create vertical profiles
+     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
+          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
+
+
+     !if(1.eq.2)then
+
+!    Run radiative transfer
+     call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
+          albedo,emis,mu0,plev,play,temp,                    &
+          tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
+          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
+          fluxabs_sw,fluxtop_dn,reffrad,tau_col,  &
+          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
+ 
+!    Iterate stratospheric temperature
+     print*,'Tstrat = ',Tstrat
+     dTstrat = Tstrat
+     Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
+     dTstrat = dTstrat-Tstrat
+
+     !endif
+
+     !dTstrat = 0
+
+     if(abs(dTstrat).lt.1.0)then
+        print*,'dTstrat = ',dTstrat
+        exit
+     endif
+
+     iter=iter+1
+     if(iter.eq.100)then
+        print*,'Stratosphere failed to converge after'
+        print*,'100 iteration, aborting run.'
+        call abort
+     endif
+
+  end do
+
+! Calculate total atmospheric energy
+  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
+     q(:,1),muvar,Eatmtot)
+
+! ------------------------
+! Save data to ascii files
+! ------------------------
+
+  print*,'Saving profiles...'
+  open(115,file='profpres.out',form='formatted')
+  open(116,file='proftemp.out',form='formatted')
+  open(117,file='profztab.out',form='formatted')
+  open(118,file='profqvar.out',form='formatted')
+  open(119,file='profvvar.out',form='formatted')
+
+  write(115,*) psurf
+  write(116,*) tsurf
+  write(117,*) 0.0
+  write(118,*) qsurf(1)
+  write(119,*) qsurf(1)*(muvar(1)/mH2O)
+  do ilay=1,nlayer
+     vmr(ilay,1) = q(ilay,1)*(muvar(ilay+1)/mH2O)
+     write(115,*) play(ilay)
+     write(116,*) temp(ilay)
+     write(117,*) zlay(ilay)
+     write(118,*) q(ilay,1)
+     write(119,*) vmr(ilay,1)
+  enddo
+  close(115)
+  close(116)
+  close(117)
+  close(118)
+  close(119)
+
+  print*,'Saving scalars...'
+  open(116,file='surf_vals.out')
+  write(116,*) tsurf,psurf,fluxtop_dn,         &
+       fluxabs_sw,fluxtop_lw 
+  close(116)
+  open(111,file='ene_vals.out')
+  write(111,*) tsurf,psurf,Eatmtot,Tstrat
+  close(111)
+
+  print*,'Saving spectra...'
+  open(117,file='OLRnu.out')
+  do iw=1,L_NSPECTI
+     write(117,*) OLR_nu(iw)
+  enddo
+  close(117)
+  
+  open(127,file='GSRnu.out')
+  do iw=1,L_NSPECTV
+     write(127,*) GSR_nu(iw)
+  enddo
+  close(127)  
+
+end program kcm1d
Index: /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 305)
@@ -0,0 +1,363 @@
+subroutine kcmprof_fn(psurf_rcm,qsurf_rcm,Tsurf_rcm,Tstra_rcm,P_rcm,Pl_rcm,z_rcm,T_rcm,q_rcm,m_rcm)
+
+use params_h
+use watercommon_h, only : mH2O
+implicit none
+
+!     ----------------------------------------------------------------
+!     Purpose: create profiles of T, rho_v, rho_n, Pv and Pn following 
+!              Kasting 1988
+!     Authour: Adapted from a code by E. Marcq by R. Wordsworth (2011)
+!     ----------------------------------------------------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "callkeys.h"
+#include "gases.h"
+
+  integer ilay, nlay
+  parameter (nlay=10000) ! number of vertical layers
+
+  ! rcm inputs
+  real Tsurf_rcm,Tstra_rcm
+
+  ! rcm outputs
+  real psurf_rcm,qsurf_rcm
+  real P_rcm(1:nlayermx)
+  real Pl_rcm(1:nlayermx+1)
+  real z_rcm(1:nlayermx)
+  real T_rcm(1:nlayermx),q_rcm(1:nlayermx)
+  real m_rcm(1:nlayermx+1)
+
+  ! rcm for interpolation (should really use log coords?)
+  !double precision p1,p2,pnew,ilay_rcm
+
+  double precision lnp1,lnp2,lnpnew
+  real Dp_rcm, dlogp_rcm
+  integer ilay_rcm,ilev_rcm,ifinal_rcm
+
+  double precision Dz, Dp
+  double precision Ptop, dlogp, Psat_max
+  parameter (Ptop=1.0)           ! Pressure at TOA [Pa]
+  parameter (Psat_max=5000.0) ! Maximum vapour pressure [Pa]
+
+  double precision T(1:nlay)                       ! temperature [K]
+  double precision Ztab(1:nlay)                    ! altitude [m]
+  double precision Pv(1:nlay),Pn(1:nlay),P(1:nlay) ! pressure [Pa]
+  double precision rho_v(1:nlay), rho_n(1:nlay)    ! density [kg m^-3]
+  double precision a_v(1:nlay)                     ! = rho_v/rho_n [kg/kg]
+  double precision q_v(1:nlay)                     ! = rho_v/rho_tot [kg/kg]
+  double precision mtot(1:nlay)                    ! = (rho_v+rho_n)/(n_v+n_n) [g/mol]
+
+  integer profil_flag(1:nlay) ! 0 = dry, 1 = moist, 2 = isothermal
+
+  ! inputs
+  double precision Tsurf   ! surface temperature [K]
+  double precision Psurf_v ! surface par. pressure (variable species) [Pa]
+  double precision Psurf_n ! surface par. pressure (incondensible species)[Pa]
+  double precision Ttop    ! stratospheric temperature [K]
+
+  double precision dTdp        ! [K/Pa]
+  double precision dPvdp,dPndp ! [Pa/Pa]
+  double precision psat_v      ! local Psat_H2O value
+  double precision Tcrit       ! Critical temperature [K]
+!  double precision Tsat        ! local Tsat_H2O value
+  double precision rho_vTEMP,rho_nTEMP
+
+  ! variables necessary for steam.f90
+  double precision rhol,rhov,nul
+
+  ! for output
+  double precision vmr
+
+  logical verbose
+  parameter(verbose=.false.)
+
+  ! initialise flags
+  profil_flag(:) = 0
+
+  !-------------------------------
+  ! assign input variables
+  m_n     = dble(mugaz/1000.)
+  ! modify/generalise later??
+
+  if(gnom(2).eq.'H2O')then
+     m_v   = dble(mH2O/1000.)
+     tcrit = 6.47d2
+  elseif(gnom(2).eq.'NH3')then
+     m_v   = 17.031/1000.
+     tcrit = 4.06d2
+  elseif(gnom(2).eq.'CH4')then
+     m_v   = 16.04/1000.
+     tcrit = 1.91d2
+     stop
+  else
+     print*,'Variable gas not recognised!'
+     call abort
+  endif
+
+  rmn     = rc/m_n
+  
+  Psurf_n = dble(psurf_rcm)
+  Ttop    = dble(Tstra_rcm)
+  Tsurf   = dble(Tsurf_rcm)
+
+
+  ! assume vapour saturation at surface (for now)
+!  if (tsurf > tcrit) then
+!     print*,'Above critical temperature for ',gnom(2),&
+!          ', at surface, assuming 10 bar pressure instead'
+!     Psurf_v = 10*1d5
+!     profil_flag(1) = 0
+!i!  endif
+
+  if(gnom(2).eq.'H2O')then
+     call Psat_H2O(tsurf,psat_v)
+  elseif(gnom(2).eq.'NH3')then
+     call Psat_NH3(tsurf,psat_v)
+  endif
+
+  ! Moist adiabat unless greater than psat_max
+  if(psat_v*1d6.lt.psat_max)then
+    Psurf_v = Psat_v*1d6
+    profil_flag(1) = 1
+  else 
+    Psurf_v = psat_max
+    profil_flag(1) = 0
+  endif
+
+
+
+  if(verbose)then
+     print*,'Psat_v  =',psat_v*1d6
+     print*,'Tsurf   =',Tsurf,' K'
+     print*,'Ttop    =',Ttop,' K'
+     print*,'Psurf_v =',Psurf_v,' Pa'
+     print*,'Psurf_n =',Psurf_n,' Pa'
+     print*,'m_n     =',m_n,' kg/mol'
+     print*,'m_v     =',m_v,' kg/mol'
+     print*,'rc      =',rc
+  endif
+
+  ! define fine pressure grid
+  psurf_rcm = real(Psurf_n+Psurf_v)
+  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayermx
+
+  P_rcm(1)  = psurf_rcm*exp(dlogp_rcm)
+  do ilay_rcm=1,nlayermx-1
+     P_rcm(ilay_rcm+1) = P_rcm(ilay_rcm)*exp(dlogp_rcm)
+  enddo
+
+  Pl_rcm(1) = psurf_rcm
+  do ilev_rcm=2,nlayermx
+     ! log-linear interpolation
+     Pl_rcm(ilev_rcm) = exp( log( P_rcm(ilev_rcm)*P_rcm(ilev_rcm-1) )/2 )
+  enddo
+
+  !-------------------------------
+  ! Layer 1
+  T(1)     = Tsurf
+  Pv(1)    = Psurf_v
+  Pn(1)    = Psurf_n
+  rho_n(1) = m_n*Pn(1)/(Rc*Tsurf)
+  rho_v(1) = m_v*Pv(1)/(Rc*Tsurf)
+  a_v(1)   = rho_v(1)/rho_n(1)
+
+  ! log pressure grid spacing (constant)
+  dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1)
+
+  call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp)
+  if(verbose)then
+     print*, 'dT/dp ground [K/Pa] =',dTdp
+  endif
+
+  ! initial delta p, delta z
+  Dp = (Pn(1) + Pv(1))*(exp(dlogp) - 1d0)
+  Dz = -Dp/(  g*(rho_n(1) + rho_v(1))  )
+
+  !-------------------------------
+  ! Layer 2
+  T(2)     = tsurf + dTdp*Dp 
+  Pv(2)    = Pv(1) + dPvdp*Dp
+  Pn(2)    = Pn(1) + dPndp*Dp
+  rho_n(2) = m_n*Pn(2)/(Rc*T(2))
+  rho_v(2) = m_v*Pv(2)/(Rc*T(2))
+  a_v(2)   = rho_v(2)/rho_n(2)
+
+  !-------------------------------
+  ! start vertical ascent
+  Ztab(1) = 0.
+  do ilay=2,nlay-1
+
+
+
+
+     ! calculate altitude levels (for diagnostic only)
+     Dz         = -Dp/(  g*(rho_n(ilay) + rho_v(ilay))  )
+     Ztab(ilay) = Dz + Ztab(ilay-1)
+
+     ! 1st assume next layer same as last one
+     profil_flag(ilay) = profil_flag(ilay-1)
+  
+     ! update delta p
+     Dp = (Pn(ilay)+Pv(ilay))*(exp(dlogp) - 1d0)
+
+     ! intial gradients call to calculate temperature at next level
+     call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),&
+                        T(ilay),dTdp,dPvdp,dPndp)
+
+     T(ilay+1) = T(ilay) + dTdp*Dp
+
+     ! test for moist adiabat at next level
+     if(gnom(2).eq.'H2O')then
+        call Psat_H2O(T(ilay+1),psat_v)
+     elseif(gnom(2).eq.'NH3')then
+        call Psat_NH3(T(ilay+1),psat_v)
+     endif
+
+     if (psat_v*1d6 .lt. Pv(ilay)+dPvdp*Dp) then
+        profil_flag(ilay)=1
+        call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),&
+                         T(ilay),dTdp,dPvdp,dPndp)
+     endif
+
+     ! test for stratosphere at next level
+     if (T(ilay+1) .le. Ttop) then
+        profil_flag(ilay)=2
+        T(ilay+1)=Ttop
+     endif
+
+     ! calculate pressures at next level
+     Pn(ilay+1) = Pn(ilay) + dPndp*Dp
+     Pv(ilay+1) = Pv(ilay) + dPvdp*Dp
+
+     if(profil_flag(ilay) .eq. 1)then
+        
+        if(gnom(2).eq.'H2O')then
+           call Psat_H2O(T(ilay+1),psat_v)
+        elseif(gnom(2).eq.'NH3')then
+           call Psat_NH3(T(ilay+1),psat_v)
+        endif
+
+        if(Pv(ilay+1) .lt. psat_v*1e6)then
+           Pv(ilay+1)=psat_v*1d6
+        endif
+
+     endif
+
+     ! calculate gas densities at next level (assume ideal)
+     rho_n(ilay+1) = m_n*Pn(ilay+1)/(rc*T(ilay+1)) 
+     select case(profil_flag(ilay)) 
+     case(2) ! isothermal
+        rho_v(ilay+1) = rho_v(ilay)/rho_n(ilay)*rho_n(ilay+1)
+     case(1) ! moist
+
+        ! dont think this is necessary
+        !call psat_est(T(ilay+1),psat_v)
+        ! modify for ammonia!!!
+
+        rho_v(ilay+1) = m_v*psat_v*1d6/(rc*T(ilay+1))
+     case(0) ! dry
+        rho_v(ilay+1) = m_v*Pv(ilay+1)/(rc*T(ilay+1)) 
+     end select
+
+
+
+!         print*,'profil_flag=',profil_flag(ilay)
+!         print*,'T=',T(ilay)
+!         print*,'a=',rho_v(ilay)/rho_n(ilay)
+!      if(profil_flag(ilay).eq.1)then
+!           stop
+!     endif
+
+
+
+
+
+  enddo
+
+  Ztab(nlay)=Ztab(nlay-1)+Dz
+
+  !-------------------------------
+  ! save to kcm1d variables 
+
+  ! surface quantities
+  psurf_rcm = Pn(1) + Pv(1)
+  qsurf_rcm = rho_v(1)/(rho_v(1) + rho_n(1))
+
+  ! create q_v, mtot for saving
+  do ilay=1,nlay
+     mtot(ilay) = 1d3*(rho_v(ilay) + rho_n(ilay)) / &
+                  (rho_v(ilay)/m_v + rho_n(ilay)/m_n)
+     q_v(ilay)  = rho_v(ilay)/(rho_v(ilay) + rho_n(ilay))
+     ! CHECK THIS
+  enddo
+
+
+  ! convert to rcm lower-res grid
+  z_rcm(:) = 0.0
+  T_rcm(:) = 0.0
+  q_rcm(:) = 0.0
+  m_rcm(:) = 0.0
+
+  m_rcm(1) = real( 1d3*(rho_v(1) + rho_n(1)) / &
+                  (rho_v(1)/m_v + rho_n(1)/m_n) )
+
+  ilay_rcm=1
+  do ilay=2,nlay
+
+     if(ilay_rcm.le.nlayermx)then
+     ! interpolate rcm variables
+
+        if(Pn(ilay)+Pv(ilay) .lt. P_rcm(ilay_rcm))then
+
+           if(ilay.eq.1)then
+              print*,'Error in create_profils: Psurf here less than Psurf in RCM!'
+              call abort
+           endif
+
+           lnp1   = log(Pn(ilay-1)+Pv(ilay-1))
+           lnp2   = log(Pn(ilay)+Pv(ilay))
+           lnpnew = dble(log(P_rcm(ilay_rcm)))
+
+           z_rcm(ilay_rcm) = real(Ztab(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
+                             + Ztab(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
+           T_rcm(ilay_rcm) = real(T(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
+                             + T(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
+           q_rcm(ilay_rcm) = real(q_v(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
+                             + q_v(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
+
+           m_rcm(ilay_rcm+1) = real(mtot(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
+                             + mtot(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
+
+           ilay_rcm = ilay_rcm+1
+        endif
+
+     endif
+  enddo
+
+  ifinal_rcm=ilay_rcm-1
+  if(ifinal_rcm.lt.nlayermx)then
+     if(verbose)then
+        print*,'Interpolation in kcmprof stopped at layer',ilay_rcm,'!'
+     endif
+
+     do ilay_rcm=ifinal_rcm+1,nlayermx
+
+        z_rcm(ilay_rcm) = z_rcm(ilay_rcm-1)
+        T_rcm(ilay_rcm) = T_rcm(ilay_rcm-1)
+        q_rcm(ilay_rcm) = q_rcm(ilay_rcm-1)
+        m_rcm(ilay_rcm+1) = m_rcm(ilay_rcm)
+
+     enddo
+  endif
+
+  do ilay=2,nlayermx
+     if(T_rcm(ilay).lt.Ttop)then
+        T_rcm(ilay)=Ttop
+     endif
+  enddo
+
+  return
+end subroutine kcmprof_fn
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 305)
@@ -1,5 +1,5 @@
       subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
            QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
-           TMID,PMID,TAUGSURF,QVAR)
+           TMID,PMID,TAUGSURF,QVAR,MUVAR)
 
       use radinc_h
@@ -59,8 +59,8 @@
       real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
       real*8 DCONT
-      double precision wn_cont, p_cont, T_cont
-
-!     Water mixing ratio variables
-      real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS)
+      double precision wn_cont, p_cont, p_air, T_cont, dtemp
+
+!     Variable species mixing ratio variables
+      real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
       real*8  KCOEF(4)
       integer NVAR(L_LEVELS)
@@ -97,5 +97,4 @@
       do K=2,L_LEVELS
          DPR(k) = PLEV(K)-PLEV(K-1)
-         U(k)   = Cmk*DPR(k)    ! only Cmk line in optci.F
 
          !--- Kasting's CIA ----------------------------------------
@@ -110,5 +109,11 @@
 
          ! if we have continuum opacities, we need dz
-         dz(k)=dpr(k)*R*TMID(K)/(g*PMID(K))
+         if(kastprof)then
+            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
+            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+         else
+            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
+            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+         endif
 
          call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
@@ -132,19 +137,30 @@
             DCONT = 0.0 ! continuum absorption
 
-            ! include H2 continuum
+            ! include continua if necessary
+            wn_cont = dble(wnoi(nw))
+            T_cont  = dble(TMID(k))
             do igas=1,ngasmx
-               if(gnom(igas).eq.'H2_')then 
-
-                  wn_cont = dble(wnoi(nw))
-                  T_cont  = dble(TMID(k))
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas))
-
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.)
-
-                  DCONT=DCONT*dz(k)
-
-
+
+               if(gfrac(igas).eq.-1)then ! variable
+                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+               else
+                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
                endif
+
+               dtemp=0.0
+               if(gnom(igas).eq.'H2_')then
+                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+
+               endif
+
+               DCONT = DCONT + dtemp
+
             enddo
+
+
+            DCONT = DCONT*dz(k) 
 
             !--- Kasting's CIA ----------------------------------------
@@ -154,8 +170,4 @@
             !----------------------------------------------------------
 
-            ! Water continuum currently inactive!
-            !if(varactive)then
-            !   call water_cont(PMID(K),TMID(K),WRATIO(K),WNOI(nw),DCO2)
-            !endif
 
             do ng=1,L_NGAUSS-1
@@ -197,4 +209,9 @@
                
                TAUGAS  = U(k)*ANS
+
+               if(graybody)then
+                  TAUGAS = 0.0
+                  DCONT  = U(k)*3.3e-26
+               endif
 
                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
@@ -311,27 +328,9 @@
       END DO
 
-
       return
+
+
     end subroutine optci
 
 
-!-------------------------------------------------------------------------
-    subroutine water_cont(p,T,wratio,nu,DCONT)
-!   Calculates the continuum opacity for H2O
-!   NOT CURRENTLY USED
-
-      implicit none
-
-      real p, T, wratio, nu, DCONT
-      real h1, h2
-
-      h1 = exp(1800.*(1./T - 0.0034))
-      h2 = 1.25e-22 + 1.67e-19*exp(-2.62e-13*nu)
-
-!      DCONT = h1*h2*p*(p*wratio)**2/(R*T)
-!      DCONT=0.0 ! to be implemented...
-
-      return
-
-    end subroutine water_cont
-
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 305)
@@ -1,5 +1,5 @@
       SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
           QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
-          TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR)
+          TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
 
       use radinc_h
@@ -30,5 +30,5 @@
 !     TLEV(L) - Temperature at the layer boundary
 !     PLEV(L) - Pressure at the layer boundary (i.e. level)
-!     GASV(NT,NPS,NW,NG) - Visual CO2 k-coefficients 
+!     GASV(NT,NPS,NW,NG) - Visible k-coefficients 
 !     
 !-------------------------------------------------------------------
@@ -65,6 +65,6 @@
       real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
 
-!     Water mixing ratio variables
-      real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
+!     Variable species mixing ratio variables
+      real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
       real*8 KCOEF(4)
       integer NVAR(L_LEVELS)
@@ -74,5 +74,5 @@
 
 !     variables for k in units m^-1
-      double precision wn_cont, p_cont, T_cont
+      double precision wn_cont, p_cont, p_air, T_cont, dtemp
       real*8 dz(L_LEVELS), DCONT
 
@@ -92,8 +92,12 @@
          DPR(k) = PLEV(K)-PLEV(K-1)
 
-         U(k)   = Cmk*DPR(k)
-
          ! if we have continuum opacities, we need dz
-         dz(k)  = dpr(k)*R*TMID(K)/(g*PMID(K))
+         if(kastprof)then
+            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
+            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+         else
+            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
+            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+         endif
 
          call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
@@ -128,22 +132,33 @@
             DCONT = 0.0 ! continuum absorption
 
-            ! include H2 continuum
+            ! include continua if necessary
+            wn_cont = dble(wnov(nw))
+            T_cont  = dble(TMID(k))
             do igas=1,ngasmx
+
+               if(gfrac(igas).eq.-1)then ! variable
+                  p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
+               else
+                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+               endif
+
+               dtemp=0.0
                if(gnom(igas).eq.'H2_')then 
-
-                  wn_cont = dble(wnov(nw))
-                  T_cont  = dble(TMID(k))
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas))
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.)
-
-                  DCONT=DCONT*dz(k) 
-
-                  if((.not.callgasvis))then
-                     DCONT=0.0
-                  endif
-
-
-               endif
+                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+               endif
+
+               DCONT = DCONT + dtemp
+
             enddo
+
+            DCONT = DCONT*dz(k) 
+
+               if((.not.callgasvis))then
+                  DCONT=0.0
+               endif
+
 
             do NG=1,L_NGAUSS-1
@@ -184,4 +199,11 @@
 
                TAUGAS          = U(k)*ANS
+
+               if(graybody)then
+                  TAUGAS = 0.0
+                  DCONT  = 0.0
+               endif
+
+
 
                !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
Index: /trunk/LMDZ.GENERIC/libf/phystd/params_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/params_h.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/params_h.F90	(revision 305)
@@ -0,0 +1,16 @@
+      module params_h
+
+!  for use with kcm subroutines
+!  RDW 22/09/11
+
+      implicit none
+
+      double precision rc,m_n,m_v,rmn
+      save m_n,m_v,rmn
+
+      logical ideal_v
+      save ideal_v
+
+      parameter (rc  = 8.314462d0) ! ideal gas constant [J mol^-1 K^-1]
+
+      end module params_h
Index: /trunk/LMDZ.GENERIC/libf/phystd/phasechange_kcm.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/phasechange_kcm.F90	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/phasechange_kcm.F90	(revision 305)
@@ -0,0 +1,1286 @@
+!  ----------------------------------------------------------------
+!  Purpose: Thermodynamic data on H2O, NH3
+!  Authour: Adapted from various sources by R. Wordsworth (2011)
+
+
+! ---------------------
+! NH3
+! ---------------------
+subroutine psat_NH3 ( T, p )
+!
+!  PSAT_NH3 makes a rough estimate of the vapor pressure.
+!
+!  Interpolated from www.engineeringtoolbox.com data by RDW 21/09/11 for
+!  temperatures between 223 and 323 K.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision P, the vapor pressure, in MegaPascals.
+!
+
+  implicit none
+
+  double precision p,T
+
+  p=exp(-1.5609d-004*T**2 + 0.1236*T - 9.1530)/1d6
+
+end subroutine psat_NH3
+
+subroutine latheat_NH3 ( T, sc, sv )
+!
+!  PSAT_NH3 makes a rough estimate of the entropies of condensation /
+!  vapourisation.
+! 
+!  Interpolated from www.engineeringtoolbox.com data by RDW 21/09/11 for
+!  temperatures between 223 and 323 K.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision sc, the entropy of condensate, in J kg^-1 K^-1.
+!            double precision sv, the entropy of gas, in J kg^-1 K^-1.
+!
+
+  implicit none
+
+  double precision T,sc,sv
+
+  sv = 0.0492*T**2 - 40.4199*T + 1.2708e+004
+  sc = -0.0215*T**2 + 28.7138*T - 5.5267e+003
+
+  !L = -11.2373*T**2 + 2.5326d+03*T + 1.4099d+06
+
+end subroutine latheat_NH3
+
+
+! ---------------------
+! H2O
+! ---------------------
+
+subroutine base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
+!
+!*******************************************************************************
+!
+!! BASE calculates quantities associated with the base Helmholtz function.
+!
+!
+!  Discussion:
+!
+!    The equation for the base Helmholtz function AB(T,RHO) is:
+!
+!      AB(T,RHO) = R * T * (
+!        - ln ( 1 - y ) 
+!        - ( beta - 1 ) / ( 1 - y ) 
+!        + ( alpha + beta + 1 ) / ( 2 * ( 1 - y )**2 )
+!        + 4 * y * ( ( Bbar / b ) - gamma ) 
+!        - 0.5 * ( alpha - beta + 3 ) 
+!        + ln ( RHO * R * T / P0 ) )
+!                                                      (Equation 2)
+!   where 
+!
+!     y = b * rho / 4, 
+!     alpha = 11,
+!     beta = 133/3, 
+!     gamma = 7/2, 
+!     P0 = 0.101325 MegaPascals = 1 atm
+!
+!   and
+!
+!     b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)**j  (Equation 3)
+!
+!     Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)**j               (Equation 4).
+!
+!   where 
+!
+!     T0=647.073 K and the coefficients b(j) and B(j) are
+!    
+!     j    b(j)                         B(j)
+!    --    -----------                  ----------
+!     0    0.7478629                    1.1278334
+!     1   -0.3540782                   -0.5944001
+!     2    0                           -5.010996
+!     3    0.007159876                  0
+!     4    0                            0.63684256
+!     5   -0.003528426                  0
+!
+!  For the derived quantities, the following relations are used:
+!
+!    Pressure:                  PB      = RHO**2 * dAB/dRHO
+!    Density derivative:        DPDRB   = 2*PB/RHO + RHO**2 * d2AB/dRHO2
+!    Temperature derivative:    DPDTB   = RHO**2 * d2AB/(dRHO dT)
+!    Specific entropy:          SB      = ( UB - AB ) / T
+!    Specific internal energy:  UB      = AB + T * SB
+!    Specific enthalpy:         HB      = UB + PB / RHO
+!    Specific heat capacity
+!      at constant volume:      CVB     = - T * d2AB/dT2
+!    Specific Gibbs function:   GB      = AB + PB / RHO
+!
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Modified:
+!
+!    03 February 2002
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Input, double precision RHO, the density, in G/CM3.
+!
+!    Output, double precision AB, the base value of the Helmholtz function,
+!    in KJ/kg.
+!
+!    Output, double precision CVB, the base value of the isochoric (constant 
+!    volume) heat capacity, in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision DPDRB, the base value of the partial
+!    derivative dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
+!
+!    Output, double precision DPDTB, the base value of the partial
+!    derivative dP(T,RHO)/dT, with RHO held fixed, in 
+!    MegaPascals/degrees Kelvin.
+!
+!    Output, double precision GB, the base value of the Gibbs free energy,
+!    in KJ/kg.
+!
+!    Output, double precision HB, the base value of enthalpy, in KJ/kg.
+!
+!    Output, double precision PB, the base pressure, in MegaPascals.
+!
+!    Output, double precision SB, the base value of entropy, 
+!    in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision UB, the base value of internal energy, 
+!    in KJ/kg.
+!
+  implicit none
+!
+  double precision ab
+  double precision, parameter :: alpha = 11.0D+00
+  double precision b1
+  double precision b1t
+  double precision b1tt
+  double precision b2
+  double precision b2t
+  double precision b2tt
+  double precision, parameter :: beta = 44.333333333333D+00
+  double precision cvb
+  double precision dpdrb
+  double precision dpdtb
+  double precision dz
+  double precision dz0
+  double precision, parameter :: gamma = 3.5D+00
+  double precision gascon
+  double precision gb
+  double precision hb
+  double precision, parameter :: p_zero = 0.101325D+00
+  double precision pb
+  double precision rho
+  double precision sb
+  double precision t
+  double precision ub
+  double precision x
+  double precision y
+  double precision z
+  double precision z0
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'BASE - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+!
+!  Refuse to handle zero or negative density.
+!
+  if ( rho <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'BASE - Fatal error!'
+    write ( *, '(a)' ) '  The input density RHO must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was RHO = ', rho
+    stop
+  end if
+!
+!  Compute auxilliary quantities for Equation 2.
+!
+  call bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
+
+  y = 0.25D+00 * b1 * rho
+
+  x = 1.0D+00 - y
+!
+!  Evaluate Equation 2.
+!
+  ab =   - log ( 1.0D+00 - y ) &
+         - ( beta - 1.0D+00 ) / ( 1.0D+00 - y ) &
+         + ( alpha + beta + 1.0D+00 ) / ( 2.0D+00 * ( 1.0D+00 - y )**2 ) &
+         + 4.0D+00 * y * ( ( b2 / b1 ) - gamma ) &
+         - 0.5D+00 * ( alpha - beta + 3.0D+00 ) &
+         + log ( rho * gascon() * t / p_zero )
+!
+!  Determine quantities defined in terms of AB.
+!
+  pb = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3 &
+    + 4.0D+00 * y * ( b2 / b1 - gamma )
+
+  z0 = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3
+
+  z = z0 + 4.0D+00 * y * ( b2 / b1 - gamma )
+
+  dz0 = ( alpha + 2.0D+00 * beta * y ) / ( 1.0D+00 - y )**3 &
+    + 3.0D+00 * ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**4
+
+  dz = dz0 + 4.0D+00 * ( b2 / b1 - gamma )
+
+  gb = ab + pb
+
+  ub = - t * b1t * ( pb - 1.0D+00 - rho * b2 ) / b1 - rho * t * b2t
+
+  hb = pb + ub
+!
+!  An incorrect version of this equation began:
+!
+!    cvb = 2.0D+00 * ub + ( pb - 1.0D+00 ) &
+!
+!  and caused me no end of trouble.  My fault, JVB, 03 February 2002
+!
+  cvb = 2.0D+00 * ub + ( z0 - 1.0D+00 ) &
+    * ( ( t * b1t / b1 )**2 - t**2 * b1tt / b1 ) &
+    - rho * t**2 * ( b2tt - gamma * b1tt ) - ( t * b1t / b1 )**2 * y * dz0
+
+  dpdtb = pb / t + rho * ( 0.25D+00 * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) ) &
+    * b1t + b2t - b2 / b1 * b1t )
+
+  sb = ub - ab
+
+  dpdrb = pb + y * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) )
+!
+!  Assign dimensions.
+!
+  ab =    gascon() * t       * ab
+  cvb =   gascon()           * cvb
+  dpdrb = gascon() * t       * dpdrb
+  dpdtb = gascon() * t * rho * dpdtb
+  gb =    gascon() * t       * gb
+  hb =    gascon() * t       * hb
+  pb =    gascon() * t * rho * pb
+  sb =    gascon()           * sb
+  ub =    gascon() * t       * ub
+
+  return
+end subroutine base
+
+subroutine bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
+!
+!*******************************************************************************
+!
+!! BB calculates the B's of equations 3 and 4.  
+!
+!
+!  Discussion:
+!
+!    Here
+!
+!      b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)**j  (Equation 3)
+!
+!      Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)**j               (Equation 4).
+!
+!    where 
+!
+!      T0 = 647.073 K 
+!
+!    and the coefficients b(j) and B(j) are
+!    
+!      j    b(j)                         B(j)
+!     --    -----------                  ----------
+!      0    0.7478629                    1.1278334
+!      1   -0.3540782                   -0.5944001
+!      2    0                           -5.010996
+!      3    0.007159876                  0
+!      4    0                            0.63684256
+!      5   -0.003528426                  0
+!
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision B1, the coefficient b from equation 3, 
+!    in CM3/G.
+!
+!    Output, double precision B2, the coefficient Bbar from equation 4, 
+!    in CM3/G.
+!
+!    Output, double precision B1T, the derivative dB1/dT, 
+!    in (CM3)/(G Degrees Kelvin).
+!
+!    Output, double precision B2T, the derivative dB2/dT, 
+!    in (CM3)/(G Degrees Kelvin).
+!
+!    Output, double precision B1TT, the second derivative of B1 with 
+!    respect to T, in (CM3)/(G (Degrees Kelvin)**2 ).
+!
+!    Output, double precision B2TT, the second derivative of B2 with 
+!    respect to T, in (CM3)/(G (Degrees Kelvin)**2 ).
+!
+  implicit none
+!
+  double precision b1
+  double precision b1t
+  double precision b1tt
+  double precision b2
+  double precision b2t
+  double precision b2tt
+  double precision, parameter, dimension ( 10 ) :: bp = (/ &
+    0.7478629D+00,   -0.3540782D+00,    0.0D+00,           0.0D+00, &
+    0.007159876D+00,  0.0D+00,         -0.003528426D+00,   0.0D+00, &
+    0.0D+00,          0.0D+00 /)
+  double precision, parameter, dimension ( 10 ) :: bq = (/ &
+    1.1278334D+00,    0.0D+00,         -0.5944001D+00,   -5.010996D+00, &
+    0.0D+00,          0.63684256D+00,   0.0D+00,          0.0D+00, &
+    0.0D+00,          0.0D+00 /)
+  integer i
+  double precision t
+  double precision, parameter :: t_ref = 647.073D+00
+  double precision v(10)
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'BB - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+!
+!  Set V(I) = ( T_REF / T )**(I-1).
+!
+  v(1) = 1.0D+00
+  do i = 2, 10
+    v(i) = v(i-1) * t_ref / t
+  end do
+!
+!  Set B1, B1T, B1TT.
+!
+  b1 = bp(1) + bp(2) * log ( 1.0D+00 / v(2) )
+  b1t = bp(2) * v(2) / t_ref
+  b1tt = 0.0D+00
+  do i = 3, 10
+    b1 = b1 + bp(i) * v(i-1)
+    b1t = b1t - dble ( i - 2 ) * bp(i) * v(i-1) / t
+    b1tt = b1tt + bp(i) * dble ( i - 2 )**2 * v(i-1) / t**2
+  end do
+
+  b1tt = b1tt -  ( b1t / t )
+!
+!  Set B2, B2T, B2TT.
+!
+  b2 = bq(1)
+  b2t = 0.0D+00
+  b2tt = 0.0D+00
+  do i = 3, 10
+    b2 = b2 + bq(i) * v(i-1)
+    b2t = b2t - dble ( i - 2 ) * bq(i) * v(i-1) / t
+    b2tt = b2tt + bq(i) * dble ( i - 2 )**2 * v(i-1) / t**2
+  end do
+
+  b2tt = b2tt - ( b2t / t )
+
+  return
+end subroutine bb
+
+
+subroutine ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
+!
+!*******************************************************************************
+!
+!! IDEAL computes ideal gas thermodynamic properties of water.
+!
+!
+!  Discussion:
+!
+!    Values for thermodynamic properties of water in the ideal
+!    gas state were reported by Woolley.  The formula for the ideal gas
+!    term of the Helmholtz function approximates a term by term summation of 
+!    contributions from each of the rotation and vibration states.  
+!    The formula, equation #6 in the reference, is:
+!   
+!    A(ideal)(T) = -R * T * ( 1 + ( C(1)/Tr + C(2) ) * ln(Tr)
+!      + Sum ( 3 <= I <= 18) C(I) * Tr**(I-6)
+!   
+!    where Tr=T/100 K.  The C(i) are tabulated coefficients.  Equation
+!    6 can be used for temperatures below 3000 K, and is accurate to
+!    within the tolerance of the gas constant for 50<=T<=2000 K.
+!     
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision AI, the Helmholtz function, in KJ/kg.
+!
+!    Output, double precision CPI, the heat capacity at constant pressure,
+!    in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision CVI, the heat capacity at constant volume,
+!    in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision GI, the Gibbs free energy, in KJ/kg.
+!
+!    Output, double precision HI, the enthalpy, in KJ/kg.
+!
+!    Output, double precision SI, the entropy, in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision UI, the internal energy, in KJ/kg.
+!
+  implicit none
+!
+  double precision ai
+  double precision, parameter, dimension ( 18 ) :: c = (/ &
+   19.730271018D+00,      20.9662681977D+00,     -0.483429455355D+00, &
+    6.05743189245D+00,    22.56023885D+00,       -9.87532442D+00, &
+   -4.3135538513D+00,      0.458155781D+00,      -0.047754901883D+00, &
+    0.0041238460633D+00,  -0.00027929052852D+00,  0.14481695261D-04, &
+   -0.56473658748D-06,     0.16200446D-07,       -0.3303822796D-09, &
+    0.451916067368D-11,   -0.370734122708D-13,    0.137546068238D-15 /)
+  double precision cpi
+  double precision cvi
+  double precision gascon
+  double precision gi
+  double precision hi
+  integer i
+  double precision si
+  double precision t
+  double precision temp
+  double precision tt
+  double precision ui
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'IDEAL - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+
+  tt = t / 100.0D+00
+
+  gi = - ( c(1) / tt + c(2) ) * log ( tt )
+  do i = 3, 18
+    gi = gi - c(i) * tt**(i-6)
+  end do
+
+  hi = c(2) + c(1) * ( 1.0D+00 - log ( tt ) ) / tt
+  do i = 3, 18
+    hi = hi + dble ( i - 6 ) * c(i) * tt**(i-6)
+  end do
+
+  cpi = c(2) - c(1) / tt
+  do i = 3, 18
+    cpi = cpi + dble ( ( i - 6 ) * ( i - 5 ) ) * c(i) * tt**(i-6)
+  end do
+
+  ai = gi - 1.0D+00
+  ui = hi - 1.0D+00
+  cvi = cpi - 1.0D+00
+  si = hi - gi
+!
+!  Assign dimensions.
+!
+  ai =  gascon() * t * ai
+  cpi = gascon()     * cpi
+  cvi = gascon()     * cvi
+  gi =  gascon() * t * gi
+  hi =  gascon() * t * hi
+  si =  gascon()     * si
+  ui =  gascon() * t * ui
+
+  return
+end subroutine ideal
+
+
+subroutine resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
+!
+!*******************************************************************************
+!
+!! RESID calculates residual contributions to thermodynamic quantities.
+!
+!
+!  Discussion:
+!
+!    The residual function consists of 40 terms.  The first 36 are
+!    used in a global least squares fit to experimental data.
+!
+!    Three terms were added that contribute only in the immediate
+!    neighborhood of the critical point 
+!      (tk-5) <= T <= (tk+5) C
+!      0.20   <= rho <= 0.44 g/cm3, 
+!
+!    A single term was added for the region of high pressure and 
+!    low temperature: T < 75 C, P > 300 MPa.
+!
+!    Except in these limited regions, the residual function is
+!    given by the first 36 terms.  The equation is
+!   
+!      A(residual)(rho,T)=
+!        sum(i=1 to 36) (g(i)/k(i)) * (T0/T)**(l(i)) (1-exp(-rho))**(k(i))
+!      + sum(i=37 to 40) g(i)*delta(i)**(k(i))
+!        * exp(-alpha(i)*delta(i)**(k(i)) - beta(i)*tau(i)**2)
+!                                                     (Equation 5)
+!   
+!    where
+!
+!      g(i) are coefficients determined by fits to data,
+!      delta(i) are reduced densities (delta(i)=((rho-rho(i))/rho(i))
+!      tau(i) are reduced temperatures (tau(i)=((T-tau(i))/tau(i))
+!      rho(i) are specified densities.
+!      tau(i) are specified temperatures.
+!      The k(i) and l(i) are specified integers.
+!   
+!  Modified:
+!
+!    22 November 1998
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Input, double precision RHO, the density, in G/CM3.
+!
+!    Output, double precision AR, the residual contribution to the 
+!    Helmholtz function, in KJ/kg.
+!
+!    Output, double precision CVR, the residual contribution to the 
+!    isochoric (constant volume) heat capacity, in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision DPDRR, the residual contribution to 
+!    the partial derivative dP(T,RHO)/dRHO, with T held fixed, in 
+!    (MegaPascals CM3)/G.
+!
+!    Output, double precision DPDTR, the residual contribution to 
+!    the partial derivative dP(T,RHO)/dT, with RHO held fixed, 
+!    in MegaPascals/degrees Kelvin.
+!
+!    Output, double precision GR, the residual contribution to the Gibbs 
+!    function, in KJ/kg.
+!
+!    Output, double precision HR, the residual contribution to the 
+!    enthalpy, in KJ/kg.
+!
+!    Output, double precision PR, the residual contribution to the pressure, 
+!    in MegaPascals.
+!
+!    Output, double precision SR, the residual contribution to the entropy, 
+!    in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision UR, the residual contribution to the 
+!    internal energy, in KJ/kg.
+!
+  implicit none
+!
+  double precision, parameter, dimension ( 4 ) :: aad = (/ &
+    34.0D+00, 40.0D+00, 30.0D+00, 1050.0D+00 /)
+  double precision, parameter, dimension ( 4 ) :: aat = (/ &
+    20000.0D+00, 20000.0D+00, 40000.0D+00, 25.0D+00 /)
+  double precision, parameter, dimension ( 4 ) :: adz = (/ &
+    0.319D+00, 0.319D+00, 0.319D+00, 1.55D+00 /)
+  double precision ar
+  double precision att
+  double precision, parameter, dimension ( 4 ) ::  atz = (/ &
+    640.0D+00, 640.0D+00, 641.6D+00, 270.0D+00 /)
+  double precision cvr
+  double precision dadt
+  double precision ddz
+  double precision del
+  double precision dex
+  double precision dfdt
+  double precision dpdrr
+  double precision dpdtr
+  double precision e
+  double precision errtol
+  double precision ex0
+  double precision ex1
+  double precision ex2
+  double precision fct
+  double precision, parameter, dimension ( 40 ) :: g = (/ &
+    -530.62968529023D+00,  0.22744901424408D+04, 0.78779333020687D+03, &
+    -69.830527374994D+00,  0.17863832875422D+05,-0.39514731563338D+05, &
+    0.33803884280753D+05, -0.13855050202703D+05,-0.25637436613260D+06, &
+    0.48212575981415D+06, -0.34183016969660D+06, 0.12223156417448D+06, &
+    0.11797433655832D+07, -0.21734810110373D+07, 0.10829952168620D+07, &
+   -0.25441998064049D+06, -0.31377774947767D+07, 0.52911910757704D+07, &
+   -0.13802577177877D+07, -0.25109914369001D+06, 0.46561826115608D+07, &
+   -0.72752773275387D+07,  0.41774246148294D+06, 0.14016358244614D+07, &
+   -0.31555231392127D+07,  0.47929666384584D+07, 0.40912664781209D+06, &
+   -0.13626369388386D+07,  0.69625220862664D+06,-0.10834900096447D+07, &
+   -0.22722827401688D+06,  0.38365486000660D+06, 0.68833257944332D+04, &
+    0.21757245522644D+05, -0.26627944829770D+04,-0.70730418082074D+05, &
+   -0.225D+00, -1.68D+00, 0.055D+00, -93.0D+00 /)
+  double precision gascon
+  double precision gr
+  double precision hr
+  integer i
+  integer, parameter, dimension ( 40 ) :: ii = (/ &
+    0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6, &
+    8,8,8,8,2,2,0,4,2,2,2,4 /)
+  integer j
+  integer, parameter, dimension ( 40 ) :: jj = (/ &
+    2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,&
+    2,3,5,7,1,4,4,4,0,2,0,0 /)
+  integer k
+  integer l
+  integer nc
+  double precision pr
+  double precision q10
+  double precision q20
+  double precision q2a
+  double precision q5t
+  double precision qm
+  double precision qp
+  double precision qr(11)
+  double precision qt(10)
+  double precision rho
+  double precision sr
+  double precision, parameter :: s_ref = 7.6180720166752D+00
+  double precision t
+  double precision, parameter :: t_ref = 647.073D+00
+  double precision tau
+  double precision tx
+  double precision, parameter :: u_ref = - 4328.4549774261D+00
+  double precision ur
+  double precision v
+!
+  errtol = sqrt ( epsilon ( errtol ) )
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'RESID - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+!
+!  Refuse to handle zero or negative density.
+!
+  if ( rho <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'RESID - Fatal error!'
+    write ( *, '(a)' ) '  The input density RHO must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was RHO = ', rho
+    stop
+  end if
+
+  nc = 36
+  dpdrr = 0.0D+00
+  pr = 0.0D+00
+  ar = 0.0D+00
+  dadt = 0.0D+00
+  cvr = 0.0D+00
+  dpdtr = 0.0D+00
+
+  ex0 = - rho
+! ex0 = max ( ex0, - 225.0D+00 )
+! ex0 = min ( ex0, 225.0D+00 )
+  e = exp ( ex0 )
+
+  q10 = rho * rho * e
+  q20 = 1.0D+00 - e
+
+  qr(1) = 0.0D+00
+  qr(2) = q10
+  do i = 2, 10
+    qr(i+1) = qr(i) * q20
+  end do
+
+  v = t_ref / t
+  qt(1) = t / t_ref
+  do i = 2, 10
+    qt(i) = qt(i-1) * v
+  end do
+
+  do i = 1, nc
+
+    k = ii(i) + 1
+    l = jj(i)
+    qp = g(i) * qr(k+1) * qt(l+1)
+    pr = pr + qp
+
+    dpdrr = dpdrr + ( 2.0D+00 / rho - ( 1.0D+00 - e * dble ( k - 1 ) / &
+      ( 1.0D+00 - e ) ) ) * qp
+
+    ar = ar + g(i) * qr(k+2) * qt(l+1) / ( rho**2 * e * dble ( k ) &
+      * gascon ( ) * t )
+
+    dfdt = ( 1.0D+00 - e )**k * dble ( 1 - l ) * qt(l+2) / t_ref / dble ( k )
+
+    dadt = dadt + g(i) * dfdt
+
+    dpdtr = dpdtr + g(i) * dfdt * rho**2 * e * dble ( k ) / ( 1.0D+00 - e )
+
+    cvr = cvr + g(i) * dble ( l ) * dfdt / gascon()
+
+  end do
+
+  qp = 0.0D+00
+  q2a = 0.0D+00
+
+  do j = 37, 40
+
+    k = ii(j)
+    ddz = adz(j-36)
+    del = rho / ddz - 1.0D+00
+
+    if ( abs ( del ) < errtol ) then
+      del = errtol
+    end if
+
+    ex1 = - aad(j-36) * del**k
+!   ex1 = max ( ex1, - 225.0D+00 )
+!   ex1 = min ( ex1, 225.0D+00 )
+    dex = exp ( ex1 ) * del**jj(j)
+
+    att = aat(j-36)
+    tx = atz(j-36)
+    tau = ( t / tx ) - 1.0D+00
+
+    ex2 = - att * tau**2
+!   ex2 = max ( ex2, - 225.0D+00 )
+!   ex2 = min ( ex2, 225.0D+00 )
+    q10 = dex * exp ( ex2 )
+
+    qm = dble ( jj(j) ) / del - dble ( k ) * aad(j-36) * del**(k-1)
+    fct = qm * rho**2 * q10 / ddz
+
+    q5t = fct * ( 2.0D+00 / rho + qm / ddz ) - ( rho / ddz )**2 * q10 * &
+      ( dble ( jj(j) ) / del**2 + dble ( k * ( k - 1 ) ) * aad(j-36) * &
+      del**(k-2) )
+
+    dpdrr = dpdrr + q5t * g(j)
+    qp = qp + g(j) * fct
+    dadt = dadt - 2.0D+00 * g(j) * att * tau * q10 / tx
+    dpdtr = dpdtr - 2.0D+00 * g(j) * att * tau * fct / tx
+
+    q2a = q2a + t * g(j) * att * ( 4.0D+00 * ex2 + 2.0D+00 ) * q10 / tx**2
+
+    ar = ar + q10 * g(j) / ( gascon() * t )
+
+  end do
+
+  cvr = cvr + q2a / gascon()
+  pr = pr + qp
+  sr = - dadt / gascon()
+  ur = ar + sr
+!
+!  Assign dimensions.
+!
+  ar =  gascon() * t *  ar
+  cvr = gascon() *     cvr
+  sr =  gascon() *      sr
+  ur =  gascon() * t *  ur
+!
+!  Adjust energies.
+!
+  ar = ar + gascon ( ) * t * s_ref - gascon ( ) * u_ref
+  sr = sr - gascon ( ) * s_ref
+  ur = ur - gascon ( ) * u_ref
+
+  gr = ar + pr / rho
+  hr = ur + pr / rho
+
+  return
+end subroutine resid
+
+subroutine psat_H2O ( t, p )
+!
+!*******************************************************************************
+!
+!! PSAT_H2O makes a rough estimate of the vapor pressure.
+!
+!
+!  Discussion:
+!
+!    The calculation agrees with tabulated data to within
+!    0.02% for temperature to within a degree or so of the critical
+!    temperature.  The approximate vapor pressure can be refined
+!    by imposing the condition that the Gibbs functions of the vapor
+!    and liquid phases be equal.
+!
+!  Modified:
+!
+!    21 November 1998
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision P, the vapor pressure, in MegaPascals.
+!
+  implicit none
+!
+  double precision, parameter, dimension ( 8 ) :: a = (/ &
+    -7.8889166D+00,   2.5514255D+00,   -6.716169D+00, 33.239495D+00, &
+    -105.38479D+00,   174.35319D+00,  -148.39348D+00, 48.631602D+00 /)
+  double precision b
+  integer i
+  double precision p
+  double precision q
+  double precision t
+  double precision, parameter :: t_ref = 647.25D+00
+  double precision v
+  double precision w
+  double precision z
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'PSAT_H2O - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+
+  if ( t <= 314.0D+00 ) then
+
+    p = 0.1D+00 * exp ( 6.3573118D+00 - 8858.843D+00 / t &
+      + 607.56335D+00 * t**( -0.6D+00 ) )
+
+  else
+
+    v = t / t_ref
+    w = abs ( 1.0D+00 - v )
+    b = 0.0D+00
+    do i = 1, 8
+      z = i
+      b = b + a(i) * w**( ( z + 1.0D+00 ) / 2.0D+00 )
+    end do
+
+    q = b / v
+    p = 22.093D+00 * exp ( q )
+
+  end if
+
+  return
+end subroutine psat_H2O
+
+
+subroutine tdpsdt ( t, dp )
+!
+!*******************************************************************************
+!
+!! TDPSDT computes the quantity T * dP(Sat)/dT.
+!
+!
+!  Discussion:
+!
+!    Here T is the temperature and P(Sat) is the vapor pressure.  
+!    It is used by TSAT_EST and TSAT.
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Output, double precision DP, the value T*(dP(Sat)/dT), 
+!    in MegaPascals.
+!
+  implicit none
+!
+  double precision, parameter, dimension ( 8 ) :: a = (/ &
+      -7.8889166D+00,   2.5514255D+00,   -6.716169D+00, 33.239495D+00, &
+    -105.38479D+00,   174.35319D+00,   -148.39348D+00,  48.631602D+00 /)
+  double precision b
+  double precision c
+  double precision dp
+  integer i
+  double precision q
+  double precision t
+  double precision, parameter :: t_ref = 647.25D+00
+  double precision v
+  double precision w
+  double precision y
+  double precision z
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'TDPSDT - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+
+  v = t / t_ref
+  w = 1.0D+00 - v
+  b = 0.0D+00
+  c = 0.0D+00
+  do i = 1, 8
+    z = dble ( i + 1 ) / 2.0D+00
+    y = a(i) * w**z
+    c = c + ( y / w ) * ( 0.5D+00 - 0.5D+00 * dble ( i ) - 1.0D+00 / v )
+    b = b + y
+  end do
+
+  q = b / v
+  dp = 22.093D+00 * exp ( q ) * c
+
+  return
+end subroutine tdpsdt
+
+
+
+subroutine therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
+!
+!*******************************************************************************
+!
+!! THERM calculates thermodynamic functions given temperature and density.
+!
+!
+!  Discussion:
+!
+!    Thermodynamic values were calculated from an analytic equation
+!    that approximates the Helmholtz function (specific Helmholtz
+!    energy) for ordinary water and steam, of the form A=A(rho,T)
+!    where A is the Helmholtz function, rho the density, and T
+!    the absolute (thermodynamic) temperature.  Any thermodynamic
+!    value for any state, liquid, vapor or metastable, may be
+!    calculated by differentiation of this equation in accord with
+!    the first and second laws of thermodynamics.
+! 
+!    The International Association for the Properties of Steam
+!    has provisionally accepted this formulation for the range
+!    273.15 <= T <= 1273.15 degrees Kelvin, where, for 423.15 <= T,
+!    the maximum pressure is Pmax = 1500 MPa = 15000 bar, and for
+!    273.15 <= T < 423.15, the maximum pressure is
+!    Pmax = 100 * (5 + (T-273.15)/15) MPa.
+! 
+!    Close to the critical point, a small region is excluded:
+!    Abs(T-Tk) < 1, abs((rho-rhok)/rhok) < 0.3.
+! 
+!    The equation has a wider useful range, namely, fluid states
+!    of pure, undissociated water and steam defined by
+!    260 <= T <= 2500 K and 0 <= P <= 3000 MPa.
+! 
+!    Thermodynamic property values for specific volume, density,
+!    specific internal energy, specific enthalpy, and specific
+!    entropy of water and steam were tabulated over the range
+!    0 <= t <= 2000 C, 0 <= P <= 3000 MPa.  The reference
+!    state is the liquid at the triple point, for which the
+!    internal energy and entropy have been assigned the value zero.
+!     
+!    Thermodynamic quantities are determined from the Helmholtz function
+!    A(rho,T), which is computed as the sum of three terms:
+!
+!      A(rho,T) = A(base)(rho,T) + A(residual)(rho,T) + A(ideal)(T)
+!                                                       (Equation 1)
+!
+!    Because A(rho,T) is everywhere single valued and analytic,
+!    we can derive closed form relations for all other properties.
+!    In the following, unless otherwise indicated, the independent
+!    variables are temperature T and density RHO, and differentiation
+!    with respect to one variable is to imply that the other is fixed.
+! 
+!    Pressure:                  P       = RHO**2 * dA/dRHO
+!    Density derivative:        dP/dRHO = 2*P/RHO + RHO**2 * d2A/dRHO2
+!    Temperature derivative:    dP/dT   = RHO**2 * d2A/(dRHO dT)
+!    Specific entropy:          S       = - dA/dT
+!    Specific internal energy:  U       = A + T*S
+!    Specific enthalpy:         H       = U + P/RHO
+!    Specific heat capacity
+!      at constant volume:      Cv      = - T * d2A/dT2
+!    Specific Gibbs function:   G       = A + P/RHO
+!    Specific heat capacity
+!      at constant pressure:    Cp      = Cv + (T*(dP/dT)**2)/(RHO**2*dP/dRHO)
+!    Speed of sound:            Omega   = Sqrt ((Cp/Cv) * dP/dRHO)
+!    Second virial coefficient: B       = 1/(2*R*T) * (d2P/dRHO2) (at RHO=0)
+!    Isothermal Joule-Thomson
+!      coefficient:             DeltaT  = (dH/dP) (fixed T) =
+!                                         (1/RHO)-(T*dP/dT)/(RHO**2*dP/dRHO)
+!    Joule-Thomson coefficient: Mu      = (dT/dP) (fixed H) = DeltaT/Cp
+!    Isentropic temperature-
+!      pressure coefficient:    BetaS   = (dT/dP) (fixed S) =
+!                                         (DeltaT - 1/RHO)/Cp
+!   
+!  Modified:
+!
+!    19 November 1998
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Input, double precision T, the temperature, in degrees Kelvin.
+!
+!    Input, double precision RHO, the fluid density, in G/CM3.
+!
+!    Output, double precision A, the Helmholtz function, in KJ/kg.
+!
+!    Output, double precision CJTH, the Joule-Thomson coefficient,
+!    in K/MegaPascals.
+!
+!    Output, double precision CJTT, the isothermal Joule-Thomson coefficient,
+!    in CM3/G.
+!
+!    Output, double precision CP, the isobaric (constant pressure) heat
+!    capacity, in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision CV, the isochoric (constant volume) heat capacity,
+!    in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision DPDR, the partial derivative 
+!    dP(T,RHO)/dRHO, with T held fixed, in MegaPascals*CM3/G.
+!
+!    Output, double precision DPDT, the partial derivative 
+!    dP(T,RHO)/dT, with RHO held fixed, in MegaPascals/degrees Kelvin.
+!
+!    Output, double precision G, the Gibbs free energy, in KJ/kg.
+!
+!    Output, double precision H, the enthalpy, in KJ/kg.
+!
+!    Output, double precision P, the pressure, in MegaPascals.
+!
+!    Output, double precision S, the entropy, in KJ/(kg degrees Kelvin).
+!
+!    Output, double precision U, the internal energy, in KJ/kg.
+!
+  implicit none
+!
+  double precision a
+  double precision ab
+  double precision ai
+  double precision ar
+  double precision cjth
+  double precision cjtt
+  double precision cp
+  double precision cpi
+  double precision cv
+  double precision cvb
+  double precision cvi
+  double precision cvr
+  logical, parameter :: debug = .false.
+!  logical, parameter :: debug = .true.
+  double precision dpdr
+  double precision dpdrb
+  double precision dpdrr
+  double precision dpdt
+  double precision dpdtb
+  double precision dpdtr
+  double precision g
+  double precision gb
+  double precision gi
+  double precision gr
+  double precision h
+  double precision hb
+  double precision hi
+  double precision hr
+  double precision p
+  double precision pb
+  double precision pr
+  double precision rho
+  double precision s
+  double precision sb
+  double precision si
+  double precision sr
+  double precision t
+  double precision u
+  double precision ub
+  double precision ui
+  double precision ur
+!
+!  Refuse to handle zero or negative temperatures.
+!
+  if ( t <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'THERM - Fatal error!'
+    write ( *, '(a)' ) '  The input temperature T must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was T = ', t
+    stop
+  end if
+!
+!  Refuse to handle zero or negative density.
+!
+  if ( rho <= 0.0D+00 ) then
+    write ( *, '(a)' ) ' '
+    write ( *, '(a)' ) 'THERM - Fatal error!'
+    write ( *, '(a)' ) '  The input density RHO must be positive.'
+    write ( *, '(a,g14.6)' ) '  Input value was RHO = ', rho
+    stop
+  end if
+
+  call ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
+
+  call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
+
+  call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
+
+  a =       ab +    ar +  ai
+  cv =     cvb +   cvr + cvi
+
+  if ( debug ) then
+    write ( *, * ) ' '
+    write ( *, * ) 'THERM:'
+    write ( *, * ) '  CVB = ', cvb
+    write ( *, * ) '  CVR = ', cvr
+    write ( *, * ) '  CVI = ', cvi
+    write ( *, * ) '  CV  = ', cv
+  end if
+
+
+  dpdr = dpdrb + dpdrr
+  dpdt = dpdtb + dpdtr
+  p =       pb +    pr
+  s =       sb +    sr +  si
+  u =       ub +    ur +  ui
+
+  if ( debug ) then
+    write ( *, * ) ' '
+    write ( *, * ) 'THERM:'
+    write ( *, * ) '  UB = ', ub
+    write ( *, * ) '  UR = ', ur
+    write ( *, * ) '  UI = ', ui
+  end if
+
+  g = a + p / rho
+  h = u + p / rho
+  cp = cv + t * dpdt**2 / ( dpdr * rho**2 )
+  cjtt = 1.0D+00 / rho - t * dpdt / ( dpdr * rho**2 )
+  cjth = - cjtt / cp
+
+  return
+end subroutine therm
+
+function gascon ( )
+!
+!*******************************************************************************
+!
+!! GASCON returns the value of the specific gas constant.
+!
+!
+!  Note:
+!
+!    The specific gas constant R is related to the universal gas
+!    constant R-bar = 8.31441 J/(mol degrees Kelvin) by the molar mass 
+!    M = 18.0152 g/mol:
+!
+!      R = R-bar / M.
+!
+!  Reference:
+!
+!    Lester Haar, John Gallagher and George Kell,
+!    NBS/NRC Steam Tables:
+!    Thermodynamic and Transport Properties and Computer Programs
+!      for Vapor and Liquid States of Water in SI Units,
+!    Hemisphere Publishing Corporation, Washington, 1984,
+!    TJ270.H3
+!
+!    C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
+!    ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
+!    American Society of Mechanical Engineers, 1967.
+!
+!  Parameters:
+!
+!    Output, double precision GASCON, the value of the specific gas 
+!    constant, in J/(g degrees Kelvin).
+!
+  implicit none
+!
+  double precision gascon
+!
+  gascon = 0.461522D+00
+ 
+  return
+end function gascon
Index: /trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F	(revision 305)
@@ -534,5 +534,5 @@
       IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>'
-         CALL abort
+         !CALL abort
       ENDIF
       xmin = 1.0E+20
@@ -541,5 +541,4 @@
       xmax = MAXVAL(totcloudfrac)
       PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax
-
 
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 305)
@@ -8,5 +8,5 @@
  
       use radinc_h, only : naerkind
-      use watercommon_h, only : mx_eau_sol, To, RLVTT
+      use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
      
       implicit none
@@ -392,6 +392,6 @@
       integer igtest
 
-!     included by RW for pressure broadening broadening study
-      integer gascut
+!     included by RW for runway greenhouse 1D study
+      real muvar(ngridmx,nlayermx+1)
 
 !     included by RW for variable H2O particle sizes
@@ -404,4 +404,5 @@
 !     included by RW for sourceevol
       real ice_initial(ngridmx)!, isoil
+      real delta_ice,ice_tot
       real ice_min(ngridmx)
       save ice_min
@@ -502,5 +503,5 @@
             close(128)
         
-            if(num_run.ne.0.and.mod(num_run,6).eq.0)then
+            if(num_run.ne.0.and.mod(num_run,3).eq.0)then
                print*,'Updating ice at end of this year!'
                ice_update=.true.
@@ -713,17 +714,8 @@
 
               if(kastprof)then
-                 ISR=fract(1)*Fat1AU/dist_star**2/2.0
-
-                 call simpleprof_fn(tsurf,pt,zpopsk,ISR,pplev(1,1))
-                 open(116,file='profpres.out',form='formatted')
-                 open(117,file='proftemp.out',form='formatted')
-                 do l=1,nlayer
-                    write(116,*) pplay(1,l)
-                    write(117,*) pt(1,l)
-                 enddo
-                 close(116)
-                 close(117)
+                 print*,'kastprof should not = true here'
+                 call abort
               endif
-
+              muvar(:,:)=0.0 ! only used for climate evolution for now
 
 !             Option to call scheme once for clear regions, once for cloudy regions (BC)
@@ -731,17 +723,17 @@
 
                  clearsky=.true.
-                 call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,           &
+                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
                       albedo,emis,mu0,pplev,pplay,pt,                      &
-                      tsurf,fract,dist_star,aerosol,cpp3D,                 &
+                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
                       zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
-                      fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,ptime,pday,  &
+                      fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,             &
                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
 
                  clearsky=.false.
-                 call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,           &
+                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
                       albedo,emis,mu0,pplev,pplay,pt,                      &
-                      tsurf,fract,dist_star,aerosol,cpp3D,                 &
+                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
                       zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
-                      fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,ptime,pday,  &
+                      fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,             &
                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
 
@@ -767,9 +759,9 @@
 
                  clearsky=.false.
-                 call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,         &
+                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
                       albedo,emis,mu0,pplev,pplay,pt,                    &
-                      tsurf,fract,dist_star,aerosol,cpp3D,               &
+                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
                       zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                      fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,  &
+                      fluxabs_sw,fluxtop_dn,reffrad,tau_col,             &
                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
 
@@ -881,5 +873,4 @@
       endif ! of if (callrad)
 
-	
 !-----------------------------------------------------------------------
 !    4. Vertical diffusion (turbulent mixing):
@@ -906,5 +897,5 @@
               zdum1,zdum2,zdh,pdq,zflubid,           &
               zdudif,zdvdif,zdhdif,zdtsdif,q2,       &
-              zdqdif,zdqsdif)
+              zdqdif,zdqsdif,lastcall)
 
          do l=1,nlayer
@@ -1143,5 +1134,5 @@
          endif
 
-         call condens_co2cloud(ngrid,nlayer,nq,ptimestep,   &
+         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
               capcal,pplay,pplev,tsurf,pt,                  &
               pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
@@ -1213,5 +1204,5 @@
                call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
                DO l = 1, nlayer  
-                  DO ig = 1, ngrid        	
+                  DO ig = 1, ngrid 
                      pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
@@ -1626,9 +1617,7 @@
 
 !     Increment surface temperature
-      if(.not.kastprof)then
-         do ig=1,ngrid
-            tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
-         enddo
-      endif
+      do ig=1,ngrid
+         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+      enddo
 
 !     Compute soil temperatures and subsurface heat flux
@@ -1783,17 +1772,4 @@
          endif
 
-         if(kastprof)then ! for f(p,T) computation (not needed in universal model)
-            open(15,file='ene_bal.out',form='formatted',access='append')
-            write(15,*) pplev(1,1), Tsurf(1), ISR, OLR, ASR
-            close(15)
-            open(16,file='Seff.out',form='formatted')
-            write(16,*) OLR/ASR
-            close(16)
-            open(17,file='alb.out',form='formatted')
-            write(17,*) 1.0-ASR/ISR
-            close(17)
-            call abort
-         endif
-
       endif
 
@@ -1830,10 +1806,13 @@
          do ig=1,ngrid
             do l=1,nlayer
-               reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * reffrad(ig,l,1) *    &
-                    (pplev(ig,l) - pplev(ig,l+1)) / g
+               if(co2cond)then
+                  reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
+                                  reffrad(ig,l,1) *    &
+                                  (pplev(ig,l) - pplev(ig,l+1)) / g
+               endif
                if(water)then
-               reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * reffrad(ig,l,2) *    &
-                    (pplev(ig,l) - pplev(ig,l+1)) / g
-
+                  reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
+                                  reffrad(ig,l,2) *    &
+                                  (pplev(ig,l) - pplev(ig,l+1)) / g
                endif
             enddo
@@ -1923,23 +1902,38 @@
 !           Update surface ice distribution to iterate to steady state if requested
             if(ice_update)then
+
                do ig = 1, ngrid
 
-                  ! if ice amount is declining, set it to zero
-                  if(qsurf_hist(ig,igcm_h2o_ice).lt.ice_initial(ig))then
-                     qsurf_hist(ig,igcm_h2o_ice)=0.0 
-                     qsurf_hist(ig,igcm_h2o_vap)=0.0
+                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
+
+                  ! add one hundred years/timesteps of evolution
+                  qsurf_hist(ig,igcm_h2o_ice) = &
+                     qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0 
+
+                  ! if ice has gone -ve, set to zero
+                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
+                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+                     qsurf_hist(ig,igcm_h2o_vap) = 0.0
                   endif
- 
-                  ! if ice amount is increasing and was never less than 0.5 kg/m2,
-                  ! set value in gridbox to 1x10^3 kg/m2
-                  !if(qsurf_hist(ig,igcm_h2o_ice).gt.0.0 .and. ice_min(ig).gt.0.5)then
-                  if(qsurf_hist(ig,igcm_h2o_ice).gt.ice_initial(ig))then
-                     if(ice_min(ig).gt.0.5)then ! ignore seasonal ice
-                        qsurf_hist(ig,igcm_h2o_ice)=1.e3
-                        qsurf_hist(ig,igcm_h2o_vap)=0.0
-                     endif
+
+                  ! set maximum to ice
+                  if(qsurf_hist(ig,igcm_h2o_ice).gt.1.e3)then
+                     qsurf_hist(ig,igcm_h2o_ice) = 1.e3 
+                     qsurf_hist(ig,igcm_h2o_vap) = 0.0
                   endif
 
                enddo
+
+               ! enforce ice conservation
+               ice_tot=0.0
+               do ig = 1, ngrid
+                  ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
+               enddo
+               do ig = 1, ngrid
+                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
+               enddo
+
+
+
             endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 305)
@@ -74,4 +74,8 @@
 
 !      integer, parameter :: L_NPREF   = 8
+!      integer, parameter :: L_NTREF   = 11
+!      integer, parameter :: L_REFVAR  = 8   ! N2_H2Ovar
+
+!      integer, parameter :: L_NPREF   = 8
 !      integer, parameter :: L_NTREF   = 5
 !      integer, parameter :: L_REFVAR  = 1   ! therm_test2
@@ -109,7 +113,6 @@
       ! equivalent temperatures are 1/10 of these values
       integer, parameter :: NTstar = 500
-!      integer, parameter :: NTstop = 9000 ! for hotstuff runs
-      integer, parameter :: NTstop = 6000 ! for GJ581d / earlymars runs
-!      integer, parameter :: NTstop = 10000 ! for H2 warming runs
+      integer, parameter :: NTstop = 9000 ! new default for all non hot Jupiter runs
+      !integer, parameter :: NTstop = 6000 ! for GJ581d / earlymars runs
 
 ! Maximum number of grain size classes
Index: /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 304)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 305)
@@ -493,17 +493,12 @@
 
       if(autozlevs)then
-         if(kastprof)then
-            Hscale=10.0
-            Hmax=400.0
-         else
             
-            open(91,file="z2sig.def",form='formatted')
-            read(91,*) Hscale
-            DO ilayer=1,nlayer-2
-               read(91,*) Hmax
-            enddo
-            close(91)
- 
-         endif
+         open(91,file="z2sig.def",form='formatted')
+         read(91,*) Hscale
+         DO ilayer=1,nlayer-2
+            read(91,*) Hmax
+         enddo
+         close(91)
+ 
             
          print*,'Hmax = ',Hmax,' km'
Index: unk/LMDZ.GENERIC/libf/phystd/simpleprof_fn.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/simpleprof_fn.F90	(revision 304)
+++ 	(revision )
@@ -1,58 +1,0 @@
-      subroutine simpleprof_fn(tsurf,temp,popsk,ISR,psurf)
-!     creates a kasting temperature profile
-use radcommon_h, only: sigma
-      implicit none
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "callkeys.h"
-
-!      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
-!      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
-      REAL psurf,tsurf      
-      REAL temp(nlayermx)   ! temperature at the middle of the layers
-      real ISR, ALB
-
-      real popsk(nlayermx)  ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx)
-                            ! but kastprof requires 1D, so its ok I guess
-      integer l,iter
-
-
-      if(startype.eq.1)then
-         ALB=0.2+0.1*log10(psurf/100000.)
-      else
-         ALB=0.15
-      end if
-
-      Tstrat=((1-ALB)*ISR/(2*sigma))**0.25 ! skin temperature (gray approx.)
-      !Tstrat=60.                           ! minimum possible given H2-H2 CIA
-      print*,'psurf',psurf/100000.
-      print*,'ALB',ALB
-      print*,'Tstrat=',Tstrat
-
-      if(ngridmx.gt.1)then
-         print*,'DO NOT RUN KASTPROF IN 3D!!!'
-         call abort
-      end if
-
-      ! 1st, we make pot. temperature constant everywhere
-      do l=1,nlayermx
-            temp(l)=tsurf*popsk(l)
-      enddo
-
-      if(1.eq.1)then
-
-      do l=1,nlayermx-1
-            if (temp(l) .le. Tstrat)then
-               temp(l) = Tstrat            
-            endif
-      enddo
-
-      temp(nlayermx) = Tstrat
-
-      endif
-
-   
-   return
- end subroutine simpleprof_fn
Index: unk/LMDZ.GENERIC/libf/phystd/writediagspec.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspec.F	(revision 304)
+++ 	(revision )
@@ -1,279 +1,0 @@
-      subroutine writediagspec(ngrid,nom,titre,unite,dim,px)
-
-!  Ecriture de variables diagnostiques au choix dans la physique 
-!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
-!  3d (ex : temperature), 2d (ex : temperature de surface), ou
-!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
-!  solaire)
-!  Dans la version 2000, la periode d'ecriture est celle de 
-!  "ecritphy " regle dans le fichier de controle de run :  run.def
-!
-!    writediagfi peut etre appele de n'importe quelle subroutine
-!    de la physique, plusieurs fois. L'initialisation et la creation du
-!    fichier se fait au tout premier appel.
-!
-! WARNING : les variables dynamique (u,v,t,q,ps)
-!  sauvees par writediagfi avec une
-! date donnee sont legerement differentes que dans le fichier histoire car 
-! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
-! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
-! avant l'ecriture dans diagfi (cf. physiq.F)
-!  
-!
-!  parametres (input) :
-!  ----------
-!      ngrid : nombres de point ou est calcule la physique
-!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
-!                 (= nlon ou klon dans la physique terrestre)
-!      
-!      unit : unite logique du fichier de sortie (toujours la meme)
-!      nom  : nom de la variable a sortir (chaine de caracteres)
-!      titre: titre de la variable (chaine de caracteres)
-!      unite : unite de la variable (chaine de caracteres)
-!      px : variable a sortir (real 0, 2, ou 3d)
-!      dim : dimension de px : 0, 2, ou 3 dimensions
-!
-!=================================================================
-!
-!      This is a modified version that accepts spectrally varying input
-!      RW (2010)
-!
-!=================================================================
- 
-! Addition by RW (2010) to allow OLR to be saved in .nc format
-      use radinc_h, only : L_NSPECTI
-
-      implicit none
-
-! Commons
-#include "dimensions.h"
-#include "dimphys.h"
-#include "paramet.h"
-#include "control.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "description.h"
-#include "netcdf.inc"
-#include "temps.h"
-#include "surfdat.h"
-
-! Arguments on input:
-      integer ngrid
-      character (len=*) :: nom,titre,unite
-      integer dim
-      real px(ngrid,L_NSPECTI)
-
-! Local variables:
-
-!      real dx3(iip1,jjp1,llm) ! to store a 3D data set
-!      real dx2(iip1,jjp1)     ! to store a 2D (surface) data set
-!      real dx0
-
-      real date
-
-      REAL phis(ip1jmp1)
-
-      integer irythme
-      integer ierr
-      integer iq
-      integer i,j,l,zmax , ig0
-
-      integer zitau
-      character firstnom*20
-      SAVE firstnom
-      SAVE zitau
-      SAVE date
-      data firstnom /'1234567890'/
-      data zitau /0/
-
-! Ajouts
-      integer, save :: ntime=0
-      integer :: idim,varid
-      integer :: nid
-      character (len =50):: fichnom
-      integer, dimension(4) :: id
-      integer, dimension(4) :: edges,corner
-
-! added by RDW for OLR output
-       real dx3(iip1,jjp1,L_NSPECTI) ! to store the data set
-
-
-!***************************************************************
-!Sortie des variables au rythme voulu
-
-      irythme = int(ecritphy) ! sortie au rythme de ecritphy
-!     irythme = iconser  ! sortie au rythme des variables de controle
-!     irythme = iphysiq  ! sortie a tous les pas physique
-!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
-!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
-
-!***************************************************************
-
-! The following test is here to enforce that writediagfi is not used with the
-! 1D version of the GCM
-      if (ngrid.eq.1) return
-      
-c     nom=trim((nom))
-c     unite=trim((unite))
-c     titre=trim((titre))
-
-! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
-! ------------------------------------------------------------------------
-! (Au tout premier appel de la subroutine durant le run.)
-
-      fichnom="diagspec.nc"
-
-      if (firstnom.eq.'1234567890') then ! .true. for the very first call
-      !  to this subroutine; now set 'firstnom'
-         firstnom = nom
-         ! just to be sure, check that firstnom is large enough to hold nom
-         if (len_trim(firstnom).lt.len_trim(nom)) then
-           write(*,*) "writediagfi: Error !!!"
-           write(*,*) "   firstnom string not long enough!!"
-           write(*,*) "   increase its size to at least ",len_trim(nom)
-           stop
-         endif
-
-         ! Create the NetCDF file
-         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
-         ! Define the 'Time' dimension
-         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
-         ! Define the 'Time' variable
-#ifdef NC_DOUBLE
-         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
-#else
-         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
-#endif
-         ! Add a long_name attribute
-         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
-     .          4,"Time")
-         ! Add a units attribute
-         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
-     .          "days since 0000-00-0 00:00:00")
-         ! Switch out of NetCDF Define mode
-         ierr = NF_ENDDEF(nid)
-
-         ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
-!         call iniwrite(nid,day_ini,phis)
-         call iniwrite_spec(nid,day_ini,phis)
-
-         zitau = -1 ! initialize zitau
-      else
-         ! Open the NetCDF file
-         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
-      endif ! if (firstnom.eq.'1234567890')
-
-! Increment time index 'zitau' if it is the "firstcall" (at given time level)
-! to writediagfi
-!------------------------------------------------------------------------
-      if (nom.eq.firstnom) then
-          zitau = zitau + iphysiq
-      end if
-
-!--------------------------------------------------------
-! Write the variables to output file if it's time to do so
-!--------------------------------------------------------
-
-      if ( MOD(zitau+1,irythme) .eq.0.) then
-
-! Compute/write/extend 'Time' coordinate (date given in days)
-! (done every "first call" (at given time level) to writediagfi)
-! Note: date is incremented as 1 step ahead of physics time
-!       (like the 'histoire' outputs)
-!--------------------------------------------------------
-
-        if (nom.eq.firstnom) then
-
-        ! We have identified a "first call" (at given date)
-           ntime=ntime+1 ! increment # of stored time steps
-           ! compute corresponding date (in days and fractions thereof)
-           date= float (zitau +1)/float (day_step)
-           ! Get NetCDF ID of 'Time' variable
-           ierr= NF_INQ_VARID(nid,"Time",varid)
-
-           !print*,'in writediagfi_spec.F, time=',varid
-
-           ! Write (append) the new date to the 'Time' array
-#ifdef NC_DOUBLE
-           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
-#else
-           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
-#endif
-           if (ierr.ne.NF_NOERR) then
-              write(*,*) "***** PUT_VAR matter in writediagspec_nc"
-              write(*,*) "***** with time"
-              write(*,*) 'ierr=', ierr   
-c             call abort
-           endif
-
-           write(6,*)'WRITEDIAGSPEC: date= ', date
-        end if ! of if (nom.eq.firstnom)
-
-
- 
-!Case of a 3D variable
-!---------------------
-        if (dim.eq.3) then
-
-!         recast (copy) variable from physics grid to dynamics grid
-           DO l=1,L_NSPECTI
-             DO i=1,iip1
-                dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
-             ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
-                   dx3(i,j,l)=px(ig0+i,l)
-                ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
-             ENDDO
-           ENDDO
-
-! name of the variable
-           ierr= NF_INQ_VARID(nid,nom,varid)
-           if (ierr /= NF_NOERR) then
-! corresponding dimensions
-              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
-              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
-              ierr= NF_INQ_DIMID(nid,"bandsIR",id(3))
-              ierr= NF_INQ_DIMID(nid,"Time",id(4))
-
-! Create the variable if it doesn't exist yet
-
-              write (*,*) "=========================="
-              write (*,*) "DIAGSPEC: creating variable ",nom
-              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
-
-           endif
-
-           corner(1)=1
-           corner(2)=1
-           corner(3)=1
-           corner(4)=ntime
-
-           edges(1)=iip1
-           edges(2)=jjp1
-           edges(3)=L_NSPECTI
-           edges(4)=1
-#ifdef NC_DOUBLE
-           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
-#else
-           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
-#endif
-
-           if (ierr.ne.NF_NOERR) then
-              write(*,*) "***** PUT_VAR problem in writediagspec"
-              write(*,*) "***** with ",nom
-              write(*,*) 'ierr=', ierr
-             call abort
-           endif 
-
-        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
-
-      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
-
-      ierr= NF_CLOSE(nid)
-
-      end
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 305)
@@ -0,0 +1,279 @@
+      subroutine writediagspecIR(ngrid,nom,titre,unite,dim,px)
+
+!  Ecriture de variables diagnostiques au choix dans la physique 
+!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
+!  3d (ex : temperature), 2d (ex : temperature de surface), ou
+!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
+!  solaire)
+!  Dans la version 2000, la periode d'ecriture est celle de 
+!  "ecritphy " regle dans le fichier de controle de run :  run.def
+!
+!    writediagfi peut etre appele de n'importe quelle subroutine
+!    de la physique, plusieurs fois. L'initialisation et la creation du
+!    fichier se fait au tout premier appel.
+!
+! WARNING : les variables dynamique (u,v,t,q,ps)
+!  sauvees par writediagfi avec une
+! date donnee sont legerement differentes que dans le fichier histoire car 
+! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
+! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
+! avant l'ecriture dans diagfi (cf. physiq.F)
+!  
+!
+!  parametres (input) :
+!  ----------
+!      ngrid : nombres de point ou est calcule la physique
+!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
+!                 (= nlon ou klon dans la physique terrestre)
+!      
+!      unit : unite logique du fichier de sortie (toujours la meme)
+!      nom  : nom de la variable a sortir (chaine de caracteres)
+!      titre: titre de la variable (chaine de caracteres)
+!      unite : unite de la variable (chaine de caracteres)
+!      px : variable a sortir (real 0, 2, ou 3d)
+!      dim : dimension de px : 0, 2, ou 3 dimensions
+!
+!=================================================================
+!
+!      This is a modified version that accepts spectrally varying input
+!      RW (2010)
+!
+!=================================================================
+ 
+! Addition by RW (2010) to allow OLR to be saved in .nc format
+      use radinc_h, only : L_NSPECTI
+
+      implicit none
+
+! Commons
+#include "dimensions.h"
+#include "dimphys.h"
+#include "paramet.h"
+#include "control.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "description.h"
+#include "netcdf.inc"
+#include "temps.h"
+#include "surfdat.h"
+
+! Arguments on input:
+      integer ngrid
+      character (len=*) :: nom,titre,unite
+      integer dim
+      real px(ngrid,L_NSPECTI)
+
+! Local variables:
+
+!      real dx3(iip1,jjp1,llm) ! to store a 3D data set
+!      real dx2(iip1,jjp1)     ! to store a 2D (surface) data set
+!      real dx0
+
+      real date
+
+      REAL phis(ip1jmp1)
+
+      integer irythme
+      integer ierr
+      integer iq
+      integer i,j,l,zmax , ig0
+
+      integer zitau
+      character firstnom*20
+      SAVE firstnom
+      SAVE zitau
+      SAVE date
+      data firstnom /'1234567890'/
+      data zitau /0/
+
+! Ajouts
+      integer, save :: ntime=0
+      integer :: idim,varid
+      integer :: nid
+      character (len =50):: fichnom
+      integer, dimension(4) :: id
+      integer, dimension(4) :: edges,corner
+
+! added by RDW for OLR output
+       real dx3(iip1,jjp1,L_NSPECTI) ! to store the data set
+
+
+!***************************************************************
+!Sortie des variables au rythme voulu
+
+      irythme = int(ecritphy) ! sortie au rythme de ecritphy
+!     irythme = iconser  ! sortie au rythme des variables de controle
+!     irythme = iphysiq  ! sortie a tous les pas physique
+!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
+!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
+
+!***************************************************************
+
+! The following test is here to enforce that writediagfi is not used with the
+! 1D version of the GCM
+      if (ngrid.eq.1) return
+      
+c     nom=trim((nom))
+c     unite=trim((unite))
+c     titre=trim((titre))
+
+! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
+! ------------------------------------------------------------------------
+! (Au tout premier appel de la subroutine durant le run.)
+
+      fichnom="diagspecIR.nc"
+
+      if (firstnom.eq.'1234567890') then ! .true. for the very first call
+      !  to this subroutine; now set 'firstnom'
+         firstnom = nom
+         ! just to be sure, check that firstnom is large enough to hold nom
+         if (len_trim(firstnom).lt.len_trim(nom)) then
+           write(*,*) "writediagfi: Error !!!"
+           write(*,*) "   firstnom string not long enough!!"
+           write(*,*) "   increase its size to at least ",len_trim(nom)
+           stop
+         endif
+
+         ! Create the NetCDF file
+         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
+         ! Define the 'Time' dimension
+         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
+         ! Define the 'Time' variable
+#ifdef NC_DOUBLE
+         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
+#else
+         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
+#endif
+         ! Add a long_name attribute
+         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
+     .          4,"Time")
+         ! Add a units attribute
+         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
+     .          "days since 0000-00-0 00:00:00")
+         ! Switch out of NetCDF Define mode
+         ierr = NF_ENDDEF(nid)
+
+         ! write "header" of file (longitudes, latitudes, geopotential, ...)
+         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
+!         call iniwrite(nid,day_ini,phis)
+         call iniwrite_specIR(nid,day_ini,phis)
+
+         zitau = -1 ! initialize zitau
+      else
+         ! Open the NetCDF file
+         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+      endif ! if (firstnom.eq.'1234567890')
+
+! Increment time index 'zitau' if it is the "firstcall" (at given time level)
+! to writediagfi
+!------------------------------------------------------------------------
+      if (nom.eq.firstnom) then
+          zitau = zitau + iphysiq
+      end if
+
+!--------------------------------------------------------
+! Write the variables to output file if it's time to do so
+!--------------------------------------------------------
+
+      if ( MOD(zitau+1,irythme) .eq.0.) then
+
+! Compute/write/extend 'Time' coordinate (date given in days)
+! (done every "first call" (at given time level) to writediagfi)
+! Note: date is incremented as 1 step ahead of physics time
+!       (like the 'histoire' outputs)
+!--------------------------------------------------------
+
+        if (nom.eq.firstnom) then
+
+        ! We have identified a "first call" (at given date)
+           ntime=ntime+1 ! increment # of stored time steps
+           ! compute corresponding date (in days and fractions thereof)
+           date= float (zitau +1)/float (day_step)
+           ! Get NetCDF ID of 'Time' variable
+           ierr= NF_INQ_VARID(nid,"Time",varid)
+
+           !print*,'in writediagfi_spec.F, time=',varid
+
+           ! Write (append) the new date to the 'Time' array
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
+#endif
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in writediagspec_nc"
+              write(*,*) "***** with time"
+              write(*,*) 'ierr=', ierr   
+c             call abort
+           endif
+
+           write(6,*)'WRITEDIAGSPEC: date= ', date
+        end if ! of if (nom.eq.firstnom)
+
+
+ 
+!Case of a 3D variable
+!---------------------
+        if (dim.eq.3) then
+
+!         recast (copy) variable from physics grid to dynamics grid
+           DO l=1,L_NSPECTI
+             DO i=1,iip1
+                dx3(i,1,l)=px(1,l)
+                dx3(i,jjp1,l)=px(ngrid,l)
+             ENDDO
+             DO j=2,jjm
+                ig0= 1+(j-2)*iim
+                DO i=1,iim
+                   dx3(i,j,l)=px(ig0+i,l)
+                ENDDO
+                dx3(iip1,j,l)=dx3(1,j,l)
+             ENDDO
+           ENDDO
+
+! name of the variable
+           ierr= NF_INQ_VARID(nid,nom,varid)
+           if (ierr /= NF_NOERR) then
+! corresponding dimensions
+              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+              ierr= NF_INQ_DIMID(nid,"bandsIR",id(3))
+              ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+! Create the variable if it doesn't exist yet
+
+              write (*,*) "=========================="
+              write (*,*) "DIAGSPEC: creating variable ",nom
+              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+           endif
+
+           corner(1)=1
+           corner(2)=1
+           corner(3)=1
+           corner(4)=ntime
+
+           edges(1)=iip1
+           edges(2)=jjp1
+           edges(3)=L_NSPECTI
+           edges(4)=1
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
+#endif
+
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR problem in writediagspec"
+              write(*,*) "***** with ",nom
+              write(*,*) 'ierr=', ierr
+             call abort
+           endif 
+
+        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
+
+      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
+
+      ierr= NF_CLOSE(nid)
+
+      end
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 305)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 305)
@@ -0,0 +1,279 @@
+      subroutine writediagspecVI(ngrid,nom,titre,unite,dim,px)
+
+!  Ecriture de variables diagnostiques au choix dans la physique 
+!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
+!  3d (ex : temperature), 2d (ex : temperature de surface), ou
+!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
+!  solaire)
+!  Dans la version 2000, la periode d'ecriture est celle de 
+!  "ecritphy " regle dans le fichier de controle de run :  run.def
+!
+!    writediagfi peut etre appele de n'importe quelle subroutine
+!    de la physique, plusieurs fois. L'initialisation et la creation du
+!    fichier se fait au tout premier appel.
+!
+! WARNING : les variables dynamique (u,v,t,q,ps)
+!  sauvees par writediagfi avec une
+! date donnee sont legerement differentes que dans le fichier histoire car 
+! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
+! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
+! avant l'ecriture dans diagfi (cf. physiq.F)
+!  
+!
+!  parametres (input) :
+!  ----------
+!      ngrid : nombres de point ou est calcule la physique
+!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
+!                 (= nlon ou klon dans la physique terrestre)
+!      
+!      unit : unite logique du fichier de sortie (toujours la meme)
+!      nom  : nom de la variable a sortir (chaine de caracteres)
+!      titre: titre de la variable (chaine de caracteres)
+!      unite : unite de la variable (chaine de caracteres)
+!      px : variable a sortir (real 0, 2, ou 3d)
+!      dim : dimension de px : 0, 2, ou 3 dimensions
+!
+!=================================================================
+!
+!      This is a modified version that accepts spectrally varying input
+!      RW (2010)
+!
+!=================================================================
+ 
+! Addition by RW (2010) to allow OSR to be saved in .nc format
+      use radinc_h, only : L_NSPECTV
+
+      implicit none
+
+! Commons
+#include "dimensions.h"
+#include "dimphys.h"
+#include "paramet.h"
+#include "control.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "description.h"
+#include "netcdf.inc"
+#include "temps.h"
+#include "surfdat.h"
+
+! Arguments on input:
+      integer ngrid
+      character (len=*) :: nom,titre,unite
+      integer dim
+      real px(ngrid,L_NSPECTV)
+
+! Local variables:
+
+!      real dx3(iip1,jjp1,llm) ! to store a 3D data set
+!      real dx2(iip1,jjp1)     ! to store a 2D (surface) data set
+!      real dx0
+
+      real date
+
+      REAL phis(ip1jmp1)
+
+      integer irythme
+      integer ierr
+      integer iq
+      integer i,j,l,zmax , ig0
+
+      integer zitau
+      character firstnom*20
+      SAVE firstnom
+      SAVE zitau
+      SAVE date
+      data firstnom /'1234567890'/
+      data zitau /0/
+
+! Ajouts
+      integer, save :: ntime=0
+      integer :: idim,varid
+      integer :: nid
+      character (len =50):: fichnom
+      integer, dimension(4) :: id
+      integer, dimension(4) :: edges,corner
+
+! added by RDW for OSR output
+       real dx3(iip1,jjp1,L_NSPECTV) ! to store the data set
+
+
+!***************************************************************
+!Sortie des variables au rythme voulu
+
+      irythme = int(ecritphy) ! sortie au rythme de ecritphy
+!     irythme = iconser  ! sortie au rythme des variables de controle
+!     irythme = iphysiq  ! sortie a tous les pas physique
+!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
+!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
+
+!***************************************************************
+
+! The following test is here to enforce that writediagfi is not used with the
+! 1D version of the GCM
+      if (ngrid.eq.1) return
+      
+c     nom=trim((nom))
+c     unite=trim((unite))
+c     titre=trim((titre))
+
+! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
+! ------------------------------------------------------------------------
+! (Au tout premier appel de la subroutine durant le run.)
+
+      fichnom="diagspecVI.nc"
+
+      if (firstnom.eq.'1234567890') then ! .true. for the very first call
+      !  to this subroutine; now set 'firstnom'
+         firstnom = nom
+         ! just to be sure, check that firstnom is large enough to hold nom
+         if (len_trim(firstnom).lt.len_trim(nom)) then
+           write(*,*) "writediagfi: Error !!!"
+           write(*,*) "   firstnom string not long enough!!"
+           write(*,*) "   increase its size to at least ",len_trim(nom)
+           stop
+         endif
+
+         ! Create the NetCDF file
+         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
+         ! Define the 'Time' dimension
+         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
+         ! Define the 'Time' variable
+#ifdef NC_DOUBLE
+         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
+#else
+         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
+#endif
+         ! Add a long_name attribute
+         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
+     .          4,"Time")
+         ! Add a units attribute
+         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
+     .          "days since 0000-00-0 00:00:00")
+         ! Switch out of NetCDF Define mode
+         ierr = NF_ENDDEF(nid)
+
+         ! write "header" of file (longitudes, latitudes, geopotential, ...)
+         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
+!         call iniwrite(nid,day_ini,phis)
+         call iniwrite_specVI(nid,day_ini,phis)
+
+         zitau = -1 ! initialize zitau
+      else
+         ! Open the NetCDF file
+         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+      endif ! if (firstnom.eq.'1234567890')
+
+! Increment time index 'zitau' if it is the "firstcall" (at given time level)
+! to writediagfi
+!------------------------------------------------------------------------
+      if (nom.eq.firstnom) then
+          zitau = zitau + iphysiq
+      end if
+
+!--------------------------------------------------------
+! Write the variables to output file if it's time to do so
+!--------------------------------------------------------
+
+      if ( MOD(zitau+1,irythme) .eq.0.) then
+
+! Compute/write/extend 'Time' coordinate (date given in days)
+! (done every "first call" (at given time level) to writediagfi)
+! Note: date is incremented as 1 step ahead of physics time
+!       (like the 'histoire' outputs)
+!--------------------------------------------------------
+
+        if (nom.eq.firstnom) then
+
+        ! We have identified a "first call" (at given date)
+           ntime=ntime+1 ! increment # of stored time steps
+           ! compute corresponding date (in days and fractions thereof)
+           date= float (zitau +1)/float (day_step)
+           ! Get NetCDF ID of 'Time' variable
+           ierr= NF_INQ_VARID(nid,"Time",varid)
+
+           !print*,'in writediagfi_spec.F, time=',varid
+
+           ! Write (append) the new date to the 'Time' array
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
+#endif
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in writediagspec_nc"
+              write(*,*) "***** with time"
+              write(*,*) 'ierr=', ierr   
+c             call abort
+           endif
+
+           write(6,*)'WRITEDIAGSPEC: date= ', date
+        end if ! of if (nom.eq.firstnom)
+
+
+ 
+!Case of a 3D variable
+!---------------------
+        if (dim.eq.3) then
+
+!         recast (copy) variable from physics grid to dynamics grid
+           DO l=1,L_NSPECTV
+             DO i=1,iip1
+                dx3(i,1,l)=px(1,l)
+                dx3(i,jjp1,l)=px(ngrid,l)
+             ENDDO
+             DO j=2,jjm
+                ig0= 1+(j-2)*iim
+                DO i=1,iim
+                   dx3(i,j,l)=px(ig0+i,l)
+                ENDDO
+                dx3(iip1,j,l)=dx3(1,j,l)
+             ENDDO
+           ENDDO
+
+! name of the variable
+           ierr= NF_INQ_VARID(nid,nom,varid)
+           if (ierr /= NF_NOERR) then
+! corresponding dimensions
+              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+              ierr= NF_INQ_DIMID(nid,"bandsVI",id(3))
+              ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+! Create the variable if it doesn't exist yet
+
+              write (*,*) "=========================="
+              write (*,*) "DIAGSPEC: creating variable ",nom
+              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+           endif
+
+           corner(1)=1
+           corner(2)=1
+           corner(3)=1
+           corner(4)=ntime
+
+           edges(1)=iip1
+           edges(2)=jjp1
+           edges(3)=L_NSPECTV
+           edges(4)=1
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
+#endif
+
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR problem in writediagspec"
+              write(*,*) "***** with ",nom
+              write(*,*) 'ierr=', ierr
+             call abort
+           endif 
+
+        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
+
+      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
+
+      ierr= NF_CLOSE(nid)
+
+      end
