Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 710)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 712)
@@ -45,7 +45,6 @@
       INTEGER     imd,jmd
       PARAMETER   (imd=360,jmd=180)
-      REAL        zdata(imd*jmd)
+      REAL        zdata(imd,jmd)
       REAL        zelat,zelon 
-! on a data(i,j) = zdata((j-1)*jmd + i) avec i longitude de 1 a imd (E vers O) et j lat de 1 a jmd (N vers S)
 
       INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
@@ -190,20 +189,20 @@
       do ig=2,ngridmx-1
       
-        ! loop over the surface file grid
-        do i=1,imd*jmd
-         zelat = 90. - 90./real(jmd) - (i-1)/imd      ! surface file latitude
-         zelon = mod(i-1,imd) - (180.-180./real(imd)) ! surface file longitude
-                  
-          if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
+        ! loop over the surface file grid      
+        do i=1,imd
+          do j=1,jmd
+            zelon = i - 180.
+            zelat = 90. - j 
+            if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
      &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
               ! count all points in that GCM grid point
               nb_ice(ig,1) = nb_ice(ig,1) + 1
-              if (zdata(i) > min_icevalue)
+              if (zdata(i,j) > min_icevalue)
                  ! count all detected points in that GCM grid point
      &           nb_ice(ig,2) = nb_ice(ig,2) + 1
-          endif
-          
-        enddo ! of do i=1,imd*jmd
-        
+             endif
+          enddo
+        enddo  
+
         ! projection of nb_ice on GCM lat and lon axes
         latice(1+(ig-2)/iim,:) =
@@ -224,15 +223,15 @@
       
      
-  !    print*, 'latice TOT', latice(:,1)
-  !    print*, 'latice FOUND', latice(:,2)
-  !    print*, 'lonice TOT', lonice(:,1)
-  !    print*, 'lonice FOUND', lonice(:,2)
-      
-  !    print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
-  !    print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
-      
-  !    print*,''
-  !    print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
-  !    print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
+!      print*, 'latice TOT', latice(:,1)
+!      print*, 'latice FOUND', latice(:,2)
+!      print*, 'lonice TOT', lonice(:,1)
+!      print*, 'lonice FOUND', lonice(:,2)
+      
+!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
+!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
+      
+!      print*,''
+!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
+!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
       
     
@@ -480,4 +479,22 @@
            if(noms(iq).eq."h2o_ice") then
              do ig=1,ngrid
+             
+              if ((watercaptag(ig).eqv..false.) 
+     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
+              print*, ''
+              print*, '!!! PROBLEM in SURFINI !!!!'
+              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
+     & WATERCAPTAG IS FALSE'
+              print*, ''
+              print*, 'ig,qsurf,threshold' , 
+     &         ig, qsurf(ig,iq), -frost_albedo_threshold
+              print*, ''
+              print*, '1) Check h2o_ice in startfi and ice
+     & distribution in surfini'
+              print*, '2) Use ini_h2osurf option in newstart'
+              print*, ''
+              CALL ABORT
+            endif
+             
               if ((piceco2(ig) .eq. 0.).and.
      &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then