Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 588)
+++ /trunk/LMDZ.GENERIC/README	(revision 589)
@@ -610,2 +610,7 @@
   in interp_horiz (which "fixes" an odd behaviour which fills some arrays with
   zeros, but only when using some versions of ifort!)
+
+== 19/03/2012 == AS
+- Cleaned rcm1d.F and made it truly generic by asking for planetary constants without any default values.
+- Settings are now needed in a rcm1d.def file. Unbeknown to the user, we create a minimal run.def file, read parameters, then remove this dummy run.def.
+- Introduced a keyword force_cpp if the user wants to give values for cpp and mugaz in def files.
Index: /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 588)
+++ /trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 589)
@@ -9,10 +9,9 @@
 !     unless option 'check_cpp_match' is set to false in 
 !     callphys.def.
-!     NOTE: for now, in 1D we do as before. Jeremy, if you're 
-!     re-writing rcm1d you may want to alter this.
 !
 !     Authors
 !     ------- 
 !     Robin Wordsworth (2009)
+!     A. Spiga: make the routine OK with latest changes in rcm1d
 !
 !==================================================================
@@ -26,5 +25,5 @@
 #include "callkeys.h"
 
-      real cpp_c
+      real cpp_c   
       real mugaz_c
 
@@ -74,22 +73,27 @@
       print*,'Predefined Mg in physics is ',mugaz,'amu'
 
-      if(ngridmx.eq.1)then
-         print*,'Automatically setting cpp & mugaz to calculated values in calc_cpp_mugaz' 
-         cpp   = cpp_c
-         mugaz = mugaz_c
-         R     = 8.314511E+0 *1000.E+0/mugaz
-         rcp   = R/cpp
-!      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.  &
+      if (check_cpp_match) then
+         print*,'REQUEST TO CHECK cpp_match :'
+         if((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!'
-            print*,'Either adjust cpp / mugaz via newstart to calculated values,'
-            print*,'or set check_cpp_match to .false. in callphys.def.'
+            ! Ehouarn: tolerate a small mismatch between computed/stored values
+            print*,'--> Values do not match!'
+            print*,'    Either adjust cpp / mugaz via newstart to calculated values,'
+            print*,'    or set check_cpp_match to .false. in callphys.def.'
             stop
+         else
+            print*,'--> OK. Settings match composition.'
          endif
       endif
 
+      if (.not.force_cpp) then
+          print*,'*** Setting cpp & mugaz to computations in calc_cpp_mugaz.'
+          mugaz = mugaz_c
+          cpp = cpp_c
+      else
+          print*,'*** Setting cpp & mugaz to predefined values.'
+      endif
+
+
       return
     end subroutine calc_cpp_mugaz
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 588)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 589)
@@ -21,4 +21,5 @@
      &   , newtonian                                                    &
      &   , check_cpp_match                                              &
+     &   , force_cpp                                                    &
      &   , tau_relax                                                    &
      &   , testradtimes                                                 &
@@ -63,4 +64,5 @@
       logical newtonian
       logical check_cpp_match
+      logical force_cpp
       logical testradtimes
       logical rayleigh
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 588)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 589)
@@ -182,5 +182,4 @@
          write(*,*) " check_cpp_match = ",check_cpp_match
 
-
          write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value
@@ -308,8 +307,8 @@
            call abort
          endif
-         if (noradsurf.and.calldifv) then
-           print*,'noradsurf not compatible with a boundary layer!'
-           call abort
-         endif
+         !if (noradsurf.and.calldifv) then
+         !  print*,'noradsurf not compatible with a boundary layer!'
+         !  call abort
+         !endif
 
          write(*,*)"Use Newtonian cooling for radiative transfer?"
@@ -508,5 +507,31 @@
          endif
 
-         mugaz=8.314*1000./pr
+         write(*,*) "Does user want to force cpp and mugaz?"
+         force_cpp=.false. ! default value
+         call getin("force_cpp",force_cpp)
+         write(*,*) " force_cpp = ",force_cpp
+
+         if (force_cpp) then
+           mugaz = -99999.
+           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
+           call getin("mugaz",mugaz)
+           IF (mugaz.eq.-99999.) THEN
+               PRINT *, "mugaz must be set if force_cpp = T"
+               STOP
+           ELSE
+               write(*,*) "mugaz=",mugaz
+           ENDIF
+           cpp = -99999.
+           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
+           call getin("cpp",cpp)
+           IF (cpp.eq.-99999.) THEN
+               PRINT *, "cpp must be set if force_cpp = T"
+               STOP
+           ELSE
+               write(*,*) "cpp=",cpp
+           ENDIF
+         else
+           mugaz=8.314*1000./pr
+         endif
          call su_gases
          call calc_cpp_mugaz
Index: /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 588)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 589)
@@ -25,4 +25,6 @@
 !     F. Montmessin (water ice added)
 !     R. Wordsworth
+!     B. Charnay
+!     A. Spiga
 !
 !==================================================================
@@ -116,26 +118,10 @@
       saveprofile=.false.
 
-c ------------------------------------------------------
-c  Default values for constants (corresponding to Mars)
-c ------------------------------------------------------
+c ----------------------------------------
+c  Default values  (corresponding to Mars)
+c ----------------------------------------
 
       pi=2.E+0*asin(1.E+0)
 
-c     Constante de la Planete Mars
-c     ----------------------------
-      rad=3397200.               ! rayon de mars (m)  ~3397200 m
-      daysec=88775.              ! duree du sol (s)  ~88775 s
-      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
-      g=3.72                     ! gravite (m.s-2) ~3.72  
-      mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
-      rcp=.256793                ! = r/cp  ~0.256793
-      r= 8.314511E+0 *1000.E+0/mugaz
-      cpp= r/rcp
-      year_day = 669           ! duree de l'annee (sols) ~668.6
-      periastr = 206.66          ! min. dist. star-planet (Mkm) ~206.66
-      apoastr = 249.22           ! max. dist. star-planet (Mkm) ~249.22
-      peri_day =  485.           ! date du periastron (sols depuis printemps)
-      obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
- 
 c     Parametres Couche limite et Turbulence 
 c     --------------------------------------
@@ -159,17 +145,13 @@
 
 c ------------------------------------------------------
-c  Load parameters from "run.def" 
+c  Load parameters from "run.def" and "gases.def" 
 c ------------------------------------------------------
 
-! also 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',
+! check if 'rcm1d.def' file is around 
+      open(90,file='rcm1d.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(*,*) 'Cannot find required file "rcm1d.def"'
+        write(*,*) 'which should contain some input parameters'
         write(*,*) ' ... might as well stop here ...'
         stop
@@ -178,4 +160,22 @@
       endif
 
+! now, run.def is needed anyway. so we create a dummy temporary one
+! for ioipsl to work. if a run.def is already here, stop the
+! program and ask the user to do a bit of cleaning
+      open(90,file='run.def',status='old',form='formatted',
+     &     iostat=ierr)
+      if (ierr.eq.0) then
+        close(90)
+        write(*,*) 'There is already a run.def file.'
+        write(*,*) 'This is not compatible with 1D runs.'
+        write(*,*) 'Please remove the file and restart the run.'
+        write(*,*) 'Runtime parameters are supposed to be in rcm1d.def'
+        stop
+      else
+        call system('touch run.def')
+        call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
+        call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
+      endif
+
 ! check if we are going to run with or without tracers
       write(*,*) "Run with or without tracer transport ?"
@@ -183,4 +183,8 @@
       call getin("tracer",tracer)
       write(*,*) " tracer = ",tracer
+
+! OK. now that run.def has been read once -- any variable is in memory.
+! so we can dump the dummy run.def
+      call system("rm -rf run.def")
 
 ! while we're at it, check if there is a 'traceur.def' file
@@ -256,4 +260,155 @@
       endif ! of if (tracer)
 
+!!! We have to check that check_cpp_match is F for 1D computations
+!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
+      check_cpp_match = .false.
+      call getin("check_cpp_match",check_cpp_match)
+      if (check_cpp_match) then
+          print*,"In 1D modeling, check_cpp_match is supposed to be F"
+          print*,"Please correct callphys.def"
+          stop
+      endif
+
+!!! GEOGRAPHICAL INITIALIZATIONS
+     !!! AREA. useless in 1D
+      area(1)=1.E+0
+      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
+     !!! GEOPOTENTIAL. useless in 1D because control bu surface pressure
+      phisfi(1)=0.E+0
+     !!! LATITUDE. can be set. 
+      lati(1)=0 ! default value for lati(1)
+      PRINT *,'latitude (in degrees) ?'
+      call getin("latitude",lati(1))
+      write(*,*) " latitude = ",lati(1)
+      lati(1)=lati(1)*pi/180.E+0
+     !!! LONGITUDE. useless in 1D.
+      long(1)=0.E+0
+      long(1)=long(1)*pi/180.E+0
+
+!!! CALL INIFIS
+!!! - read callphys.def
+!!! - calculate sine and cosine of longitude and latitude
+!!! - calculate mugaz and cp 
+!!! NB: some operations are being done dummily in inifis in 1D
+!!! - physical constants: nevermind, things are done allright below
+!!! - physical frequency: nevermind, in inifis this is a simple print
+      CALL inifis(1,llm,day0,daysec,dtphys,
+     .            lati,long,area,rad,g,r,cpp)
+!!! We check everything is OK.
+      PRINT *,"CHECK"
+      PRINT *,"--> mugaz = ",mugaz
+      PRINT *,"--> cpp = ",cpp
+      r = 8.314511E+0 * 1000.E+0 / mugaz
+      rcp = r / cpp
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!! PLANETARY CONSTANTS !!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      g = -99999.
+      PRINT *,'GRAVITY in m s-2 ?'
+      call getin("g",g)
+      IF (g.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED g IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> g = ",g
+      ENDIF
+
+      !rad = -99999.
+      !PRINT *,'PLANETARY RADIUS in m ?'
+      !call getin("rad",rad)
+      !IF (rad.eq.-99999.) THEN
+      !    PRINT *,"STOP. I NEED rad IN RUN.DEF."
+      !    STOP
+      !ELSE
+      !    PRINT *,"--> rad = ",rad
+      !ENDIF
+
+      daysec = -99999.
+      PRINT *,'LENGTH OF A DAY in s ?'
+      call getin("daysec",daysec)
+      IF (daysec.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED daysec IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> daysec = ",daysec
+      ENDIF
+      omeg=4.*asin(1.)/(daysec)
+      PRINT *,"OK. FROM THIS I WORKED OUT:"
+      PRINT *,"--> omeg = ",omeg
+
+      year_day = -99999.
+      PRINT *,'LENGTH OF A YEAR in days ?'
+      call getin("year_day",year_day)
+      IF (year_day.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED year_day IN RUN.DEF."
+          STOP
+      ELSE 
+          PRINT *,"--> year_day = ",year_day
+      ENDIF
+
+      periastr = -99999.
+      PRINT *,'MIN DIST STAR-PLANET in AU ?'
+      call getin("periastr",periastr)
+      IF (periastr.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED periastr IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> periastr = ",periastr
+      ENDIF
+      periastr=periastr*149.598 ! AU to Mkm
+
+      apoastr = -99999.
+      PRINT *,'MIN DIST STAR-PLANET in AU ?'
+      call getin("apoastr",apoastr)
+      IF (apoastr.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED apoastr IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> apoastr = ",apoastr
+      ENDIF
+      apoastr=apoastr*149.598 ! AU to Mkm
+
+      peri_day = -99999.
+      PRINT *,'DATE OF PERIASTRON in days ?'
+      call getin("peri_day",peri_day)
+      IF (peri_day.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED peri_day IN RUN.DEF."
+          STOP
+      ELSE IF (peri_day.gt.year_day) THEN
+          PRINT *,"STOP. peri_day.gt.year_day"
+          STOP
+      ELSE  
+          PRINT *,"--> peri_day = ", peri_day
+      ENDIF 
+
+      obliquit = -99999. 
+      PRINT *,'OBLIQUITY in deg ?'
+      call getin("obliquit",obliquit)
+      IF (obliquit.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED obliquit IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> obliquit = ",obliquit
+      ENDIF 
+
+      psurf = -99999.
+      PRINT *,'SURFACE PRESSURE in Pa ?'
+      call getin("psurf",psurf)
+      IF (psurf.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED psurf IN RUN.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> psurf = ",psurf
+      ENDIF
+      !! we need reference pressures
+      pa=psurf/30.
+      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!! END PLANETARY CONSTANTS !!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 c  Date et heure locale du debut du run
@@ -298,5 +453,9 @@
       ndt=ndt*day_step     
       dtphys=daysec/day_step  
-
+      write(*,*)"-------------------------------------"
+      write(*,*)"-------------------------------------"
+      write(*,*)"--> Physical timestep is ",dtphys
+      write(*,*)"-------------------------------------"
+      write(*,*)"-------------------------------------"
 
 c output spectrum?
@@ -306,30 +465,5 @@
       write(*,*)" specOLR = ",specOLR
 
-
-c Pression de surface sur la planete
-c ------------------------------------
-c
-      psurf=100000. ! default value for psurf
-      PRINT *,'Surface pressure (Pa) ?'
-      call getin("psurf",psurf)
-      write(*,*) " psurf = ",psurf
-c Pression de reference
-      pa=psurf/30.
-      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
-
-c Gravity of planet
-c -----------------
-      g=10.0 ! default value for g
-      PRINT *,'Gravity ?'
-      call getin("g",g)
-      write(*,*) " g = ",g
-
-      mugaz=0.
-      cpp=0.
-
-      call su_gases
-      call calc_cpp_mugaz
-      
-
+!!! AS12: what shall we do with this?
 c Proprietes des poussiere aerosol
 c --------------------------------
@@ -339,49 +473,7 @@
       write(*,*) " tauvis = ",tauvis
  
-c  latitude/longitude
-c  ------------------
-      lati(1)=0 ! default value for lati(1)
-      PRINT *,'latitude (in degrees) ?'
-      call getin("latitude",lati(1))
-      write(*,*) " latitude = ",lati(1)
-      lati(1)=lati(1)*pi/180.E+0
-      long(1)=0.E+0
-      long(1)=long(1)*pi/180.E+0
-
-c  periastron/apoastron
-c  --------------------
-      periastr=227.94
-      apoastr=227.94 ! default = mean martian values
-      PRINT *,'periastron (in AU) ?'
-      call getin("periastr",periastr)
-      write(*,*) "periastron = ",periastr
-      periastr=periastr*149.598 ! AU to Mkm
-      PRINT *,'apoastron (in AU) ?'
-      call getin("apoastr",apoastr)
-      write(*,*) "apoastron = ",apoastr
-      apoastr=apoastr*149.598 ! AU to Mkm
-
-c  obliquity
-c  ---------
-      obliquit=0.0! default=zero in 1d
-      PRINT *,'obliquity (in AU) ?'
-      call getin("obliquit",obliquit)
-      write(*,*) "obliquity = ",obliquit
-
-c  Initialisation albedo / inertie du sol
-c  --------------------------------------
-      albedodat(1)=0.2 ! default value for albedodat
-      PRINT *,'Albedo of bare ground ?'
-      call getin("albedo",albedodat(1))
-      write(*,*) " albedo = ",albedodat(1)
-
-      inertiedat(1,1)=400 ! default value for inertiedat
-      PRINT *,'Soil thermal inertia (SI) ?'
-      call getin("inertia",inertiedat(1,1))
-      write(*,*) " inertia = ",inertiedat(1,1)
 c
 c  pour le schema d'ondes de gravite
 c  ---------------------------------
-c
       zmea(1)=0.E+0
       zstd(1)=0.E+0
@@ -403,22 +495,4 @@
       ENDDO
 
-c   Initialisation speciales "physiq"
-c   ---------------------------------
-c   la surface de chaque maille est inutile en 1D --->
-      area(1)=1.E+0
-      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
-
-c   le geopotentiel au sol est inutile en 1D car tout est controle
-c   par la pression de surface --->
-      phisfi(1)=0.E+0
-
-c  "inifis" reproduit un certain nombre d'initialisations deja faites
-c  + lecture des clefs de callphys.def
-c  + calcul de la frequence d'appel au rayonnement
-c  + calcul des sinus et cosinus des longitude latitude
-
-!Mars possible issue with dtphys in input and include!!!
-      CALL inifis(1,llm,day0,daysec,dtphys,
-     .            lati,long,area,rad,g,r,cpp)
 c   Initialisation pour prendre en compte les vents en 1-D
 c   ------------------------------------------------------
@@ -488,5 +562,5 @@
       write(*,*) " autozlevs = ", autozlevs
 
-      pceil=100.0 ! Pascals
+      pceil = psurf / 1000.0 ! Pascals
       PRINT *,'Ceiling pressure (Pa) ?'
       call getin("pceil",pceil)
@@ -571,4 +645,15 @@
       
 
+c  Initialisation albedo / inertie du sol
+c  --------------------------------------
+      albedodat(1)=0.2 ! default value for albedodat
+      PRINT *,'Albedo of bare ground ?'
+      call getin("albedo",albedodat(1))
+      write(*,*) " albedo = ",albedodat(1)
+
+      inertiedat(1,1)=400 ! default value for inertiedat
+      PRINT *,'Soil thermal inertia (SI) ?'
+      call getin("inertia",inertiedat(1,1))
+      write(*,*) " inertia = ",inertiedat(1,1)
 
 ! Initialize soil properties and temperature
