Index: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 3515)
+++ trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 3562)
@@ -38,4 +38,9 @@
       use phys_state_var_mod, only: phys_state_var_init
       use physiq_mod, only: physiq
+      use restart1D_mod, only: writerestart1D
+      use iostart, only: length
+      use netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE,
+     &                  nf90_strerror,NF90_INQ_VARID, NF90_GET_VAR,
+     &                  NF90_CLOSE
       implicit none
 
@@ -78,4 +83,8 @@
       INTEGER ilayer,ilevel,isoil,idt,iq
       LOGICAl firstcall,lastcall
+      LOGICAL restart
+      INTEGER nid_restart1D, nid_restartfi
+      INTEGER controleid, uid, vid, tempid, tsurfid, did, tid
+      INTEGER,ALLOCATABLE :: qid(:), qsurfid(:)
 c
       INTEGER day0          ! date initial (sol ; =0 a Ls=0)
@@ -122,4 +131,5 @@
 
       logical oldcompare, earthhack,saveprofile
+      real controle1D(length)
 
 !     added by RW for zlay computation
@@ -186,4 +196,21 @@
         call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
       endif
+
+! Check restart
+      call getin("restart",restart)
+      if (restart) then
+        ierr = NF90_OPEN('start1D.nc', NF90_NOWRITE, nid_restart1D)
+        if (ierr.NE.NF90_NOERR) then
+          write(*,*)'readstart1D: problem opening 1D start file'
+          write(*,*)trim(nf90_strerror(ierr))
+          call abort_physic("readstart1D","Cannot open file",1)
+        endif
+        ierr = NF90_OPEN('startfi.nc', NF90_NOWRITE, nid_restartfi)
+        if (ierr.NE.NF90_NOERR) then
+          write(*,*)'readstart1D: problem opening startfi file'
+          write(*,*)trim(nf90_strerror(ierr))
+          call abort_physic("readstart1D","Cannot open file",1)
+        endif
+      endif ! restart
 
       ! read nq from traceur.def
@@ -474,13 +501,29 @@
       ENDIF 
 
-      psurf = -99999.
-      PRINT *,'SURFACE PRESSURE in Pa ?'
-      call getin("psurf",psurf)
-      IF (psurf.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> psurf = ",psurf
-      ENDIF
+      if (restart) then
+        ierr=NF90_INQ_VARID(nid_restart1D,'controle',controleid)
+        if (ierr==NF90_NOERR) then
+          ierr=NF90_GET_VAR(nid_restart1D,controleid,controle1D)
+          if (ierr/=NF90_NOERR) then
+            PRINT*, 'read restart 1D: Failed loading psurf'
+            CALL abort_physic("readstart1D","Failed to read variable",1)
+          endif
+        endif
+        psurf = controle1D(1)
+        PRINT *,"In restart1D.nc, day0=",controle1D(2),
+     &          " and time=",controle1D(3)
+        PRINT*,"(These values are NOT used, ",
+     &         "those given in rcm1d.def are used)"
+      else ! restart
+        psurf = -99999.
+        PRINT *,'SURFACE PRESSURE in Pa ?'
+        call getin("psurf",psurf)
+        IF (psurf.eq.-99999.) THEN
+            PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
+            STOP
+        ELSE
+            PRINT *,"--> psurf = ",psurf
+        ENDIF
+      endif ! restart
       !! we need reference pressures
       pa=psurf/30.
@@ -494,16 +537,38 @@
 c  ------------------------------------
 c    Date (en sols depuis le solstice de printemps) du debut du run
-      day0 = 0 ! default value for day0
-      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
-      call getin("day0",day0)
-      day=float(day0)
-      write(*,*) " day0 = ",day0
+      !if (restart) then
+      !  ! day    
+      !  ierr=NF90_INQ_VARID(nid_restart1D,'day',did)
+      !  if (ierr==NF90_NOERR) then
+      !    ierr=NF90_GET_VAR(nid_restart1D,did,day)
+      !    if (ierr/=NF90_NOERR) then
+      !      PRINT*, 'read restart 1D: Failed loading day'
+      !      CALL abort_physic("readstart1D","Failed to read variable",1)
+      !    endif
+      !  endif
+
+      !  ! time
+      !  ierr=NF90_INQ_VARID(nid_restart1D,'time',tid)
+      !  if (ierr==NF90_NOERR) then
+      !    ierr=NF90_GET_VAR(nid_restart1D,tid,time)
+      !    if (ierr/=NF90_NOERR) then
+      !      PRINT*, 'read restart 1D: Failed loading time'
+      !      CALL abort_physic("readstart1D","Failed to read variable",1)
+      !    endif
+      !  endif
+
+      !else
+        day0 = 0 ! default value for day0
+        write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
+        call getin("day0",day0)
+        day=float(day0)
+        write(*,*) " day0 = ",day0
 c  Heure de demarrage
-      time=0 ! default value for time
-      write(*,*)'Initial local time (in hours, between 0 and 24)?'
-      call getin("time",time)
-      write(*,*)" time = ",time
-      time=time/24.E+0 ! convert time (hours) to fraction of sol
-
+        time=0 ! default value for time
+        write(*,*)'Initial local time (in hours, between 0 and 24)?'
+        call getin("time",time)
+        write(*,*)" time = ",time
+        time=time/24.E+0 ! convert time (hours) to fraction of sol
+      !endif
 
 c  Discretisation (Definition de la grille et des pas de temps)
@@ -629,93 +694,124 @@
 
          do iq=1,nq
+
+            if (restart) then
+              if (iq.eq.1) then
+                allocate(qid(nq))
+                allocate(qsurfid(nq))
+              endif
+              ! read q
+              ierr=NF90_INQ_VARID(nid_restart1D,noms(iq),qid(iq))
+              if (ierr==NF90_NOERR) then
+                ierr=NF90_GET_VAR(nid_restart1D,qid(iq),q(:,iq),
+     &           count=(/1,nlayer,1,1/))
+                if (ierr/=NF90_NOERR) then
+                  PRINT*, 'read restart 1D: Failed loading ', noms(iq)
+                  write(*,*)trim(nf90_strerror(ierr))
+                  CALL abort_physic("readstart1D",
+     &                              "Failed to read variable",1)
+                endif
+              endif
+              ! read qsurf
+              ierr=NF90_INQ_VARID(nid_restartfi,noms(iq),qsurfid(iq))
+              if (ierr==NF90_NOERR) then
+                ierr=NF90_GET_VAR(nid_restartfi,qsurfid(iq),qsurf(iq))
+                if (ierr/=NF90_NOERR) then
+                  PRINT*, 'read restartfi: Failed loading ', noms(iq)
+                  write(*,*)trim(nf90_strerror(ierr))
+                  CALL abort_physic("readstartfi",
+     &                              "Failed to read variable",1)
+                endif
+              endif
+            else ! restart
          
-            txt=""
-            write(txt,"(a)") tname(iq)
-            write(*,*)"  tracer:",trim(txt)
+              txt=""
+              write(txt,"(a)") tname(iq)
+              write(*,*)"  tracer:",trim(txt)
               
-            ! CO2
-            if (txt.eq."co2_ice") then
-               q(:,iq)=0.   ! kg/kg of atmosphere
-               qsurf(iq)=0. ! kg/m2 at the surface               
-               ! Look for a "profile_co2_ice" input file
-               open(91,file='profile_co2_ice',status='old',
-     &         form='formatted',iostat=ierr)
-               if (ierr.eq.0) then
-                  read(91,*) qsurf(iq)
-                  do ilayer=1,nlayer
-                     read(91,*) q(ilayer,iq)
-                  enddo
-               else
-                  write(*,*) "No profile_co2_ice file!"
-               endif
-               close(91)
-            endif ! of if (txt.eq."co2")
+              ! CO2
+              if (txt.eq."co2_ice") then
+                 q(:,iq)=0.   ! kg/kg of atmosphere
+                 qsurf(iq)=0. ! kg/m2 at the surface               
+                 ! Look for a "profile_co2_ice" input file
+                 open(91,file='profile_co2_ice',status='old',
+     &           form='formatted',iostat=ierr)
+                 if (ierr.eq.0) then
+                    read(91,*) qsurf(iq)
+                    do ilayer=1,nlayer
+                       read(91,*) q(ilayer,iq)
+                    enddo
+                 else
+                    write(*,*) "No profile_co2_ice file!"
+                 endif
+                 close(91)
+              endif ! of if (txt.eq."co2")
           
-            ! WATER VAPOUR
-            if (txt.eq."h2o_vap") then
-               q(:,iq)=0.   ! kg/kg of atmosphere
-               qsurf(iq)=0. ! kg/m2 at the surface
-               ! Look for a "profile_h2o_vap" input file   
-               open(91,file='profile_h2o_vap',status='old',
-     &         form='formatted',iostat=ierr)
-               if (ierr.eq.0) then
-                  read(91,*) qsurf(iq)
-                  do ilayer=1,nlayer
-                     read(91,*) q(ilayer,iq)
-                  enddo
-               else
-                  write(*,*) "No profile_h2o_vap file!"
-               endif
-               close(91)
-            endif ! of if (txt.eq."h2o_vap")
+              ! WATER VAPOUR
+              if (txt.eq."h2o_vap") then
+                 q(:,iq)=0.   ! kg/kg of atmosphere
+                 qsurf(iq)=0. ! kg/m2 at the surface
+                 ! Look for a "profile_h2o_vap" input file   
+                 open(91,file='profile_h2o_vap',status='old',
+     &           form='formatted',iostat=ierr)
+                 if (ierr.eq.0) then
+                    read(91,*) qsurf(iq)
+                    do ilayer=1,nlayer
+                       read(91,*) q(ilayer,iq)
+                    enddo
+                 else
+                    write(*,*) "No profile_h2o_vap file!"
+                 endif
+                 close(91)
+              endif ! of if (txt.eq."h2o_vap")
             
-            ! WATER ICE
-            if (txt.eq."h2o_ice") then
-               q(:,iq)=0.   ! kg/kg of atmosphere
-               qsurf(iq)=0. ! kg/m2 at the surface
-               ! Look for a "profile_h2o_ice" input file
-               open(91,file='profile_h2o_ice',status='old',
-     &         form='formatted',iostat=ierr)
-               if (ierr.eq.0) then
-                  read(91,*) qsurf(iq)
-                  do ilayer=1,nlayer
-                     read(91,*) q(ilayer,iq)
-                  enddo
-               else
-                  write(*,*) "No profile_h2o_ice file!"
-               endif
-               close(91)
-            endif ! of if (txt.eq."h2o_ice")
-
-            !_vap
-            if((txt .ne. 'h2o_vap') .and.
-     &                     (index(txt,'_vap'   ) .ne. 0))   then
-                  q(:,iq)=0. !kg/kg of atmosphere
-                  qsurf(iq) = 0. !kg/kg of atmosphere
-                  ! Look for a "profile_...._vap" input file
-                  tracer_profile_file_name=""
-                  tracer_profile_file_name='profile_'//txt
-                  open(91,file=tracer_profile_file_name,status='old',
-     &            form="formatted",iostat=ierr)
-                  if (ierr .eq. 0) then 
-                        read(91,*)qsurf(iq)
-                        do ilayer=1,nlayer 
-                              read(91,*)q(ilayer,iq)
-                        enddo 
-                  else 
-                        write(*,*) "No initial profile "
-                        write(*,*) " for this tracer :"
-                        write(*,*) txt
-                  endif
-                  close(91)
-            endif ! (txt .eq. "_vap") 
-            !_ice 
-            if((txt.ne."h2o_ice") .and. 
-     &                      (index(txt,'_ice'   ) /= 0)) then
-                  q(:,iq)=0. !kg/kg of atmosphere
-                  qsurf(iq) = 0. !kg/kg of atmosphere
-            endif ! we only initialize the solid at 0
-         enddo ! of do iq=1,nq
-         
+              ! WATER ICE
+              if (txt.eq."h2o_ice") then
+                 q(:,iq)=0.   ! kg/kg of atmosphere
+                 qsurf(iq)=0. ! kg/m2 at the surface
+                 ! Look for a "profile_h2o_ice" input file
+                 open(91,file='profile_h2o_ice',status='old',
+     &           form='formatted',iostat=ierr)
+                 if (ierr.eq.0) then
+                    read(91,*) qsurf(iq)
+                    do ilayer=1,nlayer
+                       read(91,*) q(ilayer,iq)
+                    enddo
+                 else
+                    write(*,*) "No profile_h2o_ice file!"
+                 endif
+                 close(91)
+              endif ! of if (txt.eq."h2o_ice")
+
+              !_vap
+              if((txt .ne. 'h2o_vap') .and.
+     &                       (index(txt,'_vap'   ) .ne. 0))   then
+                    q(:,iq)=0. !kg/kg of atmosphere
+                    qsurf(iq) = 0. !kg/kg of atmosphere
+                    ! Look for a "profile_...._vap" input file
+                    tracer_profile_file_name=""
+                    tracer_profile_file_name='profile_'//txt
+                    open(91,file=tracer_profile_file_name,status='old',
+     &              form="formatted",iostat=ierr)
+                    if (ierr .eq. 0) then 
+                          read(91,*)qsurf(iq)
+                          do ilayer=1,nlayer 
+                                read(91,*)q(ilayer,iq)
+                          enddo 
+                    else 
+                          write(*,*) "No initial profile "
+                          write(*,*) " for this tracer :"
+                          write(*,*) txt
+                    endif
+                    close(91)
+              endif ! (txt .eq. "_vap") 
+              !_ice 
+              if((txt.ne."h2o_ice") .and. 
+     &                        (index(txt,'_ice'   ) /= 0)) then
+                    q(:,iq)=0. !kg/kg of atmosphere
+                    qsurf(iq) = 0. !kg/kg of atmosphere
+              endif ! we only initialize the solid at 0
+            endif ! restart
+           enddo ! of do iq=1,nq
+        
       endif ! of tracer
 
@@ -739,9 +835,36 @@
       w(1:nlayer)=0
 
-c     Initialisation des vents  au premier pas de temps
-      DO ilayer=1,nlayer
-         u(ilayer)=gru
-         v(ilayer)=grv
-      ENDDO
+      if (restart) then
+            ierr=NF90_INQ_VARID(nid_restart1D,'u',uid)
+            if (ierr==NF90_NOERR) then
+              ierr=NF90_GET_VAR(nid_restart1D,uid,u,
+     &                          count=(/1,nlayer,1,1/))
+              if (ierr/=NF90_NOERR) then
+                PRINT*, 'read restart 1D: Failed loading u'
+                write(*,*)trim(nf90_strerror(ierr))
+                CALL abort_physic("readstart1D",
+     &                            "Failed to read variable",1)
+              endif
+            endif
+            ierr=NF90_INQ_VARID(nid_restart1D,'v',vid)
+            if (ierr==NF90_NOERR) then
+              ierr=NF90_GET_VAR(nid_restart1D,vid,v,
+     &                          count=(/1,nlayer,1,1/))
+              if (ierr/=NF90_NOERR) then
+                PRINT*, 'read restart 1D: Failed loading v'
+                write(*,*)trim(nf90_strerror(ierr))
+                CALL abort_physic("readstart1D",
+     &                            "Failed to read variable",1)
+              endif
+            endif
+      else ! restart
+
+c       Initialisation des vents  au premier pas de temps
+        DO ilayer=1,nlayer
+           u(ilayer)=gru
+           v(ilayer)=grv
+        ENDDO
+
+      endif ! restart
 
 c     energie cinetique turbulente
@@ -859,22 +982,49 @@
 !      endif
 
+
 c  profil de temperature au premier appel
 c  --------------------------------------
       pks=psurf**rcp
 
+      if (restart) then
+            ! read temp
+            ierr=NF90_IN Q_VARID(nid_restart1D,'temp',tempid)
+            if (ierr==NF90_NOERR) then
+              ierr=NF90_GET_VAR(nid_restart1D,tempid,temp,
+     &                          count=(/1,nlayer,1,1/))
+              if (ierr/=NF90_NOERR) then
+                PRINT*, 'read restart 1D: Failed loading temp'
+                CALL abort_physic("readstart1D",
+     &                            "Failed to read variable",1)
+              endif
+            endif
+            ! read tsurf
+            ierr=NF90_INQ_VARID(nid_restartfi,'tsurf',tsurfid)
+            if (ierr==NF90_NOERR) then
+              ierr=NF90_GET_VAR(nid_restartfi,tsurfid,tsurf)
+              if (ierr/=NF90_NOERR) then
+                PRINT*, 'read restartfi: Failed loading tsurf'
+                CALL abort_physic("readstartfi",
+     &                            "Failed to read variable",1)
+              endif
+            endif
+      else ! restart
+
 c altitude en km dans profile: on divise zlay par 1000
-      tmp1(0)=0.E+0
-      DO ilayer=1,nlayer
-        tmp1(ilayer)=zlay(ilayer)/1000.E+0
-      ENDDO
-      call profile(nlayer+1,tmp1,tmp2)
-
-      tsurf(1)=tmp2(0)
-      DO ilayer=1,nlayer
-        temp(ilayer)=tmp2(ilayer)
-      ENDDO
-      print*,"check"
-      PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
-      PRINT*,"INPUT TEMPERATURE PROFILE",temp
+        tmp1(0)=0.E+0
+        DO ilayer=1,nlayer
+          tmp1(ilayer)=zlay(ilayer)/1000.E+0
+        ENDDO
+        call profile(nlayer+1,tmp1,tmp2)
+
+        tsurf(1)=tmp2(0)
+        DO ilayer=1,nlayer
+          temp(ilayer)=tmp2(ilayer)
+        ENDDO
+        print*,"check"
+        PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
+        PRINT*,"INPUT TEMPERATURE PROFILE",temp
+
+      endif ! restart
 
 c  Initialisation albedo / inertie du sol
@@ -937,5 +1087,6 @@
 ! -----------------
 #ifndef MESOSCALE
-      if(tracer.and.photochem) then
+      if(tracer.and.photochem.and. .not.restart) then
+           print*, "Calling inichim_1D"
            call initracer(1,nq)
            allocate(nametmp(nq))
@@ -947,4 +1098,18 @@
 #endif
 
+      if (restart) then
+        ierr = NF90_CLOSE(nid_restart1D)
+        if (ierr/=NF90_NOERR) then
+            PRINT*, 'read restart1D: Failed closing restart1D.nc'
+            CALL abort_physic("readstart1D",
+     &                        "Failed closing file",1)
+          endif
+        ierr = NF90_CLOSE(nid_restartfi)
+        if (ierr/=NF90_NOERR) then
+            PRINT*, 'read restartfi: Failed closing restartfi.nc'
+            CALL abort_physic("readstartfi",
+     &                        "Failed closing file",1)
+          endif
+      endif
 
 c  Write a "startfi" file
@@ -953,13 +1118,14 @@
 c  It is needed to transfert physics variables to "physiq"...
 
-      call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
-     &              dtphys,real(day0),time,cell_area,
-     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
-      call physdem1("startfi.nc",nsoilmx,1,llm,nq,
-     &                dtphys,time,
-     &                tsurf,tsoil,emis,albedo,q2,qsurf,
-     &                cloudfrac,totcloudfrac,hice, 
-     &                rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice)
-
+      if (.not. restart) then
+        call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
+     &                dtphys,real(day0),time,cell_area,
+     &                albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+        call physdem1("startfi.nc",nsoilmx,1,llm,nq,
+     &                  dtphys,time,
+     &                  tsurf,tsoil,emis,albedo,q2,qsurf,
+     &                  cloudfrac,totcloudfrac,hice, 
+     &                  rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice)
+      endif
 c=======================================================================
 c  BOUCLE TEMPORELLE DU MODELE 1D 
@@ -1098,4 +1264,9 @@
       endif
 
+c     Produce the restart1D.nc file (MM)
+      if (lastcall) then
+        call writerestart1D('restart1D.nc',nlayer,nsoil,day,time,psurf,
+     &                                     temp,tsoil,u,v,nq,q)    
+      endif
 
       ENDDO   ! fin de la boucle temporelle
