Changeset 164 for trunk


Ignore:
Timestamp:
Jun 17, 2011, 10:49:17 AM (13 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
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/README

    r148 r164  
    664664   Modified phymars/dimradmars.h , added directory phymars/scatterers
    665665   with script make_scatterers , and adapted makegcm* scripts.
     666
     667== 17/06/2011 == AC
     668================================================
     669======== IMPLEMENTATION OF THERMALS ============
     670================================================
     671The main goal of this revision is to start including the thermals into the model
     672for development purposes. Users should not use the thermals yet, as
     673several major configuration changes still need to be done.
     674
     675This version includes :
     676- updraft and downdraft parametrizations
     677- velocity in the thermal, including drag
     678- plume height analysis
     679- closure equation
     680- updraft transport of heat, tracers and momentum
     681- downdraft transport of heat
     682
     683This model should not be used without upcoming developments, namely :
     684- downdraft transport of tracers and momentum
     685- updraft & downdraft transport of q2 (tke)
     686- revision of vdif_kc to compute q2 for non-stratified cases
     687
     688Thermals could also include in a later revision :
     689- momentum loss during transport (horizontal drag)
     690
     691Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90
     692
     693================================================
     694================================================
     695
     696M      libf/phymars/callkeys.h
     697M      libf/phymars/inifis.F
     698
     699Added new control flags to call the thermals :
     700- calltherm (false by default) <- to call thermals
     701- outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)
     702================================================
     703M      libf/phymars/vdifc.F
     704^------> added a temporary output for thermal-related diagnostics
     705M      libf/phymars/testphys1d.F
     706^------> added treatment for a initialization from a profile of neutral gas (ar)
     707      -> will be transformed in a decaying tracer for thermal diagnostics
     708M      libf/phymars/physiq.F
     709^------> added a section to call the thermals
     710      -> changed the call to convadj
     711      -> added thermal-related outputs for diagnostics
     712M      libf/phymars/convadj.F
     713^------> takes now into account the height of thermals to execute convective adjustment
     714      => note : convective adjustment needs to be activated when using thermals, in case of a
     715         second instable layer above the thermals
     716================================================
     717A      libf/phymars/calltherm_interface.F90
     718^------> Interface between physiq.F and the thermals
     719A      libf/phymars/calltherm_mars.F90
     720^------> Routine running the sub-timestep of the thermals
     721A      libf/phymars/thermcell_main_mars.F90
     722^------> Main thermals routine specific to Martian physics
     723A      libf/phymars/thermcell_dqupdown.F90
     724^------> Thermals subroutine computing transport of quantities by updrafts and downdrafts
     725A      libf/phymars/thermcell.F90
     726^------> Module including parameters from the Earth to Mars importation. Will disappear in future dev
     727================================================
     728================================================
     729
     730== 17/06/2011 == EM
     731>>> Updates and corrections (to enable compiling/running in debug mode with
     732    ifort)
     733- removed option "-free-form" from makegcm_ifort and set mod_loc_dir="."
     734  so that module files (produced in local directory by ifort)
     735  are moved to LIBO
     736- updated initracer.F, physdem1.F, physiq.F, inichim_newstart.F
     737  to avoid referencing out-of-bound array indexes (even if unused)
     738- cosmetic updates on inwrite.F, datareadnc.F
     739- updated newstart.F to initialize and use 'datadir' when looking for files
     740- corrected bug on interpolation of sub-surface temperatures in
     741  lect_start_archive.F
  • 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
  • trunk/LMDZ.MARS/makegcm_ifort

    r148 r164  
    2222# default LMDGCM to where makegcm script is located:
    2323#setenv LMDGCM "`dirname $0`"
     24setenv LMDGCM /san/home/millour/Planeto/test_mars_141/LMDZ.MARS
    2425# You may set LIBOGCM to something else; otherwise we default to:
    25 #setenv LIBOGCM $LMDGCM/libo
     26setenv LIBOGCM $LMDGCM/libo
    2627## NetCDF Libraries: what follows is OK on GNOME
    27 #  setenv NCDFLIB /usr/local/lib
    28 #  setenv NCDFINC /usr/local/include
     28  setenv NCDFLIB /usr/local/lib
     29  setenv NCDFINC /usr/local/include
    2930####
    3031
     
    180181   set optimtru90=" -O2 -ip -mkl=sequential -align common "
    181182#   set opt_link=" -Mfree -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Bstatic "
    182    set mod_loc_dir=$LIBOGCM
     183#   set mod_loc_dir=$LIBOGCM
     184# ifort puts mod files in local directory
     185   set mod_loc_dir="."
    183186   set mod_suffix=mod
    184187else if $NEC then
     
    709712 set optim="$optim -I${libo}"
    710713 set optim90="$optim90 -I${libo}"
    711  set optimtru90="$optimtru90 -ffree-form -I${libo}"
     714 set optimtru90="$optimtru90 -I${libo}"
    712715# Ehouarn: remove set mod_loc_dir def below; mod_loc_dir=$localdir (set above)
    713716# set mod_loc_dir=$libo
Note: See TracChangeset for help on using the changeset viewer.