Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3906)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3912)
@@ -7,5 +7,6 @@
 SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref,               &
                            nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
-                           play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
+                           play,plev,latitude,longitude,cell_area,ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,      &
+                           ctrl_h2oice,relaxtime_h2oice,qref_h2oice)
 
 use ioipsl_getincom,          only: getin ! To use 'getin'
@@ -65,24 +66,24 @@
 integer, intent(inout) :: nq
 
-real, dimension(:,:,:), allocatable, intent(out) :: q     ! tracer mixing ratio (e.g. kg/kg)
-real,                                intent(out) :: time  ! time (0<time<1; time=0.5 at noon)
-real,                                intent(out) :: psurf ! Surface pressure
-real, dimension(nlayer),             intent(out) :: u, v  ! zonal, meridional wind
-real, dimension(nlayer),             intent(out) :: temp  ! temperature at the middle of the layers
+real, dimension(:,:,:), allocatable, intent(out) :: q          ! tracer mixing ratio (e.g. kg/kg)
+real,                                intent(out) :: time       ! time (0<time<1; time=0.5 at noon)
+real,                                intent(out) :: psurf      ! Surface pressure
+real, dimension(nlayer),             intent(out) :: u, v       ! zonal, meridional wind
+real, dimension(nlayer),             intent(out) :: temp       ! temperature at the middle of the layers
 integer,                             intent(out) :: ndt
 real,                                intent(out) :: ptif, pks
-real,                                intent(out) :: dttestphys                   ! testphys1d timestep
+real,                                intent(out) :: dttestphys ! testphys1d timestep
 real, dimension(:),     allocatable, intent(out) :: zqsat
-real, dimension(:,:,:), allocatable, intent(out) :: dq, dqdyn                    ! Physical and dynamical tandencies
-integer,                             intent(out) :: day0                         ! initial (sol ; =0 at Ls=0) and final date
-real,                                intent(out) :: day                          ! date during the run
-real,                                intent(out) :: gru, grv                     ! prescribed "geostrophic" background wind
-real, dimension(nlayer),             intent(out) :: w                            ! "Dummy wind" in 1D
-real, dimension(nlayer),             intent(out) :: play                         ! Pressure at the middle of the layers (Pa)
-real, dimension(nlayer + 1),         intent(out) :: plev                         ! intermediate pressure levels (pa)
+real, dimension(:,:,:), allocatable, intent(out) :: dq, dqdyn  ! Physical and dynamical tandencies
+integer,                             intent(out) :: day0       ! initial (sol ; =0 at Ls=0) and final date
+real,                                intent(out) :: day        ! date during the run
+real,                                intent(out) :: gru, grv   ! prescribed "geostrophic" background wind
+real, dimension(nlayer),             intent(out) :: w          ! "Dummy wind" in 1D
+real, dimension(nlayer),             intent(out) :: play       ! Pressure at the middle of the layers (Pa)
+real, dimension(nlayer + 1),         intent(out) :: plev       ! intermediate pressure levels (pa)
 real, dimension(1),                  intent(out) :: latitude, longitude, cell_area
-logical,                             intent(out) :: prescribed_h2ovap            ! Prescribe atmospheric water profile
-real,                                intent(out) :: h2ovap_relax_time            ! Relaxtion time towards the atmospheric water profile
-real, dimension(nlayer),             intent(out) :: q_prescribed_h2o_vap         ! User-defined atmospheric water profile
+logical,                             intent(out) :: ctrl_h2ovap, ctrl_h2oice           ! Flags to prescribe atmospheric water vapour/ice profiles
+real,                                intent(out) :: relaxtime_h2ovap, relaxtime_h2oice ! Relaxation time towards the atmospheric water vapour/ice profiles
+real, dimension(nlayer),             intent(out) :: qref_h2ovap, qref_h2oice           ! User-defined atmospheric water vapour/ice profiles
 
 !=======================================================================
@@ -781,38 +782,66 @@
 if (.not. therestartfi .and. paleoclimate) perennial_co2ice = 0.
 
-! Prescription of atmospheric water profile
-! -----------------------------------------
-prescribed_h2ovap = .false. ! Default: free atmospheric water profile
-h2ovap_relax_time = -1. ! Default: no time relaxation
-q_prescribed_h2o_vap = 0. ! Default
+! Prescription of atmospheric water vapour/ice profiles
+! ----------------------------------------------
+ctrl_h2ovap = .false. ; ctrl_h2oice = .false. ! Default: free atmospheric water vapour/ice profiles
+relaxtime_h2ovap = -1. ; relaxtime_h2oice = -1. ! Default: no time relaxation
+qref_h2ovap = 0. ; qref_h2oice = 0. ! Default
 if (water) then
-    write(*,*)'Force atmospheric water profile (mixing ratio in kg/kg)?'
-    call getin('prescribed_h2ovap',prescribed_h2ovap)
-    write(*,*) 'prescribed_h2ovap = ', prescribed_h2ovap
-    if (prescribed_h2ovap) then
-        write(*,*) 'Atmospheric water profile is prescribed'
-        open(10,file = 'profile_prescribed_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
+    write(*,*)'Force atmospheric water vapour profile (mixing ratio in kg/kg)?'
+    call getin('ctrl_h2ovap',ctrl_h2ovap)
+    write(*,*) 'ctrl_h2ovap = ', ctrl_h2ovap
+    if (ctrl_h2ovap) then
+        write(*,*) 'Atmospheric water vapour profile is prescribed'
+        open(10,file = 'profile_ref_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
         if (ierr == 0) then
             do ilayer = 1,nlayer
-                read(10,*,iostat=ierr) q_prescribed_h2o_vap(ilayer)
-                if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "profile_prescribed_h2o_vap"!'
+                read(10,*,iostat=ierr) qref_h2ovap(ilayer)
+                if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "profile_ref_h2o_vap"!'
             enddo
         else
-            error stop 'File "profile_prescribed_h2o_vap" was not found!'
+            error stop 'File "profile_ref_h2o_vap" was not found!'
         endif 
         close(10)
 
-        write(*,*) 'Relax atmospheric water profile with time constant (s)?'
+        write(*,*) 'Relax atmospheric water vapour profile with time constant (s)?'
         
-        call getin('h2ovap_relax_time',h2ovap_relax_time)
-        write(*,*) 'h2ovap_relax_time = ', h2ovap_relax_time
-        if (h2ovap_relax_time < 0.) then
-            write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).'
+        call getin('relaxtime_h2ovap',relaxtime_h2ovap)
+        write(*,*) 'relaxtime_h2ovap = ', relaxtime_h2ovap
+        if (relaxtime_h2ovap < 0.) then
+            write(*,*) 'Atmospheric water vapour profile is not relaxed (forcing).'
         else
-            write(*,*) 'Atmospheric water profile is relaxed towards the profile in "profile_prescribed_h2o_vap"'
+            write(*,*) 'Atmospheric water vapour profile is relaxed towards the profile in "profile_ref_h2o_vap"'
         endif
     else
-        write(*,*) 'Free atmospheric water vapor profile'
-        write(*,*) 'Total water is conserved in the column'
+        write(*,*) 'Free atmospheric water vapour profile'
+    endif
+
+    write(*,*)'Force atmospheric water ice profile (mixing ratio in kg/kg)?'
+    call getin('ctrl_h2oice',ctrl_h2oice)
+    write(*,*) 'ctrl_h2oice = ', ctrl_h2oice
+    if (ctrl_h2oice) then
+        write(*,*) 'Atmospheric water ice profile is prescribed'
+        open(10,file = 'profile_ref_h2o_ice',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
+        if (ierr == 0) then
+            do ilayer = 1,nlayer
+                read(10,*,iostat=ierr) qref_h2oice(ilayer)
+                if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "profile_ref_h2o_ice"!'
+            enddo
+        else
+            error stop 'File "profile_ref_h2o_ice" was not found!'
+        endif 
+        close(10)
+
+        write(*,*) 'Relax atmospheric water ice profile with time constant (s)?'
+        
+        call getin('relaxtime_h2oice',relaxtime_h2oice)
+        write(*,*) 'relaxtime_h2oice = ', relaxtime_h2oice
+        if (relaxtime_h2oice < 0.) then
+            write(*,*) 'Atmospheric water ice profile is not relaxed (forcing).'
+        else
+            write(*,*) 'Atmospheric water ice profile is relaxed towards the profile in "profile_ref_h2o_ice"'
+        endif
+    else
+        write(*,*) 'Free atmospheric water ice profile'
     endif
 endif
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3906)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3912)
@@ -6,5 +6,5 @@
 use phyredem,            only: physdem0, physdem1
 use watersat_mod,        only: watersat
-use tracer_mod,          only: igcm_h2o_vap, noms
+use tracer_mod,          only: igcm_h2o_vap, igcm_h2o_ice, noms
 use comcstfi_h,          only: pi, g, rcp, cpp
 use time_phylmdz_mod,    only: daysec
@@ -83,8 +83,8 @@
 logical :: startfiles_1D, therestart1D, therestartfi, there
 
-! Prescription of atmospheric water profile
-logical                         :: prescribed_h2ovap
-real                            :: h2ovap_relax_time
-real, dimension(nlayer)         :: q_prescribed_h2o_vap
+! Prescription of atmospheric water/ice profiles
+logical                         :: ctrl_h2ovap, ctrl_h2oice
+real                            :: relaxtime_h2ovap, relaxtime_h2oice
+real, dimension(nlayer)         :: qref_h2ovap, qref_h2oice
 real, dimension(:), allocatable :: zqsat
 !=======================================================================
@@ -142,5 +142,6 @@
 call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
                      time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
-                     play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
+                     play,plev,latitude,longitude,cell_area,                                        &
+                     ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,ctrl_h2oice,relaxtime_h2oice,qref_h2oice)
 
 ! Write a "startfi" file
@@ -208,14 +209,23 @@
     ! Water tracer increment: specific for 1D
     ! ---------------------------------------
-    ! The physics computes the tendency on q(1,:,igcm_h2o_vap), here we can mimic an effect of the "3D dynamics"
-    ! either by forcing a profile or by relaxing towards a prescribed profile
-    if (water .and. prescribed_h2ovap) then ! If atmospheric water profile is prescribed
-        if (h2ovap_relax_time < 0.) then ! Forcing
-            ! For some tests: unless it reaches saturation (maximal value)
-            !call watersat(nlayer,temp,play,zqsat)
-            !dq(1,:,igcm_h2o_vap) = (min(zqsat(:),q_prescribed_h2o_vap(:)) - q(1,:,igcm_h2o_vap))/dttestphys
-            dq(1,:,igcm_h2o_vap) = (q_prescribed_h2o_vap(:) - q(1,:,igcm_h2o_vap))/dttestphys
-        else ! Relaxation
-            dq(1,:,igcm_h2o_vap) = dq(1,:,igcm_h2o_vap) - (q(1,:,igcm_h2o_vap) - q_prescribed_h2o_vap(:))/h2ovap_relax_time
+    ! The physics computes the tendency on q for igcm_h2o_vap/igcm_h2o_ice, here we can mimic an effect
+    ! of the "3D dynamics" either by forcing a profile or by relaxing towards a prescribed profile
+    if (water) then
+        if (ctrl_h2ovap) then ! If atmospheric water vapour profile is prescribed
+            if (relaxtime_h2ovap < 0.) then ! Forcing
+                ! For some tests: unless it reaches saturation (maximal value)
+                !call watersat(nlayer,temp,play,zqsat)
+                !dq(1,:,igcm_h2o_vap) = (min(zqsat(:),qref_h2ovap(:)) - q(1,:,igcm_h2o_vap))/dttestphys
+                dq(1,:,igcm_h2o_vap) = (qref_h2ovap(:) - q(1,:,igcm_h2o_vap))/dttestphys
+            else ! Relaxation
+                dq(1,:,igcm_h2o_vap) = dq(1,:,igcm_h2o_vap) - (q(1,:,igcm_h2o_vap) - qref_h2ovap(:))/relaxtime_h2ovap
+            endif
+        endif
+        if (ctrl_h2oice) then ! If atmospheric water ice profile is prescribed
+            if (relaxtime_h2oice < 0.) then ! Forcing
+                dq(1,:,igcm_h2o_ice) = (qref_h2oice(:) - q(1,:,igcm_h2o_ice))/dttestphys
+            else ! Relaxation
+                dq(1,:,igcm_h2o_ice) = dq(1,:,igcm_h2o_ice) - (q(1,:,igcm_h2o_ice) - qref_h2oice(:))/relaxtime_h2oice
+            endif
         endif
      endif