Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 587)
+++ /trunk/LMDZ.GENERIC/README	(revision 588)
@@ -599,2 +599,13 @@
 - Temperature grid for Planck calculations can now be refined through the parameter NTfac in radinc_h. 
   Default is NTfac = 1.0D-1, i.e. Delta T = 0.1 K
+
+== 16/03/2012 == JL
+- Removed cpp3D and nonideal stuff.
+
+== 19/03/2012 == EM
+Some cleanup and bug fixing:
+- "cloudfrac" was not well written to restartfi (wrong size).
+- missing save attribute for "reffrad" in physiq.F90.
+- cleanup recomputation of surface pressure in newstart and change loop order
+  in interp_horiz (which "fixes" an odd behaviour which fills some arrays with
+  zeros, but only when using some versions of ifort!)
Index: /trunk/LMDZ.GENERIC/libf/dyn3d/interp_horiz.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dyn3d/interp_horiz.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/dyn3d/interp_horiz.F	(revision 588)
@@ -19,15 +19,15 @@
 c  """""""""
         
-       INTEGER imo, jmo ! dimensions ancienne grille (input)
-       INTEGER imn,jmn  ! dimensions nouvelle grille (input)
+       INTEGER,INTENT(IN) :: imo, jmo ! dimensions ancienne grille (input)
+       INTEGER,INTENT(IN) :: imn,jmn  ! dimensions nouvelle grille (input)
 
-       REAL rlonuo(imo+1)     !  Latitude et
-       REAL rlatvo(jmo)       !  longitude des
-       REAL rlonun(imn+1)     !  bord des 
-       REAL rlatvn(jmn)     !  cases "scalaires" (input)
+       REAL,INTENT(IN) :: rlonuo(imo+1)     !  Latitude et
+       REAL,INTENT(IN) :: rlatvo(jmo)       !  longitude des
+       REAL,INTENT(IN) :: rlonun(imn+1)     !  bord des 
+       REAL,INTENT(IN) :: rlatvn(jmn)     !  cases "scalaires" (input)
 
-       INTEGER lm ! dimension verticale (input)
-       REAL varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
-       REAL varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
+       INTEGER,INTENT(IN) :: lm ! dimension verticale (input)
+       REAL,INTENT(IN) :: varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
+       REAL,INTENT(OUT) :: varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
 
 c Autres variables
@@ -39,5 +39,5 @@
        INTEGER ii,jj,l
 
-       REAL airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
+       REAL,SAVE :: airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
        REAL airentotn	! aire totale pole nord dans la nouvelle grille
        REAL airentots	! aire totale pole sud dans la nouvelle grille
@@ -48,26 +48,18 @@
 c + des pouiemes (cas ou une maille est a cheval sur 2 ou 4 mailles)
 c  Il y a un test dans iniinterp_h pour s'assurer que ktotal < kmax
-       INTEGER kmax, k, ktotal
+       INTEGER kmax, k
+       integer,save :: ktotal
        parameter (kmax = 360*179 + 200000)
 c      parameter (kmax = 360*179 + 40000)
 
-       INTEGER iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
-       REAL intersec(kmax)
-       REAL R
+       INTEGER,SAVE :: iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
+       REAL,SAVE :: intersec(kmax)
+       REAL r
        REAL totn, tots
-       integer prev_sumdim
-       save prev_sumdim
-       data prev_sumdim /0/
-       
+       integer,save :: prev_sumdim=0
 
-       logical firsttest, aire_ok
-       save firsttest
-       data firsttest /.true./
-       data aire_ok /.true./
+       logical,save :: firsttest=.true. , aire_ok=.true.
 
-       integer imoS,jmoS,imnS,jmnS
-       save imoS,jmoS,imnS,jmnS
-       save ktotal,iik,jjk,jk,ik,intersec,airen
-       REAL pi
+       integer,save :: imoS,jmoS,imnS,jmnS
 
 c Test dimensions imnmx2 jmnmx2
@@ -115,11 +107,6 @@
       end if
 
-      do l=1,lm
-       do jj =1 , jmn+1
-        do ii=1, imn+1
-          varn(ii,jj,l) =0.
-        end do
-       end do
-      end do 
+! initialize varn() to zero
+      varn(1:imn+1,1:jmn+1,1:lm)=0.
        
 c Interpolation horizontale
@@ -128,7 +115,9 @@
 c de l'ancienne et la  nouvelle grille
 c
-     
-      do k=1,ktotal
-        do l=1,lm
+! Ehouarn 2012: for some strange reason, with ifort v12.x,
+!               when the order of the loop below is changed
+!               values of varn(:,:,l=2...) are then sometimes remain zero!   
+      do l=1,lm
+        do k=1,ktotal
          varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l) 
      &   + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k)
@@ -161,5 +150,5 @@
          end do 
 
-      ENDDO
+      ENDDO ! of DO l=1, lm
            
 
@@ -167,5 +156,4 @@
 c  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST 
       if (firsttest) then
-      pi=2.*asin(1.)
       firsttest = .false.
       write (*,*) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
@@ -192,5 +180,5 @@
       end do
       if (aire_ok) write(*,*) 'INTERP. HORIZ. : AIRES OK'
-      endif
+      endif ! of if (firsttest)
 
 c FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST
@@ -198,4 +186,3 @@
 
 
-        return
         end
Index: /trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F	(revision 588)
@@ -46,8 +46,4 @@
       INTEGER   imold,jmold,lmold,nsoilold,nqold
 
-c et autres:
-c----------
-      INTEGER lnblnk
-      EXTERNAL lnblnk
 
 c Variables pour les lectures des fichiers "ini" 
@@ -109,5 +105,5 @@
 c     REAL year_day,periheli,aphelie,peri_day
 c     REAL obliquit,z0,emin_turb,lmixmin
-c     REAL emissiv,emisice(2),albedice(2),tauvis
+c     REAL emissiv,emisice(2),albedice(2)
 c     REAL iceradius(2) , dtemisice(2)
 
@@ -994,10 +990,17 @@
       co2icetotal = 0.
       if (igcm_co2_ice.ne.0) then
-      DO j=1,jjp1
+        ! recast surface CO2 ice on new grid
+        call interp_horiz(qsurfold(1,1,igcm_co2_ice),
+     &                  qsurfs(1,1,igcm_co2_ice),
+     &                  imold,jmold,iim,jjm,1,
+     &                  rlonuold,rlatvold,rlonu,rlatv)
+       DO j=1,jjp1
          DO i=1,iim
            !co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
            co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
          END DO
-      END DO
+       END DO
+      else
+        write(*,*) "Warning: No co2_ice surface tracer"
       endif
 
@@ -1009,5 +1012,7 @@
       write(*,*)'Ancienne grille: masse de la glace CO2:',co2icetotalold
       write(*,*)'Nouvelle grille: masse de la glace CO2:',co2icetotal
+      if (co2icetotalold.ne.0.) then
       write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
+      endif
       write(*,*)
 
Index: /trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F	(revision 588)
@@ -307,8 +307,4 @@
       pa=tab_cntrl(17)  ! reference pressure at which coord is purely pressure
       !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
-      if (preff.eq.0) then
-        preff=610
-        pa=20
-      endif
       write(*,*) "Newstart: preff=",preff," pa=",pa
       yes=' '
@@ -1671,4 +1667,6 @@
             cloudfrac(ig,l)=0.5
          ENDDO
+! Ehouarn, march 2012: also add some initialisation for hice
+         hice(ig)=0.0
       ENDDO
       
@@ -1682,49 +1680,27 @@
 !      ENDDO
 
-
 c=======================================================================
 c   Correct pressure on the new grid (menu 0)
 c=======================================================================
+
 
       if ((choix_1.eq.0).and.(.not.flagps0)) then
         r = 1000.*8.31/mugaz
 
-        phi0=0.0
         do j=1,jjp1
           do i=1,iip1
-             phi0 = phi0+phis(i,j)*aire(i,j)
-          end do
-        end do
-        phi0=phi0/airetot
-
-        do j=1,jjp1
-          do i=1,iip1
-             ps(i,j) = ps(i,j) * 
-!     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
-     .            exp((phi0-phis(i,j)) /
+             ps(i,j) = ps(i,j) *
+     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
      .                                  (t(i,j,1) * r))
           end do
         end do
-  
-!     we must renormalise pressures AGAIN, because exp(-phi) is nonlinear
-        ptot=0.0
+
+c periodicite de ps en longitude
         do j=1,jjp1
-           do i=1,iip1
-              ptot=ptot+ps(i,j)*aire(i,j)
-           enddo
-        enddo
-        do j=1,jjp1
-           do i=1,iip1
-              ps(i,j)=ps(i,j)*ptotn/ptot
-           enddo
-        enddo
-        
-!     periodicity of surface ps in longitude
-        do j=1,jjp1
-           ps(1,j) = ps(iip1,j)
+          ps(1,j) = ps(iip1,j)
         end do
-        
       end if
 
+         
 c=======================================================================
 c=======================================================================
@@ -1733,4 +1709,5 @@
 c    Initialisation de la physique / ecriture de newstartfi :
 c=======================================================================
+
 
       CALL inifilr 
@@ -1764,5 +1741,4 @@
       endif
 
-
 C Calcul intermediaire
 c
@@ -1798,8 +1774,6 @@
       CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
      *                phi,w, pbaru,pbarv,day_ini+time )
-c     CALL caldyn
-c    $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
-c    $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
-
+
+          
       CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
       CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps) 
Index: /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 588)
@@ -80,5 +80,8 @@
          R     = 8.314511E+0 *1000.E+0/mugaz
          rcp   = R/cpp
-      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
+!      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
+         ! Ehouarn: tolerate a small mismatch between computed/stored values
+       elseif((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
+              (abs(1.-mugaz/mugaz_c).gt.1.e-6)) then
          if(check_cpp_match)then
             print*,'Values do not match!'
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 588)
@@ -522,6 +522,6 @@
                temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
 
-               tauaero(2*k+2,iaer)=max(temp*pweight,0.0)
-               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.0)
+               tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
+               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
 !
             end do
Index: /trunk/LMDZ.GENERIC/libf/phystd/datareadnc.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/datareadnc.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/datareadnc.F	(revision 588)
@@ -43,4 +43,6 @@
 
       use datafile_mod, only: datadir
+! to use  'getin'
+      USE ioipsl_getincom 
       implicit none
 
@@ -77,6 +79,5 @@
       REAL        zthe(imdp1*jmdp1)
 
-      INTEGER   lnblnk, ierr
-      EXTERNAL    lnblnk
+      INTEGER   ierr
 
       INTEGER   unit,nvarid
@@ -110,5 +111,7 @@
 c    Lecture NetCDF des donnees latitude et longitude
 c-----------------------------------------------------------------------
-
+!    First get the corret value of dtadir, if not already done:
+      ! default 'datadir' is set in "datadir_mod"
+      call getin("datadir",datadir)
       write(*,*) 'Choice of surface data is:'
       filestring='ls '//trim(datadir)//'/*.nc'
Index: /trunk/LMDZ.GENERIC/libf/phystd/def_var.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/def_var.F90	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/def_var.F90	(revision 588)
@@ -27,6 +27,6 @@
 ierr=NF_REDEF(nid)
 
-print*,'in def_var.F90, dimids='
-print*,dimids
+!print*,'in def_var.F90, dimids='
+!print*,dimids
 
 ! 2. Define the variable
Index: /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F	(revision 588)
@@ -48,5 +48,5 @@
       REAL day_ini
       INTEGER nsoil,nq
-      integer ierr,idim1,idim2,idim3,idim4,idim5,nvarid
+      integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
 
       REAL phystep,time
@@ -82,6 +82,4 @@
       real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
       integer ig
-      INTEGER lnblnk
-      EXTERNAL lnblnk
 
 ! flag which identifies if we are using old tracer names (qsurf01,...)
@@ -144,5 +142,4 @@
       endif
 c
-!      ierr = NF_DEF_DIM (nid,"nlayer+1",llm+1,idim4)
       ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
       if (ierr.ne.NF_NOERR) then
@@ -157,4 +154,11 @@
         WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
         write(6,*) NF_STRERROR(ierr)
+      endif
+
+      ierr = NF_DEF_DIM (nid,"nlayer",llm,idim6)
+      if (ierr.ne.NF_NOERR) then
+        WRITE(6,*)'physdem1: Problem defining nlayer dimension'
+        write(6,*) NF_STRERROR(ierr)
+	call abort
       endif
 
@@ -448,24 +452,4 @@
 c Write the physical fields
 
-! CO2 Ice Cover
-!
-!      ierr = NF_REDEF (nid)
-!#ifdef NC_DOUBLE
-!      ierr = NF_DEF_VAR (nid, "co2ice", NF_DOUBLE, 1, idim2,nvarid)
-!#else
-!      ierr = NF_DEF_VAR (nid, "co2ice", NF_FLOAT, 1, idim2,nvarid)
-!#endif
-!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 13,
-!     .                        "CO2 ice cover")
-!      ierr = NF_ENDDEF(nid)
-!#ifdef NC_DOUBLE
-!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,co2ice)
-!#else
-!      ierr = NF_PUT_VAR_REAL (nid,nvarid,co2ice)
-!#endif
-
-
-
-
 
 ! Soil Thermal inertia
@@ -574,7 +558,9 @@
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 1, idim2,nvarid)
+      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 2,
+     & (/idim2,idim6/),nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 2,
+     & (/idim2,idim6/),nvarid)
 #endif
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 588)
@@ -209,11 +209,10 @@
       real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
 
-      real reffrad(ngridmx,nlayermx,naerkind)                 ! aerosol effective radius (m)
+      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
 
 !     Tendencies due to various processes:
       real dqsurf(ngridmx,nqmx)
-      real zdtlw(ngridmx,nlayermx)                            ! (K/s)
-      real zdtsw(ngridmx,nlayermx)                            ! (K/s)
-      save zdtlw, zdtsw
+      real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
+      real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
 
       real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
@@ -1724,8 +1723,8 @@
             if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
                ! to record global radiative balance
-               open(92,file="rad_bal.out",form='formatted',access='append')
+               open(92,file="rad_bal.out",form='formatted',position='append')
                write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
                close(92)
-               open(93,file="tem_bal.out",form='formatted',access='append')
+               open(93,file="tem_bal.out",form='formatted',position='append')
                write(93,*) zday,Ts1,Ts2,Ts3,TsS
                close(93)
@@ -1739,5 +1738,5 @@
 !     ------------------------------------------------------------------
       if(testradtimes)then
-         open(38,file="tau_phys.out",form='formatted',access='append')
+         open(38,file="tau_phys.out",form='formatted',position='append')
          ig=1
          do l=1,nlayer
@@ -1813,5 +1812,5 @@
             if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
                ! to record global water balance
-               open(98,file="h2o_bal.out",form='formatted',access='append')
+               open(98,file="h2o_bal.out",form='formatted',position='append')
                write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
                close(98)
Index: /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 587)
+++ /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 588)
@@ -176,4 +176,5 @@
       mugaz = tab_cntrl(tab0+8)
       rcp = tab_cntrl(tab0+9)
+      cpp=(8.314511/(mugaz/1000.0))/rcp
       daysec = tab_cntrl(tab0+10)
       dtphys = tab_cntrl(tab0+11)
@@ -502,5 +503,5 @@
           write(*,*) 
           write(*,*) ' mugaz (new value):',mugaz
-          r=8.314/(mugaz/1000.0)
+          r=8.314511/(mugaz/1000.0)
           write(*,*) ' R (new value):',r
 
@@ -512,4 +513,5 @@
           write(*,*) 
           write(*,*) ' rcp (new value):',rcp
+          r=8.314511/(mugaz/1000.0)
           cpp=r/rcp
           write(*,*) ' cpp (new value):',cpp
