Changeset 552


Ignore:
Timestamp:
Mar 2, 2012, 9:53:27 AM (13 years ago)
Author:
emillour
Message:

Mars GCM:

bug fix in "aeronomars/moldiff_red.F90" (wrong mmol array size in included

subroutine) + adapted extreme limit test for H to altitudes of ~4000km
(compatble with 50 layer model).

bug fix in "nirco2abs.F" index of "co2" and "o" tracers were hard coded as

1 and 3!

updates from FGG of euvheat.F, callkeys.h and inifis.F to have the

"euveff" parameter read from callphys.def

EM

Location:
trunk/LMDZ.MARS
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r544 r552  
    14261426   ... Syntax for use is to add "tke_heat_flux = ..." in callphys.def
    14271427>> Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)
     1428
     1429== 02/03/12 == EM
     1430>> bug fix in "aeronomars/moldiff_red.F90" (wrong mmol array size in included
     1431   subroutine) + adapted extreme limit test for H to altitudes of ~4000km
     1432   (compatble with 50 layer model).
     1433>> bug fix in "nirco2abs.F" index of "co2" and "o" tracers were hard coded as
     1434   1 and 3!
     1435>> updates from FGG of euvheat.F, callkeys.h and inifis.F to have the
     1436   "euveff" parameter read from callphys.def
     1437
     1438 

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

    r38 r552  
    193193
    194194        do l=1,nlayermx
    195           pdteuv(ig,l)=0.16*jtot(l)/10.
     195          pdteuv(ig,l)=euveff*jtot(l)/10.
    196196     &            /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l)))
    197197     &                  *(1.52/dist_sol)**2

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

    r460 r552  
    7272     &           1.,17.,33.,18.,34.,16.,48./   ! minors
    7373      character*3 tmp ! temporary variable
    74       integer ierr,lnblnk
    75       external lnblnk
     74      integer ierr
    7675
    7776      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
    … …  
    476475c----------------------------------------------------------------------
    477476      open(210, iostat=ierr,file=
    478      & datafile(1:lnblnk(datafile))//'/atmosfera_LMD_may.dat')
     477     & trim(datafile)//'/atmosfera_LMD_may.dat')
    479478      if (ierr.ne.0) then
    480479        write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
    481480        write(*,*)'(in aeronomars/inichim.F)'
    482         write(*,*)'It should be in :', datafile(1:lnblnk(datafile)),'/'
    483         write(*,*)'1) You can change this directory address in '
    484         write(*,*)'   file phymars/datafile.h'
     481        write(*,*)'It should be in :', trim(datafile),'/'
     482        write(*,*)'1) You can change this path in callphys.def with'
     483        write(*,*)'   datadir=/path/to/datafiles/'
    485484        write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
    486485        write(*,*)'   can be obtained online on:'
    … …  
    489488      endif
    490489      open(220, iostat=ierr,file=
    491      & datafile(1:lnblnk(datafile))//'/atmosfera_LMD_min.dat')
     490     & trim(datafile)//'/atmosfera_LMD_min.dat')
    492491      if (ierr.ne.0) then
    493492        write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
    494493        write(*,*)'(in aeronomars/inichim.F)'
    495         write(*,*)'It should be in :', datafile(1:lnblnk(datafile)),'/'
    496         write(*,*)'1) You can change this directory address in '
    497         write(*,*)'   file phymars/datafile.h'
     494        write(*,*)'It should be in :', trim(datafile),'/'
     495        write(*,*)'1) You can change this path in callphys.def with'
     496        write(*,*)'   datadir=/path/to/datafiles/'
    498497        write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
    499498        write(*,*)'   can be obtained online on:'

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

    r517 r552  
    314314
    315315
    316 ! If Zmax > 2000 km there is a problem / stop
    317 
    318         if (Zmax .gt. 2000000.) then
     316! If Zmax > 4000 km there is a problem / stop
     317
     318        if (Zmax .gt. 4000000.) then
    319319        Print*,'Zmax too high',ig,zmax,zmin
    320320        do l=1,nlayermx
    … …  
    918918        SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig)
    919919        IMPLICIT NONE
     920#include "dimensions.h"
    920921
    921922        INTEGER :: nl,nq
    922923        INTEGER :: l,iq,ig
    923924        INTEGER,dimension(nq) :: gc
    924         REAL,DIMENSION(nl,nq) :: Q1,DQ
     925        REAL,DIMENSION(nl,nqmx) :: Q1,DQ
    925926        REAL*8,DIMENSION(nl,nq) :: Q2
    926927        REAL :: dtime
    … …  
    951952        SUBROUTINE MMOY(massemoy,mmol,qq,gc,nl,nq)
    952953        IMPLICIT NONE
     954#include "dimensions.h"
    953955
    954956        INTEGER :: nl,nq,l
    … …  
    956958        REAL*8,DIMENSION(nl,nq) :: qq
    957959        REAL*8,DIMENSION(nl) :: massemoy
    958         REAL,DIMENSION(nq) :: MMOL
     960        REAL,DIMENSION(nqmx) :: MMOL
    959961
    960962
    … …  
    10341036     & Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
    10351037        IMPLICIT NONE
     1038#include "dimensions.h"
    10361039       
    10371040        INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
    … …  
    10441047        REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
    10451048        REAL*8 :: facZ,dZ,H
    1046         REAL,DIMENSION(nq) :: mmol
     1049        REAL,DIMENSION(nqmx) :: mmol
    10471050        masseU=1.660538782d-27
    10481051        kbolt=1.3806504d-23
    … …  
    13541357     &    pp,M,gc,nl,nq,nlx)
    13551358        IMPLICIT NONE
     1359#include "dimensions.h"
    13561360        INTEGER :: nl,nq,nlx,il,nn,iP
    13571361        INTEGER,DIMENSION(1) :: indP
    … …  
    13591363        REAL*8,DIMENSION(nl) :: Z,P,T
    13601364        REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
    1361         REAL,DIMENSION(nq) :: M
     1365        REAL,DIMENSION(nqmx) :: M
    13621366        REAL*8,DIMENSION(nq) :: nNew
    13631367        REAL*8,DIMENSION(nlx) :: pp,tt,tnew
    … …  
    14471451     &   ,pp,M,gc,nl,nq,nlx,facM)
    14481452        IMPLICIT NONE
     1453#include "dimensions.h"
    14491454        INTEGER :: nl,nq,nlx,il,nn,iP
    14501455        INTEGER,DIMENSION(1) :: indP
    … …  
    14521457        REAL*8,DIMENSION(nl) :: Z,P,T
    14531458        REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
    1454         REAL,DIMENSION(nq) :: M
     1459        REAL,DIMENSION(nqmx) :: M
    14551460        REAL*8,DIMENSION(nq) :: nNew
    14561461        REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r544 r552  
    1818     &   ,dustbin,nqchem_min,nltemodel,nircorr
    1919     
    20       COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,            &
     20      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff,     &
    2121     &   tke_heat_flux
    2222     
    … …  
    3737      real alphan
    3838      real solarcondate
     39      real euveff
    3940      real tke_heat_flux
    4041

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

    r551 r552  
    249249     $              "1-> new correction",
    250250     $              "(matters only if callnirco2=T)"
    251          nircorr=0
     251         nircorr=1      !default value
    252252         call getin("nircorr",nircorr)
    253253         write(*,*) " nircorr = ",nircorr
    … …  
    594594         write(*,*) " solarcondate = ",solarcondate
    595595         
     596         write(*,*) "UV heating efficiency:",
     597     &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
     598     &   "lower values may be used to compensate low 15 um cooling"
     599         euveff=0.21 !default value
     600         call getin("euveff",euveff)
     601         write(*,*) " euveff = ", euveff
    596602
    597603         if (.not.callthermos) then

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

    r498 r552  
    5050#include "comdiurn.h"
    5151#include "nirdata.h"
     52#include "tracer.h"
    5253
    5354c-----------------------------------------------------------------------
    5455c    Input/Output
    5556c    ------------
    56       INTEGER ngrid,nlayer
    57 
    58       REAL pplay(ngrid,nlayer)
    59       REAL dist_sol,mu0(ngridmx),fract(ngridmx),declin
    60 
    61       REAL pdtnirco2(ngrid,nlayer)
     57      integer,intent(in) :: ngrid ! number of (horizontal) grid points
     58      integer,intent(in) :: nlayer ! number of atmospheric layers
     59      real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
     60      real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
     61      integer,intent(in) :: nq ! number of tracers
     62      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
     63      real,intent(in) :: mu0(ngridmx) ! solar angle
     64      real,intent(in) :: fract(ngridmx) ! day fraction of the time interval
     65      real,intent(in) :: declin ! latitude of sub-solar point
     66     
     67      real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
    6268c
    6369c    Local variables :
    … …  
    7379c   local saved variables
    7480c   ---------------------
    75 
     81      logical,save :: firstcall=.true.
     82      real,save :: ico2=0 ! index of "co2" tracer
     83      real,save :: io=0 ! index of "o" tracer
    7684c     p0noonlte is a pressure below which non LTE effects are significant.
    7785c     REAL p0nonlte
    … …  
    9098      real    p2011,cociente1,merge
    9199      real    cor0,oco2gcm
    92       integer nq
    93       real    pq(ngrid,nlayer,nq)
    94100
    95101c----------------------------------------------------------------------
    … …  
    97103c     Initialisation
    98104c     --------------
     105      if (firstcall) then
     106        if (nircorr.eq.1) then
     107          ! we will need co2 and o tracers
     108          ico2=igcm_co2
     109          if (ico2==0) then
     110            write(*,*) "nirco2abs error: I need a CO2 tracer"
     111            write(*,*) "     when running with nircorr==1"
     112            stop
     113          endif
     114          io=igcm_o
     115          if (io==0) then
     116            write(*,*) "nirco2abs error: I need an O tracer"
     117            write(*,*) "     when running with nircorr==1"
     118            stop
     119          endif
     120        endif
     121        firstcall=.false.
     122      endif
     123
     124
    99125c     co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
    100126      co2heat0=n_a*(1.52/dist_sol)**2/daysec
    … …  
    121147               if(nircorr.eq.1) then
    122148                  cor0=1./(1.+n_p0/pplay(ig,l))**n_b
    123                   if(pq(ig,l,1).gt.1.e-6) then
    124                      oco2gcm=pq(ig,l,3)/pq(ig,l,1)
     149                  if(pq(ig,l,ico2).gt.1.e-6) then
     150                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
    125151                  else
    126152                     oco2gcm=1.e6
    … …  
    179205                  if(nircorr.eq.1) then
    180206                     cor0=1./(1.+n_p0/pplay(ig,l))**n_b
    181                      oco2gcm=pq(ig,l,3)/pq(ig,l,1)
     207                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
    182208                     cociente1=oco2gcm/oldoco2(l)
    183209                     merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*

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

    r224 r552  
    8181c -----------------------
    8282      INTEGER setname, cluvdb, getdat
    83 
    84       INTEGER lnblnk
    8583
    8684c-----------------------------------------------------------------------
    … …  
    296294 
    297295 
    298       do while(modif(1:1).ne.'hello')
     296      do ! neverending loop
    299297        write(*,*)
    300298        write(*,*)
    … …  
    306304 
    307305        write(*,*)
    308         write(*,*) modif(1:lnblnk(modif)) , ' : '
    309 
    310         if (modif(1:lnblnk(modif)) .eq. 'day_ini') then
     306        write(*,*) trim(modif) , ' : '
     307
     308        if (trim(modif) .eq. 'day_ini') then
    311309          write(*,*) 'current value:',day_ini
    312310          write(*,*) 'enter new value:'
    … …  
    316314          write(*,*) 'day_ini (new value):',day_ini
    317315
    318         else if (modif(1:lnblnk(modif)) .eq. 'z0') then
     316        else if (trim(modif) .eq. 'z0') then
    319317          write(*,*) 'current value (m):',z0_default
    320318          write(*,*) 'enter new value (m):'
    … …  
    324322          write(*,*) ' z0 (new value):',z0_default
    325323
    326         else if (modif(1:lnblnk(modif)) .eq. 'emin_turb') then
     324        else if (trim(modif) .eq. 'emin_turb') then
    327325          write(*,*) 'current value:',emin_turb
    328326          write(*,*) 'enter new value:'
    … …  
    332330          write(*,*) ' emin_turb (new value):',emin_turb
    333331
    334         else if (modif(1:lnblnk(modif)) .eq. 'lmixmin') then
     332        else if (trim(modif) .eq. 'lmixmin') then
    335333          write(*,*) 'current value:',lmixmin
    336334          write(*,*) 'enter new value:'
    … …  
    340338          write(*,*) ' lmixmin (new value):',lmixmin
    341339
    342         else if (modif(1:lnblnk(modif)) .eq. 'emissiv') then
     340        else if (trim(modif) .eq. 'emissiv') then
    343341          write(*,*) 'current value:',emissiv
    344342          write(*,*) 'enter new value:'
    … …  
    348346          write(*,*) ' emissiv (new value):',emissiv
    349347
    350         else if (modif(1:lnblnk(modif)) .eq. 'emisice') then
     348        else if (trim(modif) .eq. 'emisice') then
    351349          write(*,*) 'current value emisice(1) North:',emisice(1)
    352350          write(*,*) 'enter new value:'
    … …  
    364362          write(*,*) ' emisice(2) (new value):',emisice(2)
    365363
    366         else if (modif(1:lnblnk(modif)) .eq. 'albedice') then
     364        else if (trim(modif) .eq. 'albedice') then
    367365          write(*,*) 'current value albedice(1) North:',albedice(1)
    368366          write(*,*) 'enter new value:'
    … …  
    380378          write(*,*) ' albedice(2) (new value):',albedice(2)
    381379
    382         else if (modif(1:lnblnk(modif)) .eq. 'iceradius') then
     380        else if (trim(modif) .eq. 'iceradius') then
    383381          write(*,*) 'current value iceradius(1) North:',iceradius(1)
    384382          write(*,*) 'enter new value:'
    … …  
    396394          write(*,*) ' iceradius(2) (new value):',iceradius(2)
    397395
    398         else if (modif(1:lnblnk(modif)) .eq. 'dtemisice') then
     396        else if (trim(modif) .eq. 'dtemisice') then
    399397          write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
    400398          write(*,*) 'enter new value:'
    … …  
    412410          write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
    413411
    414         else if (modif(1:lnblnk(modif)) .eq. 'tauvis') then
     412        else if (trim(modif) .eq. 'tauvis') then
    415413          write(*,*) 'current value:',tauvis
    416414          write(*,*) 'enter new value:'
    … …  
    420418          write(*,*) ' tauvis (new value):',tauvis
    421419
    422         else if (modif(1:lnblnk(modif)) .eq. 'obliquit') then
     420        else if (trim(modif) .eq. 'obliquit') then
    423421          write(*,*) 'current value:',obliquit
    424422          write(*,*) 'obliquit should be 25.19 on current Mars'
    … …  
    429427          write(*,*) ' obliquit (new value):',obliquit
    430428
    431         else if (modif(1:lnblnk(modif)) .eq. 'peri_day') then
     429        else if (trim(modif) .eq. 'peri_day') then
    432430          write(*,*) 'current value:',peri_day
    433431          write(*,*) 'peri_day should be 485 on current Mars'
    … …  
    438436          write(*,*) ' peri_day (new value):',peri_day
    439437
    440         else if (modif(1:lnblnk(modif)) .eq. 'periheli') then
     438        else if (trim(modif) .eq. 'periheli') then
    441439          write(*,*) 'current value:',periheli
    442440          write(*,*) 'perihelion should be 206.66 on current Mars'
    … …  
    447445          write(*,*) ' periheli (new value):',periheli
    448446 
    449         else if (modif(1:lnblnk(modif)) .eq. 'aphelie') then
     447        else if (trim(modif) .eq. 'aphelie') then
    450448          write(*,*) 'current value:',aphelie
    451449          write(*,*) 'aphelion should be 249.22 on current Mars'
    … …  
    456454          write(*,*) ' aphelie (new value):',aphelie
    457455 
    458         else if (modif(1:lnblnk(modif)) .eq. 'volcapa') then
     456        else if (trim(modif) .eq. 'volcapa') then
    459457          write(*,*) 'current value:',volcapa
    460458          write(*,*) 'enter new value:'
    … …  
    465463 
    466464        endif
    467       enddo ! of do while(modif(1:1).ne.'hello')
     465      enddo ! of do ! neverending loop
    468466
    469467 999  continue
Note: See TracChangeset for help on using the changeset viewer.