Index: /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/datareadnc.F
===================================================================
--- /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/datareadnc.F	(revision 3381)
+++ /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/datareadnc.F	(revision 3382)
@@ -58,5 +58,5 @@
 
       INTEGER    imd,jmd,imdp1,jmdp1
-      parameter    (imd=360,jmd=179,imdp1=361,jmdp1=180)
+      parameter    (imd=1440,jmd=719,imdp1=1441,jmdp1=720)
 
       INTEGER    iimp1
@@ -85,5 +85,5 @@
 
       INTEGER klatdat,ngridmixgdat
-      PARAMETER (klatdat=180,ngridmixgdat=360)
+      PARAMETER (klatdat=720,ngridmixgdat=1440)
 
 c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
@@ -221,5 +221,5 @@
           latitude(:)=(pi/180.)*latitude(:)
 
-          call grid_noro1(360, 180, longitude, latitude, zdata,
+          call grid_noro1(1440,720, longitude, latitude, zdata,
      .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
 
Index: /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/grid_noro1.F
===================================================================
--- /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/grid_noro1.F	(revision 3381)
+++ /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/grid_noro1.F	(revision 3382)
@@ -56,5 +56,5 @@
        implicit none
        integer iusn,jusn,iext
-       parameter(iusn=360,jusn=180,iext=40)
+       parameter(iusn=1440,jusn=720,iext=40)
 c!-*-      include 'param1'
 c!-*-      include 'comcstfi.h'
Index: /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F
===================================================================
--- /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F	(revision 3381)
+++ /trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F	(revision 3382)
@@ -64,5 +64,6 @@
 c Variables pour les lectures NetCDF des fichiers "start_archive" 
 c--------------------------------------------------
-      INTEGER nid_dyn, nid_fi,nid,nvarid
+      INTEGER nid_dyn, nid_fi,nid,nvarid,nvarid_input,nvarid_inputs
+      INTEGER nid_fi_input
       INTEGER length
       parameter (length = 100)
@@ -108,4 +109,6 @@
       real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
       REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
+      REAL field_input(ngridmx)
+      REAL,ALLOCATABLE :: field_inputs(:,:)
 
       INTEGER i,j,l,isoil,ig,idum,k
@@ -137,4 +140,5 @@
 c-------------------
       real choix_1,pp,choice
+      integer randtab(ngridmx)
       character*80      fichnom
       character*250  filestring
@@ -143,5 +147,7 @@
       real z_reel(iip1,jjp1)
       real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
-      real ptoto,pcap,patm,airetot,ptotn,patmn,psea
+      real tsurf_bb,tsurf_bb2
+      real ptoto,pcap,patm,airetot,ptotn,patmn,psea,geop
+      real tempsoil(22),levsoil(22)
       character*1 yes
       logical :: flagtset=.false. ,  flagps0=.false.
@@ -230,4 +236,5 @@
       
       allocate(tsoil(ngridmx,nsoilmx))
+      allocate(field_inputs(ngridmx,nsoilmx))
       allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
       
@@ -741,4 +748,62 @@
           enddo
 
+c       ptot : Modification of the total pressure: ice + current atmosphere 
+c       -------------------------------------------------------------------
+        else if (modif(1:len_trim(modif)).eq.'ptot') then
+
+c         calcul de la pression totale glace + atm actuelle
+          patm=0.
+          airetot=0.
+          pcap=0.
+          DO j=1,jjp1
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+                patm = patm + ps(i,j)*aire(i,j)
+                airetot= airetot + aire(i,j)
+             ENDDO
+          ENDDO
+          ptoto = pcap + patm
+
+          print*,'Current total pressure at surface ',
+     &       ptoto/airetot
+
+          print*,'new value?'
+          read(*,*) ptotn
+          ptotn=ptotn*airetot
+          patmn=ptotn-pcap
+          print*,'ptoto,patm,ptotn,patmn'
+          print*,ptoto,patm,ptotn,patmn
+          print*,'Mult. factor for pressure (atm only)', patmn/patm
+          do j=1,jjp1
+             do i=1,iip1
+                ps(i,j)=ps(i,j)*patmn/patm
+             enddo
+          enddo
+
+c        Correction pour la conservation des traceurs
+         yes=' '
+         do while ((yes.ne.'y').and.(yes.ne.'n'))
+            write(*,*) 'Do you wish to conserve tracer total mass (y)',
+     &              ' or tracer mixing ratio (n) ?'
+             read(*,fmt='(a)') yes
+         end do
+
+         if (yes.eq.'y') then
+           write(*,*) 'OK : conservation of tracer total mass'
+           DO iq =1, nqtot
+             DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
+                  ENDDO
+               ENDDO
+             ENDDO
+           ENDDO
+          else
+            write(*,*) 'OK : conservation of tracer mixing ratio'
+          end if
+
 c       qname : change tracer name
 c       --------------------------
@@ -1410,5 +1475,5 @@
           if(ierr.ne.0) goto 427
           write(*,*) 'Choice additional n2'
- 428     read(*,*,iostat=ierr) val6
+ 428      read(*,*,iostat=ierr) val6
           if(ierr.ne.0) goto 428
 
@@ -1536,5 +1601,5 @@
           if(ierr.ne.0) goto 503
           write(*,*) 'Choice phisfi min max where change n2'
- 504     read(*,*,iostat=ierr) val6,val11
+ 504      read(*,*,iostat=ierr) val6,val11
           if(ierr.ne.0) goto 504
           write(*,*) 'Choice multiplicative factor'
@@ -1542,8 +1607,8 @@
           if(ierr.ne.0) goto 505
           write(*,*) 'Choice additional n2'
- 506     read(*,*,iostat=ierr) val8
+ 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
+ 507      read(*,*,iostat=ierr) iref
           if(ierr.ne.0) goto 507
 
@@ -1585,5 +1650,4 @@
           ENDDO
 
-
 c       modn2_topo : adding/remove n2 ice 
 c       ----------------------------------------------------------
@@ -1604,5 +1668,5 @@
           if(ierr.ne.0) goto 444
           write(*,*) 'Choice additional n2'
- 445     read(*,*,iostat=ierr) val6
+ 445      read(*,*,iostat=ierr) val6
           if(ierr.ne.0) goto 445
 
@@ -1627,5 +1691,5 @@
           if(ierr.ne.0) goto 471
           write(*,*) 'Choice additional n2'
- 472     read(*,*,iostat=ierr) val6
+ 472      read(*,*,iostat=ierr) val6
           if(ierr.ne.0) goto 472
 
@@ -1652,5 +1716,5 @@
           if(ierr.ne.0) goto 513
           write(*,*) 'Choice phisfi min max where change n2'
- 514     read(*,*,iostat=ierr) val6,val11
+ 514      read(*,*,iostat=ierr) val6,val11
           if(ierr.ne.0) goto 514
 
@@ -1816,5 +1880,5 @@
           enddo
           write(*,*)'OK: Soil thermal inertia has been reset to referenc
-     &e surface values'
+     &               e surface values'
           write(*,*)"inertiedat(1,1):",inertiedat(1,1)
           ithfi(:,:)=inertiedat(:,:)
@@ -1841,8 +1905,8 @@
 !       subsoilice_n: Put deep ice layer in northern hemisphere soil
 !       ------------------------------------------------------------
-       else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
+        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?'
+     &       should we put subterranean ice?'
          ierr=1
          do while (ierr.ne.0)
@@ -1917,5 +1981,5 @@
          do while (ierr.ne.0)
           write(*,*)'What is the value of subterranean ice thermal inert
-     &ia? (e.g.: 2000)'
+     &      ia? (e.g.: 2000)'
           read(*,*,iostat=ierr)iceith
          enddo ! of do while
@@ -1938,5 +2002,4 @@
          enddo ! of do j=1,jjp1
  
-
          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
 
@@ -1946,5 +2009,5 @@
 
          write(*,*)'From which latitude (in deg.), down to the south pol
-     &e, should we put subterranean ice?'
+     &     e, should we put subterranean ice?'
          ierr=1
          do while (ierr.ne.0)
@@ -2011,10 +2074,10 @@
            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)'
+     &     ia? (e.g.: 2000)'
           read(*,*,iostat=ierr)iceith
          enddo ! of do while
@@ -2051,9 +2114,88 @@
 
           ! Definition SP:
-         val3=-33.   !lat1
-         val4=50.    !lat2
-         val5=140.   ! lon1
-         val6=-155.  ! lon2
-         choice=-50. ! min phisfi
+          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       subsoil_SP : choice of thermal inertia values for SP 
+c       ----------------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'subsoil_SP') then
+
+         write(*,*) 'New value for subsoil thermal inertia in SP ?'
+ 510     read(*,*,iostat=ierr) ith_bb
+         if(ierr.ne.0) goto 510
+         write(*,*) 'thermal inertia (new value):',ith_bb
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, 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:',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
+              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
+             write(*,*)'Level selected : ',layer(isoil),' m'
+             exit
+            endif
+           enddo
+         endif 
+
+         ! Definition SP:
+         val3=-50.   !lat1
+         val4=60.    !lat2
+         val5=130.   ! lon1
+         val6=-140.  ! lon2
+         choice=-100. ! min phisfi
          write(*,*) 'def SP :'
          write(*,*) 'lat :',val3,val4
@@ -2070,12 +2212,1560 @@
      &                       (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
+                IF ((phisfi(ig).le.choice) .and.
+     &                              (qsurf(ig,igcm_n2).ge.1000)) then
+                    DO j=iref,nsoilmx
+                       ithfi(ig,j)=ith_bb
+                    ENDDO
+                ENDIF 
+           ENDIF
+         ENDDO
+
+         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+c       topo_sp : change topo of SP
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'topo_sp') then
+
+         ! def SP:
+         val2=-33.   !lat1
+         val3=50.    !lat2
+         val4=140.   ! lon1
+         val5=-155.  ! lon2
+         write(*,*) 'def SP :'
+         write(*,*) 'lat :',val2,val3
+         write(*,*) 'lon :',val4,'180 / -180',val5
+         write(*,*) 'choice phisfi min (ex: -100):'
+606      read(*,*,iostat=ierr) val6
+         if(ierr.ne.0) goto 606
+         write(*,*) 'choice coefficient for phisfi (ex: 2):'
+605      read(*,*,iostat=ierr) choice
+         if(ierr.ne.0) goto 605
+
+         write(*,*) 're scale the pressure :'
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=40.
+         write(*,*) 'r, sale height temperature =',r,temp
+                 
+         do j=1,jjp1
+           do i=1,iip1
+              phisprev= phis(i,j)
+              IF (   (rlatu(j)*180./pi.ge.val2)  .and.
+     &            (rlatu(j)*180./pi.le.val3)  .and.
+     &      (  ((rlonv(i)*180./pi.ge.-180.)  .and. 
+     &                       (rlonv(i)*180./pi.le.val5))  .or.
+     &         ((rlonv(i)*180./pi.ge.val4)  .and. 
+     &                       (rlonv(i)*180./pi.le.180.))) ) then
+
+                IF (phis(i,j).le.val6) then
+                    phis(i,j)=phis(i,j)*choice
+                    ps(i,j) = ps(i,j) *
+     .              exp((phisprev-phis(i,j))/(temp * r))
+                ENDIF 
+              ENDIF  
+           end do
+         end do
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       fill_sp inv: fill inverted SP with N2 ice and adjust phisfi (flat SP)
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'fill_sp_inv') then
+
+         ! def SP:
+         val2=-33.   !lat1
+         val3=50.    !lat2
+         val4=-40.   ! lon1
+         val5=25.  ! lon2
+         write(*,*) 'def SP :'
+         write(*,*) 'lat :',val2,val3
+         write(*,*) 'lon :',val4,val5
+         write(*,*) 'choice phisfi ref SP (ex: -800):'
+706      read(*,*,iostat=ierr) val6
+         if(ierr.ne.0) goto 706
+
+         write(*,*) 're scale the pressure :'
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=40.
+         write(*,*) 'r, sale height temperature =',r,temp,g
+         qsurf=0.        
+
+         write(*,*) 'latfi=',latfi
+         !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+         write(*,*) 'phis=',phis
+         write(*,*) 'phisfi=',phisfi
+         !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         phisold = phis
+         write(*,*) 'phisold=',phisold
+         do ig=1,ngridmx
+               
+              write(*,*) 'lat=',latfi(ig)*180./pi
+              write(*,*) 'lon=',lonfi(ig)*180./pi
+              write(*,*) 'phisfi=',phisfi(ig)
+              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
+     &            (latfi(ig)*180./pi.le.val3)  .and.
+     &           (lonfi(ig)*180./pi.ge.val4)  .and. 
+     &           (lonfi(ig)*180./pi.le.val5) ) then
+                write(*,*) 'hello SP',phisfi(ig),ig
+                IF (phisfi(ig).le.val6) then
+                    qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
+                    phisfi(ig)=val6
+                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
                 ENDIF
+              ENDIF  
+         end do
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+
+c       adjust phisfi according to the amount of N2 ice
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'adjustphi') then
+
+         phisold = phis
+         do ig=1,ngridmx
+            phisfi(ig)=phisfi(ig)+qsurf(ig,igcm_n2)*g/1000.
+         end do
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=37.
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c        periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+
+c       correct phisfi 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'correctphi') then
+
+         write(*,*) 'choice latitudes:'
+537      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 537
+         write(*,*) 'choice longitudes:'
+538      read(*,*,iostat=ierr) val3,val4
+         if(ierr.ne.0) goto 538
+         write(*,*) 'choice phi min max'
+539      read(*,*,iostat=ierr) val5,val6
+         if(ierr.ne.0) goto 539
+         write(*,*) 'choice factor phi'
+535      read(*,*,iostat=ierr) val8
+         if(ierr.ne.0) goto 535
+         write(*,*) 'choice add phi'
+536      read(*,*,iostat=ierr) val7
+         if(ierr.ne.0) goto 536
+
+         do ig=1,ngridmx
+              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
+     &            (latfi(ig)*180./pi.le.val2)  .and.
+     &           (lonfi(ig)*180./pi.ge.val3)  .and. 
+     &           (lonfi(ig)*180./pi.le.val4) ) then
+
+                IF ((phisfi(ig).le.val6).and.
+     &              (phisfi(ig).ge.val5)) then
+
+                    phisfi(ig)=phisfi(ig)*val8
+                    phisfi(ig)=phisfi(ig)+val7
+                ENDIF
+
+              ENDIF
+         enddo
+         phisold = phis
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=37.
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       correct phisfi 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'correctps') then
+
+         phisold = phis
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=37.
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c        periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       No flat topo 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'toponoise') then
+
+         write(*,*) 'choice latitudes:'
+587      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 587
+         write(*,*) 'choice longitudes:'
+588      read(*,*,iostat=ierr) val3,val4
+         if(ierr.ne.0) goto 588
+         write(*,*) 'choice phi min max'
+589      read(*,*,iostat=ierr) val5,val6
+         if(ierr.ne.0) goto 589
+         write(*,*) 'choice amplitude min max new phi'
+585      read(*,*,iostat=ierr) val7,val8
+         if(ierr.ne.0) goto 585
+c         write(*,*) 'choice nb of random cases'
+c586      read(*,*,iostat=ierr) val7
+c         if(ierr.ne.0) goto 586
+
+         do ig=1,ngridmx
+              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
+     &            (latfi(ig)*180./pi.le.val2)  .and.
+     &           (lonfi(ig)*180./pi.ge.val3)  .and. 
+     &           (lonfi(ig)*180./pi.le.val4) ) then
+
+                IF ((phisfi(ig).le.val6).and.
+     &              (phisfi(ig).ge.val5)) then
+                    CALL RANDOM_NUMBER(val9)
+                    phisfi(ig)=val7+(val8-val7)*val9
+                ENDIF
+
+              ENDIF
+         enddo
+         phisold = phis
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=37.
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       fill_sp : fill SP with N2 ice and adjust phisfi (flat SP)
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'fill_sp') then
+
+         ! def SP:
+         val2=-33.   !lat1
+         val3=50.    !lat2
+         val4=140.   ! lon1
+         val5=-155.  ! lon2
+         write(*,*) 'def SP :'
+         write(*,*) 'lat :',val2,val3
+         write(*,*) 'lon :',val4,'180 / -180',val5
+         write(*,*) 'choice phisfi ref SP (ex: -800):'
+607      read(*,*,iostat=ierr) val6
+         if(ierr.ne.0) goto 607
+
+         write(*,*) 're scale the pressure :'
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=40.
+         write(*,*) 'r, sale height temperature =',r,temp,g
+         qsurf=0.        
+
+         write(*,*) 'latfi=',latfi
+         !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+         write(*,*) 'phis=',phis
+         write(*,*) 'phisfi=',phisfi
+         !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         phisold = phis
+         write(*,*) 'phisold=',phisold
+         do ig=1,ngridmx
+               
+              write(*,*) 'lat=',latfi(ig)*180./pi
+              write(*,*) 'lon=',lonfi(ig)*180./pi
+              write(*,*) 'phisfi=',phisfi(ig)
+              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
+     &            (latfi(ig)*180./pi.le.val3)  .and.
+     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
+     &                       (lonfi(ig)*180./pi.le.val5))  .or.
+     &         ((lonfi(ig)*180./pi.ge.val4)  .and. 
+     &                    (lonfi(ig)*180./pi.le.180.))) ) then
+                write(*,*) 'hello SP',phisfi(ig),ig
+                IF (phisfi(ig).le.val6) then
+                    qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
+                    phisfi(ig)=val6
+                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
+                ENDIF
+              ENDIF  
+         end do
+c        update new phis and ps
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+         do j=1,jjp1
+           do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+           enddo
+         enddo
+c       periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+
+c       bugch4 correct bug warm ch4 due to change of resolution 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'bugch4') then
+         write(*,*) 'Ok we are going to try to solve this bug'
+         write(*,*) 'First extract a profile of tsoil and tsurf'
+         write(*,*) 'that you want to set as reference :'
+         write(*,*) 'tsoil_ref and tsurf_ref '
+         write(*,*) 'number of points to check '
+546      read(*,*,iostat=ierr) jref1
+         if(ierr.ne.0) goto 546
+
+         write(*,*) 'choice acceptable tsurf range for ch4: t1,t2:'
+547      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 547
+
+         ! Check tsurf at all locations
+         DO ig=1,ngridmx
+           IF (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
+     &                         qsurf(ig,igcm_n2).eq.0.) then
+             IF ((tsurf(ig).lt.val1)  .or.
+     &            (tsurf(ig).ge.val2))  then
+               write(*,*) 'Pb tsurf point ',ig
+
+               do jref=1,jref1
+                 IF ((ig-jref.ge.1).and.qsurf(ig,igcm_n2).eq.0..and.
+     &                 (qsurf(int(ig-jref),igcm_ch4_ice).gt.0.001).and.
+     &                  (tsurf(ig-jref).gt.val1
+     &                  .and.tsurf(ig-jref).lt.val2)) then
+                    tsurf(ig)=tsurf(int(ig-jref))
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(int(ig-jref),l)
+                    ENDDO 
+                    goto 548
+                 ELSEIF ((ig+jref.le.ngridmx).and.
+     &                 qsurf(ig,igcm_n2).eq.0..and.
+     &                 (qsurf(int(ig+jref),igcm_ch4_ice).gt.0.001).and.
+     &                  (tsurf(ig+jref).gt.val1
+     &                  .and.tsurf(ig+jref).lt.val2)) then
+                    tsurf(ig)=tsurf(int(ig+jref))
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(int(ig+jref),l)
+                    ENDDO 
+                    goto 548
+                 ENDIF
+               enddo
+548            write(*,*) 'found point with ch4'
+             ENDIF 
+           ENDIF 
+         END DO
+
+c       bugn2 correct bug warm n2 due to change of resolution 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'bugn2') then
+         write(*,*) 'Ok we are going to try to solve this bug'
+         write(*,*) 'First extract a profile of tsoil and tsurf'
+         write(*,*) 'that you want to set as reference :'
+         write(*,*) 'tsoil_ref and tsurf_ref '
+         write(*,*) 'number of points to check '
+544      read(*,*,iostat=ierr) jref1
+         if(ierr.ne.0) goto 544
+
+         write(*,*) 'choice acceptable tsurf range for n2: t1,t2:'
+540      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 540
+
+         ! Check tsurf at all locations
+         DO ig=1,ngridmx
+           IF (qsurf(ig,1).gt.0.001) then
+             IF ((tsurf(ig).lt.val1)  .or.
+     &            (tsurf(ig).ge.val2))  then
+               write(*,*) 'Pb tsurf point ',ig
+
+                ! IF low topo et peu/pas de N2: add ch4, CO, N2
+               do jref=1,jref1
+                 IF ((ig-jref.ge.1).and.
+     &                 (qsurf(int(ig-jref),igcm_n2).gt.0.001).and.
+     &                  (tsurf(ig-jref).gt.val1
+     &                  .and.tsurf(ig-jref).lt.val2)) then
+                    tsurf(ig)=tsurf(int(ig-jref))
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(int(ig-jref),l)
+                    ENDDO 
+                    goto 541
+                 ELSEIF ((ig+jref.le.ngridmx).and.
+     &                 (qsurf(int(ig+jref),igcm_n2).gt.0.001).and.
+     &                  (tsurf(ig+jref).gt.val1
+     &                  .and.tsurf(ig+jref).lt.val2)) then
+                    tsurf(ig)=tsurf(int(ig+jref))
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(int(ig+jref),l)
+                    ENDDO 
+                    goto 541
+                 ENDIF
+               enddo
+541            write(*,*) 'found point with n2'
+             ENDIF 
+           ENDIF 
+         END DO
+
+         ! second tour
+!         DO ig=1,ngridmx
+!           IF (qsurf(ig,1).gt.0.001) then
+!             IF ((tsurf(ig).lt.val1)  .or.
+!     &            (tsurf(ig).ge.val2))  then
+!                ! IF low topo et peu/pas de N2: add ch4, CO, N2
+!                  IF ((ig-1.lt.1).or.(ig+1.gt.ngridmx)) then
+!                     write(*,*) 'pole : doing nothing'
+!                  ELSEIF (qsurf(ig-1,igcm_n2).gt.0.001) then
+!                    tsurf(ig)=tsurf(ig-1)
+!                    DO l=1,nsoilmx
+!                      tsoil(ig,l)=tsoil(ig-1,l)
+!                   ENDDO 
+!                 ELSEIF (qsurf(ig+1,igcm_n2).gt.0.001) then
+!                   tsurf(ig)=tsurf(ig+1)
+!                   DO l=1,nsoilmx
+!                      tsoil(ig,l)=tsoil(ig+1,l)
+!                   ENDDO 
+!                 ENDIF
+!            ENDIF 
+!          ENDIF 
+
+!        END DO
+
+c       flatnosp : flat topo outside a specific terrain (SP) 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'flatnosp') then
+
+          write(*,*) 'Choice of topography (m) below which we keep ?'
+551       read(*,*,iostat=ierr) val
+          if(ierr.ne.0) goto 551
+          write(*,*) 'gravity g is : ',g         
+          geop=g*val
+          write(*,*) 'Choice of minimum amount of N2 ice we keep ?'
+552       read(*,*,iostat=ierr) val2
+          if(ierr.ne.0) goto 552
+       
+          write(*,*) 'phis=',phis
+          write(*,*) 'phisfi=',phisfi
+          do ig=1,ngridmx
+               
+              IF (   (qsurf(ig,1).le.val2)  .and.
+     &               (phisfi(ig).gt.geop)  ) then
+                write(*,*) 'hello SP',phisfi(ig),ig
+                phisfi(ig)=0.
+              ENDIF
+          end do
+
+          phisold = phis
+          write(*,*) 're scale the pressure :'
+          r = 8.314511E+0 *1000.E+0/mugaz
+          temp=40.
+          write(*,*) 'r, sale height temperature =',r,temp,g
+c         update new phis and ps
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+          do j=1,jjp1
+            do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+            enddo
+          enddo
+c        periodicity of surface ps in longitude
+          do j=1,jjp1
+           ps(1,j) = ps(iip1,j)
+          end do
+
+c       flatregion : flat topo specific to region 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'flatregion') then
+
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 553      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 553
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 554      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 554
+          write(*,*) 'Choice of topography range ?'
+ 555      read(*,*,iostat=ierr) val5,val6
+          if(ierr.ne.0) goto 555
+          write(*,*) 'Choice topo ?'
+ 556      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 556
+       
+          write(*,*) 'phis=',phis
+          write(*,*) 'phisfi=',phisfi
+          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)  ) then 
+               
+              IF (   (phisfi(ig).ge.val5)  .and.
+     &               (phisfi(ig).le.val6)  ) then
+                write(*,*) 'hello topo',phisfi(ig),ig
+                phisfi(ig)=val7
+              ENDIF
+             ENDIF
+          end do
+
+          r = 8.314511E+0 *1000.E+0/mugaz
+          temp=40.
+          phisold = phis
+c         update new phis and ps
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+          do j=1,jjp1
+            do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+            enddo
+          enddo
+c        periodicity of surface ps in longitude
+          do j=1,jjp1
+           ps(1,j) = ps(iip1,j)
+          end do
+
+c       changetopo  
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'changetopo') then
+
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 573      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 573
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 574      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 574
+          write(*,*) 'Choice of topography range ?'
+ 575      read(*,*,iostat=ierr) val5,val6
+          if(ierr.ne.0) goto 575
+          write(*,*) 'Choice change topo factor?'
+ 576      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 576
+          write(*,*) 'Choice change topo add?'
+ 577      read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 577
+       
+          write(*,*) 'phis=',phis
+          write(*,*) 'phisfi=',phisfi
+          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)  ) then 
+               
+              IF (   (phisfi(ig).ge.val5)  .and.
+     &               (phisfi(ig).le.val6)  ) then
+                write(*,*) 'hello topo',phisfi(ig),ig
+                phisfi(ig)=phisfi(ig)*val7
+                phisfi(ig)=phisfi(ig)+val8
+              ENDIF
+             ENDIF
+          end do
+
+          r = 8.314511E+0 *1000.E+0/mugaz
+          temp=40.
+          phisold = phis
+c         update new phis and ps
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+          do j=1,jjp1
+            do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+            enddo
+          enddo
+c        periodicity of surface ps in longitude
+          do j=1,jjp1
+           ps(1,j) = ps(iip1,j)
+          end do
+
+c       corrtopo: corr topo specific bug
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'corrtopo') then
+
+          ! get min max lon
+          print*, latfi*180/pi
+          print*, '***************************'
+          print*, '***************************'
+          print*, '***************************'
+          print*, '***************************'
+          print*, '***************************'
+          print*, lonfi*180/pi
+          print*, 'iip1=',iip1
+          do ig=2,ngridmx-1
+             IF  (lonfi(ig)*180./pi.eq.-180.) THEN
+                 print*, lonfi(ig),lonfi(ig+iip1-2)
+                 phisfi(ig)=(phisfi(ig+1)+phisfi(ig+iip1))/2
+             ENDIF
+          enddo
+          phisold = phis
+          r = 8.314511E+0 *1000.E+0/mugaz
+          temp=40.
+c         update new phis and ps
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+          do j=1,jjp1
+            do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+            enddo
+          enddo
+c         periodicity of surface ps in longitude
+          do j=1,jjp1
+           ps(1,j) = ps(iip1,j)
+          end do
+
+c       asymtopo: 
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'asymtopo') then
+
+          ! get diff altitude
+          write(*,*) 'Diff N-S altitude in km (pos = N higher)  ?'
+ 591      read(*,*,iostat=ierr) val
+          if(ierr.ne.0) goto 591
+          write(*,*) 'Coeff smoot topo (small = smoother)  ?'
+ 592      read(*,*,iostat=ierr) val2
+          if(ierr.ne.0) goto 592
+
+          do ig=1,ngridmx
+           phisfi(ig)=phisfi(ig)+val*1000.*g*tanh(latfi(ig)*180/pi*val2)
+          enddo
+          phisold = phis
+          r = 8.314511E+0 *1000.E+0/mugaz
+          temp=40.
+c         update new phis and ps
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+          do j=1,jjp1
+            do i=1,iip1
+              ps(i,j) = ps(i,j) *
+     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
+            enddo
+          enddo
+c         periodicity of surface ps in longitude
+          do j=1,jjp1
+           ps(1,j) = ps(iip1,j)
+          end do
+
+c       seticeNH : set the ices in the SP crater with NH topo
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'seticeNH') then
+
+         open(333,file='./tsoil_180_30',form='formatted')
+         do i=1,22
+           read(333,*) levsoil(i), tempsoil(i)
+         enddo
+         close(333)
+         open(334,file='./tsurf_180_30',form='formatted')
+         read(334,*) val
+         close(334)
+
+         write(*,*) 'Tsoil and tsurf ref are:'
+         write(*,*) tempsoil
+         write(*,*) 'tsurf=',val
+
+         ! def SP:
+         val2=-43.   !lat1
+         val3=51.    !lat2
+         val4=143.   ! lon1
+         val5=-158.  ! lon2
+         val6=-150   ! phisfi min
+
+         ! Rm all N2 ice outside SP
+         DO ig=1,ngridmx
+           ! IF inside SP
+           IF (   (latfi(ig)*180./pi.ge.val2)  .and.
+     &            (latfi(ig)*180./pi.le.val3)  .and.
+     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
+     &                       (lonfi(ig)*180./pi.le.val5))  .or.
+     &         ((lonfi(ig)*180./pi.ge.val4)  .and. 
+     &                       (lonfi(ig)*180./pi.le.180.))) ) then
+
+                ! IF low topo et peu/pas de N2: add ch4, CO, N2
+                IF ((phisfi(ig).le.val6).and.
+     &                              (qsurf(ig,igcm_n2).le.500)) then
+                    qsurf(ig,igcm_n2)=1000.
+                    qsurf(ig,igcm_ch4_ice)=1000.
+                    qsurf(ig,igcm_co_ice)=1000.
+                    tsurf(ig)=val
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tempsoil(l)
+                    ENDDO 
+                ENDIF 
+
+                ! IF high topo : rm N2
+                IF ( (qsurf(ig,igcm_n2).ge.20.).and.
+     &                                 (phisfi(ig).ge.val6)  ) then
+                      qsurf(ig,igcm_n2)=0.
+
+                ENDIF
+                ! Mais on garde ch4 et CO en prenant la meme quantite
+                ! qu'une autre latitude
+                IF ( (qsurf(ig,igcm_ch4_ice).ge.20.)  .and.
+     &                                 (phisfi(ig).ge.val6)  ) then
+                      jref=int(ig/jjp1)+2
+                      qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+                ENDIF
+                IF ( (qsurf(ig,igcm_co_ice).ge.20.)  .and.
+     &                                 (phisfi(ig).ge.val6)  ) then
+                      jref=int(ig/jjp1)+2
+                      qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+                ENDIF
+
+           ELSE   ! Outside SP
+
+                IF (qsurf(ig,igcm_n2).ge.20.) then
+                      qsurf(ig,igcm_n2)=0.
+                ENDIF
+
+                IF (qsurf(ig,igcm_ch4_ice).ge.20.) then
+                      jref=int(ig/jjp1)+2
+                      qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+                ENDIF
+
+                IF (qsurf(ig,igcm_co_ice).ge.20.) then
+                      jref=int(ig/jjp1)+2
+                      qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+                ENDIF
+
            ENDIF
+        
+         END DO
+
+c       'nomountain '
+c       --------------
+        else if (modif(1:len_trim(modif)) .eq. 'nomountain') then
+             do j=1,jjp1
+               do i=1,iip1
+                 if (phis(i,j).gt.0.1) then 
+                       phis(i,j) = 0.
+                       ps(i,j)=ps(iim/4,j)    ! assuming no topo here
+                 endif 
+               enddo
+             enddo
+            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       'relief'
+c       --------------
+        else if (modif(1:len_trim(modif)) .eq.'relief') then
+          write(*,*) "add a circular mountain/crater"
+          write(*,*) 'Center: lat lon  ?'
+ 707      read(*,*,iostat=ierr) lat0, lon0
+          if(ierr.ne.0) goto 707
+          if(lon0.gt.180.) lon0=lon0-360. 
+          lat0 = lat0 * pi/180.
+          lon0 = lon0 * pi/180.
+
+          DO ig=1,ngridmx
+             if(abs(latfi(ig)-lat0).lt.pi/float(jjm) ) then
+               if(abs(lonfi(ig)-lon0).lt.2*pi/float(iim) ) then
+                  ig0 = ig
+               end if
+             end if
+          END DO
+          write(*,*) "Reference Point to be used:"
+          write(*,*) 'ig0',ig0
+          write(*,*) 'lat=',latfi(ig0)*180./pi
+          write(*,*) 'lon=',lonfi(ig0)*180./pi
+
+          write(*,*) 'radius (km), height (km) negative if crater ?'
+ 508      read(*,*,iostat=ierr) radm, height
+          if(ierr.ne.0) goto 508
+
+          write(*,*) 'Sale height temperature (ex:38) ?'
+ 509      read(*,*,iostat=ierr) temp
+          if(ierr.ne.0) goto 509
+
+          r = 8.314511E+0 *1000.E+0/mugaz
+          do j=1,jjp1
+            do i=1,iip1
+               dist= dist_pluto(lat0,lon0,rlatu(j),rlonv(i))
+               if (dist.le.radm) then
+                  phisprev= phis(i,j)
+                  phis(i,j)=phisprev+1000.*g*height*(radm-dist)/radm 
+                  write(*,*) 'lat=',rlatu(j)*180./pi
+                  write(*,*) 'lon=',rlonv(i)*180./pi
+                  write(*,*) 'dist', dist
+                  write(*,*) 'z(m)=', phis(i,j)/g
+                  ps(i,j) = ps(i,j) *
+     .            exp((phis(i,j))/(temp * r))
+               end if 
+            end do
+          end do
+
+c         periodicity of surface ps in longitude
+          do j=1,jjp1
+            ps(1,j) = ps(iip1,j)
+          end do
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       subsoil_n2 : choice of thermal inertia values for N2 only
+c       ----------------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then
+
+         write(*,*) 'New value for subsoil thermal inertia  ?'
+ 102     read(*,*,iostat=ierr) ith_bb
+         if(ierr.ne.0) goto 102
+         write(*,*) 'thermal inertia (new value):',ith_bb
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, 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:',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
+              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
+             write(*,*)'Level selected : ',layer(isoil),' m'
+             exit
+            endif
+           enddo
+         endif 
+
+         DO ig=1,ngridmx
+             DO j=iref,nsoilmx
+c                if((qsurf(ig,igcm_ch4_ice).lt.1.).and.
+c     &                         (qsurf(ig,igcm_co_ice).lt.1.))then
+                if (qsurf(ig,igcm_n2).gt.0.001) then
+                   ithfi(ig,j)=ith_bb
+                endif
+             ENDDO
+         ENDDO
+
+         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+c       subsoil_ch4 : choice of thermal inertia values for CH4 only
+c       ----------------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'subsoil_ch4') then
+
+         write(*,*) 'New value for subsoil thermal inertia  ?'
+ 103     read(*,*,iostat=ierr) ith_bb
+         if(ierr.ne.0) goto 103
+         write(*,*) 'thermal inertia (new value):',ith_bb
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, 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:',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
+              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
+             write(*,*)'Level selected : ',layer(isoil),' m'
+             exit
+            endif
+           enddo
+         endif 
+
+         DO ig=1,ngridmx
+             DO j=iref,nsoilmx
+                if (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
+     &                            qsurf(ig,igcm_n2).lt.0.001) then
+                   ithfi(ig,j)=ith_bb
+                endif
+             ENDDO
+         ENDDO
+
+         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+c       subsoil_alb : choice of thermal inertia values from albedo
+c       ----------------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'subsoil_alb') then
+
+         write(*,*) 'Choice range albedo for defining TI'
+ 145     read(*,*,iostat=ierr) val,val2
+         if(ierr.ne.0) goto 145
+         write(*,*) 'New value for subsoil thermal inertia  ?'
+ 144     read(*,*,iostat=ierr) ith_bb
+         if(ierr.ne.0) goto 144
+         write(*,*) 'thermal inertia (new value):',ith_bb
+
+         write(*,*)'At which depth (in m.) does the ice layer begin?'
+         write(*,*)'(currently, 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) val3
+          if (ierr.eq.0) then ! got a value, but do a sanity check
+            if(val3.gt.layer(nsoilmx)) then
+              write(*,*)'Depth should be less than ',layer(nsoilmx)
+              ierr=1
+            endif
+            if(val3.lt.layer(1)) then
+              if(val3.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(val3.eq.0) then
+           iref=1
+           write(*,*)'Level selected is first level: ',layer(iref),' m' 
+         else
+           do isoil=1,nsoilmx-1
+            if ((val3.gt.layer(isoil)).and.(val3.lt.layer(isoil+1)))
+     &       then
+             iref=isoil+1
+             write(*,*)'Level selected : ',layer(isoil+1),' m'
+             exit
+            endif
+           enddo
+         endif 
+
+         DO ig=1,ngridmx
+             IF (   (albfi(ig).ge.val)  .and.
+     &              (albfi(ig).le.val2)  ) then
+               DO j=iref,nsoilmx
+                    ithfi(ig,j)=ith_bb
+               ENDDO
+             ENDIF
+         ENDDO
+
+         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+c       subsoil_all : choice of thermal inertia values (global)
+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(*,*)'(currently, 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:',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
+              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
+             !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)))
+     &       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),')'
+         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
+            !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)))
+     &       then
+             iref=isoil+1
+             write(*,*)'Level selected : ',layer(isoil+1),' m'
+             exit
+            endif
+         enddo
+
+         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)
+
+c       Choice surface temperature
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'initsurf') then
+
+          write(*,*) 'New value for initial surface temperature  ?'
+ 105      read(*,*,iostat=ierr) tsurf_bb
+          if(ierr.ne.0) goto 105
+          write(*,*) 'new surface temperature (K):',tsurf_bb
+          DO ig=1,ngridmx
+                   tsurf(ig)=tsurf_bb
+          ENDDO
+
+c       Modify surface temperature (for GCM start)
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modtsurf') then
+
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 455      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 455
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 456      read(*,*,iostat=ierr) val4,val5
+          if(ierr.ne.0) goto 456
+          write(*,*) 'Choice topo min max '
+ 473      read(*,*,iostat=ierr) val9,val10
+          if(ierr.ne.0) goto 473
+          write(*,*) 'Choice tsurf min max '
+ 474      read(*,*,iostat=ierr) val11,val12
+          if(ierr.ne.0) goto 474
+          write(*,*) 'Choice multiplicative factor'
+ 457      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 457
+          write(*,*) 'Choice additional tsurf'
+ 458      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 458
+          write(*,*) 'Choice multiplicative factor tsoil'
+ 459      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 459
+          write(*,*) 'Choice additional tsoil'
+ 469      read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 469
+
+          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.
+     &              (phisfi(ig).ge.val9)  .and.
+     &              (phisfi(ig).lt.val10)  .and.
+     &              (tsurf(ig).ge.val11) .and.
+     &              (tsurf(ig).lt.val12)  ) then
+
+                    tsurf(ig)=tsurf(ig)*val3
+                    tsurf(ig)=tsurf(ig)+val6
+                    DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(ig,l)*val7
+                       tsoil(ig,l)=tsoil(ig,l)+val8
+                    ENDDO
+             ENDIF
+          ENDDO
+
+c       copy latitudes tsurf / tsoil
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'copylat') then
+           
+          write(*,*) 'all latitudes : ',rlatu*180./pi
+          write(*,*) 'Choice band to be modified lat1 ?'
+ 465      read(*,*,iostat=ierr) val
+          if(ierr.ne.0) goto 465
+          write(*,*) 'Choice band copied lat2 ?'
+ 466      read(*,*,iostat=ierr) val2
+          if(ierr.ne.0) goto 466
+          write(*,*) 'Choice multiplicative factor'
+ 467      read(*,*,iostat=ierr) val3
+          if(ierr.ne.0) goto 467
+          write(*,*) 'Choice additional tsurf'
+ 468     read(*,*,iostat=ierr) val4
+          if(ierr.ne.0) goto 468
+
+          j=1
+          DO ig=1,ngridmx
+             IF ((latfi(ig)*180./pi.lt.val2+180./iip1) .and.
+     &          (latfi(ig)*180./pi.gt.val2-180./iip1))  then
+                randtab(j)=ig
+                j=j+1
+             ENDIF
+          ENDDO
+          print*,j, ' latitudes found'
+          k=1
+          DO ig=1,ngridmx
+             IF ((latfi(ig)*180./pi.lt.val+180./iip1) .and.
+     &          (latfi(ig)*180./pi.gt.val-180./iip1))  then
+                tsurf(ig)=tsurf(randtab(k))*val3
+                tsurf(ig)=tsurf(ig)+val4
+                DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(randtab(k),l)*val3
+                       tsoil(ig,l)=tsoil(ig,l)+val4
+                ENDDO 
+                k=k+1
+             ENDIF
+          ENDDO
+          print*,k, ' latitudes copied'
+          IF (j.ne.k) stop
+
+c       copy longitudes
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'copylon') then
+           
+          write(*,*) 'all longitudes : ',rlonu*180./pi
+          write(*,*) 'Choice band to be modified lon1 ?'
+ 475      read(*,*,iostat=ierr) val
+          if(ierr.ne.0) goto 475
+          write(*,*) 'Choice band copied lon2 ?'
+ 476      read(*,*,iostat=ierr) val2
+          if(ierr.ne.0) goto 476
+
+          j=1
+          DO ig=2,ngridmx-1
+             IF ((lonfi(ig)*180./pi.lt.183.) .and.
+     &          (lonfi(ig)*180./pi.gt.180.))  then
+                randtab(j)=ig
+                j=j+1
+                print*, 'lon = ',lonfi(ig)
+             ENDIF
+          ENDDO
+          print*,j, ' longitudes found'
+          k=1
+          DO ig=2,ngridmx-1
+             IF ((lonfi(ig)*180./pi.lt.180.) .and.
+     &          (lonfi(ig)*180./pi.gt.179.))  then
+                tsurf(ig)=tsurf(randtab(k))
+                DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(randtab(k),l)
+                ENDDO
+                qsurf(ig,1)=qsurf(randtab(k),1)  
+                phisfi(ig)=phisfi(randtab(k))  
+                k=k+1
+             ENDIF
+          ENDDO
+          print*,k, ' longitudes copied'
+          IF (j.ne.k) stop
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+
+c       copy latlon
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'copylatlon') then
+           
+          write(*,*) 'all longitudes : ',rlonu*180./pi
+          write(*,*) 'Choice lat/lon to be copied ?'
+ 495      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 495
+          write(*,*) 'Choice indx lon1 lon2 for band ?'
+ 496      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 496
+          write(*,*) 'Choice latitude indx range where we copy ?'
+ 497      read(*,*,iostat=ierr) val5,val6
+          if(ierr.ne.0) goto 497
+
+          ! choice coordinate
+          DO ig=2,ngridmx-1
+             IF ( (lonfi(ig)*180./pi.gt.val2-1) .and.
+     &            (lonfi(ig)*180./pi.lt.val2+1) .and.
+     &            (latfi(ig)*180./pi.lt.val+1)  .and.
+     &            (latfi(ig)*180./pi.gt.val-1) ) then
+              ig0=ig
+              print*,'lat/lon=',latfi(ig)*180./pi,lonfi(ig)*180./pi,ig0
+             ENDIF
+          ENDDO
+
+          DO ig=2,ngridmx-1
+             IF ( (lonfi(ig)*180./pi.lt.val4) .and.
+     &            (lonfi(ig)*180./pi.gt.val3) .and.
+     &            (latfi(ig)*180./pi.lt.val6) .and.
+     &            (latfi(ig)*180./pi.gt.val5) .and.
+     &            (qsurf(ig,1).lt.0.001) ) then
+                tsurf(ig)=tsurf(ig0)
+                print*,'corrd=',latfi(ig)*180./pi,lonfi(ig)*180./pi
+                DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(ig0,l)
+                ENDDO
+             ENDIF
+          ENDDO
+
+c       Choice surface temperature depending on N2 ice distribution
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'tsurfice') then
+
+          write(*,*) 'Initial soil and surface temp below n2 ?'
+ 107      read(*,*,iostat=ierr) tsurf_bb
+          if(ierr.ne.0) goto 107
+          write(*,*) 'new surface/soil temp N2 (K):',tsurf_bb
+          write(*,*) 'Initial soil and surface temp when no n2 ?'
+ 108      read(*,*,iostat=ierr) tsurf_bb2
+          if(ierr.ne.0) goto 108
+          write(*,*) 'new surface/soil temp when no n2 (K):',tsurf_bb2
+          DO ig=1,ngridmx
+                   if (qsurf(ig,igcm_n2).gt.0.001) then
+                      tsurf(ig)=tsurf_bb
+                      do l=1,nsoilmx
+                         tsoil(ig,l) = tsurf_bb
+                      end do
+                   else if (tsurf_bb2.ne.0) then
+                      tsurf(ig)=tsurf_bb2
+                      do l=1,nsoilmx
+                         tsoil(ig,l) = tsurf_bb2
+                      end do
+                   endif
+          ENDDO
+
+c       read an albedo map
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'albedomap') then
+           
+          ! Get field 2D
+          fichnom = 'albedo.nc'
+          ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
+          IF (ierr.NE.NF_NOERR) THEN
+            write(6,*)' Problem opening albedo file:',fichnom
+            write(6,*)' ierr = ', ierr
+            CALL ABORT
+          ENDIF
+
+          ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"), 
+     &                                  nvarid_input)
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not find asked variable in albedo.nc"
+           CALL abort
+          ENDIF
+
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
+     &                                 field_input)
+#else
+          ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
+#endif
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not get asked variable in albedo.nc"
+           CALL abort
+          ELSE 
+           PRINT*, "Got variable in albedo.nc"
+          ENDIF
+
+          DO ig=1,ngridmx
+                      albfi(ig)=field_input(ig)
+          ENDDO
+
+c       slope : add slope on all longitudes
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'slope') then
+
+         write(*,*) 'choice latitude alt minimum & altitude value (m):'
+603      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 603
+         write(*,*) 'choice latitude alt maximum & altitude value (m):'
+604      read(*,*,iostat=ierr) val3,val4
+         if(ierr.ne.0) goto 604
+
+         write(*,*) 're scale the pressure :'
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=40.
+         write(*,*) 'r, sale height temperature =',r,temp
+                 
+         do j=1,jjp1
+           do i=1,iip1
+              phisprev= phis(i,j)
+              IF ( (  (rlatu(j)*180./pi.ge.val1)  .and.
+     &            (rlatu(j)*180./pi.le.val3) ) .or.
+     &          (  (rlatu(j)*180./pi.le.val1)  .and.
+     &            (rlatu(j)*180./pi.ge.val3) )) then
+
+                    val5=abs(val2-val4)/abs(val1-val3)*
+     &                           abs(val1-rlatu(j)*180./pi)+val2
+                    phis(i,j)=val5*g
+                    ps(i,j) = ps(i,j) *
+     .              exp((phisprev-phis(i,j))/(temp * r))
+              ENDIF  
+           end do
+         end do
+c        periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       digsp : change topography of SP
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'digsp') then
+
+         write(*,*) 'choice latitudes:'
+517      read(*,*,iostat=ierr) val1,val2
+         if(ierr.ne.0) goto 517
+         write(*,*) 'choice longitudes:'
+518      read(*,*,iostat=ierr) val3,val4
+         if(ierr.ne.0) goto 518
+         write(*,*) 'choice phi limite'
+519      read(*,*,iostat=ierr) val5
+         if(ierr.ne.0) goto 519
+         write(*,*) 'choice change alt (m)'
+520      read(*,*,iostat=ierr) val6
+         if(ierr.ne.0) goto 520
+
+         write(*,*) 're scale the pressure :'
+         r = 8.314511E+0 *1000.E+0/mugaz
+         temp=40.
+         write(*,*) 'r, sale height temperature =',r,temp
+                 
+         do j=1,jjp1
+           do i=1,iip1
+              phisprev= phis(i,j)
+              IF ( (  (rlatu(j)*180./pi.ge.val1)  .and.
+     &            (rlatu(j)*180./pi.le.val2) ) .and.
+     &          (  (rlonv(j)*180./pi.ge.val3)  .and.
+     &            (rlonv(j)*180./pi.le.val4) ) .and.
+     &            (phis(i,j).le.val5)) then
+
+                    phis(i,j)=phis(i,j)+val6*g
+                    ps(i,j) = ps(i,j) *
+     .              exp((phisprev-phis(i,j))/(temp * r))
+              ENDIF  
+           end do
+         end do
+c        periodicity of surface ps in longitude
+         do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+         end do
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       copy field 2D
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'copytsoil') then
+           
+          ! Get field 2D
+          fichnom = 'startfi_input.nc'
+          ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
+          IF (ierr.NE.NF_NOERR) THEN
+            write(6,*)' Problem opening file:',fichnom
+            write(6,*)' ierr = ', ierr
+            CALL ABORT
+          ENDIF
+
+          ! Choice variable to be copied
+c          write(*,*) 'Choice variable to be copied ?'
+c515       read(*,fmt='(a20)',iostat=ierr) modif
+c          if(ierr.ne.0) goto 515
+c          write(*,*) 'variable asked :',trim(modif)
+
+          ierr = NF_INQ_VARID (nid_fi_input, trim("tsurf"),nvarid_input)
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not find asked variable in startfi_input.nc"
+           CALL abort
+          ENDIF
+          ierr = NF_INQ_VARID (nid_fi_input, trim("tsoil"),
+     &                                             nvarid_inputs)
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not find asked variable s in startfi_input.nc"
+           CALL abort
+          ENDIF
+
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
+     &                                              field_input)
+#else
+          ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
+#endif
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not get asked variable in startfi_input.nc"
+           CALL abort
+          ELSE 
+           PRINT*, "Got variable in startfi_input.nc"
+          ENDIF
+#ifdef NC_DOUBLE
+          ierr=NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_inputs,
+     &                        field_inputs)
+#else
+          ierr=NF_GET_VAR_REAL(nid_fi_input,nvarid_inputs,field_inputs)
+#endif
+          IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "Could not get asked variable s in startfi_input.nc"
+           CALL abort
+          ELSE 
+           PRINT*, "Got variable s in startfi_input.nc"
+          ENDIF
+
+          write(*,*) 'Choice domain lon1 lon2 ?'
+525       read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 525
+          write(*,*) 'Choice domain lat1 lat2 ?'
+526       read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 526
+          write(*,*) 'No change if N2 ice (0) / Change (1) ?'
+527       read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 527
+
+          ! Loop
+          DO ig=1,ngridmx
+             IF ( (lonfi(ig)*180./pi.ge.val) .and.
+     &            (lonfi(ig)*180./pi.le.val2) .and.
+     &            (latfi(ig)*180./pi.ge.val3)  .and.
+     &            (latfi(ig)*180./pi.le.val4) ) then
+                   if (qsurf(ig,igcm_n2).lt.0.001.or.val5.eq.1) then
+                      tsurf(ig)=field_input(ig)
+                      do l=1,nsoilmx
+                         tsoil(ig,l) = field_input(ig)
+                        !tsoil(ig,l) = field_inputs(ig,l)
+                      end do
+                   endif
+             ENDIF
+          ENDDO
+
+c       modify ice distribution depending on albedo
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'mod_distrib_ice') then
+
+          write(*,*) 'Choice range albedo for CH4 ice'
+ 220      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 220
+          write(*,*) 'Choice range albedo for N2 ice'
+ 221      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 221
+          write(*,*) 'Choice extra mass of CH4 ice'
+ 222      read(*,*,iostat=ierr) val5
+          if(ierr.ne.0) goto 222
+          write(*,*) 'Choice extra mass of N2 ice'
+ 223      read(*,*,iostat=ierr) val6
+          if(ierr.ne.0) goto 223
+
+          DO ig=1,ngridmx
+             IF (   (albfi(ig).ge.val)  .and.
+     &              (albfi(ig).le.val2) ) then
+                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
+             ENDIF
+             IF (   (albfi(ig).ge.val3)  .and.
+     &              (albfi(ig).le.val4) ) then
+                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
+             ENDIF
           ENDDO
 
@@ -2121,11 +3811,11 @@
                if (map_pluto_dat(i,j).eq.3) then
                   N2_pluto_dat(i,j)=100000.
-               elseif (map_pluto_dat(i,j).eq.4) then
+               else if (map_pluto_dat(i,j).eq.4) then
                   CH4_pluto_dat(i,j)=100000.
-               elseif (map_pluto_dat(i,j).eq.0) then
+               else if (map_pluto_dat(i,j).eq.0) then
                   alb_dat(i,j)=0.3
-               elseif (map_pluto_dat(i,j).eq.6) then
+               else if (map_pluto_dat(i,j).eq.6) then
                   alb_dat(i,j)=0.6
-               elseif (map_pluto_dat(i,j).eq.7) then
+               else if (map_pluto_dat(i,j).eq.7) then
                   alb_dat(i,j)=0.1
                endif
@@ -2352,6 +4042,6 @@
       real dlon, dlat
       real a,c
-      real rad
-      rad=1190. ! km
+      real radp
+      radp=1190. ! km
 
       dlon = lon2 - lon1
@@ -2359,5 +4049,5 @@
       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
+      dist_pluto = radp * c
       return
       end
