Index: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F	(revision 3370)
+++ trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F	(revision 3371)
@@ -6,11 +6,10 @@
 c   Auteur:   Christophe Hourdin/Francois Forget/Yann Wanherdrick
 c   ------
-c             Derniere modif : 12/03
-c
+c             Adapted to Pluto: Tanguy Bertrand
+c             Derniere modif : 06/2024
 c
 c   Objet:  Create or modify the initial state for the LMD Mars GCM
 c   -----           (fichiers NetCDF start et startfi)
 c
-c
 c=======================================================================
 
@@ -19,10 +18,9 @@
      &                              is_master
       use infotrac, only: infotrac_init, nqtot, tname
-      USE tracer_h, ONLY: igcm_n2
+      USE tracer_h, ONLY: igcm_n2,igcm_ch4_ice,igcm_co_ice,igcm_haze,
+     &                   igcm_prec_haze,igcm_co_gas,igcm_ch4_gas,noms
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
       USE surfdat_h, ONLY: phisfi, albedodat,
      &                     zmea, zstd, zsig, zgam, zthe
-! TB24      USE nonoro_gwd_ran_mod, ONLY: du_nonoro_gwd, dv_nonoro_gwd,
-!     &                              east_gwstress, west_gwstress
       use datafile_mod, only: datadir, surfdir
       use ioipsl_getin_p_mod, only: getin_p
@@ -30,6 +28,4 @@
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
-!      use slab_ice_h, only:noceanmx
-!      USE ocean_slab_mod, ONLY: nslay
       use filtreg_mod, only: inifilr
       USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -66,5 +62,4 @@
       CHARACTER        relief*3
 
-
 c Variables pour les lectures NetCDF des fichiers "start_archive" 
 c--------------------------------------------------
@@ -90,6 +85,4 @@
       REAL w(iip1,jjp1,llm+1)
       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-!      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-!      REAL dh(ip1jmp1,llm),dp(ip1jmp1)
       REAL phi(iip1,jjp1,llm)
 
@@ -106,26 +99,25 @@
       REAL tsurf(ngridmx)        ! surface temperature
       REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
-!      REAL n2ice(ngridmx)        ! N2 ice layer
       REAL emis(ngridmx)        ! surface emissivity
       real emisread             ! added by RW
       REAL,ALLOCATABLE :: qsurf(:,:)
+      REAL,ALLOCATABLE :: qsurf_tmp(:,:)
       REAL q2(ngridmx,llm+1)
-!      REAL rnaturfi(ngridmx)
       real alb(iip1,jjp1),albfi(ngridmx) ! albedos
       real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
       real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
-!      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
-
-      INTEGER i,j,l,isoil,ig,idum
+      REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
+
+      INTEGER i,j,l,isoil,ig,idum,k
       real mugaz ! molar mass of the atmosphere
 
-      integer ierr,iref 
+      integer ierr 
 
 c Variables on the new grid along scalar points 
 c------------------------------------------------------
-!      REAL p(iip1,jjp1)
       REAL t(iip1,jjp1,llm)
       REAL tset(iip1,jjp1,llm)
       real phisold_newgrid(iip1,jjp1)
+      real phisold(iip1,jjp1)
       REAL :: teta(iip1, jjp1, llm)
       REAL :: pk(iip1,jjp1,llm)
@@ -135,5 +127,4 @@
       REAL :: xpn,xps,xppn(iim),xpps(iim)
       REAL :: p3d(iip1, jjp1, llm+1)
-!      REAL dteta(ip1jmp1,llm)
 
 c Variable de l'ancienne grille 
@@ -145,5 +136,5 @@
 c variables diverses
 c-------------------
-      real choix_1,pp
+      real choix_1,pp,choice
       character*80      fichnom
       character*250  filestring
@@ -153,9 +144,11 @@
       real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
       real ptoto,pcap,patm,airetot,ptotn,patmn,psea
-!      real ssum
       character*1 yes
       logical :: flagtset=.false. ,  flagps0=.false.
-      real val, val2, val3, val4 ! to store temporary variables
+      real val, val1, val2, val3, val4, dist, ref ! to store temporary variables
+      real val5, val6, val7, val8, val9, val10,val11, val12 ! to store temporary variables
       real :: iceith=2000 ! thermal inertia of subterranean ice
+      integer :: iref,jref,iref1,jref1,iref2,jref2
+      integer :: igref,igref1,igref2,igref3
 
       INTEGER :: itau
@@ -174,17 +167,43 @@
       real fact2
 
+c Special Pluto Map from Lellouch & Grundy 
+c------------------------------------------
+       integer,parameter :: im_plu=360 ! from topo
+       integer,parameter :: jm_plu=180   
+       real latv_plu(jm_plu+1),lonu_plu(im_plu+1)
+       real map_pluto_dat(im_plu,jm_plu+1)
+
+       real N2_pluto_dat(im_plu,jm_plu+1)
+       real CH4_pluto_dat(im_plu,jm_plu+1)
+       real CO_pluto_dat(im_plu,jm_plu+1)
+       real alb_dat(im_plu,jm_plu+1)
+
+       real N2_pluto_rein(im_plu+1,jm_plu+1)
+       real CH4_pluto_rein(im_plu+1,jm_plu+1)
+       real CO_pluto_rein(im_plu+1,jm_plu+1)
+       real alb_rein(im_plu+1,jm_plu+1)
+       real N2_pluto_gcm(iip1,jjp1)
+       real CH4_pluto_gcm(iip1,jjp1)
+       real CO_pluto_gcm(iip1,jjp1)
+       real alb_gcm(iip1,jjp1)
+
+c Special Topo Map mountain
+       real lat0, lon0
+       integer ig0
+       real dist_pluto
+       real radm,height, phisprev, temp
+       real preffnew,panew
+c Special copy of terrain
+       real actualmin,angle
+       integer array_ind(ngridmx)
+       real array_dist(ngridmx)
+       real array_angle(ngridmx)
 
 c sortie visu pour les champs dynamiques
-c---------------------------------------
-!      INTEGER :: visuid
-!      real :: time_step,t_ops,t_wrt
-!      CHARACTER*80 :: visu_file
-      !TB24
       CALL conf_gcm( 99, .TRUE. )
 
-
       cpp    = 0.
-      preff  = 0.
-      pa     = 0. ! to ensure disaster rather than minor error if we don`t
+      preff  = 2.
+      pa     = 0.2 ! to ensure disaster rather than minor error if we don`t
                   ! make deliberate choice of these values elsewhere.
 
@@ -204,4 +223,5 @@
       allocate(q(iip1,jjp1,llm,nqtot))
       allocate(qsurf(ngridmx,nqtot))
+      allocate(qsurf_tmp(ngridmx,nqtot))
       
 ! get value of nsoilmx and allocate corresponding arrays
@@ -272,5 +292,4 @@
 
       endif
-
 
 c=======================================================================
@@ -354,6 +373,4 @@
      .        nqtot,day_ini,time,
      .        tsurf,tsoil,emis,q2,qsurf)    !) ! temporary modif by RDW
-!     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
-!     .        sea_ice)
 
         ! copy albedo and soil thermal inertia on (local) physics grid
@@ -432,4 +449,21 @@
       endif
 
+c Initialisation coordonnees /aires
+c -------------------------------
+! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
+!       rlatu() and rlonv() are given in radians
+      latfi(1)=rlatu(1)
+      lonfi(1)=0.
+      DO j=2,jjm
+         DO i=1,iim
+            latfi((j-2)*iim+1+i)=rlatu(j)
+            lonfi((j-2)*iim+1+i)=rlonv(i)
+         ENDDO
+      ENDDO
+      latfi(ngridmx)=rlatu(jjp1)
+      lonfi(ngridmx)=0.
+
+      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
+
       rad = p_rad
       omeg = p_omeg
@@ -438,10 +472,4 @@
       mugaz = p_mugaz
       daysec = p_daysec
-
-      if (p_omeg .eq. -9999.) then
-        p_omeg = 8.*atan(1.)/p_daysec
-        print*,"new value of omega",p_omeg
-      endif
-
       kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW
 
@@ -536,5 +564,4 @@
       endif ! of if (choix_1.eq.0)
 
-
 c=======================================================================
 c  Lecture des fichiers (start ou start_archive)
@@ -549,6 +576,4 @@
      &   surfith,nid)
 
-!     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
-! TB24     &   du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -572,19 +597,86 @@
       do ! infinite loop on list of changes
 
-      write(*,*)
-      write(*,*)
-      write(*,*) 'List of possible changes :'
-      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-      write(*,*)
-      write(*,*) 'flat : no topography ("aquaplanet")'
-      write(*,*) 'set_ps_to_preff : used if changing preff with topo'
-      write(*,*) 'qname : change tracer name'
-      write(*,*) 't=profile  : read temperature profile in profile.in'
-      write(*,*) 'q=0 : ALL tracer =zero'
-      write(*,*) 'q=x : give a specific uniform value to one tracer'
-      write(*,*) 'qs=x : give a uniform value to a surface tracer'
-      write(*,*) 'q=profile    : specify a profile for a tracer'
-      write(*,*) 'subsoil_all : set seasonal subsurface thermal inertia'
-      write(*,*) 'diurnal_TI : set diurnal subsurface thermal inertia'
+        write(*,*)
+        write(*,*)
+        write(*,*) 'List of possible changes :'
+        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+        write(*,*)
+        write(*,*) 'flat : no topography ("aquaplanet")'
+        write(*,*) 'set_ps_to_preff : used if changing preff with topo'
+        write(*,*) 'qname : change tracer name'
+        write(*,*) 't=profile  : read temperature profile in profile.in'
+        write(*,*) 'q=0 : ALL tracer =zero'
+        write(*,*) 'q=x : give a specific uniform value to one tracer'
+        write(*,*) 'qs=x : give a uniform value to a surface tracer'
+        write(*,*) 'q=profile    : specify a profile for a tracer'
+        write(*,*) 'qnogcm: initial tracer for GCM from nogcm (CH4,CO)'
+        write(*,*) 'isotherm : Isothermal Temperatures'
+        write(*,*) 'tempprof : specific Temperature profile'
+        write(*,*) 'initsoil : initialize soil Temperatures'
+        write(*,*) 'ptot : change total pressure'
+        write(*,*) 'therm_ini_s: set soil TI to reference surf. values'
+        write(*,*) 'inert3d: give uniform value of thermal inertia'
+        write(*,*) 'subsoilice_n: put deep underground ice in N. hemis'
+        write(*,*) 'subsoilice_s: put deep underground ice in S. hemis'
+        write(*,*) 'reservoir: increase/decrease reservoir of ice'
+        write(*,*) 'reservoir_SP: increase/decrease reservoir in SP'
+        write(*,*) 'plutomap: initialize surface like pluto from ...'
+        write(*,*) 'subsoil_n2: diu-sea thermal inertia for N2 only'
+        write(*,*) 'subsoil_ch4: diu-sea thermal inertia for CH4 only'
+        write(*,*) 'subsoil_all: diu-sea thermal inertia for all terr'
+        write(*,*) 'subsoil_alb: diu-sea thermal inertia from albedo'
+        write(*,*) 'diurnal_TI: diurnal thermal inertia for all terr'
+        write(*,*) 'initsurf: initial surface temperature (global)'
+        write(*,*) 'modtsurf: initial surface temperature (global)'
+        write(*,*) 'copylat: copy tsurf and tsoil latitude'
+        write(*,*) 'copylatlon: copy tsurf and tsoil lat/lon'
+        write(*,*) 'copylon: copy tsurf tsoil lat, n2surf and phisfi'
+        write(*,*) 'tsurfice: temperature depending on N2 ice'
+        write(*,*) 'icarus: option to change surface/soil temperature'
+        write(*,*) 'icarus_ch4: special option to add ch4' 
+        write(*,*) 'delfrostch4: delete ch4 frost' 
+        write(*,*) 'modch4: remove/modify amount ch4 frost' 
+        write(*,*) 'modch4n2: modify amount ch4 frost according N2' 
+        write(*,*) 'modco: remove/modify amount co frost' 
+        write(*,*) 'modn2: remove/modify amount n2 ice'
+        write(*,*) 'modcoch4: remove/modify co ch4 where no n2 '
+        write(*,*) 'checktsoil: change tsoil where no n2 '
+        write(*,*) 'non2noco: if no n2, no co '
+        write(*,*) 'modn2_topo: modify n2 ice topo according to topo' 
+        write(*,*) 'modwheren2: modify n2 where already n2' 
+        write(*,*) 'modn2HD: modify n2 for HD runs' 
+        write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP' 
+        write(*,*) 'globch4co: add/remove global amount ch4co frost' 
+        write(*,*) 'source_ch4: add source ch4' 
+        write(*,*) 'nomountain: remove pluto mountains (!)' 
+        write(*,*) 'relief: add pluto crater or mountain' 
+        write(*,*) 'seticeNH: set ice for initial start with NH topo'
+        write(*,*) 'topo_sp: change topography of SP'
+        write(*,*) 'fill_sp: fill sp with N2 ice and adjust phisfi'
+        write(*,*) 'fill_sp_inv: fill inverted sp and adjust'
+        write(*,*) 'subsoil_SP: diu-sea thermal inertia for SP '
+        write(*,*) 'bladed: add ch4 on bladed terrains'
+        write(*,*) 'cpbladed: add extension bladed terrains'
+        write(*,*) 'slope: add slope over all long. (for triton)'
+        write(*,*) 'digsp: specific to dig SP'
+        write(*,*) 'bugn2: correct bug of warm n2 due to HighRes'
+        write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes'
+        write(*,*) 'flatnosp : topo flat except SP (topo controled)'
+        write(*,*) 'flatregion: topo flat for specific region'
+        write(*,*) 'changetopo: change topo'
+        write(*,*) 'asymtopo: N-S asym topo in tanh'
+        write(*,*) 'corrtopo: correction topo bug'
+        write(*,*) 'adjustphi: adjust phisfi according to N2 mass'
+        write(*,*) 'correctphi: adjust phisfi'
+        write(*,*) 'correctps: test to correct ps'
+        write(*,*) 'toponoise: no flat topo'
+        write(*,*) 'asymtriton: asymetry in longitude for triton'
+        write(*,*) 'copytsoil: copy 2D tsoil field from nc file'
+        write(*,*) 'albedomap: read in an albedomap albedo.nc'
+        write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo'
+      
+        write(*,*)
+        write(*,*) 'Please note that emis and albedo are set in physiq'
+        write(*,*)
 
         write(*,*)
@@ -825,79 +917,976 @@
              endif
 
-
-
-c       subsoil_all : initialize subsurface thermal inertia
-c       --------------------------------------------------
-        else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then
-
-          write(*,*) 'New value for subsoil thermal inertia  ?'
- 104      read(*,*,iostat=ierr) ith_bb
-          if(ierr.ne.0) goto 104
-          write(*,*) 'thermal inertia (new value):',ith_bb
-
-          write(*,*)'At which depth (in m.) does the ice layer begin?'
-          write(*,*)'(here, the deepest soil layer extends down to:'
-     &              ,layer(1),' - ',layer(nsoilmx),')'
-          write(*,*)'write 0 for uniform value for all subsurf levels?'
-          ierr=1
-          do while (ierr.ne.0)
-           read(*,*,iostat=ierr) val2
-           write(*,*)'val2 in subsoil_all:',val2,'ierr=',ierr
-           if (ierr.eq.0) then ! got a value, but do a sanity check
-             if(val2.gt.layer(nsoilmx)) then
-               write(*,*)'Depth should be less than ',layer(nsoilmx)
+c       qnogcm : initialise tracer from nogcm (CH4, CO)  
+c       --------------------------------
+        else if (modif(1:len_trim(modif)).eq.'qnogcm') then
+             write(*,*) 'Which tracer do you want to modify ?'
+             write(*,*) 'Should be ch4_gas and co_gas'
+             do iq=1,nqtot
+               write(*,*)iq,' : ',trim(noms(iq)),' : ',q(1,1,1,iq)
+             enddo
+             write(*,*) '(choose between 1 and ',nqtot,')'
+             read(*,*) iq 
+             DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    q(i,j,l,iq)=q(i,j,1,iq)
+                  ENDDO
+               ENDDO
+             ENDDO
+
+c       isothermal temperatures 
+c       ------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
+
+          write(*,*)'Isothermal temperature of the atmosphere'
+          write(*,*) 'Value of THIS temperature ? :'
+ 203      read(*,*,iostat=ierr) Tset(1,1,1)
+          if(ierr.ne.0) goto 203
+          flagtset=.true.
+          DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    Tset(i,j,l)=Tset(1,1,1)
+                  ENDDO
+               ENDDO
+          ENDDO
+          write(*,*) 'atmospheric temp =', Tset(2,2,2)
+
+c       specific temperature profile 
+c       ------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'tempprof') then
+
+          write(*,*) 'phi='
+          write(*,*) phi(1,1,:)
+          write(*,*) 'temperature profile in the atmosphere'
+          write(*,*) 'Value of THIS temperature ? :'
+ 204      read(*,*,iostat=ierr) Tset(1,1,1)
+          if(ierr.ne.0) goto 204
+          flagtset=.true.
+          DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    Tset(i,j,l)=Tset(1,1,1)
+                  ENDDO
+               ENDDO
+          ENDDO
+          write(*,*) 'atmospheric temp =', Tset(2,2,2)
+
+c       initsoil: subsurface temperature
+c       ---------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'initsoil') then
+
+          write(*,*)'Isothermal temperature of the subsurface'
+          write(*,*) 'Value of this temperature ? :'
+ 303      read(*,*,iostat=ierr) Tiso
+          if(ierr.ne.0) goto 303
+
+          do l=1,nsoilmx
+            do ig=1, ngridmx
+              tsoil(ig,l) = Tiso
+            end do
+          end do
+
+c       icarus : changing tsoil tsurf on latitudinal bands
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'icarus') then
+
+          write(*,*) 'At which lat lon do you want to extract the
+     &                              reference tsurf/tsoil profile ?'
+ 407      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 407
+          write(*,*) 'You want lat =',val
+          write(*,*) 'You want lon =',val2
+          dist=0
+          ref=1000
+          do j=1,ngridmx-1
+                 dist=sqrt((lonfi(j)*180./pi-val2)**2+
+     &                              (latfi(j)*180./pi-val)**2)
+                 IF (dist.lt.ref) then
+                    ref=dist
+                    jref=j
+                 ENDIF
+          ENDDO
+
+          write(*,*)'Will use nearest grid point which is:',
+     &              latfi(jref)*180./pi,lonfi(jref)*180./pi
+
+          ! Extraction of the profile
+          write(*,*) 'Profile is =',tsoil(jref,:)
+          write(*,*) 'tsurf is =',tsurf(jref)
+          write(*,*) 'Choice lat for latitudinal band with this profile'
+ 408      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 408
+          write(*,*) 'Choice delta lat for this latitudinal band'
+ 409      read(*,*,iostat=ierr) val4
+          if(ierr.ne.0) goto 409
+          write(*,*) 'Choice delta temp (K) for this latitudinal band'
+ 410      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 410
+          write(*,*) 'How much N2 ice should I put on this band (kg/m2)'
+ 415      read(*,*,iostat=ierr) choice
+          if(ierr.ne.0) goto 415
+          write(*,*) 'Choice lat2'
+ 411      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 411
+          write(*,*) 'Choice delta lat for this latitudinal band'
+ 412      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 412
+          write(*,*) 'Choice delta temp (K) for this latitudinal band'
+ 413      read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 413
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.(val3-val4/2.))  .and.
+     &              (latfi(ig)*180./pi.le.(val3+val4/2.))  .and.
+     &              (qsurf(ig,igcm_n2).lt.0.001)   )      then
+                    tsurf(ig)=tsurf(jref)+val5
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(jref,l)+val5
+                    ENDDO
+                    qsurf(ig,igcm_n2)=choice
+             ENDIF
+             IF (   (latfi(ig)*180./pi.ge.(val6-val7/2.))  .and.
+     &              (latfi(ig)*180./pi.le.(val6+val7/2.))  .and.
+     &              (qsurf(ig,igcm_n2).lt.0.001)   )      then
+                    tsurf(ig)=tsurf(jref)+val8
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(jref,l)+val8
+                    ENDDO
+             ENDIF
+          ENDDO
+
+c       icarus_ch4 : changing tsoil tsurf and adding ch4
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'icarus_ch4') then
+
+          write(*,*) 'At which lat lon do you want to extract the
+     &                              reference tsurf/tsoil profile ?'
+ 416      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 416
+          write(*,*) 'You want lat =',val
+          write(*,*) 'You want lon =',val2
+          dist=0
+          ref=1000
+          do j=1,ngridmx-1
+                 dist=sqrt((lonfi(j)*180./pi-val2)**2+
+     &                              (latfi(j)*180./pi-val)**2)
+                 IF (dist.lt.ref) then
+                    ref=dist
+                    jref=j
+                 ENDIF
+          ENDDO
+
+          write(*,*)'Will use nearest grid point which is:',
+     &              latfi(jref)*180./pi,lonfi(jref)*180./pi
+
+          ! Extraction of the profile
+          write(*,*) 'Profile is =',tsoil(jref,:)
+          write(*,*) 'tsurf is =',tsurf(jref)
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 417      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 417
+          write(*,*) 'Choice temp coefficient for this latitudinal band'
+ 418      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 418
+          write(*,*) 'How much CH4 ice do I put on this band (kg/m2)'
+ 419      read(*,*,iostat=ierr) choice
+          if(ierr.ne.0) goto 419
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val3)  .and.
+     &              (latfi(ig)*180./pi.le.val4)  .and.
+     &              (qsurf(ig,igcm_ch4_ice).lt.0.001) )      then
+                    tsurf(ig)=tsurf(jref)*val5
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(jref,l)*val5
+                    ENDDO
+                    qsurf(ig,igcm_ch4_ice)=choice
+             ENDIF
+          ENDDO
+
+c       globch4co : adding/remove global ch4 co
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'globch4co') then
+
+          write(*,*) 'Adding or removing ch4 co '
+          write(*,*) 'Choice amount add/less ch4 ice (0 = rm all)'
+ 431      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 431
+          write(*,*) 'Choice amount add/less co ice (0 = rm all)'
+ 432      read(*,*,iostat=ierr) val4
+          if(ierr.ne.0) goto 432
+          IF (val3.eq.0.) then
+             DO ig=1,ngridmx
+                    qsurf(ig,igcm_ch4_ice)=0.
+             ENDDO
+          ELSE
+             DO ig=1,ngridmx
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val3
+             ENDDO
+          ENDIF
+          IF (val4.eq.0.) then
+             DO ig=1,ngridmx
+                    qsurf(ig,igcm_co_ice)=0.
+             ENDDO
+          ELSE
+             DO ig=1,ngridmx
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val4
+             ENDDO
+          ENDIF
+ 
+c       bladed : adding/remove ch4 on bladed terrains
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'bladed') then
+
+          write(*,*) 'Adding or removing ch4 on the bladed terrains'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 450      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 450
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 451      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 451
+          write(*,*) 'Choice phisfi minimum ?'
+ 454      read(*,*,iostat=ierr) choice
+          if(ierr.ne.0) goto 454
+          write(*,*) 'Choice multiplicative factor'
+ 452      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 452
+          write(*,*) 'Choice additional ch4'
+ 453      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 453
+          write(*,*) 'Choice index for tsurf tsoil'
+ 449      read(*,*,iostat=ierr) iref
+          if(ierr.ne.0) goto 449
+
+          ! shift 
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val3)  .and.
+     &              (lonfi(ig)*180./pi.le.val4)  .and.
+     &              (phisfi(ig).gt.choice) )  then     
+               qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5
+               qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
+               tsurf(ig)=tsurf(iref)
+               DO l=1,nsoilmx
+                  tsoil(ig,l)=tsoil(iref,l)
+               ENDDO
+             ENDIF
+          ENDDO
+
+c       cpbladed : add extension bladed terrains
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'cpbladed') then
+
+          write(*,*) 'Choice lat,lon, centre of band to be copied ?'
+ 560      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 560
+          write(*,*) 'Choice distance (km) from there defining the area'
+ 561      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 561
+          write(*,*) 'Nb of copies ?'
+ 562      read(*,*,iostat=ierr) l
+          if(ierr.ne.0) goto 562
+
+          ! Get index correponding to central points  
+          ! 1/ Reference point
+          igref=1
+          actualmin=1000.
+          do ig=1,ngridmx
+               val6=dist_pluto(latfi(ig),lonfi(ig),
+     &                              val*pi/180.,val2*pi/180.)
+               if (val6.lt.actualmin) then
+                   actualmin=val6
+                   igref=ig
+               endif 
+          enddo
+
+          do k=1,l 
+
+            write(*,*) 'Choice lat,lon where terrain is copied'
+ 563        read(*,*,iostat=ierr) val4,val5
+            if(ierr.ne.0) goto 563
+          
+            if (val5.gt.180.) then
+              val5=val5-360.
+            endif
+
+            ! 2/ Target point
+            igref2=1
+            actualmin=1000.
+            do ig=1,ngridmx
+               val6=dist_pluto(latfi(ig),lonfi(ig),
+     &                              val4*pi/180.,val5*pi/180.)
+               if (val6.lt.actualmin) then
+                   actualmin=val6
+                   igref2=ig
+               endif 
+            enddo
+
+            ! for each point within A1, get distance and angle
+            ! save in a table
+            i=1
+            array_ind(:)=0
+            array_dist(:)=1000.
+            array_angle(:)=0.
+            do ig=1,ngridmx
+               val6=dist_pluto(latfi(ig),lonfi(ig),
+     &                              val*pi/180.,val2*pi/180.)
+               if (val6.lt.val3) then
+                array_ind(i)=ig
+                array_dist(i)=val6
+                array_angle(i)=atan2(val-latfi(ig)*180./pi,
+     &                              val2-lonfi(ig)*180./pi)
+                i=i+1
+               endif
+            enddo
+
+          ! for each point within A2, get distance and angle 
+          ! and look where fit with previous table, and set value 
+
+            do ig=1,ngridmx
+               val6=dist_pluto(latfi(ig),lonfi(ig),
+     &                              val4*pi/180.,val5*pi/180.)
+               if (val6.lt.val3) then
+                ! dist = val6
+                ! get angle:
+                angle=atan2(val4-latfi(ig)*180./pi,
+     &                              val5-lonfi(ig)*180./pi)
+                ! Loop where min
+                actualmin=1000.
+                do j=1,i
+                   if ( (array_angle(j).lt.angle+0.52).and.
+     &                  (array_angle(j).gt.angle-0.52).and.
+     &                  (array_dist(j)-val6.lt.actualmin) ) then
+                        actualmin=array_dist(j)-val6
+                        igref3=j
+                   endif
+                enddo
+                phisfi(ig)=phisfi(array_ind(igref3))
+                qsurf(ig,igcm_ch4_ice)=
+     &                         qsurf(array_ind(igref3),igcm_ch4_ice)
+                tsurf(ig)=tsurf(array_ind(igref3))
+               endif
+            enddo
+          enddo
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+
+c       delfrostch4/ delete ch4 frost on a latitudinal band
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'delfrostch4') then
+
+          write(*,*) 'removing ch4 latitudinally'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+          read(*,*,iostat=ierr) val,val2
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 522      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 522
+          write(*,*) 'Choice amount max below whcih ch4 is removed'
+ 521      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 521
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5) .and.
+     &              (qsurf(ig,igcm_ch4_ice).lt.val3) )      then
+                    qsurf(ig,igcm_ch4_ice)=0.
+             ENDIF
+          ENDDO
+
+c       modch4 : adding/remove ch4 frost on a latitudinal band
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modch4') then
+
+          write(*,*) 'Adding or removing ch4 latitudinally'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+          read(*,*,iostat=ierr) val,val2
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 422      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 422
+          write(*,*) 'Choice multiplicative factor'
+ 421      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 421
+          write(*,*) 'Choice additional ch4'
+ 423      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 423
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5)) then
+!     &              (qsurf(ig,igcm_n2).gt.50))  then     
+!     &              (qsurf(ig,igcm_ch4_ice).lt.10) )      then
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
+             ENDIF
+          ENDDO
+
+c       modch4n2 : adding/remove ch4 frost on a latitudinal band
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modch4n2') then
+
+          write(*,*) 'Adding or removing ch4 latitudinally'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+          read(*,*,iostat=ierr) val,val2
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+          read(*,*,iostat=ierr) val4,val5
+          write(*,*) 'Choice range n2 for modif'
+          read(*,*,iostat=ierr) val7,val8
+          write(*,*) 'Choice multiplicative factor ch4'
+          read(*,*,iostat=ierr) val3
+          write(*,*) 'Choice additional ch4'
+          read(*,*,iostat=ierr) val6
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5) .and.
+     &              (qsurf(ig,igcm_n2).gt.val7)  .and.     
+     &              (qsurf(ig,igcm_n2).lt.val8) )      then
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
+             ENDIF
+          ENDDO
+
+c       non2noco : if no n2 no co
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'non2noco') then
+          DO ig=1,ngridmx
+             IF (   (qsurf(ig,igcm_n2).lt.0.001))  then     
+                    qsurf(ig,igcm_co_ice)=0.
+             ENDIF
+          ENDDO
+
+c       modco : adding/remove co frost on a latitudinal band
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modco') then
+
+          write(*,*) 'Adding or removing co where N2 is present '
+          write(*,*) 'Choice limit mini n2 pour mettre co?'
+ 460      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 460
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 461      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 461
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 462      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 462
+          write(*,*) 'Choice multiplicative factor'
+ 463      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 463
+          write(*,*) 'Choice additional co'
+ 464      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 464
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5)  .and. 
+     &              (qsurf(ig,igcm_n2).lt.val7))  then     
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
+             ENDIF
+          ENDDO
+
+c       modn2 : adding/remove n2 ice 
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modn2') then
+
+          write(*,*) 'Adding or removing n2 '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 425      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 425
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 426      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 426
+          write(*,*) 'Choice multiplicative factor'
+ 427      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 427
+          write(*,*) 'Choice additional n2'
+ 428     read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 428
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5)  ) then 
+c    &              (qsurf(ig,igcm_n2).lt.50))  then     
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
+             ENDIF
+!             IF ((phisfi(ig).gt.-1000.)) then 
+!                qsurf(ig,igcm_n2)=0.
+!             ENDIF
+          ENDDO
+
+c       modcoch4 : adding/remove co and ch4 frost where co without n2
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modcoch4') then
+
+          write(*,*) 'Adding or removing co where N2 is not present '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 491      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 491
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 492      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 492
+          write(*,*) 'Choice multiplicative factor'
+ 493      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 493
+          write(*,*) 'Choice additional co ch4'
+ 494      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 494
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.le.val5)  .and. 
+     &              (qsurf(ig,igcm_n2).lt.10.))  then     
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
+             ENDIF
+          ENDDO
+
+c       modn2HD : adding/remove n2 ice for HD runs
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modn2HD') then
+
+          write(*,*) 'Adding or removing n2 '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 480      read(*,*,iostat=ierr) val,val1
+          if(ierr.ne.0) goto 480
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 481      read(*,*,iostat=ierr) val2,val3
+          if(ierr.ne.0) goto 481
+          write(*,*) 'Choice amount N2 min max where to modify'
+ 482      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 482
+          write(*,*) 'Choice phisfi min max where change n2'
+ 483     read(*,*,iostat=ierr) val6,val11
+          if(ierr.ne.0) goto 483
+          write(*,*) 'Choice multiplicative factor'
+ 484      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 484
+          write(*,*) 'Choice additional n2'
+ 485     read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 485
+!          write(*,*) 'Choice ind lon equivalent for change tsurf tsoil'
+! 486     read(*,*,iostat=ierr) val9
+!          if(ierr.ne.0) goto 486
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val1)  .and.
+     &              (lonfi(ig)*180./pi.ge.val2)  .and.
+     &              (lonfi(ig)*180./pi.lt.val3)  .and.
+     &              (qsurf(ig,igcm_n2).ge.val4)  .and.
+     &              (qsurf(ig,igcm_n2).lt.val5)  .and.
+     &              (phisfi(ig).ge.val6) .and.
+     &              (phisfi(ig).le.val11)  ) then 
+
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8
+                    !qsurf(ig,igcm_ch4_ice)=0.
+                    qsurf(ig,igcm_co_ice)=0.
+                    ! index of ig
+                    !if (val9.eq.0.) then
+                    !   iref=2546
+                    !else
+                    !   val10=modulo((ig-1),96)
+                    !   choice=ig+int(96/2)-val10
+                    !   !choice=int(1+96*(val9-1)+val10)
+                    !endif
+                    !tsurf(ig)=tsurf(iref)
+                    !DO l=1,nsoilmx
+                    !   tsoil(ig,l)=tsoil(iref,l)
+                    !ENDDO 
+                    tsurf(ig)=tsurf(ig)+4
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(ig,l)+4
+                    ENDDO 
+             ENDIF
+!             IF ((phisfi(ig).gt.-1000.)) then 
+!                qsurf(ig,igcm_n2)=0.
+!             ENDIF
+          ENDDO
+
+c       modn2HD_SP : adding/remove n2 ice for HD runs
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modn2HD_SP') then
+
+          write(*,*) 'Adding or removing n2 '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 501      read(*,*,iostat=ierr) val,val1
+          if(ierr.ne.0) goto 501
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 502      read(*,*,iostat=ierr) val2,val3
+          if(ierr.ne.0) goto 502
+          write(*,*) 'Choice amount N2 min max where to modify'
+ 503      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 503
+          write(*,*) 'Choice phisfi min max where change n2'
+ 504     read(*,*,iostat=ierr) val6,val11
+          if(ierr.ne.0) goto 504
+          write(*,*) 'Choice multiplicative factor'
+ 505      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 505
+          write(*,*) 'Choice additional n2'
+ 506     read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 506
+          write(*,*) 'Choice ind lon equivalent for change tsurf tsoil'
+ 507     read(*,*,iostat=ierr) iref
+          if(ierr.ne.0) goto 507
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val1)  .and.
+     &              (lonfi(ig)*180./pi.ge.val2)  .and.
+     &              (lonfi(ig)*180./pi.lt.val3)  .and.
+     &              (qsurf(ig,igcm_n2).ge.val4)  .and.
+     &              (qsurf(ig,igcm_n2).lt.val5)  .and.
+     &              (phisfi(ig).ge.val6) .and.
+     &              (phisfi(ig).le.val11)  ) then 
+
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8
+                    qsurf(ig,igcm_ch4_ice)=0.
+                    qsurf(ig,igcm_co_ice)=0.
+                    ! index of ig
+                    !if (val9.eq.0.) then
+                    !   iref=2546
+                    !else
+                    !val10=modulo((ig-1),96)
+                    !choice=ig+int(96/2)-val10
+                    !choice=int(1+96*(val9-1)+val10)
+                    if (iref.ge.1) then
+                     tsurf(ig)=tsurf(iref)
+                     DO l=1,nsoilmx
+                        tsoil(ig,l)=tsoil(iref,l)
+                     ENDDO 
+                    else if (iref.eq.0) then
+                     choice=int(ig/96.)*96+84
+                     print*, 'choice=',choice
+                     tsurf(ig)=tsurf(int(choice))
+                     DO l=1,nsoilmx
+                        tsoil(ig,l)=tsoil(int(choice),l)
+                     ENDDO 
+                    endif
+             ENDIF
+          ENDDO
+
+
+c       modn2_topo : adding/remove n2 ice 
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then
+
+          write(*,*) 'Adding or removing n2 '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 441      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 441
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 442      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 442
+          write(*,*) 'Choice topo 1 et 2 (phi) where change is active'
+ 443      read(*,*,iostat=ierr) val7,val8
+          if(ierr.ne.0) goto 443
+          write(*,*) 'Choice multiplicative factor'
+ 444      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 444
+          write(*,*) 'Choice additional n2'
+ 445     read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 445
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val4)  .and.
+     &              (lonfi(ig)*180./pi.lt.val5)  .and.
+     &              (phisfi(ig).le.val8) .and.
+     &              (phisfi(ig).ge.val7)   ) then
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
+             ENDIF
+          ENDDO
+
+c       modwheren2 : adding/remove n2 ice where already n2
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modwheren2') then
+
+          write(*,*) 'Choice multiplicative factor'
+ 471      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 471
+          write(*,*) 'Choice additional n2'
+ 472     read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 472
+
+          DO ig=1,ngridmx
+             IF (   qsurf(ig,igcm_n2).gt.0.001 ) then
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
+             ENDIF
+          ENDDO
+
+c       checktsoil : changing tsoil if no n2
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then
+
+          write(*,*) 'selecting area where tsoil to be check '
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 511      read(*,*,iostat=ierr) val,val1
+          if(ierr.ne.0) goto 511
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 512      read(*,*,iostat=ierr) val2,val3
+          if(ierr.ne.0) goto 512
+          write(*,*) 'Choice temp min max in which change is made'
+ 513      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 513
+          write(*,*) 'Choice phisfi min max where change n2'
+ 514     read(*,*,iostat=ierr) val6,val11
+          if(ierr.ne.0) goto 514
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val1)  .and.
+     &              (lonfi(ig)*180./pi.ge.val2)  .and.
+     &              (lonfi(ig)*180./pi.le.val3)  .and.
+     &              (((tsurf(ig).ge.val4)  .and.
+     &              (tsurf(ig).le.val5)) .or.
+     &              ((tsoil(ig,nsoilmx).ge.val4)  .and.
+     &              (tsoil(ig,nsoilmx).le.val5))) .and.
+     &              (phisfi(ig).ge.val6) .and.
+     &              (phisfi(ig).le.val11) .and.
+     &              (qsurf(ig,igcm_n2).le.0.001) ) then 
+
+!                    DO i=1,ngridmx
+!                       IF (   (latfi(i).eq.latfi(ig))  .and.
+!     &              (lonfi(i).eq.0.) ) then
+!
+                         tsurf(ig)=50. 
+                         !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice)
+!
+                         DO l=1,nsoilmx
+                           tsoil(ig,l)=50. !tsoil(i,l)
+                         ENDDO 
+                       !ENDIF
+                    !ENDDO
+             ENDIF
+          ENDDO
+
+c       asymtriton : changing ice, tsurf and tsoil on triton
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'asymtriton') then
+
+          write(*,*) 'selecting area where tsoil to be check '
+          write(*,*) 'Choice center latitude of sinuisoid distrib?'
+ 531      read(*,*,iostat=ierr) val1
+          if(ierr.ne.0) goto 531
+          write(*,*) 'Choice maginutde in latitude?'
+ 532      read(*,*,iostat=ierr) val2
+          if(ierr.ne.0) goto 532
+          write(*,*) 'Choice center longitude?'
+ 533      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 533
+!          write(*,*) 'iip1,jjp1=',iip1,jjp1,ngridmx
+
+          DO ig=1,ngridmx
+             ! Latitude of the sinusoid:
+             val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.)
+             ! If above line and ice: remove ice
+             IF ( (latfi(ig)*180./pi.ge.val11) .and.
+     &             (latfi(ig)*180./pi.le.val1+val2) .and.
+     &             (qsurf(ig,igcm_n2).gt.0.) ) then
+               ! Look for same longitude point where no ice           
+               val5=1.
+               jref=ig
+               ! 1
+               ! ... iip1 ... x (jjp1-2)     32 x 23
+               ! 1
+               do while (val5.ge.1..and.jref.gt.iip1)
+                  ! northward point
+                  jref=jref-iip1+1
+                  ! ice at the point
+                  val5=qsurf(jref,igcm_n2)
+!                 write(*,*) jref,
+!     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5  
+               enddo  
+               if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE'  
+
+               ! Copy point in the selected area
+               tsurf(ig)=tsurf(jref) 
+               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
+               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+               qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               DO l=1,nsoilmx
+                  tsoil(ig,l)=tsoil(jref,l)
+               ENDDO 
+             ENDIF
+             ! If below line and no ice: add ice
+             IF ( (latfi(ig)*180./pi.le.val11) .and.
+     &             (latfi(ig)*180./pi.le.val1+val2) .and.
+     &             (qsurf(ig,igcm_n2).eq.0.) ) then
+               ! Look for same longitude point where ice           
+               val5=1.
+               jref=ig
+               do while (val5.le.1.and.jref.lt.ngridmx-iip1)
+                  ! southward point
+                  jref=jref+iip1-1
+                  ! ice at the point
+                  val5=qsurf(jref,igcm_n2)
+                  write(*,*) jref,
+     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5  
+               enddo  
+               if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE'  
+
+               ! Copy point in the selected area
+               tsurf(ig)=tsurf(jref) 
+               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
+               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+               qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               DO l=1,nsoilmx
+                  tsoil(ig,l)=tsoil(jref,l)
+               ENDDO 
+             ENDIF
+
+          ENDDO
+
+c       source_ch4 : adding source ch4
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then
+
+          write(*,*) 'Adding ch4 ice latitudinally ' 
+          write(*,*) 'Choice : lat1 and lat2 ?'
+ 433      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 433
+          write(*,*) 'Choice : lon1 and lon2 ?'
+ 434      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 434
+          write(*,*) 'Choice amount (kg/m2)'
+ 435      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 435
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val3)  .and.
+     &              (lonfi(ig)*180./pi.lt.val4)  ) then 
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
+             ENDIF
+          ENDDO
+
+c       source_co : adding source co
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'source_co') then
+
+          write(*,*) 'Adding co ice latitudinally ' 
+          write(*,*) 'Choice : lat1 and lat2 ?'
+ 436      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 436
+          write(*,*) 'Choice : lon1 and lon2 ?'
+ 437      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 437
+          write(*,*) 'Choice amount (kg/m2)'
+ 438      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 438
+
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val3)  .and.
+     &              (lonfi(ig)*180./pi.lt.val4)  ) then 
+                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5
+             ENDIF
+          ENDDO
+
+!       therm_ini_s: (re)-set soil thermal inertia to reference surface values
+!       ----------------------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then
+!          write(*,*)"surfithfi(1):",surfithfi(1)
+          do isoil=1,nsoilmx
+             inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
+          enddo
+          write(*,*)'OK: Soil thermal inertia has been reset to referenc
+     &e surface values'
+          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
+          ithfi(:,:)=inertiedat(:,:)
+         ! recast ithfi() onto ith()
+         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+!       inert3d: set soil thermal inertia to specific uniform value 
+!       ----------------------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'inert3d') then
+          write(*,*) 'Actual value of surf Thermal inertia at ig=1: ',
+     &                    inertiedat(1,1), ' SI'
+          write(*,*) 'Value of Thermal inertia (SI) ? '
+          read(*,*) val
+          do isoil=1,nsoilmx
+            do ig=1,ngridmx
+              inertiedat(ig,isoil)=val
+            enddo
+          enddo
+          write(*,*)'OK: Soil TI set to ',val,' SI everywhere'
+          ithfi(:,:)=inertiedat(:,:)
+         ! recast ithfi() onto ith()
+          call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+!       subsoilice_n: Put deep ice layer in northern hemisphere soil
+!       ------------------------------------------------------------
+       else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
+
+         write(*,*)'From which latitude (in deg.), up to the north pole,
+     &should we put subterranean ice?'
+         ierr=1
+         do while (ierr.ne.0)
+          read(*,*,iostat=ierr) val
+          if (ierr.eq.0) then ! got a value
+            ! do a sanity check
+            if((val.lt.0.).or.(val.gt.90)) then
+               write(*,*)'Latitude should be between 0 and 90 deg. !!!'
                ierr=1
-             endif
-             if(val2.lt.layer(1)) then
-              if(val2.eq.0) then
-               write(*,*)'Thermal inertia set for all subsurface layers'
-               ierr=0
-              else
-               write(*,*)'Depth should be more than ',layer(1)
-               ierr=1
-              endif
-             endif
-           endif
-          enddo ! of do while
-
-          ! find the reference index iref the depth corresponds to
-          if(val2.eq.0) then
-           iref=1
-           write(*,*)'Level selected is first level: ',layer(iref),' m'
-          else
-           do isoil=1,nsoilmx-1
-            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
-     &       then
-             iref=isoil+1
-             write(*,*)'Level selected : ',layer(isoil+1),' m'
-             exit
-            endif
-           enddo
-          endif
-
-          DO ig=1,ngridmx
-             DO j=iref,nsoilmx
-                   ithfi(ig,j)=ith_bb
-             ENDDO
-          ENDDO
-
-          CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
-
-c       diurnal_TI : choice of thermal inertia values (global)
-c       ----------------------------------------------------------------
-        else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then
-
-         write(*,*) 'New value for diurnal thermal inertia  ?'
- 106     read(*,*,iostat=ierr) ith_bb
-         if(ierr.ne.0) goto 106
-         write(*,*) 'Diurnal thermal inertia (new value):',ith_bb
-
-         write(*,*)'At which depth (in m.) does the ice layer ends?'
-         write(*,*)'(currently, the soil layer 1 and nsoil are:'
-     &              ,layer(1),' - ',layer(nsoilmx),')'
+            else ! find corresponding jref (nearest latitude)
+              ! note: rlatu(:) contains decreasing values of latitude
+              !       starting from PI/2 to -PI/2
+              do j=1,jjp1
+               if ((rlatu(j)*180./pi.ge.val).and.
+     &              (rlatu(j+1)*180./pi.le.val)) then
+                  ! find which grid point is nearest to val:
+                  if (abs(rlatu(j)*180./pi-val).le.
+     &                abs((rlatu(j+1)*180./pi-val))) then
+                    jref=j
+                  else
+                    jref=j+1
+                  endif
+                 
+                  write(*,*)'Will use nearest grid latitude which is:',
+     &                     rlatu(jref)*180./pi
+               endif
+              enddo ! of do j=1,jjp1
+            endif ! of if((val.lt.0.).or.(val.gt.90))
+          endif !of if (ierr.eq.0)
+         enddo ! of do while
+
+         ! Build layers() (as in soil_settings.F)
+         val2=sqrt(mlayer(0)*mlayer(1))
+         val3=mlayer(1)/mlayer(0)
+         do isoil=1,nsoilmx
+           layer(isoil)=val2*(val3**(isoil-1))
+         enddo
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, the deepest soil layer extends down to:'
+     &              ,layer(nsoilmx),')'
          ierr=1
          do while (ierr.ne.0)
           read(*,*,iostat=ierr) val2
-          write(*,*)'val2 in diurnal_TI:',val2,'ierr=',ierr
+!	  write(*,*)'val2:',val2,'ierr=',ierr
           if (ierr.eq.0) then ! got a value, but do a sanity check
             if(val2.gt.layer(nsoilmx)) then
@@ -911,47 +1900,323 @@
           endif
          enddo ! of do while
-
-           ! find the reference index iref the depth corresponds to
+ 
+         ! find the reference index iref the depth corresponds to
+!	 if (val2.lt.layer(1)) then
+!	  iref=1
+!	 else
          do isoil=1,nsoilmx-1
-            !write(*,*)'isoil, ',isoil,val2
-            !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
-            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
+           if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
      &       then
-             iref=isoil+1
-             write(*,*)'Level selected : ',layer(isoil+1),' m'
+             iref=isoil
              exit
+           endif
+         enddo
+
+         ! thermal inertia of the ice:
+         ierr=1
+         do while (ierr.ne.0)
+          write(*,*)'What is the value of subterranean ice thermal inert
+     &ia? (e.g.: 2000)'
+          read(*,*,iostat=ierr)iceith
+         enddo ! of do while
+
+         ! recast ithfi() onto ith()
+         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+     
+         do j=1,jref
+!	    write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
+            do i=1,iip1 ! loop on longitudes
+            ! Build "equivalent" thermal inertia for the mixed layer
+             ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
+     &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
+     &                      ((layer(iref+1)-val2)/(iceith)**2)))
+             ! Set thermal inertia of lower layers
+             do isoil=iref+2,nsoilmx
+              ith(i,j,isoil)=iceith ! ice
+             enddo
+            enddo ! of do i=1,iip1 
+         enddo ! of do j=1,jjp1
+ 
+
+         CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+
+!       subsoilice_s: Put deep ice layer in southern hemisphere soil
+!       ------------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
+
+         write(*,*)'From which latitude (in deg.), down to the south pol
+     &e, should we put subterranean ice?'
+         ierr=1
+         do while (ierr.ne.0)
+          read(*,*,iostat=ierr) val
+          if (ierr.eq.0) then ! got a value
+            ! do a sanity check
+            if((val.gt.0.).or.(val.lt.-90)) then
+              write(*,*)'Latitude should be between 0 and -90 deg. !!!'
+              ierr=1
+            else ! find corresponding jref (nearest latitude)
+              ! note: rlatu(:) contains decreasing values of latitude
+              !       starting from PI/2 to -PI/2
+              do j=1,jjp1
+                if ((rlatu(j)*180./pi.ge.val).and.
+     &              (rlatu(j+1)*180./pi.le.val)) then
+                ! find which grid point is nearest to val:
+                  if (abs(rlatu(j)*180./pi-val).le.
+     &                abs((rlatu(j+1)*180./pi-val))) then
+                   jref=j
+                  else
+                   jref=j+1
+                  endif
+ 
+                 write(*,*)'Will use nearest grid latitude which is:',
+     &                     rlatu(jref)*180./pi
+                endif
+              enddo ! of do j=1,jjp1
+            endif ! of if((val.lt.0.).or.(val.gt.90))
+          endif !of if (ierr.eq.0)
+         enddo ! of do while
+
+         ! Build layers() (as in soil_settings.F)
+         val2=sqrt(mlayer(0)*mlayer(1))
+         val3=mlayer(1)/mlayer(0)
+         do isoil=1,nsoilmx
+           layer(isoil)=val2*(val3**(isoil-1))
+         enddo
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, the deepest soil layer extends down to:'
+     &              ,layer(nsoilmx),')'
+         ierr=1
+         do while (ierr.ne.0)
+          read(*,*,iostat=ierr) val2
+!	  write(*,*)'val2:',val2,'ierr=',ierr
+          if (ierr.eq.0) then ! got a value, but do a sanity check
+            if(val2.gt.layer(nsoilmx)) then
+              write(*,*)'Depth should be less than ',layer(nsoilmx)
+              ierr=1
             endif
+            if(val2.lt.layer(1)) then
+              write(*,*)'Depth should be more than ',layer(1)
+              ierr=1
+            endif
+          endif
+         enddo ! of do while
+ 
+         ! find the reference index iref the depth corresponds to
+          do isoil=1,nsoilmx-1
+           if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
+     &       then
+             iref=isoil
+             exit
+           endif
+          enddo
+ 
+         ! thermal inertia of the ice:
+         ierr=1
+         do while (ierr.ne.0)
+          write(*,*)'What is the value of subterranean ice thermal inert
+     &ia? (e.g.: 2000)'
+          read(*,*,iostat=ierr)iceith
+         enddo ! of do while
+ 
+         ! recast ithfi() onto ith()
+         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+ 
+         do j=jref,jjp1
+!	    write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
+            do i=1,iip1 ! loop on longitudes
+            ! Build "equivalent" thermal inertia for the mixed layer
+             ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
+     &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
+     &                      ((layer(iref+1)-val2)/(iceith)**2)))
+             ! Set thermal inertia of lower layers
+             do isoil=iref+2,nsoilmx
+              ith(i,j,isoil)=iceith ! ice
+             enddo
+            enddo ! of do i=1,iip1 
+         enddo ! of do j=jref,jjp1
+
+         CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+
+c       ----------------------------------------------------------
+c       'reservoir_SP' : increase or decrase ices reservoir in SP by factor
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then
+          write(*,*) "Increase/Decrease N2 reservoir by factor :"
+          read(*,*) val7
+          write(*,*) "Increase/Decrease CH4 reservoir by factor :"
+          read(*,*) val8
+          write(*,*) "Increase/Decrease CO reservoir by factor :"
+          read(*,*) val9
+
+          ! Definition SP:
+         val3=-33.   !lat1
+         val4=50.    !lat2
+         val5=140.   ! lon1
+         val6=-155.  ! lon2
+         choice=-50. ! min phisfi
+         write(*,*) 'def SP :'
+         write(*,*) 'lat :',val3,val4
+         write(*,*) 'lon :',val5,'180 / -180',val6
+         write(*,*) 'choice phisfi min :',choice
+
+         DO ig=1,ngridmx
+
+           IF (   (latfi(ig)*180./pi.ge.val3)  .and.
+     &            (latfi(ig)*180./pi.le.val4)  .and.
+     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
+     &                       (lonfi(ig)*180./pi.le.val6))  .or.
+     &         ((lonfi(ig)*180./pi.ge.val5)  .and. 
+     &                       (lonfi(ig)*180./pi.le.180.))) ) then
+
+                IF ((phisfi(ig).le.choice)) then  !.and.
+c    &                              (qsurf(ig,igcm_n2).ge.50)) then
+
+                  qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
+                  qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8
+                  qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9
+                ENDIF
+           ENDIF
+          ENDDO
+
+c       ----------------------------------------------------------
+c       'reservoir' : increase or decrase ices reservoir by factor
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'reservoir') then
+          write(*,*) "Increase/Decrease N2 reservoir by factor :"
+          read(*,*) val
+          write(*,*) "Increase/Decrease CH4 reservoir by factor :"
+          read(*,*) val2
+          write(*,*) "Increase/Decrease CO reservoir by factor :"
+          read(*,*) val3
+
+          DO ig=1,ngridmx
+           qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val
+           qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2
+           qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
+          ENDDO
+
+c       --------------------------------------------------------
+c       'plutomap' : initialize pluto ices distribution
+c       --------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'plutomap') then
+
+         write(*,*) 'pluto_map.dat has to be in your simulation file !!'
+         open(27,file='pluto_map.dat')
+         N2_pluto_dat(:,:) =0.
+         CH4_pluto_dat(:,:) =0.
+         CO_pluto_dat(:,:) =0.
+
+         ! Dimension file pluto_map
+         IF (jm_plu.ne.180) then
+             write(*,*) 'Newstart : you should set jm_plu to 180'
+             stop
+         ENDIF 
+         ! Read values
+         do j=1,jm_plu+1
+           read(27,*) (map_pluto_dat(i,j),i=1,im_plu)
+           do i=1,im_plu
+
+               !!!!!! BTD_v2
+               if (map_pluto_dat(i,j).eq.3) then
+                  N2_pluto_dat(i,j)=100000.
+               elseif (map_pluto_dat(i,j).eq.4) then
+                  CH4_pluto_dat(i,j)=100000.
+               elseif (map_pluto_dat(i,j).eq.0) then
+                  alb_dat(i,j)=0.3
+               elseif (map_pluto_dat(i,j).eq.6) then
+                  alb_dat(i,j)=0.6
+               elseif (map_pluto_dat(i,j).eq.7) then
+                  alb_dat(i,j)=0.1
+               endif
+
+               !!!!!! BTD
+               !if (map_pluto_dat(i,j).eq.3) then
+               !   CH4_pluto_dat(i,j)=100000.
+               !endif
+           end do
+         end do
+         close (27)
+         ! Interpolate
+
+         !! latitude and longitude in REindexed pluto map  :
+         latv_plu(1)=90. *pi/180.
+         do j=2,jm_plu
+          latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180.
+         end do
+         latv_plu(jm_plu+1)= -90. *pi/180.
+
+         do i=1,im_plu+1
+          lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180.
          enddo
 
+         !Reindexation to shift the longitude grid like a LMD GCM grid...
+         do j=1,jm_plu+1
+            N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+
+     &                           N2_pluto_dat(im_plu,j))/2
+            CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+
+     &                          CH4_pluto_dat(im_plu,j))/2
+            CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+
+     &                           CO_pluto_dat(im_plu,j))/2
+            alb_rein(1,j)=(alb_dat(1,j)+
+     &                           alb_dat(im_plu,j))/2
+            do i=2,im_plu
+               N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+
+     &                           N2_pluto_dat(i,j))/2
+               CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+
+     &                           CH4_pluto_dat(i,j))/2
+               CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+
+     &                           CO_pluto_dat(i,j))/2
+               alb_rein(i,j)=(alb_dat(i-1,j)+
+     &                           alb_dat(i,j))/2
+            end do
+            N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j)
+            CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j)
+            CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j)
+            alb_rein(im_plu+1,j) = alb_rein(1,j)
+         end do
+
+         call interp_horiz(N2_pluto_rein,N2_pluto_gcm,
+     & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
+         call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm,
+     & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
+         call interp_horiz(CO_pluto_rein,CO_pluto_gcm,
+     & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
+         call interp_horiz(alb_rein,alb_gcm,
+     & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
+
+         ! Switch dyn grid to fi grid
+         qsurf_tmp(:,:) =0.
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
+     &         N2_pluto_gcm,qsurf_tmp(1,igcm_n2))
+         if (igcm_ch4_ice.ne.0) then
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
+     &         CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice))
+         endif
+         if (igcm_co_ice.ne.0) then
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
+     &         CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice))
+         endif
+         alb=alb_gcm
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi)
+         !print*, 'N2dat=',N2_pluto_gcm 
+         !print*, 'COdat=',CO_pluto_gcm 
+         !print*, 'CH4dat=',CH4_pluto_gcm 
+         print*, 'N2dat=',qsurf_tmp(:,igcm_n2)
+         print*, 'COdat=',qsurf_tmp(:,igcm_co_ice)
+         print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice)
+
+         ! Fill
          DO ig=1,ngridmx
-             DO j=1,iref
-                   ithfi(ig,j)=ith_bb
-             ENDDO
-         ENDDO
-
-         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+           qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2)
+           if (igcm_ch4_ice.ne.0) then
+            qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+
+     &                               qsurf_tmp(ig,igcm_ch4_ice)
+           endif
+           if (igcm_co_ice.ne.0) then
+            qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+
+     &                               qsurf_tmp(ig,igcm_co_ice)
+           endif
+         ENDDO        
 
         endif
@@ -961,33 +2226,7 @@
  999  continue
 
- 
-c=======================================================================
-c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
-c=======================================================================
-!      DO ig=1, ngridmx
-!         totalfrac(ig)=0.5
-!         DO l=1,llm
-!            cloudfrac(ig,l)=0.5
-!         ENDDO
-! Ehouarn, march 2012: also add some initialisation for hice
-!         hice(ig)=0.0
-!      ENDDO
-      
-c========================================================
-
-!      DO ig=1,ngridmx
-!         IF(tsurf(ig) .LT. 273.16-1.8) THEN
-!            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
-!            hice(ig)=min(hice(ig),1.0)
-!         ENDIF
-!      ENDDO
-
-
-
-
 c=======================================================================
 c   Correct pressure on the new grid (menu 0)
 c=======================================================================
-
 
       if ((choix_1.eq.0).and.(.not.flagps0)) then
@@ -1008,12 +2247,7 @@
       end if
 
-         
-c=======================================================================
-c=======================================================================
-
 c=======================================================================
 c    Initialisation de la physique / ecriture de newstartfi :
 c=======================================================================
-
 
       CALL inifilr 
@@ -1111,2 +2345,20 @@
       end
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      function dist_pluto(lat1,lon1,lat2,lon2)
+      implicit none
+      real dist_pluto
+      real lat1,lon1,lat2,lon2
+      real dlon, dlat
+      real a,c
+      real rad
+      rad=1190. ! km
+
+      dlon = lon2 - lon1
+      dlat = lat2 - lat1
+      a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
+      c = 2 * atan2( sqrt(a), sqrt(1-a) )
+      dist_pluto = rad * c
+      return
+      end
+
