Changeset 164 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Jun 17, 2011, 10:49:17 AM (14 years ago)
Author:
emillour
Message:

Mars GCM:

Updates and corrections (to enable compiling/running in debug mode with

ifort)

  • removed option "-free-form" from makegcm_ifort and set mod_loc_dir="." so that module files (produced in local directory by ifort) are moved to LIBO
  • updated initracer.F, physdem1.F, physiq.F, inichim_newstart.F to avoid referencing out-of-bound array indexes (even if unused)
  • cosmetic updates on inwrite.F, datareadnc.F
  • updated newstart.F to initialize and use 'datadir' when looking for files
  • corrected bug on interpolation of sub-surface temperatures in lect_start_archive.F

EM

Location:
trunk/LMDZ.MARS/libf
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F

    r38 r164  
    120120              write(*,*) 'initracer: error expected dustbin=2'
    121121            else
    122               noms(1)='dust_mass'   ! dust mass mixing ratio
    123               noms(2)='dust_number' ! dust number mixing ratio
     122!              noms(1)='dust_mass'   ! dust mass mixing ratio
     123!              noms(2)='dust_number' ! dust number mixing ratio
     124! same as above, but avoid explict possible out-of-bounds on array
     125               noms(1)='dust_mass'   ! dust mass mixing ratio
     126               do iq=2,2
     127               noms(iq)='dust_number' ! dust number mixing ratio
     128               enddo
    124129            endif
    125130          endif
    … …  
    129134          noms(nqmx)='h2o_vap'
    130135          mmol(nqmx)=18.
    131           noms(nqmx-1)='h2o_ice'
    132           mmol(nqmx-1)=18.
     136!            noms(nqmx-1)='h2o_ice'
     137!            mmol(nqmx-1)=18.
     138          !"loop" to avoid potential out-of-bounds in array
     139          do iq=nqmx-1,nqmx-1
     140            noms(iq)='h2o_ice'
     141            mmol(iq)=18.
     142          enddo
    133143        endif
    134144        ! 3. Chemistry
    … …  
    164174        write(*,*)'inichim_newstart: moving surface water ice to index '
    165175     &            ,nqmx-1
    166         qsurf(1:ngridmx,nqmx-1)=qsurf(1:ngridmx,nqmx)
     176        ! "loop" to avoid potential out-of-bounds on arrays
     177        do iq=nqmx-1,nqmx-1
     178        qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
     179        enddo
    167180        qsurf(1:ngridmx,nqmx)=0
    168181      endif
    … …  
    395408      ! as in the old days, water vapour is last and water ice,
    396409      ! if present, is just before water vapour
    397       nqchem(1)=igcm_co2
    398       nqchem(2)=igcm_co
    399       nqchem(3)=igcm_o
    400       nqchem(4)=igcm_o1d
    401       nqchem(5)=igcm_o2
    402       nqchem(6)=igcm_o3
    403       nqchem(7)=igcm_h
    404       nqchem(8)=igcm_h2
    405       nqchem(9)=igcm_oh
    406       nqchem(10)=igcm_ho2
    407       nqchem(11)=igcm_h2o2
    408       nqchem(12)=igcm_n2
    409       nqchem(13)=igcm_ar
    410       nqchem(14)=igcm_h2o_ice
    411       nqchem(15)=igcm_h2o_vap
     410      do iq=1,15 ! loop so as to avoid out-of-bounds on array
     411      if (iq==1) nqchem(i)=igcm_co2
     412      if (iq==2) nqchem(i)=igcm_co
     413      if (iq==3) nqchem(i)=igcm_o
     414      if (iq==4) nqchem(i)=igcm_o1d
     415      if (iq==5) nqchem(i)=igcm_o2
     416      if (iq==6) nqchem(i)=igcm_o3
     417      if (iq==7) nqchem(i)=igcm_h
     418      if (iq==8) nqchem(i)=igcm_h2
     419      if (iq==9) nqchem(i)=igcm_oh
     420      if (iq==10) nqchem(i)=igcm_ho2
     421      if (iq==11) nqchem(i)=igcm_h2o2
     422      if (iq==12) nqchem(i)=igcm_n2
     423      if (iq==13) nqchem(i)=igcm_ar
     424      if (iq==14) nqchem(i)=igcm_h2o_ice
     425      if (iq==15) nqchem(i)=igcm_h2o_vap
     426      enddo
    412427
    413428! Load in chemistry data for initialization:

  • trunk/LMDZ.MARS/libf/dyn3d/lect_start_archive.F

    r146 r164  
    11361136        allocate(oldval(nsoilold+1))
    11371137        allocate(newval(nsoilmx))
    1138         do i=1,imold+1
    1139          do j=1,jmold+1
     1138        do i=1,iip1
     1139         do j=1,jjp1
    11401140           ! copy values
    11411141           oldval(1)=tsurfold(i,j)
    … …  
    11681168        oldgrid(1)=0. ! ground
    11691169        oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
    1170         do i=1,imold+1
    1171          do j=1,jmold+1
     1170        do i=1,iip1
     1171         do j=1,jjp1
    11721172           ! copy values
    11731173           oldval(1)=tsurfold(i,j)

  • trunk/LMDZ.MARS/libf/dyn3d/newstart.F

    r38 r164  
    1414c
    1515c=======================================================================
     16
     17! to use  'getin'
     18      USE ioipsl_getincom
    1619
    1720      implicit none
    … …  
    4043#include"advtrac.h"
    4144#include"tracer.h"
     45#include "datafile.h"
    4246c=======================================================================
    4347c   Declarations
    … …  
    5054c et autres:
    5155c----------
    52       INTEGER lnblnk
    53       EXTERNAL lnblnk
    5456
    5557c Variables pour les lectures NetCDF des fichiers "start_archive"
    … …  
    358360        relief="mola"
    359361c     enddo
     362
     363! before using datareadnc, "datafile" must be set (normaly done in inifis)
     364      datafile="/u/forget/WWW/datagcm/datafile" ! default value
     365      call getin("datadir",datafile) ! in case user specified another path
    360366
    361367      CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
    … …  
    499505
    500506        write(*,*)
    501         write(*,*) modif(1:lnblnk(modif)) , ' : '
     507        write(*,*) trim(modif) , ' : '
    502508
    503509c       'flat : no topography ("aquaplanet")'
    504510c       -------------------------------------
    505         if (modif(1:lnblnk(modif)) .eq. 'flat') then
     511        if (trim(modif) .eq. 'flat') then
    506512c         set topo to zero
    507513          CALL initial0(ip1jmp1,z_reel)
    … …  
    537543c       bilball : albedo, inertie thermique uniforme
    538544c       --------------------------------------------
    539         else if (modif(1:lnblnk(modif)) .eq. 'bilball') then
     545        else if (trim(modif) .eq. 'bilball') then
    540546          write(*,*) 'constante albedo and iner.therm:'
    541547          write(*,*) 'New value for albedo (ex: 0.25) ?'
    … …  
    564570c       coldspole : sous-sol de la calotte sud toujours froid
    565571c       -----------------------------------------------------
    566         else if (modif(1:lnblnk(modif)) .eq. 'coldspole') then
     572        else if (trim(modif) .eq. 'coldspole') then
    567573          write(*,*)'new value for the subsurface temperature',
    568574     &   ' beneath the permanent southern polar cap ? (eg: 141 K)'
    … …  
    615621c       ptot : Modification of the total pressure: ice + current atmosphere
    616622c       -------------------------------------------------------------------
    617         else if (modif(1:lnblnk(modif)) .eq. 'ptot') then
     623        else if (trim(modif) .eq. 'ptot') then
    618624
    619625c         calcul de la pression totale glace + atm actuelle
    … …  
    699705c       q=0 : set tracers to zero
    700706c       -------------------------
    701         else if (modif(1:lnblnk(modif)) .eq. 'q=0') then
     707        else if (trim(modif) .eq. 'q=0') then
    702708c          mise a 0 des q (traceurs)
    703709          write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
    … …  
    721727c       q=x : initialise tracer manually
    722728c       --------------------------------
    723         else if (modif(1:lnblnk(modif)) .eq. 'q=x') then
     729        else if (trim(modif) .eq. 'q=x') then
    724730             write(*,*) 'Which tracer do you want to modify ?'
    725731             do iq=1,nqmx
    … …  
    747753c       ini_q : Initialize tracers for chemistry
    748754c       -----------------------------------------------
    749         else if (modif(1:lnblnk(modif)) .eq. 'ini_q') then
     755        else if (trim(modif) .eq. 'ini_q') then
    750756c         For more than 32 layers, possible to initiate thermosphere only     
    751757          thermo=0
    … …  
    779785c       ini_q-H2O : as above exept for the water vapour tracer
    780786c       ------------------------------------------------------
    781         else if (modif(1:lnblnk(modif)) .eq. 'ini_q-H2O') then
     787        else if (trim(modif) .eq. 'ini_q-H2O') then
    782788          ! for more than 32 layers, possible to initiate thermosphere only     
    783789          thermo=0
    … …  
    812818c       ini_q-iceH2O : as above exept for ice et H2O
    813819c       -----------------------------------------------
    814         else if (modif(1:lnblnk(modif)) .eq. 'ini_q-iceH2O') then
     820        else if (trim(modif) .eq. 'ini_q-iceH2O') then
    815821c         For more than 32 layers, possible to initiate thermosphere only     
    816822          thermo=0
    … …  
    846852c      wetstart : wet atmosphere with a north to south gradient
    847853c      --------------------------------------------------------
    848        else if (modif(1:lnblnk(modif)) .eq. 'wetstart') then
     854       else if (trim(modif) .eq. 'wetstart') then
    849855        ! check that there is indeed a water vapor tracer
    850856        if (igcm_h2o_vap.eq.0) then
    … …  
    867873c      noglacier : remove tropical water ice (to initialize high res sim)
    868874c      --------------------------------------------------
    869         else if (modif(1:lnblnk(modif)) .eq. 'noglacier') then
     875        else if (trim(modif) .eq. 'noglacier') then
    870876           do ig=1,ngridmx
    871877             j=(ig-2)/iim +2
    … …  
    880886c      watercapn : H20 ice on permanent northern cap
    881887c      --------------------------------------------------
    882         else if (modif(1:lnblnk(modif)) .eq. 'watercapn') then
     888        else if (trim(modif) .eq. 'watercapn') then
    883889           do ig=1,ngridmx
    884890             j=(ig-2)/iim +2
    … …  
    895901c      watercaps : H20 ice on permanent southern cap
    896902c      -------------------------------------------------
    897         else if (modif(1:lnblnk(modif)) .eq. 'watercaps') then
     903        else if (trim(modif) .eq. 'watercaps') then
    898904           do ig=1,ngridmx
    899905               j=(ig-2)/iim +2
    … …  
    910916c       isotherm : Isothermal temperatures and no winds
    911917c       ------------------------------------------------
    912         else if (modif(1:lnblnk(modif)) .eq. 'isotherm') then
     918        else if (trim(modif) .eq. 'isotherm') then
    913919
    914920          write(*,*)'Isothermal temperature of the atmosphere,
    … …  
    933939c       co2ice=0 : remove CO2 polar ice caps'
    934940c       ------------------------------------------------
    935         else if (modif(1:lnblnk(modif)) .eq. 'co2ice=0') then
     941        else if (trim(modif) .eq. 'co2ice=0') then
    936942           do ig=1,ngridmx
    937943              co2ice(ig)=0
    … …  
    942948!       ----------------------------------------------------------------------
    943949
    944         else if (modif(1:lnblnk(modif)).eq.'therm_ini_s') then
     950        else if (trim(modif).eq.'therm_ini_s') then
    945951!          write(*,*)"surfithfi(1):",surfithfi(1)
    946952          do isoil=1,nsoilmx
    … …  
    965971!       ------------------------------------------------------------
    966972
    967         else if (modif(1:lnblnk(modif)).eq.'subsoilice_n') then
     973        else if (trim(modif).eq.'subsoilice_n') then
    968974
    969975         write(*,*)'From which latitude (in deg.), up to the north pole,
    … …  
    10781084!       ------------------------------------------------------------
    10791085
    1080         else if (modif(1:lnblnk(modif)).eq.'subsoilice_s') then
     1086        else if (trim(modif).eq.'subsoilice_s') then
    10811087
    10821088         write(*,*)'From which latitude (in deg.), down to the south pol
    … …  
    11791185c       'mons_ice' : use MONS data to build subsurface ice table
    11801186c       --------------------------------------------------------
    1181         else if (modif(1:lnblnk(modif)).eq.'mons_ice') then
     1187        else if (trim(modif).eq.'mons_ice') then
    11821188       
    11831189       ! 1. Load MONS data
    … …  
    12851291        else
    12861292          write(*,*) '       Unknown (misspelled?) option!!!'
    1287         end if ! of if (modif(1:lnblnk(modif)) .eq. '...') elseif ...
     1293        end if ! of if (trim(modif) .eq. '...') elseif ...
    12881294       
    12891295       enddo ! of do ! infinite loop on liste of changes

  • trunk/LMDZ.MARS/libf/phymars/datareadnc.F

    r146 r164  
    123123        write(*,*)'1) You can set this path in the callphys.def file:'
    124124        write(*,*)'   datadir=/path/to/the/datafiles'
    125         write(*,*)'2) If necessary surface.nc (and other datafiles)'
     125        write(*,*)'   OR specify the path to use in callphys.def, as:'
     126        write(*,*)'   datadir=/path/to/the/directory'
     127        write(*,*)'2) If necessary, surface.nc (and other datafiles)'
    126128        write(*,*)'   can be obtained online on:'
    127129        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'

  • trunk/LMDZ.MARS/libf/phymars/initracer.F

    r91 r164  
    111111              write(*,*) 'initracer: error expected dustbin=2'
    112112            else
    113               noms(1)='dust_mass'   ! dust mass mixing ratio
    114               noms(2)='dust_number' ! dust number mixing ratio
     113!              noms(1)='dust_mass'   ! dust mass mixing ratio
     114!              noms(2)='dust_number' ! dust number mixing ratio
     115! same as above, but avoid explict possible out-of-bounds on array
     116               noms(1)='dust_mass'   ! dust mass mixing ratio
     117               do iq=2,2
     118               noms(iq)='dust_number' ! dust number mixing ratio
     119               enddo
    115120            endif
    116121          endif
    … …  
    120125          noms(nqmx)='h2o_vap'
    121126          mmol(nqmx)=18.
    122           noms(nqmx-1)='h2o_ice'
    123           mmol(nqmx-1)=18.
     127!            noms(nqmx-1)='h2o_ice'
     128!            mmol(nqmx-1)=18.
     129          !"loop" to avoid potential out-of-bounds in array
     130          do iq=nqmx-1,nqmx-1
     131            noms(iq)='h2o_ice'
     132            mmol(iq)=18.
     133          enddo
    124134        endif
    125135        ! 3. Chemistry
    … …  
    157167      if (oldnames.and.water) then
    158168        write(*,*)'initracer: moving surface water ice to index ',nqmx-1
    159         qsurf(1:ngridmx,nqmx-1)=qsurf(1:ngridmx,nqmx)
     169        ! "loop" to avoid potential out-of-bounds on arrays
     170        do iq=nqmx-1,nqmx-1
     171          qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
     172        enddo
    160173        qsurf(1:ngridmx,nqmx)=0
    161174      endif

  • trunk/LMDZ.MARS/libf/phymars/iniwrite.F

    r38 r164  
    3636c   ----------
    3737
    38       integer nid        ! NetCDF file ID
    39       INTEGER*4 idayref  ! date (initial date for this run)
    40       REAL phis(ip1jmp1) ! surface geopotential
     38      integer,intent(in) :: nid        ! NetCDF file ID
     39      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
     40      real,intent(in) :: phis(ip1jmp1) ! surface geopotential
    4141
    4242c   Local:

  • trunk/LMDZ.MARS/libf/phymars/physdem1.F

    r38 r164  
    565565        ! back to qsurf(nqmx)
    566566        IF (water) THEN
     567          !"loop" to avoid potential out-of-bounds on arrays
    567568          write(*,*)'physdem1: moving surface water ice to index ',nqmx
    568           qsurf(1:ngridmx,nqmx)=qsurf(1:ngridmx,nqmx-1)
    569           qsurf(1:ngridmx,nqmx-1)=0
     569          do iq=nqmx,nqmx
     570          qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq-1)
     571          qsurf(1:ngridmx,iq-1)=0
     572          enddo
    570573        ENDIF
    571574      endif ! of if (count.eq.nqmx)

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r161 r164  
    13081308c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
    13091309c    &                  emis)
    1310          call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
    1311          call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
     1310!         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
     1311!         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
    13121312         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
    13131313     &                  tsurf)
    13141314         call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
    1315 c        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
    1316 c    &                  co2ice)
     1315        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
     1316     &                  co2ice)
    13171317c         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
    13181318c     &                  "K",2,zt(1,7))
    … …  
    13311331         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
    13321332c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
    1333         call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
     1333!        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
    13341334c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
    13351335c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
    … …  
    15361536     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
    15371537         if (nqmx .gt. 1) then
     1538             iq=2 ! to avoid out-of-bounds spotting by picky compilers
     1539                  ! when gcm is compiled with only one tracer
    15381540             dummycol(ig)=dummycol(ig)+
    1539      &                  zq(ig,l,2)*(pplev(ig,l)-pplev(ig,l+1))/g
     1541     &                  zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
    15401542         endif
    15411543         enddo
Note: See TracChangeset for help on using the changeset viewer.