Index: trunk/LMDZ.MARS/libf/aeronomars/euvheat.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/euvheat.F (revision 551) +++ trunk/LMDZ.MARS/libf/aeronomars/euvheat.F (revision 552) @@ -193,5 +193,5 @@ do l=1,nlayermx - pdteuv(ig,l)=0.16*jtot(l)/10. + pdteuv(ig,l)=euveff*jtot(l)/10. & /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l))) & *(1.52/dist_sol)**2 Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F (revision 551) +++ trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F (revision 552) @@ -72,6 +72,5 @@ & 1.,17.,33.,18.,34.,16.,48./ ! minors character*3 tmp ! temporary variable - integer ierr,lnblnk - external lnblnk + integer ierr logical :: oldnames ! =.true. if old tracer naming convention (q01,...) @@ -476,11 +475,11 @@ c---------------------------------------------------------------------- open(210, iostat=ierr,file= - & datafile(1:lnblnk(datafile))//'/atmosfera_LMD_may.dat') + & trim(datafile)//'/atmosfera_LMD_may.dat') if (ierr.ne.0) then write(*,*)'Error : cannot open file atmosfera_LMD_may.dat ' write(*,*)'(in aeronomars/inichim.F)' - write(*,*)'It should be in :', datafile(1:lnblnk(datafile)),'/' - write(*,*)'1) You can change this directory address in ' - write(*,*)' file phymars/datafile.h' + write(*,*)'It should be in :', trim(datafile),'/' + write(*,*)'1) You can change this path in callphys.def with' + write(*,*)' datadir=/path/to/datafiles/' write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)' write(*,*)' can be obtained online on:' @@ -489,11 +488,11 @@ endif open(220, iostat=ierr,file= - & datafile(1:lnblnk(datafile))//'/atmosfera_LMD_min.dat') + & trim(datafile)//'/atmosfera_LMD_min.dat') if (ierr.ne.0) then write(*,*)'Error : cannot open file atmosfera_LMD_min.dat ' write(*,*)'(in aeronomars/inichim.F)' - write(*,*)'It should be in :', datafile(1:lnblnk(datafile)),'/' - write(*,*)'1) You can change this directory address in ' - write(*,*)' file phymars/datafile.h' + write(*,*)'It should be in :', trim(datafile),'/' + write(*,*)'1) You can change this path in callphys.def with' + write(*,*)' datadir=/path/to/datafiles/' write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)' write(*,*)' can be obtained online on:' Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90 (revision 551) +++ trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90 (revision 552) @@ -314,7 +314,7 @@ -! If Zmax > 2000 km there is a problem / stop - - if (Zmax .gt. 2000000.) then +! If Zmax > 4000 km there is a problem / stop + + if (Zmax .gt. 4000000.) then Print*,'Zmax too high',ig,zmax,zmin do l=1,nlayermx @@ -918,9 +918,10 @@ SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig) IMPLICIT NONE +#include "dimensions.h" INTEGER :: nl,nq INTEGER :: l,iq,ig INTEGER,dimension(nq) :: gc - REAL,DIMENSION(nl,nq) :: Q1,DQ + REAL,DIMENSION(nl,nqmx) :: Q1,DQ REAL*8,DIMENSION(nl,nq) :: Q2 REAL :: dtime @@ -951,4 +952,5 @@ SUBROUTINE MMOY(massemoy,mmol,qq,gc,nl,nq) IMPLICIT NONE +#include "dimensions.h" INTEGER :: nl,nq,l @@ -956,5 +958,5 @@ REAL*8,DIMENSION(nl,nq) :: qq REAL*8,DIMENSION(nl) :: massemoy - REAL,DIMENSION(nq) :: MMOL + REAL,DIMENSION(nqmx) :: MMOL @@ -1034,4 +1036,5 @@ & Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig) IMPLICIT NONE +#include "dimensions.h" INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig @@ -1044,5 +1047,5 @@ REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk REAL*8 :: facZ,dZ,H - REAL,DIMENSION(nq) :: mmol + REAL,DIMENSION(nqmx) :: mmol masseU=1.660538782d-27 kbolt=1.3806504d-23 @@ -1354,4 +1357,5 @@ & pp,M,gc,nl,nq,nlx) IMPLICIT NONE +#include "dimensions.h" INTEGER :: nl,nq,nlx,il,nn,iP INTEGER,DIMENSION(1) :: indP @@ -1359,5 +1363,5 @@ REAL*8,DIMENSION(nl) :: Z,P,T REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk - REAL,DIMENSION(nq) :: M + REAL,DIMENSION(nqmx) :: M REAL*8,DIMENSION(nq) :: nNew REAL*8,DIMENSION(nlx) :: pp,tt,tnew @@ -1447,4 +1451,5 @@ & ,pp,M,gc,nl,nq,nlx,facM) IMPLICIT NONE +#include "dimensions.h" INTEGER :: nl,nq,nlx,il,nn,iP INTEGER,DIMENSION(1) :: indP @@ -1452,5 +1457,5 @@ REAL*8,DIMENSION(nl) :: Z,P,T REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk - REAL,DIMENSION(nq) :: M + REAL,DIMENSION(nqmx) :: M REAL*8,DIMENSION(nq) :: nNew REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h =================================================================== --- trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 551) +++ trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 552) @@ -18,5 +18,5 @@ & ,dustbin,nqchem_min,nltemodel,nircorr - COMMON/callkeys_r/topdustref,solarcondate,semi,alphan, & + COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff, & & tke_heat_flux @@ -37,4 +37,5 @@ real alphan real solarcondate + real euveff real tke_heat_flux Index: trunk/LMDZ.MARS/libf/phymars/inifis.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 551) +++ trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 552) @@ -249,5 +249,5 @@ $ "1-> new correction", $ "(matters only if callnirco2=T)" - nircorr=0 + nircorr=1 !default value call getin("nircorr",nircorr) write(*,*) " nircorr = ",nircorr @@ -594,4 +594,10 @@ write(*,*) " solarcondate = ",solarcondate + write(*,*) "UV heating efficiency:", + & "measured values between 0.19 and 0.23 (Fox et al. 1996)", + & "lower values may be used to compensate low 15 um cooling" + euveff=0.21 !default value + call getin("euveff",euveff) + write(*,*) " euveff = ", euveff if (.not.callthermos) then Index: trunk/LMDZ.MARS/libf/phymars/nirco2abs.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/nirco2abs.F (revision 551) +++ trunk/LMDZ.MARS/libf/phymars/nirco2abs.F (revision 552) @@ -50,14 +50,20 @@ #include "comdiurn.h" #include "nirdata.h" +#include "tracer.h" c----------------------------------------------------------------------- c Input/Output c ------------ - INTEGER ngrid,nlayer - - REAL pplay(ngrid,nlayer) - REAL dist_sol,mu0(ngridmx),fract(ngridmx),declin - - REAL pdtnirco2(ngrid,nlayer) + integer,intent(in) :: ngrid ! number of (horizontal) grid points + integer,intent(in) :: nlayer ! number of atmospheric layers + real,intent(in) :: pplay(ngrid,nlayer) ! Pressure + real,intent(in) :: dist_sol ! Sun-Mars distance (in AU) + integer,intent(in) :: nq ! number of tracers + real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers + real,intent(in) :: mu0(ngridmx) ! solar angle + real,intent(in) :: fract(ngridmx) ! day fraction of the time interval + real,intent(in) :: declin ! latitude of sub-solar point + + real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s) c c Local variables : @@ -73,5 +79,7 @@ c local saved variables c --------------------- - + logical,save :: firstcall=.true. + real,save :: ico2=0 ! index of "co2" tracer + real,save :: io=0 ! index of "o" tracer c p0noonlte is a pressure below which non LTE effects are significant. c REAL p0nonlte @@ -90,6 +98,4 @@ real p2011,cociente1,merge real cor0,oco2gcm - integer nq - real pq(ngrid,nlayer,nq) c---------------------------------------------------------------------- @@ -97,4 +103,24 @@ c Initialisation c -------------- + if (firstcall) then + if (nircorr.eq.1) then + ! we will need co2 and o tracers + ico2=igcm_co2 + if (ico2==0) then + write(*,*) "nirco2abs error: I need a CO2 tracer" + write(*,*) " when running with nircorr==1" + stop + endif + io=igcm_o + if (io==0) then + write(*,*) "nirco2abs error: I need an O tracer" + write(*,*) " when running with nircorr==1" + stop + endif + endif + firstcall=.false. + endif + + c co2heat is the heating by CO2 at 700Pa for a zero zenithal angle. co2heat0=n_a*(1.52/dist_sol)**2/daysec @@ -121,6 +147,6 @@ if(nircorr.eq.1) then cor0=1./(1.+n_p0/pplay(ig,l))**n_b - if(pq(ig,l,1).gt.1.e-6) then - oco2gcm=pq(ig,l,3)/pq(ig,l,1) + if(pq(ig,l,ico2).gt.1.e-6) then + oco2gcm=pq(ig,l,io)/pq(ig,l,ico2) else oco2gcm=1.e6 @@ -179,5 +205,5 @@ if(nircorr.eq.1) then cor0=1./(1.+n_p0/pplay(ig,l))**n_b - oco2gcm=pq(ig,l,3)/pq(ig,l,1) + oco2gcm=pq(ig,l,io)/pq(ig,l,ico2) cociente1=oco2gcm/oldoco2(l) merge=alog10(cociente1)*alfa2(l)+alog10(cor0)* Index: trunk/LMDZ.MARS/libf/phymars/tabfi.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/tabfi.F (revision 551) +++ trunk/LMDZ.MARS/libf/phymars/tabfi.F (revision 552) @@ -81,6 +81,4 @@ c ----------------------- INTEGER setname, cluvdb, getdat - - INTEGER lnblnk c----------------------------------------------------------------------- @@ -296,5 +294,5 @@ - do while(modif(1:1).ne.'hello') + do ! neverending loop write(*,*) write(*,*) @@ -306,7 +304,7 @@ write(*,*) - write(*,*) modif(1:lnblnk(modif)) , ' : ' - - if (modif(1:lnblnk(modif)) .eq. 'day_ini') then + write(*,*) trim(modif) , ' : ' + + if (trim(modif) .eq. 'day_ini') then write(*,*) 'current value:',day_ini write(*,*) 'enter new value:' @@ -316,5 +314,5 @@ write(*,*) 'day_ini (new value):',day_ini - else if (modif(1:lnblnk(modif)) .eq. 'z0') then + else if (trim(modif) .eq. 'z0') then write(*,*) 'current value (m):',z0_default write(*,*) 'enter new value (m):' @@ -324,5 +322,5 @@ write(*,*) ' z0 (new value):',z0_default - else if (modif(1:lnblnk(modif)) .eq. 'emin_turb') then + else if (trim(modif) .eq. 'emin_turb') then write(*,*) 'current value:',emin_turb write(*,*) 'enter new value:' @@ -332,5 +330,5 @@ write(*,*) ' emin_turb (new value):',emin_turb - else if (modif(1:lnblnk(modif)) .eq. 'lmixmin') then + else if (trim(modif) .eq. 'lmixmin') then write(*,*) 'current value:',lmixmin write(*,*) 'enter new value:' @@ -340,5 +338,5 @@ write(*,*) ' lmixmin (new value):',lmixmin - else if (modif(1:lnblnk(modif)) .eq. 'emissiv') then + else if (trim(modif) .eq. 'emissiv') then write(*,*) 'current value:',emissiv write(*,*) 'enter new value:' @@ -348,5 +346,5 @@ write(*,*) ' emissiv (new value):',emissiv - else if (modif(1:lnblnk(modif)) .eq. 'emisice') then + else if (trim(modif) .eq. 'emisice') then write(*,*) 'current value emisice(1) North:',emisice(1) write(*,*) 'enter new value:' @@ -364,5 +362,5 @@ write(*,*) ' emisice(2) (new value):',emisice(2) - else if (modif(1:lnblnk(modif)) .eq. 'albedice') then + else if (trim(modif) .eq. 'albedice') then write(*,*) 'current value albedice(1) North:',albedice(1) write(*,*) 'enter new value:' @@ -380,5 +378,5 @@ write(*,*) ' albedice(2) (new value):',albedice(2) - else if (modif(1:lnblnk(modif)) .eq. 'iceradius') then + else if (trim(modif) .eq. 'iceradius') then write(*,*) 'current value iceradius(1) North:',iceradius(1) write(*,*) 'enter new value:' @@ -396,5 +394,5 @@ write(*,*) ' iceradius(2) (new value):',iceradius(2) - else if (modif(1:lnblnk(modif)) .eq. 'dtemisice') then + else if (trim(modif) .eq. 'dtemisice') then write(*,*) 'current value dtemisice(1) North:',dtemisice(1) write(*,*) 'enter new value:' @@ -412,5 +410,5 @@ write(*,*) ' dtemisice(2) (new value):',dtemisice(2) - else if (modif(1:lnblnk(modif)) .eq. 'tauvis') then + else if (trim(modif) .eq. 'tauvis') then write(*,*) 'current value:',tauvis write(*,*) 'enter new value:' @@ -420,5 +418,5 @@ write(*,*) ' tauvis (new value):',tauvis - else if (modif(1:lnblnk(modif)) .eq. 'obliquit') then + else if (trim(modif) .eq. 'obliquit') then write(*,*) 'current value:',obliquit write(*,*) 'obliquit should be 25.19 on current Mars' @@ -429,5 +427,5 @@ write(*,*) ' obliquit (new value):',obliquit - else if (modif(1:lnblnk(modif)) .eq. 'peri_day') then + else if (trim(modif) .eq. 'peri_day') then write(*,*) 'current value:',peri_day write(*,*) 'peri_day should be 485 on current Mars' @@ -438,5 +436,5 @@ write(*,*) ' peri_day (new value):',peri_day - else if (modif(1:lnblnk(modif)) .eq. 'periheli') then + else if (trim(modif) .eq. 'periheli') then write(*,*) 'current value:',periheli write(*,*) 'perihelion should be 206.66 on current Mars' @@ -447,5 +445,5 @@ write(*,*) ' periheli (new value):',periheli - else if (modif(1:lnblnk(modif)) .eq. 'aphelie') then + else if (trim(modif) .eq. 'aphelie') then write(*,*) 'current value:',aphelie write(*,*) 'aphelion should be 249.22 on current Mars' @@ -456,5 +454,5 @@ write(*,*) ' aphelie (new value):',aphelie - else if (modif(1:lnblnk(modif)) .eq. 'volcapa') then + else if (trim(modif) .eq. 'volcapa') then write(*,*) 'current value:',volcapa write(*,*) 'enter new value:' @@ -465,5 +463,5 @@ endif - enddo ! of do while(modif(1:1).ne.'hello') + enddo ! of do ! neverending loop 999 continue