Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 282)
+++ trunk/LMDZ.MARS/README	(revision 283)
@@ -863,2 +863,18 @@
 - Changed the definition for HFX computation in the LES in meso_inc/meso_inc_les.F. New definition yields
   very similar results to old one and follows a strict definition of what HFX should be.
+
+== 24/08/11 == TN
+Attempts to tune the water cycle by adding outliers 
+       + A few structural changes !!
+* watercap.h is now obsolete and removed -- all is in surfdat.h
+* watercaptag initialized in surfini.F (up to 7 areas defined) instead of initracer.F
+  - settings proposed by AS commented
+  - experiments by TN decommented. use with caution.
+* water ice albedo and thermal inertia in callphys.def and inifis.F
+* water ice albedo in surfini.F
+* water ice albedo computation in albedocaps.F90
+* alb_surfice is now obsolete in physiq.F, albedo_h2o_ice is used instead
+* frost_albedo_threshold defined in surfdat.h
+* water ice thermal inertia in soil.F
+TODO: * calibrate thermal inertia and ice albedo
+      * have a look at subgrid-scale ice with dryness ?
Index: trunk/LMDZ.MARS/deftank/callphys.def.outliers
===================================================================
--- trunk/LMDZ.MARS/deftank/callphys.def.outliers	(revision 283)
+++ trunk/LMDZ.MARS/deftank/callphys.def.outliers	(revision 283)
@@ -0,0 +1,127 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.true.
+
+#Diurnal cycle ?  if diurnal=False, diurnal averaged solar heating
+diurnal=.true.
+
+#Seasonal cycle ? if season=False, Ls stays constant, to value set in "start"
+season = .true.   
+
+#write some more output on the screen ?
+lwrite = .false.    
+
+#Save statistics in file "stats.nc" ?
+callstats =.true.
+
+#Save EOF profiles in file "profiles" for Climate Database?
+calleofdump = .false.
+
+## Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+## ~~~~~~~~~~~~~
+#  =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario,
+#  =4 Mars Year 24 from TES assimilation (old version of MY24; dust_tes.nc file)
+#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
+#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
+#  =26 Mars Year 26 from TES assimilation
+iaervar = 26 
+# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.2
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 3
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 55.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+# Impose polar cap surface albedos as observed by TES?
+TESicealbedo = .true.
+
+## Radiative transfer options :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# the rad.transfer is computed every "iradia" physical timestep
+iradia    = 1
+# Output of the exchange coefficient mattrix ? for diagnostic only
+callg2d  = .false.
+# Rayleigh scattering : (should be .false. for now)
+rayleigh  = .false.
+
+## Tracer (dust water, ice and/or chemical species) options (used if tracer=T):
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# DUST: Transported dust ? (if >0, use 'dustbin' dust bins)
+dustbin    = 2
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .true.
+# DUST: use mass and number mixing ratios to predict dust size ?
+doubleq   = .true.
+# DUST: use a small population of dust particules (submicron dust)?
+submicron = .false.
+# DUST: lifted by GCM surface winds ?
+lifting = .true.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by CO2 snowfall ?
+scavenging = .false.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .true.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .true.
+# WATER: Compute water cycle
+water = .true.
+# WATER: current permanent caps at both poles. True IS RECOMMENDED
+#        (with .true., North cap is a source of water and South pole
+#         is a cold trap)
+caps  = .true.
+# WATER: water ice albedo ? STILL TBD !!
+albedo_h2o_ice = 0.4
+# WATER: water ice thermal inertia (SI) ? STILL TBD !!
+# (je trouve dans les 2400 SI à 173 K)
+inert_h2o_ice = 1000.
+temptag = .false.
+# WATER: water ice frost minimal thickness (ie kg.m^-2, ie 0.92 mm) for albedo
+frost_albedo_threshold = 0.005
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+## Thermospheric options (relevant if tracer=T) :
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call thermosphere ?
+callthermos = .false.
+# WATER: included without cycle (only if water=.false.)
+thermoswater = .false.
+# call thermal conduction ? (only if callthermos=.true.)
+callconduct = .false.
+# call EUV heating ? (only if callthermos=.true.)
+calleuv=.false.
+# call molecular viscosity ? (only if callthermos=.true.)
+callmolvis = .false.
+# call molecular diffusion ? (only if callthermos=.true.)
+callmoldiff = .false.
+# call thermospheric photochemistry ? (only if callthermos=.true.)
+thermochem = .false.
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6)
+solarcondate = 1993.4
+
+       
+
+
Index: trunk/LMDZ.MARS/libf/phymars/albedocaps.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/albedocaps.F90	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/albedocaps.F90	(revision 283)
@@ -12,4 +12,5 @@
 #include"dimphys.h"
 #include"surfdat.h"
+#include"callkeys.h"
 
 ! arguments:
@@ -70,4 +71,11 @@
       psolaralb(ig,2)=albedice(icap)
     endif
+  else if (watercaptag(ig) .and. water) then
+    ! there is a water ice cap: set the surface albedo to the water ice one
+    ! to do : emissivity
+    !write(*,*) "watercaptag in albedocaps:",ig
+    emisref(ig) = 1
+    psolaralb(ig,1)=albedo_h2o_ice
+    psolaralb(ig,2)=albedo_h2o_ice
   else
     ! set emissivity of surface to be bare ground emissivity
@@ -78,5 +86,4 @@
   endif ! of if (piceco2(ig).gt.0)
 enddo ! of ig=1,ngrid
-
 end subroutine albedocaps
 
Index: trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 283)
@@ -412,4 +412,31 @@
          call getin("caps",caps)
          write(*,*) " caps = ",caps
+
+
+         write(*,*) "water ice albedo ?"
+         albedo_h2o_ice=0.45
+         call getin("albedo_h2o_ice",albedo_h2o_ice)
+         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
+         
+         write(*,*) "water ice thermal inertia ?"
+         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
+         call getin("inert_h2o_ice",inert_h2o_ice)
+         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
+         
+         write(*,*) "frost thickness threshold for albedo ?"
+         frost_albedo_threshold=0.005 ! 5.4 µm (i.e 0.005 kg/m²)
+         call getin("frost_albedo_threshold",
+     &    frost_albedo_threshold)
+         write(*,*) " frost_albedo_threshold = ",
+     &            frost_albedo_threshold
+
+!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
+         write(*,*) "temp tag for water caps ?"
+         temptag=.false. 
+         call getin("temptag",temptag)
+         write(*,*) " temptag = ",temptag
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!              
 
          write(*,*) "photochemistry: include chemical species"
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 283)
@@ -43,5 +43,4 @@
 #include "advtrac.h"
 #include "comgeomfi.h"
-#include "watercap.h"
 #include "chimiedata.h"
 
@@ -50,5 +49,5 @@
       real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g.  kg.m-2)
       real co2ice(ngridmx)           ! co2 ice mass on surface (e.g.  kg.m-2)
-      integer iq,ig,count, yeyey
+      integer iq,ig,count
       real r0_lift , reff_lift, nueff_lift
 c     Ratio of small over large dust particles (used when both 
@@ -436,5 +435,4 @@
 	
 
-
       enddo ! of do iq=1,nqmx
 !      count=count+nbqchem
@@ -587,74 +585,4 @@
          alpha_lift(igcm_h2o_vap) =0.
          alpha_devil(igcm_h2o_vap)=0.
-
-c       "Dryness coefficient" controlling the evaporation and
-c        sublimation from the ground water ice (close to 1)
-c        HERE, the goal is to correct for the fact
-c        that the simulated permanent water ice polar caps
-c        is larger than the actual cap and the atmospheric
-c        opacity not always realistic.
-
-         do ig=1,ngridmx
-           if (ngridmx.ne.1) watercaptag(ig)=.false.
-           dryness(ig) = 1.
-         enddo
-
-         IF (caps) THEN
-c Perennial H20 north cap defined by watercaptag=true (allows surface to be
-c hollowed by sublimation in vdifc).
-         yeyey = 0
-         do ig=1,ngridmx
-!          !!! TESTS TESTS outliers
-!          !!! TESTS TESTS outliers
-!          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
-!     .         ( lati(ig)*180./pi      .le.  77 ) .and.
-!     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
-!     .              ( long(ig)*180./pi .le. 120. ) ) 
-!     .             .or.
-!     .             ( ( long(ig)*180./pi .ge. -130. ) .and.
-!     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
-!             if (yeyey .eq. 0) then  !!! 1/2 en 64x48 sinon trop large en lat
-!              write(*,*) "outliers ", lati(ig)*180./pi, long(ig)*180./pi
-!              if (ngridmx.ne.1) watercaptag(ig)=.true.
-!              dryness(ig) = 1.
-!              albedodat(ig) = 0.45 !! comme alb_surfice
-!              yeyey = 1
-!             else
-!              yeyey = 0
-!             endif
-!          endif
-!          !!! TESTS TESTS outliers
-!          !!! TESTS TESTS outliers
-!
-!          !!! TESTS TESTS addcap
-!          !!! TESTS TESTS addcap     
-!          if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
-!     .         ( lati(ig)*180./pi      .le.  84 ) .and.
-!     .         ( ( long(ig)*180./pi .gt. -030. ) .and.
-!     .              ( long(ig)*180./pi .lt. 090. ) ) ) then
-!              write(*,*) "capadd ", lati(ig)*180./pi, long(ig)*180./pi
-!              if (ngridmx.ne.1) watercaptag(ig)=.true.
-!              albedodat(ig) = 0.45 !! comme alb_surfice
-!              dryness(ig) = 1.
-!          endif
-!          !!! TESTS TESTS addcap
-!          !!! TESTS TESTS addcap 
-    
-           if (lati(ig)*180./pi.gt.84) then
-             if (ngridmx.ne.1) watercaptag(ig)=.true.
-             dryness(ig) = 1.
-c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
-c	      if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
-c	        if (ngridmx.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c	      endif
-c             if (lati(ig)*180./pi.ge.85) then
-c               if (ngridmx.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c	      endif
-           endif  ! (lati>80 deg)
-         end do ! (ngridmx)
-        ENDIF ! (caps)
-
          if(water.and.(nqmx.ge.2)) then
            radius(igcm_h2o_ice)=3.e-6
Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_caps.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_caps.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_caps.F	(revision 283)
@@ -1,2 +1,10 @@
+
+      !!!! see meso_inc_ini 
+      ! alb_lim = 0.26
+      ! lat_lim = 70.
+      ! inertie_lim = 800.
+
+      !!! This has to go after initracer which change dryness and watercaptag
+
       !!!!!! MARS MESOSCALE MODELING
       !!!!!! TEST TEST TEST TEST  AS+JBM 28/02/11
@@ -14,18 +22,20 @@
       !!!!           du GCM lorsqu'ils sont consequents
       !!!! 
-      IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
+      !IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
+      IF ( caps .and. water ) THEN
           PRINT *, 'OVERWRITING watercaptag DEFINITION in INITRACER'
-          PRINT *, 'lat>70 et alb>0.26 => watercaptag=T' 
+          PRINT *, 'lat > lat_lim et alb > alb_lim => watercaptag=T'
+          PRINT *, 'ind for water ice: ', igcm_h2o_ice 
           !! Perennial H20 north cap defined by watercaptag=true (allows surface to be
           !! hollowed by sublimation in vdifc).
           do ig=1,ngridmx
             qsurf(ig,igcm_h2o_ice)=0.  !! on jette les inputs GCM
-            if ( (lati(ig)*180./pi.gt.70.) .and.
-     .           (albedodat(ig).ge.0.26) )  then
-                    watercaptag(ig)=.true.
-                    dryness(ig) = 1.
+            if ( ( lati(ig)*180./pi .gt. lat_lim ) .and.
+     .           ( albedodat(ig) .ge. alb_lim    ) )  then
+                    watercaptag(ig)  = .true.
+                    dryness(ig)      = 1.
             else
-                    watercaptag(ig)=.false.
-                    dryness(ig) = 1.
+                    watercaptag(ig)  = .false.
+                    dryness(ig)      = 1.
             endif  ! (lati, albedodat)
           end do ! (ngridmx)
Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_ini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_ini.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_ini.F	(revision 283)
@@ -68,2 +68,26 @@
 c
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      !!! see meso_inc_caps
+      !!! this is a test to change outliers' albedo and thermal inertia
+      alb_lim = 0.26
+      lat_lim = 70.
+      inertie_lim = 800. 
+      PRINT *, 'lat_lim ',lat_lim
+      PRINT *, 'alb_lim ',alb_lim
+      PRINT *, 'inertie_lim ',inertie_lim 
+      !!!
+      !!!
+      IF ( caps .and. water ) THEN
+          do ig=1,ngridmx
+            if ( lati(ig)*180./pi .gt. lat_lim ) then
+               if ( albedodat(ig) .ge. alb_lim ) then
+                    albedodat(ig) = alb_surfice
+                    inertiedat(ig,1) = inertie_lim  
+               endif
+               if (inertiedat(ig,1) .ge. inertie_lim ) then
+                    inertiedat(ig,1) = inertie_lim
+               endif
+            endif  ! (lati, albedodat)
+          end do ! (ngridmx)
+      ENDIF ! (caps)
Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_var.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_var.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_var.F	(revision 283)
@@ -18,2 +18,6 @@
       INTEGER tracerset    !!! this corresponds to config%mars
       CHARACTER (len=20) :: wtnom(nqmx) ! tracer name
+
+      REAL alb_lim 
+      REAL lat_lim 
+      REAL inertie_lim 
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 283)
@@ -128,5 +128,4 @@
 
 #include "chimiedata.h"
-#include "watercap.h"
 #include "param.h"
 #include "param_v3.h"
@@ -191,5 +190,7 @@
       REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
       REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
-      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
+      
+      REAL watercapflag(ngridmx)     ! water cap flag
 
 c     Variables used by the water ice microphysical scheme:
@@ -197,7 +198,4 @@
       REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
                                      !   of the size distribution
-c     Albedo of deposited surface ice
-      !!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
-      REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB 
 
 c     Variables used by the slope model
@@ -422,5 +420,6 @@
 
         IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
-          write(*,*)"physiq: water_param Surface ice alb:",alb_surfice
+          write(*,*)"physiq: water_param Surface water ice albedo:", 
+     .                  albedo_h2o_ice
         ENDIF
                    
@@ -1095,5 +1094,5 @@
 #endif
 c       -------------------------------------------------------------
-c       Change of surface albedo (set to 0.4) in case of ground frost
+c       Change of surface albedo in case of ground frost
 c       everywhere except on the north permanent cap and in regions
 c       covered by dry ice. 
@@ -1102,7 +1101,10 @@
          do ig=1,ngrid
            if ((co2ice(ig).eq.0).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
-              albedo(ig,1) = alb_surfice
-              albedo(ig,2) = alb_surfice
+     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
+              albedo(ig,1) = albedo_h2o_ice
+              albedo(ig,2) = albedo_h2o_ice
+c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
+c              write(*,*) "physiq.F frost :"
+c     &        ,lati(ig)*180./pi, long(ig)*180./pi
            endif
          enddo  ! of do ig=1,ngrid
@@ -1512,4 +1514,16 @@
      &                       'surface h2o_ice',
      &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+
+            if (caps) then
+             do ig=1,ngridmx
+                if (watercaptag(ig)) watercapflag(ig) = 1
+             enddo
+             CALL WRITEDIAGFI(ngridmx,'watercaptag',
+     &                         'Ice water caps',
+     &                         '',2,watercapflag)
+            endif
+            CALL WRITEDIAGFI(ngridmx,'albedo',
+     &                         'albedo',
+     &                         '',2,albedo(1:ngridmx,1))
            endif !(water)
 
Index: trunk/LMDZ.MARS/libf/phymars/soil.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 283)
@@ -18,4 +18,7 @@
 
 #include"comsoil.h"
+
+#include"surfdat.h"
+#include"callkeys.h"
 
 c-----------------------------------------------------------------------
@@ -53,8 +56,14 @@
 ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
       do ig=1,ngrid
-        do ik=0,nsoil-1
-	  mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
-!	  write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
-	enddo
+        if (watercaptag(ig)) then
+          do ik=0,nsoil-1
+! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
+              mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
+          enddo
+        else
+          do ik=0,nsoil-1
+           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
+          enddo
+        endif
       enddo
 
Index: trunk/LMDZ.MARS/libf/phymars/surfdat.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfdat.h	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/surfdat.h	(revision 283)
@@ -5,7 +5,9 @@
      &   zmea(ngridmx),zstd(ngridmx),                                   &
      &   zsig(ngridmx),zgam(ngridmx),zthe(ngridmx),                     &
-     &   z0(ngridmx),z0_default
+     &   z0(ngridmx),z0_default, albedo_h2o_ice, inert_h2o_ice,         &
+     &   frost_albedo_threshold
 
-      COMMON/surfdatl/TESicealbedo
+      COMMON/surfdatl/TESicealbedo,watercaptag,temptag
+
 
       real albedodat ! albedo of bare ground
@@ -15,4 +17,11 @@
       real emissiv ! emissivity of bare ground
       logical TESicealbedo ! use TES ice cap albedoes (if set to .true.)
+      logical watercaptag(ngridmx) ! flag for water ice surface
+      
+      logical temptag !temp tag for water caps
+      
+      real albedo_h2o_ice ! water ice albedo
+      real inert_h2o_ice ! water ice thermal inertia
+      real frost_albedo_threshold ! water frost thickness on the ground (kg.m^-2, ie mm)
       real TESice_Ncoef ! coefficient for TES ice albedo in Northern hemisphere
       real TESice_Scoef ! coefficient for TES ice albedo in Southern hemisphere
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 283)
@@ -14,6 +14,9 @@
 #include "callkeys.h"
 #include "tracer.h"
-c
-      INTEGER ngrid,ig,icap
+#include "comgeomfi.h"
+#include "comconst.h"
+
+c
+      INTEGER ngrid,ig,icap,iq,alternate
       REAL  piceco2(ngrid),psolaralb(ngrid,2)
       REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
@@ -24,28 +27,175 @@
 c=======================================================================
 
-c
-c     calcul de piceco2 (kg/m2) a l'etat initial
+
+c     water ice outliers
 c     ------------------------------------------
 
-      DO 100 ig=1,ngrid
-         psolaralb(ig,1)=albedodat(ig)
-         psolaralb(ig,2)=albedodat(ig)
-100   CONTINUE
-
-      PRINT*,'minimum des donnees albedo',
+      IF ((water) .and. (caps)) THEN
+      
+c Perennial H20 north cap defined by watercaptag=true (allows surface to be
+c hollowed by sublimation in vdifc).
+
+c We might not want albedodat to be modified because it is used to write 
+c restart files. Instead, albedo is directly modified when needed (i.e. 
+c if we have watercaptag and no co2 ice), below and in albedocaps.F90
+
+c       "Dryness coefficient" controlling the evaporation and
+c        sublimation from the ground water ice (close to 1)
+c        HERE, the goal is to correct for the fact
+c        that the simulated permanent water ice polar caps
+c        is larger than the actual cap and the atmospheric
+c        opacity not always realistic.
+
+         alternate = 0
+         
+         do ig=1,ngridmx
+         
+         dryness (ig) = 1
+         
+
+c Towards olympia planitia water caps ('relatively' low latitude ones)
+c---------------- proposition par AS --------------------
+c--------------------------------------------------------
+c          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
+c     .         ( lati(ig)*180./pi      .le.  77 ) .and.
+c     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
+c     .             ( long(ig)*180./pi .le. 120. ) ) 
+c     .           .or.
+c     .           ( ( long(ig)*180./pi .ge. -130. ) .and.
+c     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
+c---------------- proposition par TN --------------------
+c---------------- HIGHLY EXPERIMENTAL -------------------
+c--------------------------------------------------------     
+	  if ( ( ( lati(ig)*180./pi .ge. 73.  ) .and. ! cap #1
+     .           ( lati(ig)*180./pi .le. 75.1 ) .and.
+     .           ( long(ig)*180./pi .ge. 95.  ) .and.
+     .           ( long(ig)*180./pi .le. 110. ) )
+     .		.or.
+     .         ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
+     .           ( lati(ig)*180./pi .le. 80.  ) .and.
+     .           ( long(ig)*180./pi .ge. 110. ) .and.
+     .           ( long(ig)*180./pi .le. 140. ) )
+     .		.or.
+     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
+     .           ( lati(ig)*180./pi .le. 78.  ) .and.
+     .           ( long(ig)*180./pi .ge. 155. ) .and.
+     .           ( long(ig)*180./pi .le. 180. ) )
+c     .		.or.
+c     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
+c     .           ( lati(ig)*180./pi .le. 72.  ) .and.
+c     .           ( long(ig)*180./pi .ge. 163. ) .and.
+c     .           ( long(ig)*180./pi .le. 168. ) )
+     .		.or.
+     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
+     .           ( lati(ig)*180./pi .le. 78.  ) .and.
+     .           ( long(ig)*180./pi .ge. -160.) .and.
+     .           ( long(ig)*180./pi .le. -120.) ) )
+     .	        then
+     
+             if (temptag) then
+             
+               if ((alternate .eq. 0)) then  !!! 1/2 en 64x48 sinon trop large en lat
+             	 if (ngridmx.ne.1) watercaptag(ig)=.true.
+                  write(*,*) "outliers ", lati(ig)*180./pi,
+     .              long(ig)*180./pi
+              	!dryness(ig) = 1.
+                 alternate = 1
+               else
+              	 alternate = 0
+               endif !end if alternate = 0
+               
+             else
+             
+               if (ngridmx.ne.1) watercaptag(ig)=.true.
+                  write(*,*) "outliers ", lati(ig)*180./pi,
+     .              long(ig)*180./pi
+     
+             endif ! end if temptag
+             
+           endif
+
+
+c Opposite olympia planitia water cap
+c---------------- proposition par AS --------------------
+c-------------------------------------------------------- 
+c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
+c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
+c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
+c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
+c---------------- proposition par TN --------------------
+c-------------------------------------------------------- 
+           if ( ( lati(ig)*180./pi     .ge.  80 ) .and.
+     .          ( lati(ig)*180./pi     .le.  84 ) .and.
+     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
+     .            ( long(ig)*180./pi .lt.  090. ) ) ) then
+              if (ngridmx.ne.1) then
+                watercaptag(ig)=.true.
+                write(*,*) "central cap add ", lati(ig)*180./pi,
+     .            long(ig)*180./pi
+              endif
+              !dryness(ig) = 1.
+           endif
+
+c Central cap
+c---------------- proposition par TN --------------------
+c-------------------------------------------------------- 
+           if (lati(ig)*180./pi.gt.84) then
+             PRINT*,'central cap', lati(ig)*180./pi,
+     .         long(ig)*180./pi
+             if (ngridmx.ne.1) watercaptag(ig)=.true.
+             !dryness(ig) = 1.
+c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
+c	      if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
+c	        if (ngridmx.ne.1) watercaptag(ig)=.true.
+c               dryness(ig) = 1.
+c	      endif
+c             if (lati(ig)*180./pi.ge.85) then
+c               if (ngridmx.ne.1) watercaptag(ig)=.true.
+c               dryness(ig) = 1.
+c	      endif
+           endif  ! (lati>84 deg)
+           
+         end do ! (ngridmx)
+        ENDIF ! (caps & water)
+
+c ===============================================================
+c      INITIAL ALBEDO
+c ===============================================================
+
+         write(*,*)"surfini: water frost thickness",
+     s     frost_albedo_threshold
+         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
+         write(*,*)"surfini: water ice TI:", inert_h2o_ice
+
+c        To start with : Initial albedo = observed dataset 
+c        -------------------------------------------------
+         DO ig=1,ngrid
+              psolaralb(ig,1)=albedodat(ig)
+              psolaralb(ig,2)=albedodat(ig)
+         END DO
+         PRINT*,'minimum albedo sans water caps',
      s     albedodat(ISMIN(ngrid,albedodat,1))
-      PRINT*,'maximum des donnees albedo',
+         PRINT*,'maximum albedo sans water caps',
      s     albedodat(ISMAX(ngrid,albedodat,1))
-c      calcul de psolaralb 
-c      -------------------
-      DO 115 ig=1,ngrid
-
-c	 IF (water) THEN 
-c	   if (qsurf(ig,nqmx).gt.0.005) then
-c             psolaralb(ig,1) = 0.4
-c             psolaralb(ig,2) = 0.4 
-c           endif 
-c         ENDIF
-c IF there is more than 5 pr. um of h2o ice but no C02 ice, surface albedo is set to 0.4. 
+
+c        initial albedo if permanent H2O ice is present 
+c        ------------------------------------------------
+         IF ((water) .and. (caps)) THEN
+           DO ig=1,ngrid
+            IF (watercaptag(ig)) THEN
+              psolaralb(ig,1) = albedo_h2o_ice
+              psolaralb(ig,2) = albedo_h2o_ice
+            ENDIF
+           END DO
+           PRINT*,'minimum albedo avec water caps',
+     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
+           PRINT*,'maximum albedo avec water caps',
+     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
+         ENDIF
+
+c      changing initial albedo if CO2 ice is present 
+c      -------------------------------------------
+
+       DO ig=1,ngrid
          IF (piceco2(ig) .GT. 0.) THEN
              IF(ig.GT.ngrid/2+1) THEN
@@ -55,12 +205,31 @@
              ENDIF
              psolaralb(ig,1) = albedice(icap)
-             psolaralb(ig,2) =  albedice(icap)
+             psolaralb(ig,2) = albedice(icap)    
          END IF
-115   CONTINUE      
-
-      PRINT*,'minimum des donnees albedo',
+       END DO
+
+c      changing initial albedo if water ice frost is present
+c      -------------------------------------------
+       IF (water) THEN
+          do iq=1,nqmx
+c          if there is frost and surface albedo is set to albedo_h2o_ice
+           if(noms(iq).eq."h2o_ice") then
+             do ig=1,ngrid
+              if ((piceco2(ig) .eq. 0.).and.
+     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
+                     psolaralb(ig,1) = albedo_h2o_ice
+                     psolaralb(ig,2) = albedo_h2o_ice
+c                     PRINT*,'surfini.F frost',
+c     &                  lati(ig)*180./pi, long(ig)*180./pi
+               endif 
+              enddo
+           endif 
+          end do
+          PRINT*,'minimum albedo avec givre et co2',
      s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
-      PRINT*,'maximum des donnees albedo',
+          PRINT*,'maximum albedo avec givre et co2',
      s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
+       END IF
+         
 
       RETURN
Index: trunk/LMDZ.MARS/libf/phymars/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 283)
@@ -40,5 +40,4 @@
 #include "netcdf.inc"
 #include "comg1d.h"
-#include "watercap.h"
 #include "logic.h"
 #include "advtrac.h"
@@ -572,7 +571,7 @@
 
       call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
-     .              dtphys,float(day0),
-     .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
-     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+     .              dtphys,float(day0),time,tsurf,
+     .              tsoil,co2ice,emis,q2,qsurf,area,albedodat,
+     .              inertiedat,zmea,zstd,zsig,zgam,zthe)
 c=======================================================================
 c  BOUCLE TEMPORELLE DU MODELE 1D 
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 282)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 283)
@@ -36,5 +36,4 @@
 #include "tracer.h"
 
-#include "watercap.h"
 c
 c   arguments:
Index: trunk/LMDZ.MARS/libf/phymars/watercap.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercap.h	(revision 282)
+++ 	(revision )
@@ -1,6 +1,0 @@
-c-----------------------------------------------------------------------
-c INCLUDE 'watercap.h'
-
-      logical watercaptag(ngridmx)
-
-      common/permcap/watercaptag
