Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2501)
+++ trunk/LMDZ.MARS/README	(revision 2502)
@@ -3329,2 +3329,8 @@
   where the addition of a dynamical timestep is no longer needed
   
+== 22/04/2021 == JN
+- Adding a new option for icelocationmode in surfini.F, which uses the 
+  predefined 64x48 water caps for any resolution. Now default option 
+  should be icelocationmode=4.
+
+
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 2501)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 2502)
@@ -4,5 +4,6 @@
       use netcdf
       use tracer_mod, only: nqmx, noms
-      use geometry_mod, only: longitude, latitude ! in radians
+      use geometry_mod, only: longitude, latitude, ! in radians
+     &                     cell_area ! for watercaptag diagnosis
       use surfdat_h, only: watercaptag, frost_albedo_threshold,
      &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
@@ -36,10 +37,12 @@
       LOGICAL, DIMENSION(100000) :: longwatercaptag
       
-! There are 3 different modes for ice distribution:
+! There are 4 different modes for ice distribution:
 ! icelocationmode = 1 ---> based on data from surface.nc
 ! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
 ! icelocationmode = 3 ---> based on logical relations for latitude and longitude
+! icelocationmode = 4 ---> predefined 64x48 but usable with every
+! resolution, and easily adaptable for dynamico 
 ! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
-      INTEGER,SAVE :: icelocationmode = 2
+      INTEGER,SAVE :: icelocationmode = 4
        
        
@@ -430,9 +433,18 @@
          end do ! of (klon_glo)
 
+        ELSE IF (icelocationmode .eq. 4) THEN
+      
+         print*,'icelocationmode = 4'
+         print*,'Surfini: ice caps defined using manual 64x48 settings'
+         print*,'(although, it should work with any resolution)'
+         call locate_watercaptag(klon_glo,lati_glo,
+     &            long_glo,watercaptag_glo)
+
+!         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
 
         ELSE
       
          print*, 'In surfini.F, icelocationmode is ', icelocationmode
-         print*, 'It should be 1, 2 or 3.'
+         print*, 'It should be 1, 2, 3 or 4 (default is 4)'
          call abort_physic("surfini","wrong icelocationmode",1)
 
@@ -448,4 +460,10 @@
              print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
      &              long_glo(ig)*180./pi, ig
+!             write(1,*),ig, lati_glo(ig)*180./pi,
+!     &              cell_area(ig)
+!             write(2,*), lati_glo(ig)*180./pi,
+!     &              long_glo(ig)*180./pi,cell_area(ig)
+!             write(3,*), ig, lati_glo(ig)*180./pi,
+!     &              long_glo(ig)*180./pi,cell_area(ig)
           endif
         enddo
@@ -464,4 +482,122 @@
 ! end of #else of #ifndef MESOSCALE
 #endif       
-
-      END
+!      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
+      END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
+
+      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,
+     &            long_glo,watercaptag_glo)
+
+      USE comcstfi_h, ONLY: pi
+
+      integer, intent(in) :: klon_glo
+      real, intent(in) :: lati_glo(klon_glo)
+      real, intent(in) :: long_glo(klon_glo)
+      logical, intent(out) :: watercaptag_glo(klon_glo)
+      integer :: ig,i
+!      real, dimension(klon_glo,120) :: wcap
+      real, dimension(120,2) :: latedge
+      real, dimension(120,2) :: lonedge
+
+! In icelocationmode=2 there are 120 manually predefined grid points where
+! watercaptag is true (for the 64x48 resolution). The grid cells corners
+! coordinates in latitude and longitude are written below. With this
+! routine, we check if the grid cell center is in between any of those
+! points. If so, watercaptag = true. 
+
+
+
+
+      latedge(:,1)=(/ 
+     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375, 
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
+     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
+     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
+     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
+     & 80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875,
+     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125,
+     & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
+
+
+      latedge(:,2)=(/ 
+     & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,   
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
+     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
+     & 84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
+     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875,
+     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
+
+
+      lonedge(:,1)=(/
+     &-180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125,
+     &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625,
+     &-120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125,
+     & -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625,
+     & -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125,
+     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
+     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
+     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
+     &  81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875,
+     & 115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375,
+     & 149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875,
+     &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375,
+     & -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125,
+     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
+     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
+     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
+     &  81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125,
+     &  -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125,
+     & 149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875,
+     &  98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/) 
+
+      lonedge(:,2)=(/ 
+     & 180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875,
+     &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375,
+     &-115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875,
+     & -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375,
+     & -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875,
+     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
+     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
+     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
+     &  87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125,
+     & 120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625,
+     & 154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625,
+     &-120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125,
+     & -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875,
+     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
+     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
+     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
+     &  87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875,
+     &  -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375,
+     & 154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625,
+     & 104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
+
+
+      watercaptag_glo(:) = .false.
+      DO ig=1, klon_glo
+        DO i=1, 120
+           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))
+     &         .and.(long_glo(ig)*180./pi.le.lonedge(i,2))
+     &         .and.(lati_glo(ig)*180./pi.ge.latedge(i,1))
+     &         .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
+             watercaptag_glo(ig) = .true.
+           endif
+        ENDDO !i=1, 120
+      ENDDO ! ig=1, klon_glo
+
+      END SUBROUTINE locate_watercaptag
