Changeset 3631


Ignore:
Timestamp:
Feb 19, 2025, 2:26:27 PM (5 months ago)
Author:
afalco
Message:

Pluto: initialize variable to 0.
AF

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90

    r1422 r3631  
    3838  jm      = jjm
    3939  lllm    = llm
    40   imp1    = iim 
     40  imp1    = iim
    4141  jmp1    = jjm + 1
    4242  lllmm1  = llm - 1
    … …  
    4444
    4545  !-----------------------------------------------------------------------
    46 
    4746  dtphys  = iphysiq * dtvr
    4847  unsim   = 1./iim

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

    r3539 r3631  
    4343      use exner_hyb_m, only: exner_hyb
    4444      use geometry_mod, only: longitude,  ! longitudes (rad)
    45      &                         latitude,  ! latitudes (rad)                       
    46      &                         cell_area ! physics grid area (m2)                       
     45     &                         latitude,  ! latitudes (rad)
     46     &                         cell_area ! physics grid area (m2)
    4747      implicit none
    4848
    4949      include "dimensions.h"
    50       integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
     50      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
    5151      include "paramet.h"
    5252      include "comgeom2.h"
    … …  
    6262      CHARACTER        relief*3
    6363
    64 c Variables pour les lectures NetCDF des fichiers "start_archive" 
     64c Variables pour les lectures NetCDF des fichiers "start_archive"
    6565c--------------------------------------------------
    6666      INTEGER nid_dyn, nid_fi,nid,nvarid,nvarid_input,nvarid_inputs
    … …  
    7575      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
    7676
    77 c Variable histoire 
     77c Variable histoire
    7878c------------------
    7979      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
    … …  
    9191      PARAMETER (klatdat=180,klongdat=360)
    9292
    93 c Physique sur grille scalaire 
     93c Physique sur grille scalaire
    9494c----------------------------
    9595      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
    … …  
    116116      real mugaz ! molar mass of the atmosphere
    117117
    118       integer ierr 
    119 
    120 c Variables on the new grid along scalar points 
     118      integer ierr
     119
     120c Variables on the new grid along scalar points
    121121c------------------------------------------------------
    122122      REAL t(iip1,jjp1,llm)
    … …  
    132132      REAL :: p3d(iip1, jjp1, llm+1)
    133133
    134 c Variable de l'ancienne grille 
     134c Variable de l'ancienne grille
    135135c------------------------------
    136136      real time
    … …  
    160160
    161161      INTEGER :: itau
    162      
     162
    163163      character(len=20) :: txt ! to store some text
    164164      character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
    … …  
    174174      real fact2
    175175
    176 c Special Pluto Map from Lellouch & Grundy 
     176c Special Pluto Map from Lellouch & Grundy
    177177c------------------------------------------
    178178       integer,parameter :: im_plu=360 ! from topo
    179        integer,parameter :: jm_plu=180   
     179       integer,parameter :: jm_plu=180
    180180       real latv_plu(jm_plu+1),lonu_plu(im_plu+1)
    181181       real map_pluto_dat(im_plu,jm_plu+1)
    … …  
    224224      is_omp_root=.true.
    225225      is_master=.true.
    226      
     226
    227227! Load tracer number and names:
    228228      call infotrac_init
    … …  
    231231      allocate(qsurf(ngridmx,nqtot))
    232232      allocate(qsurf_tmp(ngridmx,nqtot))
    233      
     233
    234234! get value of nsoilmx and allocate corresponding arrays
    235235      ! a default value of nsoilmx is set in comsoil_h
    236236      call getin_p("nsoilmx",nsoilmx)
    237      
     237
    238238      allocate(inertiedat_simple(ngridmx,nsoilmx))
    239239      allocate(tsoil(ngridmx,nsoilmx))
    240240      allocate(field_inputs(ngridmx,nsoilmx))
    241241      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
    242      
     242
    243243c=======================================================================
    244244c   Choice of the start file(s) to use
    … …  
    248248      write(*,*) '    0 - from a file start_archive'
    249249      write(*,*) '    1 - from files start and startfi'
    250  
     250
    251251c-----------------------------------------------------------------------
    252252c   Open file(s) to modify (start or start_archive)
    … …  
    272272          CALL ABORT
    273273        ENDIF
    274         tab0 = 50 
     274        tab0 = 50
    275275        Lmodif = 1
    276276
    … …  
    288288          CALL ABORT
    289289        ENDIF
    290  
     290
    291291        fichnom = 'startfi.nc'
    292292        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
    … …  
    297297        ENDIF
    298298
    299         tab0 = 0 
     299        tab0 = 0
    300300        Lmodif = 0
    301301
    … …  
    308308! Initialize global tracer indexes (stored in tracer.h)
    309309! ... this has to be done before phyetat0
    310 ! and requires that "datadir" be correctly initialized 
     310! and requires that "datadir" be correctly initialized
    311311      call getin_p("datadir",datadir)
    312312      call initracer(ngridmx,nqtot)
    … …  
    358358          write(*,*) i,tab_cntrl(i)
    359359        enddo
    360        
     360
    361361        write(*,*) 'Reading file START'
    362362        fichnom = 'start.nc'
    … …  
    375375     &                  1)
    376376
    377         ! Lmodif set to 0 to disable modifications possibility in phyeta0                           
     377        ! Lmodif set to 0 to disable modifications possibility in phyeta0
    378378        Lmodif=0
    379379        write(*,*) 'Reading file STARTFI'
    … …  
    399399        call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    400400        call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
    401      
     401
    402402      endif
    403403c-----------------------------------------------------------------------
    … …  
    405405c-----------------------------------------------------------------------
    406406
    407       kappa = tab_cntrl(9) 
     407      kappa = tab_cntrl(9)
    408408      etot0 = tab_cntrl(12)
    409409      ptot0 = tab_cntrl(13)
    … …  
    448448         fichnom = 'startfi.nc'
    449449         call open_startphy(fichnom)
    450          Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
     450         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
    451451         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
    452452     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    … …  
    482482      daysec = p_daysec
    483483      kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW
     484      dtvr = 0 ! will be changed later but we need it initalized in iniconst
    484485
    485486c=======================================================================
    … …  
    488489! Load parameters from run.def file
    489490      CALL defrun_new( 99, .TRUE. )
    490 ! Initialize dynamics geometry and co. (which, when using start.nc may 
     491! Initialize dynamics geometry and co. (which, when using start.nc may
    491492! have changed because of modifications (tabi, preff, pa) above)
    492       CALL iniconst 
     493      CALL iniconst
    493494      CALL inigeom
    494495      idum=-1
    … …  
    509510c=======================================================================
    510511
    511       if (choix_1.eq.0) then  ! for start_archive files, 
     512      if (choix_1.eq.0) then  ! for start_archive files,
    512513                              ! where an external "surface.nc" file is needed
    513514
    … …  
    551552          zgamS(:,:)=0
    552553          ztheS(:,:)=0
    553          
     554
    554555          write(*,*) "Enter value of albedo of the bare ground:"
    555556          read(*,*) alb(1,1)
    556557          alb(:,:)=alb(1,1)
    557          
     558
    558559          write(*,*) "Enter value of thermal inertia of soil:"
    559560          read(*,*) surfith(1,1)
    560561          surfith(:,:)=surfith(1,1)
    561        
     562
    562563        endif ! of if (surfacefile.ne."none")
    563564
    … …  
    588589        ! copy soil thermal inertia
    589590        ithfi(:,:)=inertiedat(:,:)
    590        
     591
    591592        ierr= NF_CLOSE(nid)
    592593
    593       else if (choix_1.eq.1) then 
     594      else if (choix_1.eq.1) then
    594595         !do nothing, start and startfi have already been read
    595       else 
     596      else
    596597        CALL exit(1)
    597598      endif
    … …  
    601602
    602603c=======================================================================
    603 c 
     604c
    604605c=======================================================================
    605606
    … …  
    643644        write(*,*) 'tsurfice: temperature depending on N2 ice'
    644645        write(*,*) 'icarus: option to change surface/soil temperature'
    645         write(*,*) 'icarus_ch4: special option to add ch4' 
    646         write(*,*) 'delfrostch4: delete ch4 frost' 
    647         write(*,*) 'modch4: remove/modify amount ch4 frost' 
    648         write(*,*) 'modch4n2: modify amount ch4 frost according N2' 
    649         write(*,*) 'modco: remove/modify amount co frost' 
     646        write(*,*) 'icarus_ch4: special option to add ch4'
     647        write(*,*) 'delfrostch4: delete ch4 frost'
     648        write(*,*) 'modch4: remove/modify amount ch4 frost'
     649        write(*,*) 'modch4n2: modify amount ch4 frost according N2'
     650        write(*,*) 'modco: remove/modify amount co frost'
    650651        write(*,*) 'modn2: remove/modify amount n2 ice'
    651652        write(*,*) 'modcoch4: remove/modify co ch4 where no n2 '
    652653        write(*,*) 'checktsoil: change tsoil where no n2 '
    653654        write(*,*) 'non2noco: if no n2, no co '
    654         write(*,*) 'modn2_topo: modify n2 ice topo according to topo' 
    655         write(*,*) 'modwheren2: modify n2 where already n2' 
    656         write(*,*) 'modn2HD: modify n2 for HD runs' 
    657         write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP' 
    658         write(*,*) 'globch4co: add/remove global amount ch4co frost' 
    659         write(*,*) 'source_ch4: add source ch4' 
    660         write(*,*) 'nomountain: remove pluto mountains (!)' 
    661         write(*,*) 'relief: add pluto crater or mountain' 
     655        write(*,*) 'modn2_topo: modify n2 ice topo according to topo'
     656        write(*,*) 'modwheren2: modify n2 where already n2'
     657        write(*,*) 'modn2HD: modify n2 for HD runs'
     658        write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP'
     659        write(*,*) 'globch4co: add/remove global amount ch4co frost'
     660        write(*,*) 'source_ch4: add source ch4'
     661        write(*,*) 'nomountain: remove pluto mountains (!)'
     662        write(*,*) 'relief: add pluto crater or mountain'
    662663        write(*,*) 'seticeNH: set ice for initial start with NH topo'
    663664        write(*,*) 'topo_sp: change topography of SP'
    … …  
    684685        write(*,*) 'albedomap: read in an albedomap albedo.nc'
    685686        write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo'
    686      
     687
    687688        write(*,*)
    688689        write(*,*) 'Please note that emis and albedo are set in physiq'
    … …  
    705706c       -------------------------------------
    706707        if (trim(modif) .eq. 'flat') then
    707 c         set topo to zero 
     708c         set topo to zero
    708709          z_reel(1:iip1,1:jjp1)=0
    709710          phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
    … …  
    711712          write(*,*) 'topography set to zero.'
    712713          write(*,*) 'WARNING : the subgrid topography parameters',
    713      &    ' were not set to zero ! => set calllott to F'                   
     714     &    ' were not set to zero ! => set calllott to F'
    714715
    715716c        Choice of surface pressure
    … …  
    741742         end if
    742743
    743 c       'set_ps_to_preff' : used if changing preff with topo 
     744c       'set_ps_to_preff' : used if changing preff with topo
    744745c       ----------------------------------------------------
    745746        else if (trim(modif) .eq. 'set_ps_to_preff') then
    … …  
    750751          enddo
    751752
    752 c       ptot : Modification of the total pressure: ice + current atmosphere 
     753c       ptot : Modification of the total pressure: ice + current atmosphere
    753754c       -------------------------------------------------------------------
    754755        else if (modif(1:len_trim(modif)).eq.'ptot') then
    … …  
    826827            tname(iq)=txt
    827828            write(*,*)'Do you want to change another tracer name (y/n)?'
    828             read(*,'(a)') yes 
     829            read(*,'(a)') yes
    829830          else
    830831! inapropiate value of iq; quit this option
    … …  
    855856           ENDDO
    856857
    857 c       q=x : initialise tracer manually 
     858c       q=x : initialise tracer manually
    858859c       --------------------------------
    859860        else if (trim(modif).eq.'q=x') then
    … …  
    863864             enddo
    864865             write(*,*) '(choose between 1 and ',nqtot,')'
    865              read(*,*) iq 
     866             read(*,*) iq
    866867             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
    867868     &                 ' ? (kg/kg)'
    … …  
    880881c                 qsurf(ig,iq)=val
    881882c             ENDDO
    882              
    883 c       qs=x : initialise surface tracer manually 
     883
     884c       qs=x : initialise surface tracer manually
    884885c       --------------------------------
    885886        else if (trim(modif).eq.'qs=x') then
    … …  
    889890             enddo
    890891             write(*,*) '(choose between 1 and ',nqtot,')'
    891              read(*,*) iq 
     892             read(*,*) iq
    892893             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
    893894     &                   ' ? (kg/m2)'
    … …  
    949950             enddo
    950951             write(*,*) '(choose between 1 and ',nqtot,')'
    951              read(*,*) iq 
     952             read(*,*) iq
    952953             if ((iq.lt.1).or.(iq.gt.nqtot)) then
    953954               ! wrong value for iq, go back to menu
    … …  
    984985             endif
    985986
    986 c       qnogcm : initialise tracer from nogcm (CH4, CO) 
     987c       qnogcm : initialise tracer from nogcm (CH4, CO)
    987988c       --------------------------------
    988989        else if (modif(1:len_trim(modif)).eq.'qnogcm') then
    … …  
    993994             enddo
    994995             write(*,*) '(choose between 1 and ',nqtot,')'
    995              read(*,*) iq 
     996             read(*,*) iq
    996997             DO l=1,llm
    997998               DO j=1,jjp1
    … …  
    10021003             ENDDO
    10031004
    1004 c       isothermal temperatures 
     1005c       isothermal temperatures
    10051006c       ------------------------------------------------
    10061007        else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
    … …  
    10201021          write(*,*) 'atmospheric temp =', Tset(2,2,2)
    10211022
    1022 c       specific temperature profile 
     1023c       specific temperature profile
    10231024c       ------------------------------------------------
    10241025        else if (modif(1:len_trim(modif)) .eq. 'tempprof') then
    … …  
    12021203             ENDDO
    12031204          ENDIF
    1204  
     1205
    12051206c       bladed : adding/remove ch4 on bladed terrains
    12061207c       ----------------------------------------------------------
    … …  
    12271228          if(ierr.ne.0) goto 449
    12281229
    1229           ! shift 
     1230          ! shift
    12301231          DO ig=1,ngridmx
    12311232             IF (   (latfi(ig)*180./pi.ge.val)  .and.
    … …  
    12331234     &              (lonfi(ig)*180./pi.ge.val3)  .and.
    12341235     &              (lonfi(ig)*180./pi.le.val4)  .and.
    1235      &              (phisfi(ig).gt.choice) )  then     
     1236     &              (phisfi(ig).gt.choice) )  then
    12361237               qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5
    12371238               qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
    … …  
    12571258          if(ierr.ne.0) goto 562
    12581259
    1259           ! Get index correponding to central points 
     1260          ! Get index correponding to central points
    12601261          ! 1/ Reference point
    12611262          igref=1
    … …  
    12671268                   actualmin=val6
    12681269                   igref=ig
    1269                endif 
     1270               endif
    12701271          enddo
    12711272
    1272           do k=1,l 
     1273          do k=1,l
    12731274
    12741275            write(*,*) 'Choice lat,lon where terrain is copied'
    12751276 563        read(*,*,iostat=ierr) val4,val5
    12761277            if(ierr.ne.0) goto 563
    1277          
     1278
    12781279            if (val5.gt.180.) then
    12791280              val5=val5-360.
    … …  
    12891290                   actualmin=val6
    12901291                   igref2=ig
    1291                endif 
     1292               endif
    12921293            enddo
    12931294
    … …  
    13101311            enddo
    13111312
    1312           ! for each point within A2, get distance and angle 
    1313           ! and look where fit with previous table, and set value 
     1313          ! for each point within A2, get distance and angle
     1314          ! and look where fit with previous table, and set value
    13141315
    13151316            do ig=1,ngridmx
    … …  
    13861387     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    13871388     &              (lonfi(ig)*180./pi.lt.val5)) then
    1388 !     &              (qsurf(ig,igcm_n2).gt.50))  then     
     1389!     &              (qsurf(ig,igcm_n2).gt.50))  then
    13891390!     &              (qsurf(ig,igcm_ch4_ice).lt.10) )      then
    13901391                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
    … …  
    14141415     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    14151416     &              (lonfi(ig)*180./pi.lt.val5) .and.
    1416      &              (qsurf(ig,igcm_n2).gt.val7)  .and.     
     1417     &              (qsurf(ig,igcm_n2).gt.val7)  .and.
    14171418     &              (qsurf(ig,igcm_n2).lt.val8) )      then
    14181419                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
    … …  
    14251426        else if (modif(1:len_trim(modif)) .eq. 'non2noco') then
    14261427          DO ig=1,ngridmx
    1427              IF (   (qsurf(ig,igcm_n2).lt.0.001))  then     
     1428             IF (   (qsurf(ig,igcm_n2).lt.0.001))  then
    14281429                    qsurf(ig,igcm_co_ice)=0.
    14291430             ENDIF
    … …  
    14551456     &              (latfi(ig)*180./pi.le.val2)  .and.
    14561457     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    1457      &              (lonfi(ig)*180./pi.lt.val5)  .and. 
    1458      &              (qsurf(ig,igcm_n2).gt.val7))  then     
     1458     &              (lonfi(ig)*180./pi.lt.val5)  .and.
     1459     &              (qsurf(ig,igcm_n2).gt.val7))  then
    14591460                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
    14601461                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
    … …  
    14621463          ENDDO
    14631464
    1464 c       modn2 : adding/remove n2 ice 
     1465c       modn2 : adding/remove n2 ice
    14651466c       ----------------------------------------------------------
    14661467        else if (modif(1:len_trim(modif)) .eq. 'modn2') then
    … …  
    14841485     &              (latfi(ig)*180./pi.le.val2)  .and.
    14851486     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    1486      &              (lonfi(ig)*180./pi.lt.val5)  ) then 
    1487 c    &              (qsurf(ig,igcm_n2).lt.50))  then     
     1487     &              (lonfi(ig)*180./pi.lt.val5)  ) then
     1488c    &              (qsurf(ig,igcm_n2).lt.50))  then
    14881489                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
    14891490                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
    14901491             ENDIF
    1491 !             IF ((phisfi(ig).gt.-1000.)) then 
     1492!             IF ((phisfi(ig).gt.-1000.)) then
    14921493!                qsurf(ig,igcm_n2)=0.
    14931494!             ENDIF
    … …  
    15161517     &              (latfi(ig)*180./pi.le.val2)  .and.
    15171518     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    1518      &              (lonfi(ig)*180./pi.le.val5)  .and. 
    1519      &              (qsurf(ig,igcm_n2).lt.10.))  then     
     1519     &              (lonfi(ig)*180./pi.le.val5)  .and.
     1520     &              (qsurf(ig,igcm_n2).lt.10.))  then
    15201521                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
    15211522                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
    … …  
    15601561     &              (qsurf(ig,igcm_n2).lt.val5)  .and.
    15611562     &              (phisfi(ig).ge.val6) .and.
    1562      &              (phisfi(ig).le.val11)  ) then 
     1563     &              (phisfi(ig).le.val11)  ) then
    15631564
    15641565                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
    … …  
    15771578                    !DO l=1,nsoilmx
    15781579                    !   tsoil(ig,l)=tsoil(iref,l)
    1579                     !ENDDO 
     1580                    !ENDDO
    15801581                    tsurf(ig)=tsurf(ig)+4
    15811582                    DO l=1,nsoilmx
    15821583                       tsoil(ig,l)=tsoil(ig,l)+4
    1583                     ENDDO 
     1584                    ENDDO
    15841585             ENDIF
    1585 !             IF ((phisfi(ig).gt.-1000.)) then 
     1586!             IF ((phisfi(ig).gt.-1000.)) then
    15861587!                qsurf(ig,igcm_n2)=0.
    15871588!             ENDIF
    … …  
    16231624     &              (qsurf(ig,igcm_n2).lt.val5)  .and.
    16241625     &              (phisfi(ig).ge.val6) .and.
    1625      &              (phisfi(ig).le.val11)  ) then 
     1626     &              (phisfi(ig).le.val11)  ) then
    16261627
    16271628                    qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
    … …  
    16401641                     DO l=1,nsoilmx
    16411642                        tsoil(ig,l)=tsoil(iref,l)
    1642                      ENDDO 
     1643                     ENDDO
    16431644                    else if (iref.eq.0) then
    16441645                     choice=int(ig/96.)*96+84
    … …  
    16471648                     DO l=1,nsoilmx
    16481649                        tsoil(ig,l)=tsoil(int(choice),l)
    1649                      ENDDO 
     1650                     ENDDO
    16501651                    endif
    16511652             ENDIF
    16521653          ENDDO
    16531654
    1654 c       modn2_topo : adding/remove n2 ice 
     1655c       modn2_topo : adding/remove n2 ice
    16551656c       ----------------------------------------------------------
    16561657        else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then
    … …  
    17321733     &              (phisfi(ig).ge.val6) .and.
    17331734     &              (phisfi(ig).le.val11) .and.
    1734      &              (qsurf(ig,igcm_n2).gt.0.001) ) then 
     1735     &              (qsurf(ig,igcm_n2).gt.0.001) ) then
    17351736
    17361737!                    DO i=1,ngridmx
    … …  
    17381739!     &              (lonfi(i).eq.0.) ) then
    17391740!
    1740                          tsurf(ig)=34.7 
     1741                         tsurf(ig)=34.7
    17411742                         !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice)
    17421743!
    17431744                         DO l=1,nsoilmx
    17441745                           tsoil(ig,l)=34.7 !tsoil(i,l)
    1745                          ENDDO 
     1746                         ENDDO
    17461747                       !ENDIF
    17471748                    !ENDDO
    … …  
    17721773     &             (latfi(ig)*180./pi.le.val1+val2) .and.
    17731774     &             (qsurf(ig,igcm_n2).gt.0.) ) then
    1774                ! Look for same longitude point where no ice           
     1775               ! Look for same longitude point where no ice
    17751776               val5=1.
    17761777               jref=ig
    … …  
    17841785                  val5=qsurf(jref,igcm_n2)
    17851786!                 write(*,*) jref,
    1786 !     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 
    1787                enddo 
    1788                if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE' 
     1787!     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5
     1788               enddo
     1789               if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE'
    17891790
    17901791               ! Copy point in the selected area
    1791                tsurf(ig)=tsurf(jref) 
     1792               tsurf(ig)=tsurf(jref)
    17921793               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
    17931794               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
    … …  
    17951796               DO l=1,nsoilmx
    17961797                  tsoil(ig,l)=tsoil(jref,l)
    1797                ENDDO 
     1798               ENDDO
    17981799             ENDIF
    17991800             ! If below line and no ice: add ice
    … …  
    18011802     &             (latfi(ig)*180./pi.le.val1+val2) .and.
    18021803     &             (qsurf(ig,igcm_n2).eq.0.) ) then
    1803                ! Look for same longitude point where ice           
     1804               ! Look for same longitude point where ice
    18041805               val5=1.
    18051806               jref=ig
    … …  
    18101811                  val5=qsurf(jref,igcm_n2)
    18111812                  write(*,*) jref,
    1812      &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 
    1813                enddo 
    1814                if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE' 
     1813     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5
     1814               enddo
     1815               if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE'
    18151816
    18161817               ! Copy point in the selected area
    1817                tsurf(ig)=tsurf(jref) 
     1818               tsurf(ig)=tsurf(jref)
    18181819               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
    18191820               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
    … …  
    18211822               DO l=1,nsoilmx
    18221823                  tsoil(ig,l)=tsoil(jref,l)
    1823                ENDDO 
     1824               ENDDO
    18241825             ENDIF
    18251826
    … …  
    18301831        else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then
    18311832
    1832           write(*,*) 'Adding ch4 ice latitudinally ' 
     1833          write(*,*) 'Adding ch4 ice latitudinally '
    18331834          write(*,*) 'Choice : lat1 and lat2 ?'
    18341835 433      read(*,*,iostat=ierr) val,val2
    … …  
    18451846     &              (latfi(ig)*180./pi.le.val2)  .and.
    18461847     &              (lonfi(ig)*180./pi.ge.val3)  .and.
    1847      &              (lonfi(ig)*180./pi.lt.val4)  ) then 
     1848     &              (lonfi(ig)*180./pi.lt.val4)  ) then
    18481849                    qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
    18491850             ENDIF
    … …  
    18541855        else if (modif(1:len_trim(modif)) .eq. 'source_co') then
    18551856
    1856           write(*,*) 'Adding co ice latitudinally ' 
     1857          write(*,*) 'Adding co ice latitudinally '
    18571858          write(*,*) 'Choice : lat1 and lat2 ?'
    18581859 436      read(*,*,iostat=ierr) val,val2
    … …  
    18691870     &              (latfi(ig)*180./pi.le.val2)  .and.
    18701871     &              (lonfi(ig)*180./pi.ge.val3)  .and.
    1871      &              (lonfi(ig)*180./pi.lt.val4)  ) then 
     1872     &              (lonfi(ig)*180./pi.lt.val4)  ) then
    18721873                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5
    18731874             ENDIF
    … …  
    18881889         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    18891890
    1890 !       inert3d: set soil thermal inertia to specific uniform value 
     1891!       inert3d: set soil thermal inertia to specific uniform value
    18911892!       ----------------------------------------------------------------------
    18921893        else if (modif(1:len_trim(modif)).eq.'inert3d') then
    … …  
    19321933                    jref=j+1
    19331934                  endif
    1934                  
     1935
    19351936                  write(*,*)'Will use nearest grid latitude which is:',
    19361937     &                     rlatu(jref)*180./pi
    … …  
    19661967          endif
    19671968         enddo ! of do while
    1968  
     1969
    19691970         ! find the reference index iref the depth corresponds to
    19701971!        if (val2.lt.layer(1)) then
    … …  
    19891990         ! recast ithfi() onto ith()
    19901991         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    1991      
     1992
    19921993         do j=1,jref
    19931994!           write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
    … …  
    20012002              ith(i,j,isoil)=iceith ! ice
    20022003             enddo
    2003             enddo ! of do i=1,iip1 
     2004            enddo ! of do i=1,iip1
    20042005         enddo ! of do j=1,jjp1
    2005  
     2006
    20062007         CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    20072008
    … …  
    20332034                   jref=j+1
    20342035                  endif
    2035  
     2036
    20362037                 write(*,*)'Will use nearest grid latitude which is:',
    20372038     &                     rlatu(jref)*180./pi
    … …  
    20672068          endif
    20682069         enddo ! of do while
    2069  
     2070
    20702071         ! find the reference index iref the depth corresponds to
    20712072          do isoil=1,nsoilmx-1
    … …  
    20762077           endif
    20772078          enddo
    2078  
     2079
    20792080         ! thermal inertia of the ice:
    20802081         ierr=1
    … …  
    20842085          read(*,*,iostat=ierr)iceith
    20852086         enddo ! of do while
    2086  
     2087
    20872088         ! recast ithfi() onto ith()
    20882089         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    2089  
     2090
    20902091         do j=jref,jjp1
    20912092!           write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
    … …  
    20992100              ith(i,j,isoil)=iceith ! ice
    21002101             enddo
    2101             enddo ! of do i=1,iip1 
     2102            enddo ! of do i=1,iip1
    21022103         enddo ! of do j=jref,jjp1
    21032104
    … …  
    21302131            IF (   (latfi(ig)*180./pi.ge.val3)  .and.
    21312132     &            (latfi(ig)*180./pi.le.val4)  .and.
    2132      &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
     2133     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
    21332134     &                       (lonfi(ig)*180./pi.le.val6))  .or.
    2134      &         ((lonfi(ig)*180./pi.ge.val5)  .and. 
     2135     &         ((lonfi(ig)*180./pi.ge.val5)  .and.
    21352136     &                       (lonfi(ig)*180./pi.le.180.))) ) then
    21362137
    … …  
    21452146          ENDDO
    21462147
    2147 c       subsoil_SP : choice of thermal inertia values for SP 
     2148c       subsoil_SP : choice of thermal inertia values for SP
    21482149c       ----------------------------------------------------------------
    21492150        else if (modif(1:len_trim(modif)) .eq. 'subsoil_SP') then
    … …  
    21712172               write(*,*)'Thermal inertia set for all subsurface layers'
    21722173               ierr=0
    2173               else 
     2174              else
    21742175               write(*,*)'Depth should be more than ',layer(1)
    21752176               ierr=1
    … …  
    21822183         if(val2.eq.0) then
    21832184           iref=1
    2184            write(*,*)'Level selected is first level: ',layer(iref),' m' 
     2185           write(*,*)'Level selected is first level: ',layer(iref),' m'
    21852186         else
    21862187           do isoil=1,nsoilmx-1
    … …  
    21922193            endif
    21932194           enddo
    2194          endif 
     2195         endif
    21952196
    21962197         ! Definition SP:
    … …  
    22092210           IF (   (latfi(ig)*180./pi.ge.val3)  .and.
    22102211     &            (latfi(ig)*180./pi.le.val4)  .and.
    2211      &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
     2212     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
    22122213     &                       (lonfi(ig)*180./pi.le.val6))  .or.
    2213      &         ((lonfi(ig)*180./pi.ge.val5)  .and. 
     2214     &         ((lonfi(ig)*180./pi.ge.val5)  .and.
    22142215     &                       (lonfi(ig)*180./pi.le.180.))) ) then
    22152216
    … …  
    22192220                       ithfi(ig,j)=ith_bb
    22202221                    ENDDO
    2221                 ENDIF 
     2222                ENDIF
    22222223           ENDIF
    22232224         ENDDO
    … …  
    22482249         temp=40.
    22492250         write(*,*) 'r, sale height temperature =',r,temp
    2250                  
     2251
    22512252         do j=1,jjp1
    22522253           do i=1,iip1
    … …  
    22542255              IF (   (rlatu(j)*180./pi.ge.val2)  .and.
    22552256     &            (rlatu(j)*180./pi.le.val3)  .and.
    2256      &      (  ((rlonv(i)*180./pi.ge.-180.)  .and. 
     2257     &      (  ((rlonv(i)*180./pi.ge.-180.)  .and.
    22572258     &                       (rlonv(i)*180./pi.le.val5))  .or.
    2258      &         ((rlonv(i)*180./pi.ge.val4)  .and. 
     2259     &         ((rlonv(i)*180./pi.ge.val4)  .and.
    22592260     &                       (rlonv(i)*180./pi.le.180.))) ) then
    22602261
    … …  
    22632264                    ps(i,j) = ps(i,j) *
    22642265     .              exp((phisprev-phis(i,j))/(temp * r))
    2265                 ENDIF 
    2266               ENDIF 
     2266                ENDIF
     2267              ENDIF
    22672268           end do
    22682269         end do
    … …  
    22932294         temp=40.
    22942295         write(*,*) 'r, sale height temperature =',r,temp,g
    2295          qsurf=0.       
     2296         qsurf=0.
    22962297
    22972298         write(*,*) 'latfi=',latfi
    … …  
    23032304         write(*,*) 'phisold=',phisold
    23042305         do ig=1,ngridmx
    2305                
     2306
    23062307              write(*,*) 'lat=',latfi(ig)*180./pi
    23072308              write(*,*) 'lon=',lonfi(ig)*180./pi
    … …  
    23092310              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
    23102311     &            (latfi(ig)*180./pi.le.val3)  .and.
    2311      &           (lonfi(ig)*180./pi.ge.val4)  .and. 
     2312     &           (lonfi(ig)*180./pi.ge.val4)  .and.
    23122313     &           (lonfi(ig)*180./pi.le.val5) ) then
    23132314                write(*,*) 'hello SP',phisfi(ig),ig
    … …  
    23172318                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
    23182319                ENDIF
    2319               ENDIF 
     2320              ENDIF
    23202321         end do
    23212322c        update new phis and ps
    … …  
    23552356         end do
    23562357
    2357 c       correct phisfi 
     2358c       correct phisfi
    23582359c       -----------------------------------------------------
    23592360        else if (modif(1:len_trim(modif)) .eq. 'correctphi') then
    … …  
    23782379              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
    23792380     &            (latfi(ig)*180./pi.le.val2)  .and.
    2380      &           (lonfi(ig)*180./pi.ge.val3)  .and. 
     2381     &           (lonfi(ig)*180./pi.ge.val3)  .and.
    23812382     &           (lonfi(ig)*180./pi.le.val4) ) then
    23822383
    … …  
    24142415         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    24152416
    2416 c       correct phisfi 
     2417c       correct phisfi
    24172418c       -----------------------------------------------------
    24182419        else if (modif(1:len_trim(modif)) .eq. 'correctps') then
    … …  
    24422443         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    24432444
    2444 c       No flat topo 
     2445c       No flat topo
    24452446c       -----------------------------------------------------
    24462447        else if (modif(1:len_trim(modif)) .eq. 'toponoise') then
    … …  
    24652466              IF (   (latfi(ig)*180./pi.ge.val1)  .and.
    24662467     &            (latfi(ig)*180./pi.le.val2)  .and.
    2467      &           (lonfi(ig)*180./pi.ge.val3)  .and. 
     2468     &           (lonfi(ig)*180./pi.ge.val3)  .and.
    24682469     &           (lonfi(ig)*180./pi.le.val4) ) then
    24692470
    … …  
    25202521         temp=40.
    25212522         write(*,*) 'r, sale height temperature =',r,temp,g
    2522          qsurf=0.       
     2523         qsurf=0.
    25232524
    25242525         write(*,*) 'latfi=',latfi
    … …  
    25302531         write(*,*) 'phisold=',phisold
    25312532         do ig=1,ngridmx
    2532                
     2533
    25332534              write(*,*) 'lat=',latfi(ig)*180./pi
    25342535              write(*,*) 'lon=',lonfi(ig)*180./pi
    … …  
    25362537              IF (   (latfi(ig)*180./pi.ge.val2)  .and.
    25372538     &            (latfi(ig)*180./pi.le.val3)  .and.
    2538      &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
     2539     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
    25392540     &                       (lonfi(ig)*180./pi.le.val5))  .or.
    2540      &         ((lonfi(ig)*180./pi.ge.val4)  .and. 
     2541     &         ((lonfi(ig)*180./pi.ge.val4)  .and.
    25412542     &                    (lonfi(ig)*180./pi.le.180.))) ) then
    25422543                write(*,*) 'hello SP',phisfi(ig),ig
    … …  
    25462547                    write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
    25472548                ENDIF
    2548               ENDIF 
     2549              ENDIF
    25492550         end do
    25502551c        update new phis and ps
    … …  
    25612562         end do
    25622563
    2563 c       bugch4 correct bug warm ch4 due to change of resolution 
     2564c       bugch4 correct bug warm ch4 due to change of resolution
    25642565c       -----------------------------------------------------
    25652566        else if (modif(1:len_trim(modif)) .eq. 'bugch4') then
    … …  
    25922593                    DO l=1,nsoilmx
    25932594                       tsoil(ig,l)=tsoil(int(ig-jref),l)
    2594                     ENDDO 
     2595                    ENDDO
    25952596                    goto 548
    25962597                 ELSEIF ((ig+jref.le.ngridmx).and.
    … …  
    26022603                    DO l=1,nsoilmx
    26032604                       tsoil(ig,l)=tsoil(int(ig+jref),l)
    2604                     ENDDO 
     2605                    ENDDO
    26052606                    goto 548
    26062607                 ENDIF
    26072608               enddo
    26082609548            write(*,*) 'found point with ch4'
    2609              ENDIF 
    2610            ENDIF 
     2610             ENDIF
     2611           ENDIF
    26112612         END DO
    26122613
    2613 c       bugn2 correct bug warm n2 due to change of resolution 
     2614c       bugn2 correct bug warm n2 due to change of resolution
    26142615c       -----------------------------------------------------
    26152616        else if (modif(1:len_trim(modif)) .eq. 'bugn2') then
    … …  
    26422643                    DO l=1,nsoilmx
    26432644                       tsoil(ig,l)=tsoil(int(ig-jref),l)
    2644                     ENDDO 
     2645                    ENDDO
    26452646                    goto 541
    26462647                 ELSEIF ((ig+jref.le.ngridmx).and.
    … …  
    26512652                    DO l=1,nsoilmx
    26522653                       tsoil(ig,l)=tsoil(int(ig+jref),l)
    2653                     ENDDO 
     2654                    ENDDO
    26542655                    goto 541
    26552656                 ENDIF
    26562657               enddo
    26572658541            write(*,*) 'found point with n2'
    2658              ENDIF 
    2659            ENDIF 
     2659             ENDIF
     2660           ENDIF
    26602661         END DO
    26612662
    … …  
    26722673!                    DO l=1,nsoilmx
    26732674!                      tsoil(ig,l)=tsoil(ig-1,l)
    2674 !                   ENDDO 
     2675!                   ENDDO
    26752676!                 ELSEIF (qsurf(ig+1,igcm_n2).gt.0.001) then
    26762677!                   tsurf(ig)=tsurf(ig+1)
    26772678!                   DO l=1,nsoilmx
    26782679!                      tsoil(ig,l)=tsoil(ig+1,l)
    2679 !                   ENDDO 
     2680!                   ENDDO
    26802681!                 ENDIF
    2681 !            ENDIF 
    2682 !          ENDIF 
     2682!            ENDIF
     2683!          ENDIF
    26832684
    26842685!        END DO
    26852686
    2686 c       flatnosp : flat topo outside a specific terrain (SP) 
     2687c       flatnosp : flat topo outside a specific terrain (SP)
    26872688c       -----------------------------------------------------
    26882689        else if (modif(1:len_trim(modif)) .eq. 'flatnosp') then
    … …  
    26912692551       read(*,*,iostat=ierr) val
    26922693          if(ierr.ne.0) goto 551
    2693           write(*,*) 'gravity g is : ',g         
     2694          write(*,*) 'gravity g is : ',g
    26942695          geop=g*val
    26952696          write(*,*) 'Choice of minimum amount of N2 ice we keep ?'
    26962697552       read(*,*,iostat=ierr) val2
    26972698          if(ierr.ne.0) goto 552
    2698        
     2699
    26992700          write(*,*) 'phis=',phis
    27002701          write(*,*) 'phisfi=',phisfi
    27012702          do ig=1,ngridmx
    2702                
     2703
    27032704              IF (   (qsurf(ig,1).le.val2)  .and.
    27042705     &               (phisfi(ig).gt.geop)  ) then
    … …  
    27262727          end do
    27272728
    2728 c       flatregion : flat topo specific to region 
     2729c       flatregion : flat topo specific to region
    27292730c       -----------------------------------------------------
    27302731        else if (modif(1:len_trim(modif)) .eq. 'flatregion') then
    … …  
    27422743 556      read(*,*,iostat=ierr) val7
    27432744          if(ierr.ne.0) goto 556
    2744        
     2745
    27452746          write(*,*) 'phis=',phis
    27462747          write(*,*) 'phisfi=',phisfi
    … …  
    27492750     &              (latfi(ig)*180./pi.le.val2)  .and.
    27502751     &              (lonfi(ig)*180./pi.ge.val3)  .and.
    2751      &              (lonfi(ig)*180./pi.le.val4)  ) then 
    2752                
     2752     &              (lonfi(ig)*180./pi.le.val4)  ) then
     2753
    27532754              IF (   (phisfi(ig).ge.val5)  .and.
    27542755     &               (phisfi(ig).le.val6)  ) then
    … …  
    27752776          end do
    27762777
    2777 c       changetopo 
     2778c       changetopo
    27782779c       -----------------------------------------------------
    27792780        else if (modif(1:len_trim(modif)) .eq. 'changetopo') then
    … …  
    27942795 577      read(*,*,iostat=ierr) val8
    27952796          if(ierr.ne.0) goto 577
    2796        
     2797
    27972798          write(*,*) 'phis=',phis
    27982799          write(*,*) 'phisfi=',phisfi
    … …  
    28012802     &              (latfi(ig)*180./pi.le.val2)  .and.
    28022803     &              (lonfi(ig)*180./pi.ge.val3)  .and.
    2803      &              (lonfi(ig)*180./pi.le.val4)  ) then 
    2804                
     2804     &              (lonfi(ig)*180./pi.le.val4)  ) then
     2805
    28052806              IF (   (phisfi(ig).ge.val5)  .and.
    28062807     &               (phisfi(ig).le.val6)  ) then
    … …  
    28632864          end do
    28642865
    2865 c       asymtopo: 
     2866c       asymtopo:
    28662867c       -----------------------------------------------------
    28672868        else if (modif(1:len_trim(modif)) .eq. 'asymtopo') then
    … …  
    29232924           IF (   (latfi(ig)*180./pi.ge.val2)  .and.
    29242925     &            (latfi(ig)*180./pi.le.val3)  .and.
    2925      &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and. 
     2926     &      (  ((lonfi(ig)*180./pi.ge.-180.)  .and.
    29262927     &                       (lonfi(ig)*180./pi.le.val5))  .or.
    2927      &         ((lonfi(ig)*180./pi.ge.val4)  .and. 
     2928     &         ((lonfi(ig)*180./pi.ge.val4)  .and.
    29282929     &                       (lonfi(ig)*180./pi.le.180.))) ) then
    29292930
    … …  
    29372938                    DO l=1,nsoilmx
    29382939                       tsoil(ig,l)=tempsoil(l)
    2939                     ENDDO 
    2940                 ENDIF 
     2940                    ENDDO
     2941                ENDIF
    29412942
    29422943                ! IF high topo : rm N2
    … …  
    29762977
    29772978           ENDIF
    2978        
     2979
    29792980         END DO
    29802981
    … …  
    29842985             do j=1,jjp1
    29852986               do i=1,iip1
    2986                  if (phis(i,j).gt.0.1) then 
     2987                 if (phis(i,j).gt.0.1) then
    29872988                       phis(i,j) = 0.
    29882989                       ps(i,j)=ps(iim/4,j)    ! assuming no topo here
    2989                  endif 
     2990                 endif
    29902991               enddo
    29912992             enddo
    … …  
    29993000 707      read(*,*,iostat=ierr) lat0, lon0
    30003001          if(ierr.ne.0) goto 707
    3001           if(lon0.gt.180.) lon0=lon0-360. 
     3002          if(lon0.gt.180.) lon0=lon0-360.
    30023003          lat0 = lat0 * pi/180.
    30033004          lon0 = lon0 * pi/180.
    … …  
    30293030               if (dist.le.radm) then
    30303031                  phisprev= phis(i,j)
    3031                   phis(i,j)=phisprev+1000.*g*height*(radm-dist)/radm 
     3032                  phis(i,j)=phisprev+1000.*g*height*(radm-dist)/radm
    30323033                  write(*,*) 'lat=',rlatu(j)*180./pi
    30333034                  write(*,*) 'lon=',rlonv(i)*180./pi
    … …  
    30363037                  ps(i,j) = ps(i,j) *
    30373038     .            exp((phis(i,j))/(temp * r))
    3038                end if 
     3039               end if
    30393040            end do
    30403041          end do
    … …  
    30723073               write(*,*)'Thermal inertia set for all subsurface layers'
    30733074               ierr=0
    3074               else 
     3075              else
    30753076               write(*,*)'Depth should be more than ',layer(1)
    30763077               ierr=1
    … …  
    30833084         if(val2.eq.0) then
    30843085           iref=1
    3085            write(*,*)'Level selected is first level: ',layer(iref),' m' 
     3086           write(*,*)'Level selected is first level: ',layer(iref),' m'
    30863087         else
    30873088           do isoil=1,nsoilmx-1
    … …  
    30933094            endif
    30943095           enddo
    3095          endif 
     3096         endif
    30963097
    30973098         DO ig=1,ngridmx
    … …  
    31333134               write(*,*)'Thermal inertia set for all subsurface layers'
    31343135               ierr=0
    3135               else 
     3136              else
    31363137               write(*,*)'Depth should be more than ',layer(1)
    31373138               ierr=1
    … …  
    31443145         if(val2.eq.0) then
    31453146           iref=1
    3146            write(*,*)'Level selected is first level: ',layer(iref),' m' 
     3147           write(*,*)'Level selected is first level: ',layer(iref),' m'
    31473148         else
    31483149           do isoil=1,nsoilmx-1
    … …  
    31543155            endif
    31553156           enddo
    3156          endif 
     3157         endif
    31573158
    31583159         DO ig=1,ngridmx
    … …  
    31953196               write(*,*)'Thermal inertia set for all subsurface layers'
    31963197               ierr=0
    3197               else 
     3198              else
    31983199               write(*,*)'Depth should be more than ',layer(1)
    31993200               ierr=1
    … …  
    32063207         if(val3.eq.0) then
    32073208           iref=1
    3208            write(*,*)'Level selected is first level: ',layer(iref),' m' 
     3209           write(*,*)'Level selected is first level: ',layer(iref),' m'
    32093210         else
    32103211           do isoil=1,nsoilmx-1
    … …  
    32163217            endif
    32173218           enddo
    3218          endif 
     3219         endif
    32193220
    32203221         DO ig=1,ngridmx
    … …  
    32553256               write(*,*)'Thermal inertia set for all subsurface layers'
    32563257               ierr=0
    3257               else 
     3258              else
    32583259               write(*,*)'Depth should be more than ',layer(1)
    32593260               ierr=1
    … …  
    32663267         if(val2.eq.0) then
    32673268           iref=1
    3268            write(*,*)'Level selected is first level: ',layer(iref),' m' 
     3269           write(*,*)'Level selected is first level: ',layer(iref),' m'
    32693270         else
    32703271           do isoil=1,nsoilmx-1
    … …  
    32783279            endif
    32793280           enddo
    3280          endif 
     3281         endif
    32813282
    32823283         DO ig=1,ngridmx
    … …  
    33993400c       -----------------------------------------------------
    34003401        else if (modif(1:len_trim(modif)) .eq. 'copylat') then
    3401            
     3402
    34023403          write(*,*) 'all latitudes : ',rlatu*180./pi
    34033404          write(*,*) 'Choice band to be modified lat1 ?'
    … …  
    34323433                       tsoil(ig,l)=tsoil(randtab(k),l)*val3
    34333434                       tsoil(ig,l)=tsoil(ig,l)+val4
    3434                 ENDDO 
     3435                ENDDO
    34353436                k=k+1
    34363437             ENDIF
    … …  
    34423443c       -----------------------------------------------------
    34433444        else if (modif(1:len_trim(modif)) .eq. 'copylon') then
    3444            
     3445
    34453446          write(*,*) 'all longitudes : ',rlonu*180./pi
    34463447          write(*,*) 'Choice band to be modified lon1 ?'
    … …  
    34693470                       tsoil(ig,l)=tsoil(randtab(k),l)
    34703471                ENDDO
    3471                 qsurf(ig,1)=qsurf(randtab(k),1) 
    3472                 phisfi(ig)=phisfi(randtab(k)) 
     3472                qsurf(ig,1)=qsurf(randtab(k),1)
     3473                phisfi(ig)=phisfi(randtab(k))
    34733474                k=k+1
    34743475             ENDIF
    … …  
    34813482c       -----------------------------------------------------
    34823483        else if (modif(1:len_trim(modif)) .eq. 'copylatlon') then
    3483            
     3484
    34843485          write(*,*) 'all longitudes : ',rlonu*180./pi
    34853486          write(*,*) 'Choice lat/lon to be copied ?'
    … …  
    35473548c       -----------------------------------------------------
    35483549        else if (modif(1:len_trim(modif)) .eq. 'albedomap') then
    3549            
     3550
    35503551          ! Get field 2D
    35513552          fichnom = 'albedo.nc'
    … …  
    35573558          ENDIF
    35583559
    3559           ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"), 
     3560          ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"),
    35603561     &                                  nvarid_input)
    35613562          IF (ierr .NE. NF_NOERR) THEN
    … …  
    35733574           PRINT*, "Could not get asked variable in albedo.nc"
    35743575           CALL abort
    3575           ELSE 
     3576          ELSE
    35763577           PRINT*, "Got variable in albedo.nc"
    35773578          ENDIF
    … …  
    35963597         temp=40.
    35973598         write(*,*) 'r, sale height temperature =',r,temp
    3598                  
     3599
    35993600         do j=1,jjp1
    36003601           do i=1,iip1
    … …  
    36103611                    ps(i,j) = ps(i,j) *
    36113612     .              exp((phisprev-phis(i,j))/(temp * r))
    3612               ENDIF 
     3613              ENDIF
    36133614           end do
    36143615         end do
    … …  
    36403641         temp=40.
    36413642         write(*,*) 'r, sale height temperature =',r,temp
    3642                  
     3643
    36433644         do j=1,jjp1
    36443645           do i=1,iip1
    … …  
    36533654                    ps(i,j) = ps(i,j) *
    36543655     .              exp((phisprev-phis(i,j))/(temp * r))
    3655               ENDIF 
     3656              ENDIF
    36563657           end do
    36573658         end do
    … …  
    36653666c       -----------------------------------------------------
    36663667        else if (modif(1:len_trim(modif)) .eq. 'copytsoil') then
    3667            
     3668
    36683669          ! Get field 2D
    36693670          fichnom = 'startfi_input.nc'
    … …  
    37023703           PRINT*, "Could not get asked variable in startfi_input.nc"
    37033704           CALL abort
    3704           ELSE 
     3705          ELSE
    37053706           PRINT*, "Got variable in startfi_input.nc"
    37063707          ENDIF
    … …  
    37143715           PRINT*, "Could not get asked variable s in startfi_input.nc"
    37153716           CALL abort
    3716           ELSE 
     3717          ELSE
    37173718           PRINT*, "Got variable s in startfi_input.nc"
    37183719          ENDIF
    … …  
    38043805             write(*,*) 'Newstart : you should set jm_plu to 180'
    38053806             stop
    3806          ENDIF 
     3807         ENDIF
    38073808         ! Read values
    38083809         do j=1,jm_plu+1
    … …  
    38923893         alb=alb_gcm
    38933894         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi)
    3894          !print*, 'N2dat=',N2_pluto_gcm 
    3895          !print*, 'COdat=',CO_pluto_gcm 
    3896          !print*, 'CH4dat=',CH4_pluto_gcm 
     3895         !print*, 'N2dat=',N2_pluto_gcm
     3896         !print*, 'COdat=',CO_pluto_gcm
     3897         !print*, 'CH4dat=',CH4_pluto_gcm
    38973898         print*, 'N2dat=',qsurf_tmp(:,igcm_n2)
    38983899         print*, 'COdat=',qsurf_tmp(:,igcm_co_ice)
    … …  
    39103911     &                               qsurf_tmp(ig,igcm_co_ice)
    39113912           endif
    3912          ENDDO       
     3913         ENDDO
    39133914
    39143915        endif
    3915              
     3916
    39163917       enddo ! of do ! infinite loop on liste of changes
    39173918
    … …  
    39433944c=======================================================================
    39443945
    3945       CALL inifilr 
     3946      CALL inifilr
    39463947      CALL pression(ip1jmp1, ap, bp, ps, p3d)
    39473948
    … …  
    40054006     *                phi,w, pbaru,pbarv,day_ini+time )
    40064007
    4007          
     4008
    40084009      CALL dynredem0("restart.nc",day_ini,phis)
    4009       CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps) 
     4010      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
    40104011C
    40114012C Ecriture etat initial physique
    … …  
    40234024
    40244025c=======================================================================
    4025 c        Formats 
     4026c        Formats
    40264027c=======================================================================
    40274028
Note: See TracChangeset for help on using the changeset viewer.