Changeset 3382


Ignore:
Timestamp:
Jun 17, 2024, 9:55:26 AM (6 months ago)
Author:
tbertrand
Message:

LMDZ.PLUTO
Adding more options in newstart from Pluto.old
TB

Location:
trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/datareadnc.F

    r3184 r3382  
    5858
    5959      INTEGER    imd,jmd,imdp1,jmdp1
    60       parameter    (imd=360,jmd=179,imdp1=361,jmdp1=180)
     60      parameter    (imd=1440,jmd=719,imdp1=1441,jmdp1=720)
    6161
    6262      INTEGER    iimp1
    … …  
    8585
    8686      INTEGER klatdat,ngridmixgdat
    87       PARAMETER (klatdat=180,ngridmixgdat=360)
     87      PARAMETER (klatdat=720,ngridmixgdat=1440)
    8888
    8989c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
    … …  
    221221          latitude(:)=(pi/180.)*latitude(:)
    222222
    223           call grid_noro1(360, 180, longitude, latitude, zdata,
     223          call grid_noro1(1440,720, longitude, latitude, zdata,
    224224     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
    225225

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/grid_noro1.F

    r3184 r3382  
    5656       implicit none
    5757       integer iusn,jusn,iext
    58        parameter(iusn=360,jusn=180,iext=40)
     58       parameter(iusn=1440,jusn=720,iext=40)
    5959c!-*-      include 'param1'
    6060c!-*-      include 'comcstfi.h'

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

    r3371 r3382  
    6464c Variables pour les lectures NetCDF des fichiers "start_archive"
    6565c--------------------------------------------------
    66       INTEGER nid_dyn, nid_fi,nid,nvarid
     66      INTEGER nid_dyn, nid_fi,nid,nvarid,nvarid_input,nvarid_inputs
     67      INTEGER nid_fi_input
    6768      INTEGER length
    6869      parameter (length = 100)
    … …  
    108109      real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
    109110      REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
     111      REAL field_input(ngridmx)
     112      REAL,ALLOCATABLE :: field_inputs(:,:)
    110113
    111114      INTEGER i,j,l,isoil,ig,idum,k
    … …  
    137140c-------------------
    138141      real choix_1,pp,choice
     142      integer randtab(ngridmx)
    139143      character*80      fichnom
    140144      character*250  filestring
    … …  
    143147      real z_reel(iip1,jjp1)
    144148      real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
    145       real ptoto,pcap,patm,airetot,ptotn,patmn,psea
     149      real tsurf_bb,tsurf_bb2
     150      real ptoto,pcap,patm,airetot,ptotn,patmn,psea,geop
     151      real tempsoil(22),levsoil(22)
    146152      character*1 yes
    147153      logical :: flagtset=.false. ,  flagps0=.false.
    … …  
    230236     
    231237      allocate(tsoil(ngridmx,nsoilmx))
     238      allocate(field_inputs(ngridmx,nsoilmx))
    232239      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
    233240     
    … …  
    741748          enddo
    742749
     750c       ptot : Modification of the total pressure: ice + current atmosphere
     751c       -------------------------------------------------------------------
     752        else if (modif(1:len_trim(modif)).eq.'ptot') then
     753
     754c         calcul de la pression totale glace + atm actuelle
     755          patm=0.
     756          airetot=0.
     757          pcap=0.
     758          DO j=1,jjp1
     759             DO i=1,iim
     760                ig=1+(j-2)*iim +i
     761                if(j.eq.1) ig=1
     762                if(j.eq.jjp1) ig=ngridmx
     763                patm = patm + ps(i,j)*aire(i,j)
     764                airetot= airetot + aire(i,j)
     765             ENDDO
     766          ENDDO
     767          ptoto = pcap + patm
     768
     769          print*,'Current total pressure at surface ',
     770     &       ptoto/airetot
     771
     772          print*,'new value?'
     773          read(*,*) ptotn
     774          ptotn=ptotn*airetot
     775          patmn=ptotn-pcap
     776          print*,'ptoto,patm,ptotn,patmn'
     777          print*,ptoto,patm,ptotn,patmn
     778          print*,'Mult. factor for pressure (atm only)', patmn/patm
     779          do j=1,jjp1
     780             do i=1,iip1
     781                ps(i,j)=ps(i,j)*patmn/patm
     782             enddo
     783          enddo
     784
     785c        Correction pour la conservation des traceurs
     786         yes=' '
     787         do while ((yes.ne.'y').and.(yes.ne.'n'))
     788            write(*,*) 'Do you wish to conserve tracer total mass (y)',
     789     &              ' or tracer mixing ratio (n) ?'
     790             read(*,fmt='(a)') yes
     791         end do
     792
     793         if (yes.eq.'y') then
     794           write(*,*) 'OK : conservation of tracer total mass'
     795           DO iq =1, nqtot
     796             DO l=1,llm
     797               DO j=1,jjp1
     798                  DO i=1,iip1
     799                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
     800                  ENDDO
     801               ENDDO
     802             ENDDO
     803           ENDDO
     804          else
     805            write(*,*) 'OK : conservation of tracer mixing ratio'
     806          end if
     807
    743808c       qname : change tracer name
    744809c       --------------------------
    … …  
    14101475          if(ierr.ne.0) goto 427
    14111476          write(*,*) 'Choice additional n2'
    1412  428     read(*,*,iostat=ierr) val6
     1477 428      read(*,*,iostat=ierr) val6
    14131478          if(ierr.ne.0) goto 428
    14141479
    … …  
    15361601          if(ierr.ne.0) goto 503
    15371602          write(*,*) 'Choice phisfi min max where change n2'
    1538  504     read(*,*,iostat=ierr) val6,val11
     1603 504      read(*,*,iostat=ierr) val6,val11
    15391604          if(ierr.ne.0) goto 504
    15401605          write(*,*) 'Choice multiplicative factor'
    … …  
    15421607          if(ierr.ne.0) goto 505
    15431608          write(*,*) 'Choice additional n2'
    1544  506     read(*,*,iostat=ierr) val8
     1609 506      read(*,*,iostat=ierr) val8
    15451610          if(ierr.ne.0) goto 506
    15461611          write(*,*) 'Choice ind lon equivalent for change tsurf tsoil'
    1547  507     read(*,*,iostat=ierr) iref
     1612 507      read(*,*,iostat=ierr) iref
    15481613          if(ierr.ne.0) goto 507
    15491614
    … …  
    15851650          ENDDO
    15861651
    1587 
    15881652c       modn2_topo : adding/remove n2 ice
    15891653c       ----------------------------------------------------------
    … …  
    16041668          if(ierr.ne.0) goto 444
    16051669          write(*,*) 'Choice additional n2'
    1606  445     read(*,*,iostat=ierr) val6
     1670 445      read(*,*,iostat=ierr) val6
    16071671          if(ierr.ne.0) goto 445
    16081672
    … …  
    16271691          if(ierr.ne.0) goto 471
    16281692          write(*,*) 'Choice additional n2'
    1629  472     read(*,*,iostat=ierr) val6
     1693 472      read(*,*,iostat=ierr) val6
    16301694          if(ierr.ne.0) goto 472
    16311695
    … …  
    16521716          if(ierr.ne.0) goto 513
    16531717          write(*,*) 'Choice phisfi min max where change n2'
    1654  514     read(*,*,iostat=ierr) val6,val11
     1718 514      read(*,*,iostat=ierr) val6,val11
    16551719          if(ierr.ne.0) goto 514
    16561720
    … …  
    18161880          enddo
    18171881          write(*,*)'OK: Soil thermal inertia has been reset to referenc
    1818      &e surface values'
     1882     &               e surface values'
    18191883          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
    18201884          ithfi(:,:)=inertiedat(:,:)
    … …  
    18411905!       subsoilice_n: Put deep ice layer in northern hemisphere soil
    18421906!       ------------------------------------------------------------
    1843        else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
     1907        else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
    18441908
    18451909         write(*,*)'From which latitude (in deg.), up to the north pole,
    1846      &should we put subterranean ice?'
     1910     &       should we put subterranean ice?'
    18471911         ierr=1
    18481912         do while (ierr.ne.0)
    … …  
    19171981         do while (ierr.ne.0)
    19181982          write(*,*)'What is the value of subterranean ice thermal inert
    1919      &ia? (e.g.: 2000)'
     1983     &      ia? (e.g.: 2000)'
    19201984          read(*,*,iostat=ierr)iceith
    19211985         enddo ! of do while
    … …  
    19382002         enddo ! of do j=1,jjp1
    19392003 
    1940 
    19412004         CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    19422005
    … …  
    19462009
    19472010         write(*,*)'From which latitude (in deg.), down to the south pol
    1948      &e, should we put subterranean ice?'
     2011     &     e, should we put subterranean ice?'
    19492012         ierr=1
    19502013         do while (ierr.ne.0)
    … …  
    20112074           endif
    20122075          enddo
    2013  
     2076  
    20142077         ! thermal inertia of the ice:
    20152078         ierr=1
    20162079         do while (ierr.ne.0)
    20172080          write(*,*)'What is the value of subterranean ice thermal inert
    2018      &ia? (e.g.: 2000)'
     2081     &     ia? (e.g.: 2000)'
    20192082          read(*,*,iostat=ierr)iceith
    20202083         enddo ! of do while
    … …  
    20512114
    20522115          ! Definition SP:
    2053          val3=-33.   !lat1
    2054          val4=50.    !lat2
    2055          val5=140.   ! lon1
    2056          val6=-155.  ! lon2
    2057          choice=-50. ! min phisfi
     2116          val3=-33.   !lat1
     2117          val4=50.    !lat2
     2118          val5=140.   ! lon1
     2119          val6=-155.  ! lon2
     2120          choice=-50. ! min phisfi
     2121          write(*,*) 'def SP :'
     2122          write(*,*) 'lat :',val3,val4
     2123          write(*,*) 'lon :',val5,'180 / -180',val6
     2124          write(*,*) 'choice phisfi min :',choice
     2125
     2126          DO ig=1,ngridmx
     2127
     2128            IF (   (latfi(ig)*180./pi.ge.val3)  .and.
     2129     &            (latfi(ig)*180./pi.le.val4)  .and.
     2130     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
     2131     &                       (lonfi(ig)*180./pi.le.val6))  .or.
     2132     &         ((lonfi(ig)*180./pi.ge.val5)  .and.
     2133     &                       (lonfi(ig)*180./pi.le.180.))) ) then
     2134
     2135                IF ((phisfi(ig).le.choice)) then  !.and.
     2136c    &                              (qsurf(ig,igcm_n2).ge.50)) then
     2137
     2138                  qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
     2139                  qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8
     2140                  qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9
     2141                ENDIF
     2142            ENDIF
     2143          ENDDO
     2144
     2145c       subsoil_SP : choice of thermal inertia values for SP
     2146c       ----------------------------------------------------------------
     2147        else if (modif(1:len_trim(modif)) .eq. 'subsoil_SP') then
     2148
     2149         write(*,*) 'New value for subsoil thermal inertia in SP ?'
     2150 510     read(*,*,iostat=ierr) ith_bb
     2151         if(ierr.ne.0) goto 510
     2152         write(*,*) 'thermal inertia (new value):',ith_bb
     2153
     2154         write(*,*)'At which depth (in m.) does the ice layer begin?'
     2155         write(*,*)'(currently, the deepest soil layer extends down to:'
     2156     &              ,layer(1),' - ',layer(nsoilmx),')'
     2157         write(*,*)'write 0 for uniform value for all subsurf levels?'
     2158         ierr=1
     2159         do while (ierr.ne.0)
     2160          read(*,*,iostat=ierr) val2
     2161          write(*,*)'val2:',val2,'ierr=',ierr
     2162          if (ierr.eq.0) then ! got a value, but do a sanity check
     2163            if(val2.gt.layer(nsoilmx)) then
     2164              write(*,*)'Depth should be less than ',layer(nsoilmx)
     2165              ierr=1
     2166            endif
     2167            if(val2.lt.layer(1)) then
     2168              if(val2.eq.0) then
     2169               write(*,*)'Thermal inertia set for all subsurface layers'
     2170               ierr=0
     2171              else 
     2172               write(*,*)'Depth should be more than ',layer(1)
     2173               ierr=1
     2174              endif
     2175           endif
     2176          endif
     2177         enddo ! of do while
     2178
     2179         ! find the reference index iref the depth corresponds to
     2180         if(val2.eq.0) then
     2181           iref=1
     2182           write(*,*)'Level selected is first level: ',layer(iref),' m'
     2183         else
     2184           do isoil=1,nsoilmx-1
     2185            if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     2186     &       then
     2187             iref=isoil
     2188             write(*,*)'Level selected : ',layer(isoil),' m'
     2189             exit
     2190            endif
     2191           enddo
     2192         endif
     2193
     2194         ! Definition SP:
     2195         val3=-50.   !lat1
     2196         val4=60.    !lat2
     2197         val5=130.   ! lon1
     2198         val6=-140.  ! lon2
     2199         choice=-100. ! min phisfi
    20582200         write(*,*) 'def SP :'
    20592201         write(*,*) 'lat :',val3,val4
    … …  
    20702212     &                       (lonfi(ig)*180./pi.le.180.))) ) then
    20712213
    2072                 IF ((phisfi(ig).le.choice)) then  !.and.
    2073 c    &                              (qsurf(ig,igcm_n2).ge.50)) then
    2074 
    2075                   qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
    2076                   qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8
    2077                   qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9
     2214                IF ((phisfi(ig).le.choice) .and.
     2215     &                              (qsurf(ig,igcm_n2).ge.1000)) then
     2216                    DO j=iref,nsoilmx
     2217                       ithfi(ig,j)=ith_bb
     2218                    ENDDO
     2219                ENDIF
     2220           ENDIF
     2221         ENDDO
     2222
     2223         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     2224
     2225c       topo_sp : change topo of SP
     2226c       -----------------------------------------------------
     2227        else if (modif(1:len_trim(modif)) .eq. 'topo_sp') then
     2228
     2229         ! def SP:
     2230         val2=-33.   !lat1
     2231         val3=50.    !lat2
     2232         val4=140.   ! lon1
     2233         val5=-155.  ! lon2
     2234         write(*,*) 'def SP :'
     2235         write(*,*) 'lat :',val2,val3
     2236         write(*,*) 'lon :',val4,'180 / -180',val5
     2237         write(*,*) 'choice phisfi min (ex: -100):'
     2238606      read(*,*,iostat=ierr) val6
     2239         if(ierr.ne.0) goto 606
     2240         write(*,*) 'choice coefficient for phisfi (ex: 2):'
     2241605      read(*,*,iostat=ierr) choice
     2242         if(ierr.ne.0) goto 605
     2243
     2244         write(*,*) 're scale the pressure :'
     2245         r = 8.314511E+0 *1000.E+0/mugaz
     2246         temp=40.
     2247         write(*,*) 'r, sale height temperature =',r,temp
     2248                 
     2249         do j=1,jjp1
     2250           do i=1,iip1
     2251              phisprev= phis(i,j)
     2252              IF (   (rlatu(j)*180./pi.ge.val2)  .and.
     2253     &            (rlatu(j)*180./pi.le.val3)  .and.
     2254     &      (  ((rlonv(i)*180./pi.ge.-180.)  .and.
     2255     &                       (rlonv(i)*180./pi.le.val5))  .or.
     2256     &         ((rlonv(i)*180./pi.ge.val4)  .and.
     2257     &                       (rlonv(i)*180./pi.le.180.))) ) then
     2258
     2259                IF (phis(i,j).le.val6) then
     2260                    phis(i,j)=phis(i,j)*choice
     2261                    ps(i,j) = ps(i,j) *
     2262     .              exp((phisprev-phis(i,j))/(temp * r))
     2263                ENDIF
     2264              ENDIF 
     2265           end do
     2266         end do
     2267c       periodicity of surface ps in longitude
     2268         do j=1,jjp1
     2269          ps(1,j) = ps(iip1,j)
     2270         end do
     2271         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2272
     2273c       fill_sp inv: fill inverted SP with N2 ice and adjust phisfi (flat SP)
     2274c       -----------------------------------------------------
     2275        else if (modif(1:len_trim(modif)) .eq. 'fill_sp_inv') then
     2276
     2277         ! def SP:
     2278         val2=-33.   !lat1
     2279         val3=50.    !lat2
     2280         val4=-40.   ! lon1
     2281         val5=25.  ! lon2
     2282         write(*,*) 'def SP :'
     2283         write(*,*) 'lat :',val2,val3
     2284         write(*,*) 'lon :',val4,val5
     2285         write(*,*) 'choice phisfi ref SP (ex: -800):'
     2286706      read(*,*,iostat=ierr) val6
     2287         if(ierr.ne.0) goto 706
     2288
     2289         write(*,*) 're scale the pressure :'
     2290         r = 8.314511E+0 *1000.E+0/mugaz
     2291         temp=40.
     2292         write(*,*) 'r, sale height temperature =',r,temp,g
     2293         qsurf=0.       
     2294
     2295         write(*,*) 'latfi=',latfi
     2296         !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2297         write(*,*) 'phis=',phis
     2298         write(*,*) 'phisfi=',phisfi
     2299         !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2300         phisold = phis
     2301         write(*,*) 'phisold=',phisold
     2302         do ig=1,ngridmx
     2303               
     2304              write(*,*) 'lat=',latfi(ig)*180./pi
     2305              write(*,*) 'lon=',lonfi(ig)*180./pi
     2306              write(*,*) 'phisfi=',phisfi(ig)
     2307              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
     2308     &            (latfi(ig)*180./pi.le.val3)  .and.
     2309     &           (lonfi(ig)*180./pi.ge.val4)  .and.
     2310     &           (lonfi(ig)*180./pi.le.val5) ) then
     2311                write(*,*) 'hello SP',phisfi(ig),ig
     2312                IF (phisfi(ig).le.val6) then
     2313                    qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
     2314                    phisfi(ig)=val6
     2315                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
    20782316                ENDIF
     2317              ENDIF 
     2318         end do
     2319c        update new phis and ps
     2320         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2321         do j=1,jjp1
     2322           do i=1,iip1
     2323              ps(i,j) = ps(i,j) *
     2324     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2325           enddo
     2326         enddo
     2327c       periodicity of surface ps in longitude
     2328         do j=1,jjp1
     2329          ps(1,j) = ps(iip1,j)
     2330         end do
     2331
     2332c       adjust phisfi according to the amount of N2 ice
     2333c       -----------------------------------------------------
     2334        else if (modif(1:len_trim(modif)) .eq. 'adjustphi') then
     2335
     2336         phisold = phis
     2337         do ig=1,ngridmx
     2338            phisfi(ig)=phisfi(ig)+qsurf(ig,igcm_n2)*g/1000.
     2339         end do
     2340c        update new phis and ps
     2341         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2342         r = 8.314511E+0 *1000.E+0/mugaz
     2343         temp=37.
     2344         do j=1,jjp1
     2345           do i=1,iip1
     2346              ps(i,j) = ps(i,j) *
     2347     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2348           enddo
     2349         enddo
     2350c        periodicity of surface ps in longitude
     2351         do j=1,jjp1
     2352          ps(1,j) = ps(iip1,j)
     2353         end do
     2354
     2355c       correct phisfi
     2356c       -----------------------------------------------------
     2357        else if (modif(1:len_trim(modif)) .eq. 'correctphi') then
     2358
     2359         write(*,*) 'choice latitudes:'
     2360537      read(*,*,iostat=ierr) val1,val2
     2361         if(ierr.ne.0) goto 537
     2362         write(*,*) 'choice longitudes:'
     2363538      read(*,*,iostat=ierr) val3,val4
     2364         if(ierr.ne.0) goto 538
     2365         write(*,*) 'choice phi min max'
     2366539      read(*,*,iostat=ierr) val5,val6
     2367         if(ierr.ne.0) goto 539
     2368         write(*,*) 'choice factor phi'
     2369535      read(*,*,iostat=ierr) val8
     2370         if(ierr.ne.0) goto 535
     2371         write(*,*) 'choice add phi'
     2372536      read(*,*,iostat=ierr) val7
     2373         if(ierr.ne.0) goto 536
     2374
     2375         do ig=1,ngridmx
     2376              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
     2377     &            (latfi(ig)*180./pi.le.val2)  .and.
     2378     &           (lonfi(ig)*180./pi.ge.val3)  .and.
     2379     &           (lonfi(ig)*180./pi.le.val4) ) then
     2380
     2381                IF ((phisfi(ig).le.val6).and.
     2382     &              (phisfi(ig).ge.val5)) then
     2383
     2384                    phisfi(ig)=phisfi(ig)*val8
     2385                    phisfi(ig)=phisfi(ig)+val7
     2386                ENDIF
     2387
     2388              ENDIF
     2389         enddo
     2390         phisold = phis
     2391c        update new phis and ps
     2392         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2393         r = 8.314511E+0 *1000.E+0/mugaz
     2394         temp=37.
     2395         do j=1,jjp1
     2396           do i=1,iip1
     2397              ps(i,j) = ps(i,j) *
     2398     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2399           enddo
     2400         enddo
     2401c       periodicity of surface ps in longitude
     2402         do j=1,jjp1
     2403           do i=1,iip1
     2404              ps(i,j) = ps(i,j) *
     2405     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2406           enddo
     2407         enddo
     2408c       periodicity of surface ps in longitude
     2409         do j=1,jjp1
     2410          ps(1,j) = ps(iip1,j)
     2411         end do
     2412         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2413
     2414c       correct phisfi
     2415c       -----------------------------------------------------
     2416        else if (modif(1:len_trim(modif)) .eq. 'correctps') then
     2417
     2418         phisold = phis
     2419c        update new phis and ps
     2420         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2421         r = 8.314511E+0 *1000.E+0/mugaz
     2422         temp=37.
     2423         do j=1,jjp1
     2424           do i=1,iip1
     2425              ps(i,j) = ps(i,j) *
     2426     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2427           enddo
     2428         enddo
     2429c       periodicity of surface ps in longitude
     2430         do j=1,jjp1
     2431           do i=1,iip1
     2432              ps(i,j) = ps(i,j) *
     2433     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2434           enddo
     2435         enddo
     2436c        periodicity of surface ps in longitude
     2437         do j=1,jjp1
     2438          ps(1,j) = ps(iip1,j)
     2439         end do
     2440         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2441
     2442c       No flat topo
     2443c       -----------------------------------------------------
     2444        else if (modif(1:len_trim(modif)) .eq. 'toponoise') then
     2445
     2446         write(*,*) 'choice latitudes:'
     2447587      read(*,*,iostat=ierr) val1,val2
     2448         if(ierr.ne.0) goto 587
     2449         write(*,*) 'choice longitudes:'
     2450588      read(*,*,iostat=ierr) val3,val4
     2451         if(ierr.ne.0) goto 588
     2452         write(*,*) 'choice phi min max'
     2453589      read(*,*,iostat=ierr) val5,val6
     2454         if(ierr.ne.0) goto 589
     2455         write(*,*) 'choice amplitude min max new phi'
     2456585      read(*,*,iostat=ierr) val7,val8
     2457         if(ierr.ne.0) goto 585
     2458c         write(*,*) 'choice nb of random cases'
     2459c586      read(*,*,iostat=ierr) val7
     2460c         if(ierr.ne.0) goto 586
     2461
     2462         do ig=1,ngridmx
     2463              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
     2464     &            (latfi(ig)*180./pi.le.val2)  .and.
     2465     &           (lonfi(ig)*180./pi.ge.val3)  .and.
     2466     &           (lonfi(ig)*180./pi.le.val4) ) then
     2467
     2468                IF ((phisfi(ig).le.val6).and.
     2469     &              (phisfi(ig).ge.val5)) then
     2470                    CALL RANDOM_NUMBER(val9)
     2471                    phisfi(ig)=val7+(val8-val7)*val9
     2472                ENDIF
     2473
     2474              ENDIF
     2475         enddo
     2476         phisold = phis
     2477c        update new phis and ps
     2478         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2479         r = 8.314511E+0 *1000.E+0/mugaz
     2480         temp=37.
     2481         do j=1,jjp1
     2482           do i=1,iip1
     2483              ps(i,j) = ps(i,j) *
     2484     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2485           enddo
     2486         enddo
     2487c       periodicity of surface ps in longitude
     2488         do j=1,jjp1
     2489           do i=1,iip1
     2490              ps(i,j) = ps(i,j) *
     2491     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2492           enddo
     2493         enddo
     2494c       periodicity of surface ps in longitude
     2495         do j=1,jjp1
     2496          ps(1,j) = ps(iip1,j)
     2497         end do
     2498         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2499
     2500c       fill_sp : fill SP with N2 ice and adjust phisfi (flat SP)
     2501c       -----------------------------------------------------
     2502        else if (modif(1:len_trim(modif)) .eq. 'fill_sp') then
     2503
     2504         ! def SP:
     2505         val2=-33.   !lat1
     2506         val3=50.    !lat2
     2507         val4=140.   ! lon1
     2508         val5=-155.  ! lon2
     2509         write(*,*) 'def SP :'
     2510         write(*,*) 'lat :',val2,val3
     2511         write(*,*) 'lon :',val4,'180 / -180',val5
     2512         write(*,*) 'choice phisfi ref SP (ex: -800):'
     2513607      read(*,*,iostat=ierr) val6
     2514         if(ierr.ne.0) goto 607
     2515
     2516         write(*,*) 're scale the pressure :'
     2517         r = 8.314511E+0 *1000.E+0/mugaz
     2518         temp=40.
     2519         write(*,*) 'r, sale height temperature =',r,temp,g
     2520         qsurf=0.       
     2521
     2522         write(*,*) 'latfi=',latfi
     2523         !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2524         write(*,*) 'phis=',phis
     2525         write(*,*) 'phisfi=',phisfi
     2526         !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2527         phisold = phis
     2528         write(*,*) 'phisold=',phisold
     2529         do ig=1,ngridmx
     2530               
     2531              write(*,*) 'lat=',latfi(ig)*180./pi
     2532              write(*,*) 'lon=',lonfi(ig)*180./pi
     2533              write(*,*) 'phisfi=',phisfi(ig)
     2534              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
     2535     &            (latfi(ig)*180./pi.le.val3)  .and.
     2536     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
     2537     &                       (lonfi(ig)*180./pi.le.val5))  .or.
     2538     &         ((lonfi(ig)*180./pi.ge.val4)  .and.
     2539     &                    (lonfi(ig)*180./pi.le.180.))) ) then
     2540                write(*,*) 'hello SP',phisfi(ig),ig
     2541                IF (phisfi(ig).le.val6) then
     2542                    qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
     2543                    phisfi(ig)=val6
     2544                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
     2545                ENDIF
     2546              ENDIF 
     2547         end do
     2548c        update new phis and ps
     2549         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2550         do j=1,jjp1
     2551           do i=1,iip1
     2552              ps(i,j) = ps(i,j) *
     2553     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2554           enddo
     2555         enddo
     2556c       periodicity of surface ps in longitude
     2557         do j=1,jjp1
     2558          ps(1,j) = ps(iip1,j)
     2559         end do
     2560
     2561c       bugch4 correct bug warm ch4 due to change of resolution
     2562c       -----------------------------------------------------
     2563        else if (modif(1:len_trim(modif)) .eq. 'bugch4') then
     2564         write(*,*) 'Ok we are going to try to solve this bug'
     2565         write(*,*) 'First extract a profile of tsoil and tsurf'
     2566         write(*,*) 'that you want to set as reference :'
     2567         write(*,*) 'tsoil_ref and tsurf_ref '
     2568         write(*,*) 'number of points to check '
     2569546      read(*,*,iostat=ierr) jref1
     2570         if(ierr.ne.0) goto 546
     2571
     2572         write(*,*) 'choice acceptable tsurf range for ch4: t1,t2:'
     2573547      read(*,*,iostat=ierr) val1,val2
     2574         if(ierr.ne.0) goto 547
     2575
     2576         ! Check tsurf at all locations
     2577         DO ig=1,ngridmx
     2578           IF (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
     2579     &                         qsurf(ig,igcm_n2).eq.0.) then
     2580             IF ((tsurf(ig).lt.val1)  .or.
     2581     &            (tsurf(ig).ge.val2))  then
     2582               write(*,*) 'Pb tsurf point ',ig
     2583
     2584               do jref=1,jref1
     2585                 IF ((ig-jref.ge.1).and.qsurf(ig,igcm_n2).eq.0..and.
     2586     &                 (qsurf(int(ig-jref),igcm_ch4_ice).gt.0.001).and.
     2587     &                  (tsurf(ig-jref).gt.val1
     2588     &                  .and.tsurf(ig-jref).lt.val2)) then
     2589                    tsurf(ig)=tsurf(int(ig-jref))
     2590                    DO l=1,nsoilmx
     2591                       tsoil(ig,l)=tsoil(int(ig-jref),l)
     2592                    ENDDO
     2593                    goto 548
     2594                 ELSEIF ((ig+jref.le.ngridmx).and.
     2595     &                 qsurf(ig,igcm_n2).eq.0..and.
     2596     &                 (qsurf(int(ig+jref),igcm_ch4_ice).gt.0.001).and.
     2597     &                  (tsurf(ig+jref).gt.val1
     2598     &                  .and.tsurf(ig+jref).lt.val2)) then
     2599                    tsurf(ig)=tsurf(int(ig+jref))
     2600                    DO l=1,nsoilmx
     2601                       tsoil(ig,l)=tsoil(int(ig+jref),l)
     2602                    ENDDO
     2603                    goto 548
     2604                 ENDIF
     2605               enddo
     2606548            write(*,*) 'found point with ch4'
     2607             ENDIF
     2608           ENDIF
     2609         END DO
     2610
     2611c       bugn2 correct bug warm n2 due to change of resolution
     2612c       -----------------------------------------------------
     2613        else if (modif(1:len_trim(modif)) .eq. 'bugn2') then
     2614         write(*,*) 'Ok we are going to try to solve this bug'
     2615         write(*,*) 'First extract a profile of tsoil and tsurf'
     2616         write(*,*) 'that you want to set as reference :'
     2617         write(*,*) 'tsoil_ref and tsurf_ref '
     2618         write(*,*) 'number of points to check '
     2619544      read(*,*,iostat=ierr) jref1
     2620         if(ierr.ne.0) goto 544
     2621
     2622         write(*,*) 'choice acceptable tsurf range for n2: t1,t2:'
     2623540      read(*,*,iostat=ierr) val1,val2
     2624         if(ierr.ne.0) goto 540
     2625
     2626         ! Check tsurf at all locations
     2627         DO ig=1,ngridmx
     2628           IF (qsurf(ig,1).gt.0.001) then
     2629             IF ((tsurf(ig).lt.val1)  .or.
     2630     &            (tsurf(ig).ge.val2))  then
     2631               write(*,*) 'Pb tsurf point ',ig
     2632
     2633                ! IF low topo et peu/pas de N2: add ch4, CO, N2
     2634               do jref=1,jref1
     2635                 IF ((ig-jref.ge.1).and.
     2636     &                 (qsurf(int(ig-jref),igcm_n2).gt.0.001).and.
     2637     &                  (tsurf(ig-jref).gt.val1
     2638     &                  .and.tsurf(ig-jref).lt.val2)) then
     2639                    tsurf(ig)=tsurf(int(ig-jref))
     2640                    DO l=1,nsoilmx
     2641                       tsoil(ig,l)=tsoil(int(ig-jref),l)
     2642                    ENDDO
     2643                    goto 541
     2644                 ELSEIF ((ig+jref.le.ngridmx).and.
     2645     &                 (qsurf(int(ig+jref),igcm_n2).gt.0.001).and.
     2646     &                  (tsurf(ig+jref).gt.val1
     2647     &                  .and.tsurf(ig+jref).lt.val2)) then
     2648                    tsurf(ig)=tsurf(int(ig+jref))
     2649                    DO l=1,nsoilmx
     2650                       tsoil(ig,l)=tsoil(int(ig+jref),l)
     2651                    ENDDO
     2652                    goto 541
     2653                 ENDIF
     2654               enddo
     2655541            write(*,*) 'found point with n2'
     2656             ENDIF
     2657           ENDIF
     2658         END DO
     2659
     2660         ! second tour
     2661!         DO ig=1,ngridmx
     2662!           IF (qsurf(ig,1).gt.0.001) then
     2663!             IF ((tsurf(ig).lt.val1)  .or.
     2664!     &            (tsurf(ig).ge.val2))  then
     2665!                ! IF low topo et peu/pas de N2: add ch4, CO, N2
     2666!                  IF ((ig-1.lt.1).or.(ig+1.gt.ngridmx)) then
     2667!                     write(*,*) 'pole : doing nothing'
     2668!                  ELSEIF (qsurf(ig-1,igcm_n2).gt.0.001) then
     2669!                    tsurf(ig)=tsurf(ig-1)
     2670!                    DO l=1,nsoilmx
     2671!                      tsoil(ig,l)=tsoil(ig-1,l)
     2672!                   ENDDO
     2673!                 ELSEIF (qsurf(ig+1,igcm_n2).gt.0.001) then
     2674!                   tsurf(ig)=tsurf(ig+1)
     2675!                   DO l=1,nsoilmx
     2676!                      tsoil(ig,l)=tsoil(ig+1,l)
     2677!                   ENDDO
     2678!                 ENDIF
     2679!            ENDIF
     2680!          ENDIF
     2681
     2682!        END DO
     2683
     2684c       flatnosp : flat topo outside a specific terrain (SP)
     2685c       -----------------------------------------------------
     2686        else if (modif(1:len_trim(modif)) .eq. 'flatnosp') then
     2687
     2688          write(*,*) 'Choice of topography (m) below which we keep ?'
     2689551       read(*,*,iostat=ierr) val
     2690          if(ierr.ne.0) goto 551
     2691          write(*,*) 'gravity g is : ',g         
     2692          geop=g*val
     2693          write(*,*) 'Choice of minimum amount of N2 ice we keep ?'
     2694552       read(*,*,iostat=ierr) val2
     2695          if(ierr.ne.0) goto 552
     2696       
     2697          write(*,*) 'phis=',phis
     2698          write(*,*) 'phisfi=',phisfi
     2699          do ig=1,ngridmx
     2700               
     2701              IF (   (qsurf(ig,1).le.val2)  .and.
     2702     &               (phisfi(ig).gt.geop)  ) then
     2703                write(*,*) 'hello SP',phisfi(ig),ig
     2704                phisfi(ig)=0.
     2705              ENDIF
     2706          end do
     2707
     2708          phisold = phis
     2709          write(*,*) 're scale the pressure :'
     2710          r = 8.314511E+0 *1000.E+0/mugaz
     2711          temp=40.
     2712          write(*,*) 'r, sale height temperature =',r,temp,g
     2713c         update new phis and ps
     2714          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2715          do j=1,jjp1
     2716            do i=1,iip1
     2717              ps(i,j) = ps(i,j) *
     2718     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2719            enddo
     2720          enddo
     2721c        periodicity of surface ps in longitude
     2722          do j=1,jjp1
     2723           ps(1,j) = ps(iip1,j)
     2724          end do
     2725
     2726c       flatregion : flat topo specific to region
     2727c       -----------------------------------------------------
     2728        else if (modif(1:len_trim(modif)) .eq. 'flatregion') then
     2729
     2730          write(*,*) 'Choice band : lat1 and lat2 ?'
     2731 553      read(*,*,iostat=ierr) val,val2
     2732          if(ierr.ne.0) goto 553
     2733          write(*,*) 'Choice band : lon1 and lon2 ?'
     2734 554      read(*,*,iostat=ierr) val3,val4
     2735          if(ierr.ne.0) goto 554
     2736          write(*,*) 'Choice of topography range ?'
     2737 555      read(*,*,iostat=ierr) val5,val6
     2738          if(ierr.ne.0) goto 555
     2739          write(*,*) 'Choice topo ?'
     2740 556      read(*,*,iostat=ierr) val7
     2741          if(ierr.ne.0) goto 556
     2742       
     2743          write(*,*) 'phis=',phis
     2744          write(*,*) 'phisfi=',phisfi
     2745          do ig=1,ngridmx
     2746             IF (   (latfi(ig)*180./pi.ge.val)  .and.
     2747     &              (latfi(ig)*180./pi.le.val2)  .and.
     2748     &              (lonfi(ig)*180./pi.ge.val3)  .and.
     2749     &              (lonfi(ig)*180./pi.le.val4)  ) then
     2750               
     2751              IF (   (phisfi(ig).ge.val5)  .and.
     2752     &               (phisfi(ig).le.val6)  ) then
     2753                write(*,*) 'hello topo',phisfi(ig),ig
     2754                phisfi(ig)=val7
     2755              ENDIF
     2756             ENDIF
     2757          end do
     2758
     2759          r = 8.314511E+0 *1000.E+0/mugaz
     2760          temp=40.
     2761          phisold = phis
     2762c         update new phis and ps
     2763          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2764          do j=1,jjp1
     2765            do i=1,iip1
     2766              ps(i,j) = ps(i,j) *
     2767     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2768            enddo
     2769          enddo
     2770c        periodicity of surface ps in longitude
     2771          do j=1,jjp1
     2772           ps(1,j) = ps(iip1,j)
     2773          end do
     2774
     2775c       changetopo 
     2776c       -----------------------------------------------------
     2777        else if (modif(1:len_trim(modif)) .eq. 'changetopo') then
     2778
     2779          write(*,*) 'Choice band : lat1 and lat2 ?'
     2780 573      read(*,*,iostat=ierr) val,val2
     2781          if(ierr.ne.0) goto 573
     2782          write(*,*) 'Choice band : lon1 and lon2 ?'
     2783 574      read(*,*,iostat=ierr) val3,val4
     2784          if(ierr.ne.0) goto 574
     2785          write(*,*) 'Choice of topography range ?'
     2786 575      read(*,*,iostat=ierr) val5,val6
     2787          if(ierr.ne.0) goto 575
     2788          write(*,*) 'Choice change topo factor?'
     2789 576      read(*,*,iostat=ierr) val7
     2790          if(ierr.ne.0) goto 576
     2791          write(*,*) 'Choice change topo add?'
     2792 577      read(*,*,iostat=ierr) val8
     2793          if(ierr.ne.0) goto 577
     2794       
     2795          write(*,*) 'phis=',phis
     2796          write(*,*) 'phisfi=',phisfi
     2797          do ig=1,ngridmx
     2798             IF (   (latfi(ig)*180./pi.ge.val)  .and.
     2799     &              (latfi(ig)*180./pi.le.val2)  .and.
     2800     &              (lonfi(ig)*180./pi.ge.val3)  .and.
     2801     &              (lonfi(ig)*180./pi.le.val4)  ) then
     2802               
     2803              IF (   (phisfi(ig).ge.val5)  .and.
     2804     &               (phisfi(ig).le.val6)  ) then
     2805                write(*,*) 'hello topo',phisfi(ig),ig
     2806                phisfi(ig)=phisfi(ig)*val7
     2807                phisfi(ig)=phisfi(ig)+val8
     2808              ENDIF
     2809             ENDIF
     2810          end do
     2811
     2812          r = 8.314511E+0 *1000.E+0/mugaz
     2813          temp=40.
     2814          phisold = phis
     2815c         update new phis and ps
     2816          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2817          do j=1,jjp1
     2818            do i=1,iip1
     2819              ps(i,j) = ps(i,j) *
     2820     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2821            enddo
     2822          enddo
     2823c        periodicity of surface ps in longitude
     2824          do j=1,jjp1
     2825           ps(1,j) = ps(iip1,j)
     2826          end do
     2827
     2828c       corrtopo: corr topo specific bug
     2829c       -----------------------------------------------------
     2830        else if (modif(1:len_trim(modif)) .eq. 'corrtopo') then
     2831
     2832          ! get min max lon
     2833          print*, latfi*180/pi
     2834          print*, '***************************'
     2835          print*, '***************************'
     2836          print*, '***************************'
     2837          print*, '***************************'
     2838          print*, '***************************'
     2839          print*, lonfi*180/pi
     2840          print*, 'iip1=',iip1
     2841          do ig=2,ngridmx-1
     2842             IF  (lonfi(ig)*180./pi.eq.-180.) THEN
     2843                 print*, lonfi(ig),lonfi(ig+iip1-2)
     2844                 phisfi(ig)=(phisfi(ig+1)+phisfi(ig+iip1))/2
     2845             ENDIF
     2846          enddo
     2847          phisold = phis
     2848          r = 8.314511E+0 *1000.E+0/mugaz
     2849          temp=40.
     2850c         update new phis and ps
     2851          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2852          do j=1,jjp1
     2853            do i=1,iip1
     2854              ps(i,j) = ps(i,j) *
     2855     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2856            enddo
     2857          enddo
     2858c         periodicity of surface ps in longitude
     2859          do j=1,jjp1
     2860           ps(1,j) = ps(iip1,j)
     2861          end do
     2862
     2863c       asymtopo:
     2864c       -----------------------------------------------------
     2865        else if (modif(1:len_trim(modif)) .eq. 'asymtopo') then
     2866
     2867          ! get diff altitude
     2868          write(*,*) 'Diff N-S altitude in km (pos = N higher)  ?'
     2869 591      read(*,*,iostat=ierr) val
     2870          if(ierr.ne.0) goto 591
     2871          write(*,*) 'Coeff smoot topo (small = smoother)  ?'
     2872 592      read(*,*,iostat=ierr) val2
     2873          if(ierr.ne.0) goto 592
     2874
     2875          do ig=1,ngridmx
     2876           phisfi(ig)=phisfi(ig)+val*1000.*g*tanh(latfi(ig)*180/pi*val2)
     2877          enddo
     2878          phisold = phis
     2879          r = 8.314511E+0 *1000.E+0/mugaz
     2880          temp=40.
     2881c         update new phis and ps
     2882          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     2883          do j=1,jjp1
     2884            do i=1,iip1
     2885              ps(i,j) = ps(i,j) *
     2886     .          exp((phisold(i,j)-phis(i,j))/(temp * r))
     2887            enddo
     2888          enddo
     2889c         periodicity of surface ps in longitude
     2890          do j=1,jjp1
     2891           ps(1,j) = ps(iip1,j)
     2892          end do
     2893
     2894c       seticeNH : set the ices in the SP crater with NH topo
     2895c       -----------------------------------------------------
     2896        else if (modif(1:len_trim(modif)) .eq. 'seticeNH') then
     2897
     2898         open(333,file='./tsoil_180_30',form='formatted')
     2899         do i=1,22
     2900           read(333,*) levsoil(i), tempsoil(i)
     2901         enddo
     2902         close(333)
     2903         open(334,file='./tsurf_180_30',form='formatted')
     2904         read(334,*) val
     2905         close(334)
     2906
     2907         write(*,*) 'Tsoil and tsurf ref are:'
     2908         write(*,*) tempsoil
     2909         write(*,*) 'tsurf=',val
     2910
     2911         ! def SP:
     2912         val2=-43.   !lat1
     2913         val3=51.    !lat2
     2914         val4=143.   ! lon1
     2915         val5=-158.  ! lon2
     2916         val6=-150   ! phisfi min
     2917
     2918         ! Rm all N2 ice outside SP
     2919         DO ig=1,ngridmx
     2920           ! IF inside SP
     2921           IF (   (latfi(ig)*180./pi.ge.val2)  .and.
     2922     &            (latfi(ig)*180./pi.le.val3)  .and.
     2923     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
     2924     &                       (lonfi(ig)*180./pi.le.val5))  .or.
     2925     &         ((lonfi(ig)*180./pi.ge.val4)  .and.
     2926     &                       (lonfi(ig)*180./pi.le.180.))) ) then
     2927
     2928                ! IF low topo et peu/pas de N2: add ch4, CO, N2
     2929                IF ((phisfi(ig).le.val6).and.
     2930     &                              (qsurf(ig,igcm_n2).le.500)) then
     2931                    qsurf(ig,igcm_n2)=1000.
     2932                    qsurf(ig,igcm_ch4_ice)=1000.
     2933                    qsurf(ig,igcm_co_ice)=1000.
     2934                    tsurf(ig)=val
     2935                    DO l=1,nsoilmx
     2936                       tsoil(ig,l)=tempsoil(l)
     2937                    ENDDO
     2938                ENDIF
     2939
     2940                ! IF high topo : rm N2
     2941                IF ( (qsurf(ig,igcm_n2).ge.20.).and.
     2942     &                                 (phisfi(ig).ge.val6)  ) then
     2943                      qsurf(ig,igcm_n2)=0.
     2944
     2945                ENDIF
     2946                ! Mais on garde ch4 et CO en prenant la meme quantite
     2947                ! qu'une autre latitude
     2948                IF ( (qsurf(ig,igcm_ch4_ice).ge.20.)  .and.
     2949     &                                 (phisfi(ig).ge.val6)  ) then
     2950                      jref=int(ig/jjp1)+2
     2951                      qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
     2952                ENDIF
     2953                IF ( (qsurf(ig,igcm_co_ice).ge.20.)  .and.
     2954     &                                 (phisfi(ig).ge.val6)  ) then
     2955                      jref=int(ig/jjp1)+2
     2956                      qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
     2957                ENDIF
     2958
     2959           ELSE   ! Outside SP
     2960
     2961                IF (qsurf(ig,igcm_n2).ge.20.) then
     2962                      qsurf(ig,igcm_n2)=0.
     2963                ENDIF
     2964
     2965                IF (qsurf(ig,igcm_ch4_ice).ge.20.) then
     2966                      jref=int(ig/jjp1)+2
     2967                      qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
     2968                ENDIF
     2969
     2970                IF (qsurf(ig,igcm_co_ice).ge.20.) then
     2971                      jref=int(ig/jjp1)+2
     2972                      qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
     2973                ENDIF
     2974
    20792975           ENDIF
     2976       
     2977         END DO
     2978
     2979c       'nomountain '
     2980c       --------------
     2981        else if (modif(1:len_trim(modif)) .eq. 'nomountain') then
     2982             do j=1,jjp1
     2983               do i=1,iip1
     2984                 if (phis(i,j).gt.0.1) then
     2985                       phis(i,j) = 0.
     2986                       ps(i,j)=ps(iim/4,j)    ! assuming no topo here
     2987                 endif
     2988               enddo
     2989             enddo
     2990            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     2991
     2992c       'relief'
     2993c       --------------
     2994        else if (modif(1:len_trim(modif)) .eq.'relief') then
     2995          write(*,*) "add a circular mountain/crater"
     2996          write(*,*) 'Center: lat lon  ?'
     2997 707      read(*,*,iostat=ierr) lat0, lon0
     2998          if(ierr.ne.0) goto 707
     2999          if(lon0.gt.180.) lon0=lon0-360.
     3000          lat0 = lat0 * pi/180.
     3001          lon0 = lon0 * pi/180.
     3002
     3003          DO ig=1,ngridmx
     3004             if(abs(latfi(ig)-lat0).lt.pi/float(jjm) ) then
     3005               if(abs(lonfi(ig)-lon0).lt.2*pi/float(iim) ) then
     3006                  ig0 = ig
     3007               end if
     3008             end if
     3009          END DO
     3010          write(*,*) "Reference Point to be used:"
     3011          write(*,*) 'ig0',ig0
     3012          write(*,*) 'lat=',latfi(ig0)*180./pi
     3013          write(*,*) 'lon=',lonfi(ig0)*180./pi
     3014
     3015          write(*,*) 'radius (km), height (km) negative if crater ?'
     3016 508      read(*,*,iostat=ierr) radm, height
     3017          if(ierr.ne.0) goto 508
     3018
     3019          write(*,*) 'Sale height temperature (ex:38) ?'
     3020 509      read(*,*,iostat=ierr) temp
     3021          if(ierr.ne.0) goto 509
     3022
     3023          r = 8.314511E+0 *1000.E+0/mugaz
     3024          do j=1,jjp1
     3025            do i=1,iip1
     3026               dist= dist_pluto(lat0,lon0,rlatu(j),rlonv(i))
     3027               if (dist.le.radm) then
     3028                  phisprev= phis(i,j)
     3029                  phis(i,j)=phisprev+1000.*g*height*(radm-dist)/radm
     3030                  write(*,*) 'lat=',rlatu(j)*180./pi
     3031                  write(*,*) 'lon=',rlonv(i)*180./pi
     3032                  write(*,*) 'dist', dist
     3033                  write(*,*) 'z(m)=', phis(i,j)/g
     3034                  ps(i,j) = ps(i,j) *
     3035     .            exp((phis(i,j))/(temp * r))
     3036               end if
     3037            end do
     3038          end do
     3039
     3040c         periodicity of surface ps in longitude
     3041          do j=1,jjp1
     3042            ps(1,j) = ps(iip1,j)
     3043          end do
     3044          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     3045
     3046c       subsoil_n2 : choice of thermal inertia values for N2 only
     3047c       ----------------------------------------------------------------
     3048        else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then
     3049
     3050         write(*,*) 'New value for subsoil thermal inertia  ?'
     3051 102     read(*,*,iostat=ierr) ith_bb
     3052         if(ierr.ne.0) goto 102
     3053         write(*,*) 'thermal inertia (new value):',ith_bb
     3054
     3055         write(*,*)'At which depth (in m.) does the ice layer begin?'
     3056         write(*,*)'(currently, the deepest soil layer extends down to:'
     3057     &              ,layer(1),' - ',layer(nsoilmx),')'
     3058         write(*,*)'write 0 for uniform value for all subsurf levels?'
     3059         ierr=1
     3060         do while (ierr.ne.0)
     3061          read(*,*,iostat=ierr) val2
     3062          write(*,*)'val2:',val2,'ierr=',ierr
     3063          if (ierr.eq.0) then ! got a value, but do a sanity check
     3064            if(val2.gt.layer(nsoilmx)) then
     3065             write(*,*)'Depth should be less than ',layer(nsoilmx)
     3066             ierr=1
     3067            endif
     3068            if(val2.lt.layer(1)) then
     3069              if(val2.eq.0) then
     3070               write(*,*)'Thermal inertia set for all subsurface layers'
     3071               ierr=0
     3072              else 
     3073               write(*,*)'Depth should be more than ',layer(1)
     3074               ierr=1
     3075              endif
     3076            endif
     3077          endif
     3078         enddo ! of do while
     3079
     3080         ! find the reference index iref the depth corresponds to
     3081         if(val2.eq.0) then
     3082           iref=1
     3083           write(*,*)'Level selected is first level: ',layer(iref),' m'
     3084         else
     3085           do isoil=1,nsoilmx-1
     3086            if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     3087     &       then
     3088             iref=isoil
     3089             write(*,*)'Level selected : ',layer(isoil),' m'
     3090             exit
     3091            endif
     3092           enddo
     3093         endif
     3094
     3095         DO ig=1,ngridmx
     3096             DO j=iref,nsoilmx
     3097c                if((qsurf(ig,igcm_ch4_ice).lt.1.).and.
     3098c     &                         (qsurf(ig,igcm_co_ice).lt.1.))then
     3099                if (qsurf(ig,igcm_n2).gt.0.001) then
     3100                   ithfi(ig,j)=ith_bb
     3101                endif
     3102             ENDDO
     3103         ENDDO
     3104
     3105         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     3106
     3107c       subsoil_ch4 : choice of thermal inertia values for CH4 only
     3108c       ----------------------------------------------------------------
     3109        else if (modif(1:len_trim(modif)) .eq. 'subsoil_ch4') then
     3110
     3111         write(*,*) 'New value for subsoil thermal inertia  ?'
     3112 103     read(*,*,iostat=ierr) ith_bb
     3113         if(ierr.ne.0) goto 103
     3114         write(*,*) 'thermal inertia (new value):',ith_bb
     3115
     3116         write(*,*)'At which depth (in m.) does the ice layer begin?'
     3117         write(*,*)'(currently, the deepest soil layer extends down to:'
     3118     &              ,layer(1),' - ',layer(nsoilmx),')'
     3119         write(*,*)'write 0 for uniform value for all subsurf levels?'
     3120         ierr=1
     3121         do while (ierr.ne.0)
     3122          read(*,*,iostat=ierr) val2
     3123          write(*,*)'val2:',val2,'ierr=',ierr
     3124          if (ierr.eq.0) then ! got a value, but do a sanity check
     3125            if(val2.gt.layer(nsoilmx)) then
     3126              write(*,*)'Depth should be less than ',layer(nsoilmx)
     3127              ierr=1
     3128            endif
     3129            if(val2.lt.layer(1)) then
     3130              if(val2.eq.0) then
     3131               write(*,*)'Thermal inertia set for all subsurface layers'
     3132               ierr=0
     3133              else 
     3134               write(*,*)'Depth should be more than ',layer(1)
     3135               ierr=1
     3136              endif
     3137            endif
     3138          endif
     3139         enddo ! of do while
     3140
     3141         ! find the reference index iref the depth corresponds to
     3142         if(val2.eq.0) then
     3143           iref=1
     3144           write(*,*)'Level selected is first level: ',layer(iref),' m'
     3145         else
     3146           do isoil=1,nsoilmx-1
     3147            if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     3148     &       then
     3149             iref=isoil
     3150             write(*,*)'Level selected : ',layer(isoil),' m'
     3151             exit
     3152            endif
     3153           enddo
     3154         endif
     3155
     3156         DO ig=1,ngridmx
     3157             DO j=iref,nsoilmx
     3158                if (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
     3159     &                            qsurf(ig,igcm_n2).lt.0.001) then
     3160                   ithfi(ig,j)=ith_bb
     3161                endif
     3162             ENDDO
     3163         ENDDO
     3164
     3165         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     3166
     3167c       subsoil_alb : choice of thermal inertia values from albedo
     3168c       ----------------------------------------------------------------
     3169        else if (modif(1:len_trim(modif)) .eq. 'subsoil_alb') then
     3170
     3171         write(*,*) 'Choice range albedo for defining TI'
     3172 145     read(*,*,iostat=ierr) val,val2
     3173         if(ierr.ne.0) goto 145
     3174         write(*,*) 'New value for subsoil thermal inertia  ?'
     3175 144     read(*,*,iostat=ierr) ith_bb
     3176         if(ierr.ne.0) goto 144
     3177         write(*,*) 'thermal inertia (new value):',ith_bb
     3178
     3179         write(*,*)'At which depth (in m.) does the ice layer begin?'
     3180         write(*,*)'(currently, the deepest soil layer extends down to:'
     3181     &              ,layer(1),' - ',layer(nsoilmx),')'
     3182         write(*,*)'write 0 for uniform value for all subsurf levels?'
     3183         ierr=1
     3184         do while (ierr.ne.0)
     3185          read(*,*,iostat=ierr) val3
     3186          if (ierr.eq.0) then ! got a value, but do a sanity check
     3187            if(val3.gt.layer(nsoilmx)) then
     3188              write(*,*)'Depth should be less than ',layer(nsoilmx)
     3189              ierr=1
     3190            endif
     3191            if(val3.lt.layer(1)) then
     3192              if(val3.eq.0) then
     3193               write(*,*)'Thermal inertia set for all subsurface layers'
     3194               ierr=0
     3195              else 
     3196               write(*,*)'Depth should be more than ',layer(1)
     3197               ierr=1
     3198              endif
     3199            endif
     3200          endif
     3201         enddo ! of do while
     3202
     3203         ! find the reference index iref the depth corresponds to
     3204         if(val3.eq.0) then
     3205           iref=1
     3206           write(*,*)'Level selected is first level: ',layer(iref),' m'
     3207         else
     3208           do isoil=1,nsoilmx-1
     3209            if ((val3.gt.layer(isoil)).and.(val3.lt.layer(isoil+1)))
     3210     &       then
     3211             iref=isoil+1
     3212             write(*,*)'Level selected : ',layer(isoil+1),' m'
     3213             exit
     3214            endif
     3215           enddo
     3216         endif
     3217
     3218         DO ig=1,ngridmx
     3219             IF (   (albfi(ig).ge.val)  .and.
     3220     &              (albfi(ig).le.val2)  ) then
     3221               DO j=iref,nsoilmx
     3222                    ithfi(ig,j)=ith_bb
     3223               ENDDO
     3224             ENDIF
     3225         ENDDO
     3226
     3227         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     3228
     3229c       subsoil_all : choice of thermal inertia values (global)
     3230c       ----------------------------------------------------------------
     3231        else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then
     3232
     3233         write(*,*) 'New value for subsoil thermal inertia  ?'
     3234 104     read(*,*,iostat=ierr) ith_bb
     3235         if(ierr.ne.0) goto 104
     3236         write(*,*) 'thermal inertia (new value):',ith_bb
     3237
     3238         write(*,*)'At which depth (in m.) does the ice layer begin?'
     3239         write(*,*)'(currently, the deepest soil layer extends down to:'
     3240     &              ,layer(1),' - ',layer(nsoilmx),')'
     3241         write(*,*)'write 0 for uniform value for all subsurf levels?'
     3242         ierr=1
     3243         do while (ierr.ne.0)
     3244          read(*,*,iostat=ierr) val2
     3245          write(*,*)'val2:',val2,'ierr=',ierr
     3246          if (ierr.eq.0) then ! got a value, but do a sanity check
     3247            if(val2.gt.layer(nsoilmx)) then
     3248              write(*,*)'Depth should be less than ',layer(nsoilmx)
     3249              ierr=1
     3250            endif
     3251            if(val2.lt.layer(1)) then
     3252              if(val2.eq.0) then
     3253               write(*,*)'Thermal inertia set for all subsurface layers'
     3254               ierr=0
     3255              else 
     3256               write(*,*)'Depth should be more than ',layer(1)
     3257               ierr=1
     3258              endif
     3259            endif
     3260          endif
     3261         enddo ! of do while
     3262
     3263         ! find the reference index iref the depth corresponds to
     3264         if(val2.eq.0) then
     3265           iref=1
     3266           write(*,*)'Level selected is first level: ',layer(iref),' m'
     3267         else
     3268           do isoil=1,nsoilmx-1
     3269             !write(*,*)'isoil, ',isoil,val2
     3270             !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
     3271            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     3272     &       then
     3273             iref=isoil+1
     3274             write(*,*)'Level selected : ',layer(isoil+1),' m'
     3275             exit
     3276            endif
     3277           enddo
     3278         endif
     3279
     3280         DO ig=1,ngridmx
     3281             DO j=iref,nsoilmx
     3282                   ithfi(ig,j)=ith_bb
     3283             ENDDO
     3284         ENDDO
     3285
     3286         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     3287
     3288c       diurnal_TI : choice of thermal inertia values (global)
     3289c       ----------------------------------------------------------------
     3290        else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then
     3291
     3292         write(*,*) 'New value for diurnal thermal inertia  ?'
     3293 106     read(*,*,iostat=ierr) ith_bb
     3294         if(ierr.ne.0) goto 106
     3295         write(*,*) 'Diurnal thermal inertia (new value):',ith_bb
     3296
     3297         write(*,*)'At which depth (in m.) does the ice layer ends?'
     3298         write(*,*)'(currently, the soil layer 1 and nsoil are:'
     3299     &              ,layer(1),' - ',layer(nsoilmx),')'
     3300         ierr=1
     3301         do while (ierr.ne.0)
     3302          read(*,*,iostat=ierr) val2
     3303          write(*,*)'val2:',val2,'ierr=',ierr
     3304          if (ierr.eq.0) then ! got a value, but do a sanity check
     3305            if(val2.gt.layer(nsoilmx)) then
     3306              write(*,*)'Depth should be less than ',layer(nsoilmx)
     3307              ierr=1
     3308            endif
     3309            if(val2.lt.layer(1)) then
     3310              write(*,*)'Depth should be more than ',layer(1)
     3311              ierr=1
     3312            endif
     3313          endif
     3314         enddo ! of do while
     3315
     3316         ! find the reference index iref the depth corresponds to
     3317         do isoil=1,nsoilmx-1
     3318            !write(*,*)'isoil, ',isoil,val2
     3319            !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
     3320            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     3321     &       then
     3322             iref=isoil+1
     3323             write(*,*)'Level selected : ',layer(isoil+1),' m'
     3324             exit
     3325            endif
     3326         enddo
     3327
     3328         DO ig=1,ngridmx
     3329             DO j=1,iref
     3330                   ithfi(ig,j)=ith_bb
     3331             ENDDO
     3332         ENDDO
     3333
     3334         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     3335
     3336c       Choice surface temperature
     3337c       -----------------------------------------------------
     3338        else if (modif(1:len_trim(modif)) .eq. 'initsurf') then
     3339
     3340          write(*,*) 'New value for initial surface temperature  ?'
     3341 105      read(*,*,iostat=ierr) tsurf_bb
     3342          if(ierr.ne.0) goto 105
     3343          write(*,*) 'new surface temperature (K):',tsurf_bb
     3344          DO ig=1,ngridmx
     3345                   tsurf(ig)=tsurf_bb
     3346          ENDDO
     3347
     3348c       Modify surface temperature (for GCM start)
     3349c       -----------------------------------------------------
     3350        else if (modif(1:len_trim(modif)) .eq. 'modtsurf') then
     3351
     3352          write(*,*) 'Choice band : lat1 and lat2 ?'
     3353 455      read(*,*,iostat=ierr) val,val2
     3354          if(ierr.ne.0) goto 455
     3355          write(*,*) 'Choice band : lon1 and lon2 ?'
     3356 456      read(*,*,iostat=ierr) val4,val5
     3357          if(ierr.ne.0) goto 456
     3358          write(*,*) 'Choice topo min max '
     3359 473      read(*,*,iostat=ierr) val9,val10
     3360          if(ierr.ne.0) goto 473
     3361          write(*,*) 'Choice tsurf min max '
     3362 474      read(*,*,iostat=ierr) val11,val12
     3363          if(ierr.ne.0) goto 474
     3364          write(*,*) 'Choice multiplicative factor'
     3365 457      read(*,*,iostat=ierr) val3
     3366          if(ierr.ne.0) goto 457
     3367          write(*,*) 'Choice additional tsurf'
     3368 458      read(*,*,iostat=ierr) val6
     3369          if(ierr.ne.0) goto 458
     3370          write(*,*) 'Choice multiplicative factor tsoil'
     3371 459      read(*,*,iostat=ierr) val7
     3372          if(ierr.ne.0) goto 459
     3373          write(*,*) 'Choice additional tsoil'
     3374 469      read(*,*,iostat=ierr) val8
     3375          if(ierr.ne.0) goto 469
     3376
     3377          DO ig=1,ngridmx
     3378             IF (   (latfi(ig)*180./pi.ge.val)  .and.
     3379     &              (latfi(ig)*180./pi.le.val2)  .and.
     3380     &              (lonfi(ig)*180./pi.ge.val4)  .and.
     3381     &              (lonfi(ig)*180./pi.le.val5)  .and.
     3382     &              (phisfi(ig).ge.val9)  .and.
     3383     &              (phisfi(ig).lt.val10)  .and.
     3384     &              (tsurf(ig).ge.val11) .and.
     3385     &              (tsurf(ig).lt.val12)  ) then
     3386
     3387                    tsurf(ig)=tsurf(ig)*val3
     3388                    tsurf(ig)=tsurf(ig)+val6
     3389                    DO l=1,nsoilmx
     3390                       tsoil(ig,l)=tsoil(ig,l)*val7
     3391                       tsoil(ig,l)=tsoil(ig,l)+val8
     3392                    ENDDO
     3393             ENDIF
     3394          ENDDO
     3395
     3396c       copy latitudes tsurf / tsoil
     3397c       -----------------------------------------------------
     3398        else if (modif(1:len_trim(modif)) .eq. 'copylat') then
     3399           
     3400          write(*,*) 'all latitudes : ',rlatu*180./pi
     3401          write(*,*) 'Choice band to be modified lat1 ?'
     3402 465      read(*,*,iostat=ierr) val
     3403          if(ierr.ne.0) goto 465
     3404          write(*,*) 'Choice band copied lat2 ?'
     3405 466      read(*,*,iostat=ierr) val2
     3406          if(ierr.ne.0) goto 466
     3407          write(*,*) 'Choice multiplicative factor'
     3408 467      read(*,*,iostat=ierr) val3
     3409          if(ierr.ne.0) goto 467
     3410          write(*,*) 'Choice additional tsurf'
     3411 468     read(*,*,iostat=ierr) val4
     3412          if(ierr.ne.0) goto 468
     3413
     3414          j=1
     3415          DO ig=1,ngridmx
     3416             IF ((latfi(ig)*180./pi.lt.val2+180./iip1) .and.
     3417     &          (latfi(ig)*180./pi.gt.val2-180./iip1))  then
     3418                randtab(j)=ig
     3419                j=j+1
     3420             ENDIF
     3421          ENDDO
     3422          print*,j, ' latitudes found'
     3423          k=1
     3424          DO ig=1,ngridmx
     3425             IF ((latfi(ig)*180./pi.lt.val+180./iip1) .and.
     3426     &          (latfi(ig)*180./pi.gt.val-180./iip1))  then
     3427                tsurf(ig)=tsurf(randtab(k))*val3
     3428                tsurf(ig)=tsurf(ig)+val4
     3429                DO l=1,nsoilmx
     3430                       tsoil(ig,l)=tsoil(randtab(k),l)*val3
     3431                       tsoil(ig,l)=tsoil(ig,l)+val4
     3432                ENDDO
     3433                k=k+1
     3434             ENDIF
     3435          ENDDO
     3436          print*,k, ' latitudes copied'
     3437          IF (j.ne.k) stop
     3438
     3439c       copy longitudes
     3440c       -----------------------------------------------------
     3441        else if (modif(1:len_trim(modif)) .eq. 'copylon') then
     3442           
     3443          write(*,*) 'all longitudes : ',rlonu*180./pi
     3444          write(*,*) 'Choice band to be modified lon1 ?'
     3445 475      read(*,*,iostat=ierr) val
     3446          if(ierr.ne.0) goto 475
     3447          write(*,*) 'Choice band copied lon2 ?'
     3448 476      read(*,*,iostat=ierr) val2
     3449          if(ierr.ne.0) goto 476
     3450
     3451          j=1
     3452          DO ig=2,ngridmx-1
     3453             IF ((lonfi(ig)*180./pi.lt.183.) .and.
     3454     &          (lonfi(ig)*180./pi.gt.180.))  then
     3455                randtab(j)=ig
     3456                j=j+1
     3457                print*, 'lon = ',lonfi(ig)
     3458             ENDIF
     3459          ENDDO
     3460          print*,j, ' longitudes found'
     3461          k=1
     3462          DO ig=2,ngridmx-1
     3463             IF ((lonfi(ig)*180./pi.lt.180.) .and.
     3464     &          (lonfi(ig)*180./pi.gt.179.))  then
     3465                tsurf(ig)=tsurf(randtab(k))
     3466                DO l=1,nsoilmx
     3467                       tsoil(ig,l)=tsoil(randtab(k),l)
     3468                ENDDO
     3469                qsurf(ig,1)=qsurf(randtab(k),1) 
     3470                phisfi(ig)=phisfi(randtab(k)) 
     3471                k=k+1
     3472             ENDIF
     3473          ENDDO
     3474          print*,k, ' longitudes copied'
     3475          IF (j.ne.k) stop
     3476          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
     3477
     3478c       copy latlon
     3479c       -----------------------------------------------------
     3480        else if (modif(1:len_trim(modif)) .eq. 'copylatlon') then
     3481           
     3482          write(*,*) 'all longitudes : ',rlonu*180./pi
     3483          write(*,*) 'Choice lat/lon to be copied ?'
     3484 495      read(*,*,iostat=ierr) val,val2
     3485          if(ierr.ne.0) goto 495
     3486          write(*,*) 'Choice indx lon1 lon2 for band ?'
     3487 496      read(*,*,iostat=ierr) val3,val4
     3488          if(ierr.ne.0) goto 496
     3489          write(*,*) 'Choice latitude indx range where we copy ?'
     3490 497      read(*,*,iostat=ierr) val5,val6
     3491          if(ierr.ne.0) goto 497
     3492
     3493          ! choice coordinate
     3494          DO ig=2,ngridmx-1
     3495             IF ( (lonfi(ig)*180./pi.gt.val2-1) .and.
     3496     &            (lonfi(ig)*180./pi.lt.val2+1) .and.
     3497     &            (latfi(ig)*180./pi.lt.val+1)  .and.
     3498     &            (latfi(ig)*180./pi.gt.val-1) ) then
     3499              ig0=ig
     3500              print*,'lat/lon=',latfi(ig)*180./pi,lonfi(ig)*180./pi,ig0
     3501             ENDIF
     3502          ENDDO
     3503
     3504          DO ig=2,ngridmx-1
     3505             IF ( (lonfi(ig)*180./pi.lt.val4) .and.
     3506     &            (lonfi(ig)*180./pi.gt.val3) .and.
     3507     &            (latfi(ig)*180./pi.lt.val6) .and.
     3508     &            (latfi(ig)*180./pi.gt.val5) .and.
     3509     &            (qsurf(ig,1).lt.0.001) ) then
     3510                tsurf(ig)=tsurf(ig0)
     3511                print*,'corrd=',latfi(ig)*180./pi,lonfi(ig)*180./pi
     3512                DO l=1,nsoilmx
     3513                       tsoil(ig,l)=tsoil(ig0,l)
     3514                ENDDO
     3515             ENDIF
     3516          ENDDO
     3517
     3518c       Choice surface temperature depending on N2 ice distribution
     3519c       -----------------------------------------------------
     3520        else if (modif(1:len_trim(modif)) .eq. 'tsurfice') then
     3521
     3522          write(*,*) 'Initial soil and surface temp below n2 ?'
     3523 107      read(*,*,iostat=ierr) tsurf_bb
     3524          if(ierr.ne.0) goto 107
     3525          write(*,*) 'new surface/soil temp N2 (K):',tsurf_bb
     3526          write(*,*) 'Initial soil and surface temp when no n2 ?'
     3527 108      read(*,*,iostat=ierr) tsurf_bb2
     3528          if(ierr.ne.0) goto 108
     3529          write(*,*) 'new surface/soil temp when no n2 (K):',tsurf_bb2
     3530          DO ig=1,ngridmx
     3531                   if (qsurf(ig,igcm_n2).gt.0.001) then
     3532                      tsurf(ig)=tsurf_bb
     3533                      do l=1,nsoilmx
     3534                         tsoil(ig,l) = tsurf_bb
     3535                      end do
     3536                   else if (tsurf_bb2.ne.0) then
     3537                      tsurf(ig)=tsurf_bb2
     3538                      do l=1,nsoilmx
     3539                         tsoil(ig,l) = tsurf_bb2
     3540                      end do
     3541                   endif
     3542          ENDDO
     3543
     3544c       read an albedo map
     3545c       -----------------------------------------------------
     3546        else if (modif(1:len_trim(modif)) .eq. 'albedomap') then
     3547           
     3548          ! Get field 2D
     3549          fichnom = 'albedo.nc'
     3550          ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
     3551          IF (ierr.NE.NF_NOERR) THEN
     3552            write(6,*)' Problem opening albedo file:',fichnom
     3553            write(6,*)' ierr = ', ierr
     3554            CALL ABORT
     3555          ENDIF
     3556
     3557          ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"),
     3558     &                                  nvarid_input)
     3559          IF (ierr .NE. NF_NOERR) THEN
     3560           PRINT*, "Could not find asked variable in albedo.nc"
     3561           CALL abort
     3562          ENDIF
     3563
     3564#ifdef NC_DOUBLE
     3565          ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
     3566     &                                 field_input)
     3567#else
     3568          ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
     3569#endif
     3570          IF (ierr .NE. NF_NOERR) THEN
     3571           PRINT*, "Could not get asked variable in albedo.nc"
     3572           CALL abort
     3573          ELSE
     3574           PRINT*, "Got variable in albedo.nc"
     3575          ENDIF
     3576
     3577          DO ig=1,ngridmx
     3578                      albfi(ig)=field_input(ig)
     3579          ENDDO
     3580
     3581c       slope : add slope on all longitudes
     3582c       -----------------------------------------------------
     3583        else if (modif(1:len_trim(modif)) .eq. 'slope') then
     3584
     3585         write(*,*) 'choice latitude alt minimum & altitude value (m):'
     3586603      read(*,*,iostat=ierr) val1,val2
     3587         if(ierr.ne.0) goto 603
     3588         write(*,*) 'choice latitude alt maximum & altitude value (m):'
     3589604      read(*,*,iostat=ierr) val3,val4
     3590         if(ierr.ne.0) goto 604
     3591
     3592         write(*,*) 're scale the pressure :'
     3593         r = 8.314511E+0 *1000.E+0/mugaz
     3594         temp=40.
     3595         write(*,*) 'r, sale height temperature =',r,temp
     3596                 
     3597         do j=1,jjp1
     3598           do i=1,iip1
     3599              phisprev= phis(i,j)
     3600              IF ( (  (rlatu(j)*180./pi.ge.val1)  .and.
     3601     &            (rlatu(j)*180./pi.le.val3) ) .or.
     3602     &          (  (rlatu(j)*180./pi.le.val1)  .and.
     3603     &            (rlatu(j)*180./pi.ge.val3) )) then
     3604
     3605                    val5=abs(val2-val4)/abs(val1-val3)*
     3606     &                           abs(val1-rlatu(j)*180./pi)+val2
     3607                    phis(i,j)=val5*g
     3608                    ps(i,j) = ps(i,j) *
     3609     .              exp((phisprev-phis(i,j))/(temp * r))
     3610              ENDIF 
     3611           end do
     3612         end do
     3613c        periodicity of surface ps in longitude
     3614         do j=1,jjp1
     3615          ps(1,j) = ps(iip1,j)
     3616         end do
     3617         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     3618
     3619c       digsp : change topography of SP
     3620c       -----------------------------------------------------
     3621        else if (modif(1:len_trim(modif)) .eq. 'digsp') then
     3622
     3623         write(*,*) 'choice latitudes:'
     3624517      read(*,*,iostat=ierr) val1,val2
     3625         if(ierr.ne.0) goto 517
     3626         write(*,*) 'choice longitudes:'
     3627518      read(*,*,iostat=ierr) val3,val4
     3628         if(ierr.ne.0) goto 518
     3629         write(*,*) 'choice phi limite'
     3630519      read(*,*,iostat=ierr) val5
     3631         if(ierr.ne.0) goto 519
     3632         write(*,*) 'choice change alt (m)'
     3633520      read(*,*,iostat=ierr) val6
     3634         if(ierr.ne.0) goto 520
     3635
     3636         write(*,*) 're scale the pressure :'
     3637         r = 8.314511E+0 *1000.E+0/mugaz
     3638         temp=40.
     3639         write(*,*) 'r, sale height temperature =',r,temp
     3640                 
     3641         do j=1,jjp1
     3642           do i=1,iip1
     3643              phisprev= phis(i,j)
     3644              IF ( (  (rlatu(j)*180./pi.ge.val1)  .and.
     3645     &            (rlatu(j)*180./pi.le.val2) ) .and.
     3646     &          (  (rlonv(j)*180./pi.ge.val3)  .and.
     3647     &            (rlonv(j)*180./pi.le.val4) ) .and.
     3648     &            (phis(i,j).le.val5)) then
     3649
     3650                    phis(i,j)=phis(i,j)+val6*g
     3651                    ps(i,j) = ps(i,j) *
     3652     .              exp((phisprev-phis(i,j))/(temp * r))
     3653              ENDIF 
     3654           end do
     3655         end do
     3656c        periodicity of surface ps in longitude
     3657         do j=1,jjp1
     3658          ps(1,j) = ps(iip1,j)
     3659         end do
     3660         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     3661
     3662c       copy field 2D
     3663c       -----------------------------------------------------
     3664        else if (modif(1:len_trim(modif)) .eq. 'copytsoil') then
     3665           
     3666          ! Get field 2D
     3667          fichnom = 'startfi_input.nc'
     3668          ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
     3669          IF (ierr.NE.NF_NOERR) THEN
     3670            write(6,*)' Problem opening file:',fichnom
     3671            write(6,*)' ierr = ', ierr
     3672            CALL ABORT
     3673          ENDIF
     3674
     3675          ! Choice variable to be copied
     3676c          write(*,*) 'Choice variable to be copied ?'
     3677c515       read(*,fmt='(a20)',iostat=ierr) modif
     3678c          if(ierr.ne.0) goto 515
     3679c          write(*,*) 'variable asked :',trim(modif)
     3680
     3681          ierr = NF_INQ_VARID (nid_fi_input, trim("tsurf"),nvarid_input)
     3682          IF (ierr .NE. NF_NOERR) THEN
     3683           PRINT*, "Could not find asked variable in startfi_input.nc"
     3684           CALL abort
     3685          ENDIF
     3686          ierr = NF_INQ_VARID (nid_fi_input, trim("tsoil"),
     3687     &                                             nvarid_inputs)
     3688          IF (ierr .NE. NF_NOERR) THEN
     3689           PRINT*, "Could not find asked variable s in startfi_input.nc"
     3690           CALL abort
     3691          ENDIF
     3692
     3693#ifdef NC_DOUBLE
     3694          ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
     3695     &                                              field_input)
     3696#else
     3697          ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
     3698#endif
     3699          IF (ierr .NE. NF_NOERR) THEN
     3700           PRINT*, "Could not get asked variable in startfi_input.nc"
     3701           CALL abort
     3702          ELSE
     3703           PRINT*, "Got variable in startfi_input.nc"
     3704          ENDIF
     3705#ifdef NC_DOUBLE
     3706          ierr=NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_inputs,
     3707     &                        field_inputs)
     3708#else
     3709          ierr=NF_GET_VAR_REAL(nid_fi_input,nvarid_inputs,field_inputs)
     3710#endif
     3711          IF (ierr .NE. NF_NOERR) THEN
     3712           PRINT*, "Could not get asked variable s in startfi_input.nc"
     3713           CALL abort
     3714          ELSE
     3715           PRINT*, "Got variable s in startfi_input.nc"
     3716          ENDIF
     3717
     3718          write(*,*) 'Choice domain lon1 lon2 ?'
     3719525       read(*,*,iostat=ierr) val,val2
     3720          if(ierr.ne.0) goto 525
     3721          write(*,*) 'Choice domain lat1 lat2 ?'
     3722526       read(*,*,iostat=ierr) val3,val4
     3723          if(ierr.ne.0) goto 526
     3724          write(*,*) 'No change if N2 ice (0) / Change (1) ?'
     3725527       read(*,*,iostat=ierr) val5
     3726          if(ierr.ne.0) goto 527
     3727
     3728          ! Loop
     3729          DO ig=1,ngridmx
     3730             IF ( (lonfi(ig)*180./pi.ge.val) .and.
     3731     &            (lonfi(ig)*180./pi.le.val2) .and.
     3732     &            (latfi(ig)*180./pi.ge.val3)  .and.
     3733     &            (latfi(ig)*180./pi.le.val4) ) then
     3734                   if (qsurf(ig,igcm_n2).lt.0.001.or.val5.eq.1) then
     3735                      tsurf(ig)=field_input(ig)
     3736                      do l=1,nsoilmx
     3737                         tsoil(ig,l) = field_input(ig)
     3738                        !tsoil(ig,l) = field_inputs(ig,l)
     3739                      end do
     3740                   endif
     3741             ENDIF
     3742          ENDDO
     3743
     3744c       modify ice distribution depending on albedo
     3745c       -----------------------------------------------------
     3746        else if (modif(1:len_trim(modif)) .eq. 'mod_distrib_ice') then
     3747
     3748          write(*,*) 'Choice range albedo for CH4 ice'
     3749 220      read(*,*,iostat=ierr) val,val2
     3750          if(ierr.ne.0) goto 220
     3751          write(*,*) 'Choice range albedo for N2 ice'
     3752 221      read(*,*,iostat=ierr) val3,val4
     3753          if(ierr.ne.0) goto 221
     3754          write(*,*) 'Choice extra mass of CH4 ice'
     3755 222      read(*,*,iostat=ierr) val5
     3756          if(ierr.ne.0) goto 222
     3757          write(*,*) 'Choice extra mass of N2 ice'
     3758 223      read(*,*,iostat=ierr) val6
     3759          if(ierr.ne.0) goto 223
     3760
     3761          DO ig=1,ngridmx
     3762             IF (   (albfi(ig).ge.val)  .and.
     3763     &              (albfi(ig).le.val2) ) then
     3764                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
     3765             ENDIF
     3766             IF (   (albfi(ig).ge.val3)  .and.
     3767     &              (albfi(ig).le.val4) ) then
     3768                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
     3769             ENDIF
    20803770          ENDDO
    20813771
    … …  
    21213811               if (map_pluto_dat(i,j).eq.3) then
    21223812                  N2_pluto_dat(i,j)=100000.
    2123                elseif (map_pluto_dat(i,j).eq.4) then
     3813               else if (map_pluto_dat(i,j).eq.4) then
    21243814                  CH4_pluto_dat(i,j)=100000.
    2125                elseif (map_pluto_dat(i,j).eq.0) then
     3815               else if (map_pluto_dat(i,j).eq.0) then
    21263816                  alb_dat(i,j)=0.3
    2127                elseif (map_pluto_dat(i,j).eq.6) then
     3817               else if (map_pluto_dat(i,j).eq.6) then
    21283818                  alb_dat(i,j)=0.6
    2129                elseif (map_pluto_dat(i,j).eq.7) then
     3819               else if (map_pluto_dat(i,j).eq.7) then
    21303820                  alb_dat(i,j)=0.1
    21313821               endif
    … …  
    23524042      real dlon, dlat
    23534043      real a,c
    2354       real rad
    2355       rad=1190. ! km
     4044      real radp
     4045      radp=1190. ! km
    23564046
    23574047      dlon = lon2 - lon1
    … …  
    23594049      a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
    23604050      c = 2 * atan2( sqrt(a), sqrt(1-a) )
    2361       dist_pluto = rad * c
     4051      dist_pluto = radp * c
    23624052      return
    23634053      end
Note: See TracChangeset for help on using the changeset viewer.