== 17/09/08 == >>> Build a version with new soil but old radiative transfer, but keeping possibility of switching back to new radiative transfer), which incorporates changes & improvements currently included in the 'reference version' GCM (see /u/emlmd/LMDZ.MARS.mixdyn) >>> start by modifying makegcm as in /u/emlmd/LMDZ.MARS.mixdyn, so that it runs without environment variables and set LIBOGCM to /tmp15/emlmd/libo >>> directory contents of 'aeronomars', 'grid' and 'filtrez' are simillar to those in /u/emlmd/LMDZ.MARS.mixdyn >>> in bibio , only file mxva.F needed be upgraded >>> get phymars and dyn3d contents from /u/emlmd/LMDZ.MARS.mixdyn (and remove all *old files) >>> check differences between dyn3d and /u/emlmd/LMDZ.MARS.170908/libf/dyn3d and upgrade when necessary: - removed 'netcdf.inc' file (has nothing to do there!) - comgeom.h and comgeom.h : made fortran90 compliant - control.h : made fortran90 compliant - dynredem.F : more read/write controls + comments in english - ini_archive.F : new soil/thermal inertia changes - integrd.F : added additional information to output when crashing - lect_start_archive.F : new soil/thermal inertia changes - newstart.F : new soil/thermal inertia changes + comments in english - start2archive.F : new soil/thermal inertia changes - vanleer.F : removed inapropriate 'external' statement - write_archive.F : enable writting a subterranean field >>> check differences between phymars and /u/emlmd/LMDZ.MARS.170908/libf/phymars and upgrade when necessary: New soil stuff: - added comsoil.h - iniwrite.F : new soil changes - added iniwritesoil.F90 and writediagsoil.F90 for subterranean fields - added interp_line.F (for subterranean grid interpolation) - adapted phyetat0.F for subterranean temperature & inertia - adapted physdem1.F to include new soil stuff - physiq.F : added calls to writediagsoil - soil.F : new routine (fixed vertical grid + variable thermal inertia) - added soil_settings.F (to read/initialize/interpolate soil properties) - updated surfdat.h (since thermal inertia is now in comsoil.h) - updated tabfi.F : include new soil properties - updated testphys1d.F - updated dimphys.h (set nsoil=18 as default) == 18/09/08 == >>> add the possibility of easily switching to Tran radiative transfert - updated aerdust.h.ocke97 (changed some variables name) so it can replace aerdust.h (which is currently the same as aerdust.h.clan91). - imported Tran's 'gfluxv.F' routine - imported Trans' version of 'swr.F' routine, saved it as 'swr.F.toon' >>> Backup of 'old' Morcrette swr.F is 'swr.F.morc' NB: to switch from one radiative transfer to the other, just copy swr.F.morc or swr.F.toon to swr.F (and eventually 'touch swr.F' so that makegcm recompiles swr.F) No other dependencies (swr.F.toon uses 'gfluxv.F' and swr.F.morc uses 'dedd.F'). >>> Changed the latter, so that users can switch from one to the other - modified swr.F.toon to become swr_toon.F (and to include gfluxv.F) - modified swr.F.morc to become swr_fouquart.F (and to include dedd.F) - added a flag in callkeys.h, swrtype (parameter to be set/changed by the user 1=Fouquart and 2=Toon) - update readtesassim so that the coefficient by which opacity is multiplied is set according to the 'swrtype' parameter == 25/09/08 == >>> Implement the use of tracer-by-name in physics - in phymars/tracer.h set 'noms' length to 20 (instead of 10) - in phymars/callsedim2q.F and phymars/callsedim.F, use tracers by name - in phymars/dustopacity.F, use tracers by name - in phymars/vdifc.F, use tracers by name == 26/09/08 == >>> Change implementation strategy (for now); don't move surface tracer around i.e.: surface ice remains equivalent to qsurf(nqmx)=qsurf(i_h2o_vap) and likewise for surface tendencies ... - modified vdifc.F and callsedim.F back - modified initracer.F (so that water names are h2o_vap & h2o_ice) == 29/09/08 == - modified aeronomars/init_chimie_B (cosmetics) - corrected aeronomars/moldiff.F internal routine tridag; changed "pause" error messages to 'stop' messages - modified phymars/watercloud.F to use tracers by name - corrected aeronomars/molvis.F (undefined 'fac' and 'Akk' written to output at first call) == 30/09/08 == - modified aeronomars/calchim.F to use tracers by name - adapted aeronomars/photochemist_B.F to use tracers by name - adapted aeronomars/chemtermos.F to use tracers by name - adapted aeronomars/concentrations.F to use tracers by name - corrected aeronomars/conduction.F (undefined 'Akk' written to output at first call) - adapted aeronomars/euvheat.F to use tracers by name - adapted aeronomars/moldiff.F and moldiffcoeff.F to use tracers by name == 01/10/08 == - For more compatibility with LMDZ4; mimic reading a 'traceur.def' file in the dynamics via a call to a routine 'iniadvtrac.F' and saving tracers names in 'advtrac.h' -> created 'iniadvtrac.F', 'advtrac.h' and modified gcm.F - modified 'dynetat0.F' so that tracers are loaded from 'start.nc' by name - modified 'dynredem.F' so that tracers are written to 'restart.nc' by name - modified 'initracer.F' tu use tracers by name == 02/10/08 == - removed use of 'nqchem_min' everywhere: adapted 'euvheat.F','inifis.F','physiq.F' (leave 'inichim.F' for later) - updated 'phyetat0' and 'physdem1.F' to read/write surface tracers by name - modify things so that surface water ice index is the same as atmospheric water ice (except when running without water ice; then simply set i_h2o_ice=i_h2o_vap). NB: the easiest is to have global storage of tracer names/indexes in tracer.h => changed initracer.F & tracer.h to have global igcm_something indexes == 03/10/08 == - adaptations for surface ice index, modified files: phyetat0.F : if there is a dynamical tracer 'h2o_vap' then load surface tracer called 'h2o_ice' instead initracer.F : in 'old' tracer name case: move qsurf(nqmx)->qsurf(nqmx-1) and set i_h2o_ice=i_h2o_vap if iceparty=.false. physdem1.F : if old tracer names: move qsurf(nqmx-1)->qsurf(nqmx) if iceparty=.false., write surface tracer 'h2o_ice' (and not 'h2o_vap') to file. adapted vdifc.F, callsedim.F & watercloud.F & physiq.F so that surface ice is now identified as qsurf(i_h2o_ice) - updated aeronomars/perosat.F (cosmetics) == 06/10/08 == - modify newstart.F and lect_start_archive.F to use tracers by name == 07/10/08 == - adapted inichim_newstart.F (added qsurf to arguments) and inichim_readcallphys.F == 08/10/08 == - implement reading traceur.def in dyn3d/iniadvtrac.F == 16/10/08== -small change in inifis.F (only warn if too many tracers, compared to the expected number, not stop). - corrected bug in initracer.F == 21/10/08 == - modified newstart.F to load B.Diez subsurface ice maps. - corrected small bug (uninitialized variable) in interp_horiz.F == 22/10/08 == - updated iniwritediagsoil.F so that thermal inertia is written to diagsoil.nc == 31/10/08 == - changed xvik.F program so it works even if we don't have atmospheric temperature at hand (then it uses a 10km reference scale height) and so that it does surface pressure interpolation log-wise. == 03/11/08 == - modified physiq to compute (and output) co2 column. - added improvement by Francois in newcondens.F about computing CO2 partial pressure. This behavior is turned on by setting internal logical flag 'improved_ztcond' to '.true.' (and running with a co2 tracer) - updated 'start2archive' to work with 'new' gcm output (soil, tracers ...) == 04/11/08 == - upgraded xvik program to look for temperature in 7th layer variable if there is no global atmospheric temperature field at hand. == 05/11/08 == - more modifs to newcondens.F: added another internal flag 'bound_qco2' to enforce (if set to .true.) that co2 mass mixing ratio remains bounded. == 07/11/08 == - corrected 'writediagfi' & 'writediagsoil' so that an error message is issued if called with a variable name which is too long. == 18/12/08 == - corrected bug in dyn3d 'addfi.F', (dimensions of local array p()) == 23/02/09 == - modified "aeronomars/param_read.F" to do strictly fortran data initialization (otherwise xlf compiler complains) - changed a few '1.e-30' to '1.d-30' in aeronomars/photochemist_B.F so that max functions has 2 doubles as arguments (otherwise xlf compiler complains) ==07/04/09 == -cosmetic changes/minor improvements in the handling of tracers in: aeronomars/photochemist_B.F aeronomars/perosat.F aeronomoars/euvheat.F aeronomars/moldiffcoeff.F aeronomars/moldiff.F aeronomars/cocentrations.F aeronomars/chemtermos.F aeronomars/calchim.F --> NB: still there are differences in outputs when order of tracers is changed == 09/04/09 == >>> fixed problem in 'vdifc.F' which lead to different results when moving tracers around. == 10/04/09 == >>> corrected small bug in diagnostic outputs of 'watercloud.F' (tendencies were not added to tracer values). == 21/04/09 == >>> corrected small bug in "physdem1.F" about writing water ice surface tracer to file == 07/05/09 == >>> very minor correction (firstcall not set to true after first call if no tracers) in convadj.F == 30/06/09 == >>> Implement reading *def files with IOIPSL ersatz 'getin' function - import "ioipsl_errioipsl.F90","ioipsl_getincom.F90","ioipsl_stringop.F90" in bibio - adapted 'dyn3d/defrun_new.F' to use "getin" function - adapted 'phymars/inifis.F' to use "getin" function == 01/07/09 == >>> Adapted 'create_make_gcm' so that the "use" in *.F files is identified and corresponding dependencies included in the makefile rules. >>> Added the 3D scattering from aerosols by JB Madeleine: - minor changes in aerave.F - added the calls to aeropacity.F, and aeroptproperties.F in callradite.F - changed the calls to lwu.F and swr.F in lwmain.F and swmain.F, respectively - added 3D scattering properties in lwu.F and swr.F - added the new aeroptproperties.F, aeropacity.F and suaer.F90 routines (removed dustopacity.F) - updated aeropacity.F with new tracer names - changed the call to callradite.F in physiq.F, added the initialization of reffrad and nueffrad (aerosol effective radius and variance) - removed all the lines relative to the old "activice" option, including temperature variation due to latent heat release (now in comments) - renamed nsize into naersize in watercloud.F, watersat.F and newsedim.F, to avoid conflicts with another "nsize" variable in the radiative transfer - added the statement of nuice in watercloud.F, which is the effective variance of the log-normal distribution for ice - updated yomaer.h and removed aerice.h (and corresponding "includes") == 02/07/09 == >>> Adapted 'aeronomars/inichim_readcallphys.F' (called by newstart) to use "getin" routine. + minor correction in 'inifis.F' (close 'iradia.def' file) >>> Minor correction in 'dyn3d/dynetat0.F' and 'phymars/phyetat0.F'; do not attempt to reindex tracers if none were found. >>> in 'deftank' added examples of 'traceur.def' files (traceur.def.co2 : 1 co2 tracer; traceur.def.watercycle : 2 traceurs, water vapour and water ice tracer.def.chemistry : all 15 species) == 06/07/09 == >>> Modified 'makegcm' and makegcm_g95' so that modules files are put with libraries (and not in current directory) == 21/07/09 == >>> Added in "testphys1d.F" a check that the (required) 'run.def' file is around (that file should contain the "INCLUDEDEF=callphys.def" instruction otherwise getin() calls won't work. Also added reding of 'traceur.def' (or initialisation of tracer names to dummy values q01,q02, ...) in testphys1d.F == 22-24/07/09 == JBM >>> Removed "iceparty" everywhere (calchim.F, inichim_readcallphys.F, chimie_data.h, inichim_newstart.F, photochemist_B.F, callkeys.h, callradite.F, callsedim.F, aeropacity.F, inifis.F, initracer.F, physdem1.F, physiq.F, watercloud.F). Water = .true. now implies the use of two tracers, i.e. water vapor and water ice. >>> Removed fisice.h and stated the corresponding variables in the right places (physiq.F, callsedim.F, watercloud.F, newcondens.F, calchim.F); also removed unused cltop and clsurf variables in physiq.F. >>> Removed the variable's splitting, which is now obsolete, in callradite.F (and its subroutines lwi.F, lwxb.F, lwxn.F, lwflux.F, lwmain.F, lwxd.F). Also removed the variable's splitting in calldrag_noro.F. Finally removed ndomain from dimradmars.h. >>> Removed useless tests in aeropacity.F. == 27-30/07/09 == JBM >>> Cleaned the WRITEDIAGFI section in physiq.F, and moved the "mtot", "icetot" and "tauTES" variables from watercloud.F to physiq.F. Also cleaned the albedo change due to water ice deposition. >>> Renamed rdust into rnuclei in callsedim.F, physiq.F and watercloud.F. >>> Added a logical test for (water.and..not.tracer) in inifis.F. >>> Removed "qsurf","zls" and "icount" from the list of inputs in watercloud.F (these variables were not used by the subroutine). >>> Added a call to watercloud.F at firstcall, using typical dust optical depth (taufirstcall), in order to correctly initialize the ice particle size distribution for the radiative transfer scheme. Also created a new subroutine to load the effective radius and variance of the aerosols used by the radiative transfer scheme. Its name is updatereffrad.F, and it is called before aeroptproperties.F in callradite.F. == 05/08/09 == JBM >>> TES water-ice opacity is now fully computed using the radiatively active aerosol scheme. Absorption coefficient is calculated using Qext and omega at the IR reference wavelength. Omega was not computed before; it was only computed in the GCM channels, not at the reference wavelength. Thus it has been added in suaer.F90, aerave.F, yomaer.h, callradite.F, aeropacity.F and aeroptproperties.F. If "activice" is false, the TES opacity is computed using the old method (fit of the Qabs as a function of reff curve). == 07/08/09 == EM >>> Removed "iceparty" option from callphys.def >>> modified physiq.F and initracer.F so that building of array niqchem() which contains the indexes of all chemistry tracers + water vapour and water ice is done in initracer.F (array niqchem(:) is now a common in tracer.h ). Also had to adapt inichim_newstart.F to behave similarly. == 24/08/09 == EM included corrections by JBM >>> Added declaration of nqchem(nqmx) local array in aeronomars/inichim_newstart.F >>> Encapsulated calls to writediagfi & wstats in if (ngrid.ne.1) clause in phymars/aeropacity.F (otherwise it crashes in 1D). == 26/08/09 == EM >>> modified tracer.h, initracer.F, inichim_newstart.F and physiq.F to not use an niqchem() array (added in 07/08/09 changes) >>> modified phymars/readtesassim.F90 and phymars/aeropacity.F so that assimilated dust for MY24 or MY25 or MY26 may be used (with iaervar= 24,25 or 26); we keep iaervar=4 to also read MY24 dust for compatibility with older versions of code. Modified deftank/callphys.def : added comments about new iaervar values == 25/09/09 == EM >>> modified phymars/testphys1D.F : added incrementation of tracer values after call to physiq(). == 27/11/09 == EM >>> updated comments in makegcm (translated help to english) >> shifted to reading file traceur.def (dyn3d/iniadvtrac.F) in an Earth-LMDZ4 like way: first line == number of tracers and then tracer name (1 per line; later we'll add advection scheme type) >> also updated example 'traceur.def' files in deftank accordingly == 02/12/09 == EM >>> upgraded testphys1d.F and profile.F to load run parameters with getin() function (from run.def ; no need for a "testphys1d.def" any more) added an example 'run.def.1d' file to 'deftank' ==10/12/09 == EM >>> minor correction in testphys1d.F (was still checking if there is a testphys1d.def file around ; which is not used anymore) ==15/01/10 == JBM >>> aeropacity.F: implemented a weighting of the dust opacity profile by using the dust size vertical profile defined in updatereffrad.F. >>> aeroptproperties.F: changed the integration scheme (Gauss-Legendre) of the scattering parameters. >>> suaer.F90: removed the use of an ad-hoc "solsir" factor. It is now directly computed from the scattering properties read in the ASCII files. Consequently, the IR extinction coefficient has been divided by solsir=2 in the ASCII file (called optprop_dustir_x0.5.dat instead of optprop_dustir.dat to allow compatibility with the older version and avoid chaos). >>> updatereffrad.F: changed the dust size vertical profile. >>> yomaer.h: the particle radius variable is now in simple precision, because the new scattering property integration scheme has changed. >>> aeronomars/inichim_readcallphys.F (small) bug correction load value of 'water' before testing its value... == 18/01/10 == EM >> added possibility to read (in inifis.F and aeropacity.F) the value of dust opacity tauvis from callphys.def file == 01/02/10 == EM >> added JBM updates of "callradite.F" (coments) "aeroptproperties.F" (bug fix of bad array bounds) and "aeropacity.F" (encapsulate wstats calls in if (callstats) ) >> added implementation of TES Cap albedos: "albedocaps.F90" and adapted "newcondens.F" (and physiq.F, because added 'zls' argument to newcondens) and surfdat.h >> changed default settings for dust : set a 1.3 factor in readtesassim.F90 when using Toon radiative transfer, use M.Wolff-T-Matrix files in suaer.F90 == 03/02/10 == EM >> Updated newstart.F in dyn3d, so that sub-surface thermal inertia values may be different in North and South hemispheres. >> Updated "makegcm" and "makegcm_g95" scripts (cosmetic + default compilation option changes) >> Minor changes in aeronomars/init_chimie_B.F (do not use lnblnk(); F90 trim() intrinsic is much safer and better), and in initracer.F (better control over a possible array bound underflow). Also, in dyn3d/iniadvtrac.F, close input file properly, and in infis.F, more verbose message to output. == 26/02/10 == EM >> Updated makegcm and makegcm_g95 : default usage is now to set everything ("environment variables") in the script. Changed some default compilation options. >>> removed 'float()' instructions in tabfi.F and iniwrite.F use "real()" to be compliant with standards. >> Corrected small bug in testphys1d.F (look for file traceur.def) also added initialisation of tracers >> Cleaned up inifis.F and initracer.F (some sanity checks were obsolete and/or wrong) >> Improved writediagfi.F so that 1D (individual column in the GCM, or fields in testphys1d) data can be written in the diagfi.nc file >> Minor changes/improvements in calls to writediagfi from physiq.F for dust == 08/03/10 == EM >> Minor update of "makegcm" and "makegcm_g95": use instruction "./makdim" (instead of "makdim"; in case "." is not in user's path) >> put "real()" instructions instead of "float()" in dyn3d routines: disvert.F , dynredem.F , fluxstoke.F , fxhyp.F , fxy.F , fxysinus.F , fyhyp.F , gcm.F , grid_atob.F , grid_noro.F , grid_noro1.F , ini_archive.F , inigeom.F , newstart.F , ran1.F , sortvarc.F , sortvarc0.F >> Minor update of aeropacity.F (added if (callstats) around call to wstats) == 28/04/10 == EM >> Put the splitting in radiative transfer back in the model (JB): updated calldrag_noro.F callradite.F dimradmars.h lwflux.F lwi.F lwmain.F lwxb.F lwxd.F lwxn.F swmain.F >> Fix bug (AS) in callradite.F (wrong loop boundaries line 332) >> Fix bug (AS+JB) in "swr_toon.F" to enable running with more than 100 levels... >> Fix bug (JBM) in callsedim2q.F about setting pdqs_sed(:niq(iq)) to zero == 25/08/10 == EM >> Add a 'makegcm_gfortran' for compiling with gfortran and a 'makegcm_ifort' for compiling with ifort (on Gnome) == 03/09/10 == EM >> Modifications to enable running in double precision (using starts in r4 or r8); just add options '-r8 -DNC_DOUBLE' to compile GCM in double precision -> adapted dyn3d/dynetat0.F, physmars/physdem1.F, phymars/soil_settings.F, phymars/readtesassim.F90, phymars/writediagfi.F, phymars/def_var.F90, phymars/writediagsoil.F90, phymars/wstats.F90, phymars/inistats.F >> Added dyn3d/writediagdyn.F90 routine (to output scalar dynamical fields), adapted 'comconst.h' and 'comvert.h' to be Fortran77/Fortran 90 compatible. == 14/12/10 == EM >> Add -f option to #!/bin/csh in makegcm* scripts (to make sure that it is the bash environment compiler that is used as a default) >> Update convadj.F with RW's version (fixes bug of non conservation of tracers in cases where convection stops at one level and starts at the next). == 13/12/10 == EM >> Update testphys1d.F so that initial tracer profiles may be loaded at initialization == 24/01/11 == JBM(+ some cleanup by EM) >> Reactivated the "doubleq" method (two-moment scheme for dust transport) and connected it with the radiative transfer code. The opacity is set constant below a level indicated by the variable cstdustlevel in aeropacity.F to remove the thick layer of dust near the surface created by the constant lifting rate. The "density scaled opacity" used by the MCS team is computed and saved in dsodust. Updated routines: callradite.F, aeropacity.F, updatereffrad.F, callsedim.F, newsedim.F, initracer.F, vdifc.F, suaer.F90. >> Added the use of named scatterers (same method as tracer-by-name) in the radiative transfer code. Scatterers are declared at firstcall in callradite.F (which is the equivalent of traceur.def) and the corresponding indices are saved in the common called aerkind.h (the equivalent of tracer.h). Updated routines: callradite.F, updatereffrad.F, aeropacity.F, suaer.F90, aerkind.h. EM: aerdust.h, aerdata.h and aerice.h are not used any more >> Merged callsedim2q.F and callsedim.F in one single routine callsedim.F to allow doubleq to be used with other tracers. >> Added an input parameter called beta to newsedim.F, that allows to account for the shape of the particles in the computation of the sedimentation velocity. >> Added the ability to transport a radiatively active population of submicron dust particles (flag "submicron" in callphys.def). Updated routines: callradite.F, updatereffrad.F, aeropacity.F, initracer.F, vdifc.F, suaer.F90, inifis.F, callkeys.h. >> Connected the predicted size and amount of dust to the water cycle. The size of the particles is now given by rdust in updatereffrad.F, and the amount of cloud condensation nuclei (CCN) is given by ccn in aeropacity.F. The calculation of ccn is done in aeropacity.F because it's deduced from the opacity when doubleq is not used. rnuclei and dustcores are removed from watercloud.F and replaced by rdust and ccn. Updated routines: physiq.F, callradite.F, updatereffrad.F, aeropacity.F, watercloud.F. >> Removed the "fake" call to watercloud.F in physiq.F which was used to give the size of the ice particles to the radiative transfer code at firstcall. Instead, rice is computed in updatereffrad.F using a simple equation and a typical amount of dust nuclei (ccn0). >> Increased the number of Gauss integration points in aeroptproperties.F. >> Added the ability to write the 3D scattering parameters of a given aerosol in the outputs using the out_qwg flag in aeroptproperties. >> Changed "DO iir=1,4" into "DO iir=1,nir-1" in suaer.F90 (in case the number of infrared channels is changed). >> Added nuice_ref in tracer.h and initracer.F, which is the effective variance of the log-normal distribution used for water-ice particles in the radiative transfer code. >> Updated the computation of rice, reffice and rsedice in updatereffrad.F and callsedim.F. >> Added nuice_sed in callsedim.F, which is the effective variance of the lognormal distribution used for the sedimentation of water-ice particles. >> Added ccn_factor in watercloud.F, which is the ratio of the total number of dust particles over the number of condensation nuclei. >> The variable beta is not saved anymore in newsedim.F. >> Corrected the BIG bug in lwu.F that was responsible for unstabilities when clouds were radiatively active (FF+EM+JBM!) >> Turned ilwd, ilwn and ilwb to 1 in inifis.F. >> Added dust and ice visible opacities in the outputs. Modified routine: aeropacity.F. >> Named water cycle tuneable parameters (nuice_sed, nuice_ref, alb_surfice, ccn_factor) which are mentioned at firstcall by the GCM (flag "water_param"). Modified routines: callsedim.F, physiq.F, watercloud.F. == 25/01/2011 == EM >> updated testphys1d.F: removed #include "aerdust.h" >> cleanup in suaer.F90: no more calls to zerophys; added test to check there is no overflow of isize index (in visible domain averaged properties case) >> minor cleanup in newsedim.F (make it more F90 oriented) >> added flag 'TESicealbedo' (set to .true. to impose polar surface albedos as observed by TES) in callphys.def == 28/01/2011 == EM >> added additional tests (to check correct reading of input files) in suaer.F90 >> updated inifis.F so that the directory where external data files are to be found (e.g. TES opacities, dust properties, etc.) can be specified in run.def (or callphys.def) as " datadir = /path/to/the/directory " == 29/01/2011 == AS >> added updated mesoscale-related routines in phymars: ---------------------------------------------------------------------------- NAME CHANGES compared to GCM counterpart ---------------------------------------------------------------------------- meso_callkeys.h --> one variable is added [consider merging w/ GCM?] meso_dustlift.F --> stress + alpha default, or read in a file stress.def if here [consider merging w/ GCM?] meso_newcondens.F --> correction on U V T tendencies is switched off (unstable in mesoscale) meso_physiq.F --> major modifications mainly related to I/O meso_slope.h --> additional common for slope scheme meso_dimphys.h_ref --> reference common serving as a basis for a compilation script (makemeso) meso_inifis.F --> major modifications mainly related to I/O meso_param_slope.F90 --> slope scheme by Spiga and Forget GRL 2008 [consider adding to GCM?] meso_readtesassim.F90 --> an old version because the new F90-compliant version needs the new makegcm scripts [TBD] meso_testphys1d.F --> similar to GCM except for routine names ---------------------------------------------------------------------------- NB: in meso_dustlift.F and meso_readtesassim.F90, the subroutines have the same name as in the GCM. this is because those files are supposed to be copied in specific temporary folders for compilation >> any future change in the following GCM routines in phymars: - callkeys.h - dustlift.F - newcondens.F - physiq.F - inifis.F - readtesassim.F90 - testphys1d.F will be in need to be impacted to the corresponding meso_ routines [hence it is important to document this README file] >> any change in any other GCM routines than the ones listed will have an effect in mesoscale simulations as well: -- the two models are being kept updated at the same time :) -- the two models would be possibly broken at the same time :( == 15/02/2011 == EM >> updated dissipation coefficients in indissip.F == 25/02/2011 == EM >> corrected bug in 'inifis.F' and 'datafile.h' to really be able to specify (in callphys.def; using datadir=/whatever/path/to/use ) the path to external datafiles (topography, surface properties, etc.) == 28/02/2011 == JBM + AS >> used settings reached by JBM to obtain his PhD results alb_surfice = 0.45 --- in physiq.F and meso_physiq.F ccn_factor = 4.5 --- in watercloud.F nuice_sed = 0.45 --- in callsedim.F >> NB: this is supposed to be further refined in the future == 01/03/2011 == AS + JBM >> nasty bug in the water cycle when iradia != 1 [no problem when iradia = 1] --> mesoscale runs w/ water cycle had strange 5-hour fluctuations in RICE, from 80mic to 5mic >> PB: calculation of ccn [condensation nuclei] is done in callradite.F * ccn must be saved --> corrected in physiq.F and meso_physiq.F * ccn must not be modified elsewhere [e.g. in watercloud, when divided by ccn_factor] --> all calculations on ccn are now moved in callradite == 03/03/2011 == AS >> Added a pre-compilation flag MESOSCALE so that the LMDZ.MARS GCM will compile without stating errors because of mesoscale routines. [meso_physiq.F, meso_inifis.F] >> Now, this MESOSCALE precompilation flag can be used to lower the number of meso_* routines when adaptations for mesoscale applications are not very extended. --> meso_testphys1d.F, meso_testphys1d.F, meso_dustlift.F routines were deleted and changes are now moved under the MESOSCALE flag in the original GCM routines --> Completely transparent for GCM compilation since it is devoid of the -DMESOSCALE option --> Very good for syncing because changes in dustlift, newcondens will be directly available in the mesoscale model == 04/03/2011 == AS + JBM >> new version version for aeroptproperties.F in phymars to limit uncertainties and be able to play with ngau >> this was coded by JBM in his personal reference version but not transmitted to the team reference version == 10/05/2011 == AS + JF >> in newsedim.F used for mesoscale computations, spurious values close to the surface --> this was related to 1 - exp(-x) calculated as zero in w(ig,l) if x is very small --> fix: when this happens, replace exp(-x) by 1 - x since x ~ 0 >> in newsedim.F, "if (dztop.gt.epaisseur(ig,l)) then" was closed too soon by an "endif" --> hence basically the simple method was never used and useless calculations with the complex method were carried out --> fix by moving the closing "endif" [in addition to corrections mentioned in the previous point] == 17/05/2011 == EM >> set internal computations using double precision in growthrate.F and watercloud.F (otherwise we sometimes end up with Nans). >> add extra checks in newcondens.F to avoid possibility of out of bounds evaluation of array masse() == 25/05/2011 == AS >> found that the 10/05/2011 bug fix in newsedim.F is also useful for GCM runs. >> no more need to modify callradite.F prior to compilation [but still dimradmars.h must be modified] --> in callradite.F we now have -- DEFAULT name_iaer(1) is "dust_conrath" -- IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" -- IF (water.AND.activice) name_iaer(2) = "h2o_ice" == 27/05/2011 == EM >> minor bug correction in writeg1d (JYC) >> add "-check" to debug option in makegcm_ifort == 08/06/2011 == EM >> minor bug fix in lect_start_archive.F (using wrong surface temp. array). + swiched output messages to english and added that tracers not found in file must be initialized by user. >> minor bug fix in datareadnc.F : 'datafile' path must be initialized. == 08/06/2011 == EM >> Significant update on how the number of scatterers is managed: Instead of having to manualy change 'nearkind' in dimradmars.h, the number of scatterers must now be set when compiling, using makegcm "makegcm -s 1" for one scatterer (dust) or "makegcm -s 2" (e.g. dust and water ice), default behaviour (ie not specifying -s #) is -s 1 Modified phymars/dimradmars.h , added directory phymars/scatterers with script make_scatterers , and adapted makegcm* scripts. == 17/06/2011 == AC ================================================ ======== IMPLEMENTATION OF THERMALS ============ ================================================ The main goal of this revision is to start including the thermals into the model for development purposes. Users should not use the thermals yet, as several major configuration changes still need to be done. This version includes : - updraft and downdraft parametrizations - velocity in the thermal, including drag - plume height analysis - closure equation - updraft transport of heat, tracers and momentum - downdraft transport of heat This model should not be used without upcoming developments, namely : - downdraft transport of tracers and momentum - updraft & downdraft transport of q2 (tke) - revision of vdif_kc to compute q2 for non-stratified cases Thermals could also include in a later revision : - momentum loss during transport (horizontal drag) Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90 ================================================ ================================================ M libf/phymars/callkeys.h M libf/phymars/inifis.F Added new control flags to call the thermals : - calltherm (false by default) <- to call thermals - outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes) ================================================ M libf/phymars/vdifc.F ^------> added a temporary output for thermal-related diagnostics M libf/phymars/testphys1d.F ^------> added treatment for a initialization from a profile of neutral gas (ar) -> will be transformed in a decaying tracer for thermal diagnostics M libf/phymars/physiq.F ^------> added a section to call the thermals -> changed the call to convadj -> added thermal-related outputs for diagnostics M libf/phymars/convadj.F ^------> takes now into account the height of thermals to execute convective adjustment => note : convective adjustment needs to be activated when using thermals, in case of a second instable layer above the thermals ================================================ A libf/phymars/calltherm_interface.F90 ^------> Interface between physiq.F and the thermals A libf/phymars/calltherm_mars.F90 ^------> Routine running the sub-timestep of the thermals A libf/phymars/thermcell_main_mars.F90 ^------> Main thermals routine specific to Martian physics A libf/phymars/thermcell_dqupdown.F90 ^------> Thermals subroutine computing transport of quantities by updrafts and downdrafts A libf/phymars/thermcell.F90 ^------> Module including parameters from the Earth to Mars importation. Will disappear in future dev ================================================ ================================================ == 17/06/2011 == EM >>> 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 == 17/06/2011 == AC >>> Important updates to thermals parameters - Tuned aspect ratio of thermals to suit Buoyancy estimations from LES in CLOSURE relation - Renormalization of alim_star after plume - Removed alimentation mixing of estimated Teta in plume >>> Minor change in makegcm_ifort == 22/06/2011 == EM - added modifications (from JYC & FGG) to tracer.h & initracer.F for ions - minor improvement to newstart.F (q=x option, check that tracer index provided by user is valid). - minor correction to callradite.F (to enable compilation in debug mode with ifort when there is only one tracer). == 17/06/2011 == AC - Added maximum vertical velocity and heat flux output from thermals - Added buoyancy diagnostics - Minor modifications in thermals routines == 23/06/2011 == EM - correct bug (introduced previously) in lect_start_archive.F on loop boundaries for soil temperature. == 01/07/2011 == AC - Added new settings for the Martian thermals from new LES observations - Revamped thermcell's module variables to allow it's removal - Minor changes in physiq and meso_physiq for the call to thermals - Switched from dynamic to static memory allocation for all thermals variable to gain computation speed == 04/07/2011 == AC - Minor setting modification to thermals - Added new flux optimization in thermcell_main_mars.F90, to run properly in 3D == 14/07/2011 == JBM - Tidying up dust properties in DATAFILE for better consistency (cf. suaer.F90) - Cosmetic changes in aeropacity.F (changed comment and put a print inside a water flag) == 15/07/2011 == EM >> Implemented using 'z0' roughness length map (important: 'z0' reference field is in datafile surface.nc, which has also been updated). - made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc). - adapted 'tabfi.F' to use 'z0_default'. - adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used. - modified 'physdem1.F' to write 'z0' field to restart.nc - adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'. - adapted 'dustdevil.F' to use 'z0_default'. - 'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant"). - modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file. == 15/07/2011 == AS - Modified the mesoscale part so that the previous change by EM does not imply an error in the mesoscale case. More development is needed though to get the "varying z0" capability in the mesoscale model. - Worked on versions of meso_physiq and meso_inifis as close as possible to physiq and inifis for more continuity in the process of impacting changes (and even possibly to reach a common version of physiq and inifis). >>> The main point is to make the mesoscale significant specific parts coded into include files in meso_inc so that meso_physiq and meso_inifis looks very close to physiq and inifis. >>> This is completely transparent for GCM users who does not need the contents of meso_inc. - Slight cosmetic changes to physiq.f and inifis.F --- some of them e.g. to prepare convergence between meso_physiq and physiq - Dropped the meso_ prefix from the slope routines (currently those are only used in the mesoscale model, but one could imagine using those in hi-res GCM runs) - Added a MESOSCALE flag to writediagfi so that mesoscale simulations are not bothered with calls put in routines by GCM people == 19/07/2011 == AS - Finished converging meso_physiq.F and meso_inifis.F towards physiq.F and inifis.F --> see previous point 15/07/2011 --> meso_ routines no longer exist (everything is in meso_inc and transparent to GCM users) --> GCM routines include mesoscale parts within MESOSCALE precompiler commands --> MESOSCALE parts are as hidden as possible not to mess up with GCM users/developers - Cleaned inelegant or useless #ifdef [or] #ifndef MESOSCALE in physiq and inifis so that a minimum amount of such precompiler commands is now reached [mainly related to I/O] - Added the SF08 slope insolation model in the general physics parameterizations. Added a callslope keyword in inifis.F and callkeys.h --> This keyword is False by default / True if you use -DMESOSCALE - Removed the obsolete call to Viking Lander 1 diagnostic Replaced it with a diagnostic for opacity at the domain center [valid for GCM and mesoscale] == 01/08/2011 == EM - Fixed bug in readtesassim (timeflag was not always set before being used) == 03/08/2011 == AC - Modified physiq.F to interface new SL parametrization - New SL parametrization based on a bulk Richardson Monin-Obukhov theory formulation Stability functions are taken from D.E. England et al. (95) Similarity functions coefficients based on Dyer and Hicks (70). Includes thermal roughness length computation, heat and momentum drag coefficient computation Can be used to output hydrodynamic-related SL quantities (bulk Richardson, turbulent Prandtl number estimation, Reynolds number) - Minor modification in thermcell_dqupdown.F to suit picky compilers - vdifc.F Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity) - Minor modification in meso_inc_les.F: u* is now taken from the new vdifc and not recomputed from a simple law == 08/08/2011 == AC - Formula correction for the Reynolds in vdif_cd : will affect the thermal roughness length. Results obtained with the precedent formula (previous commit) probably underestimate by about 1K the temperature in the mixing layer. - Added a new subroutine to libf/phymars called surflayer_interpol.F. This routine can be called as a diagnostic for interpolation of horizontal velocity norm and temperature in the surface layer, in 1D or 3D. The output height is controlled by a parameter z_out, to be set in physiq.F. This might be later added to a .def file. z_out is set to 0. by default, resulting in no call to this subroutine. == 16/08/2011 == AC - Minor modification in physiq.F to pass Ch from vdifc to meso_inc_les - Major modification to the formulation of integrals in surflayer_interpol.F Now stable for most cases. Some cases with highly negative Monin Obukhov length remain to be explored. - Added gustiness to the Richardson computation in vdif_cd.F. Gustiness factor is for now of beta=1., after several comparisons with LES aerodynamic conductances. May be subject to a minor change (+/- 0.1) in the near future. (almost transparent for the user) - Minor modifications relative to variables in vdifc.F. - Esthetic change to calltherm_mars.F - Changed the definition for HFX computation in the LES in meso_inc/meso_inc_les.F. New definition yields very similar results to old one and follows a strict definition of what HFX should be. == 24/08/11 == TN Attempts to tune the water cycle by adding outliers + A few structural changes !! * watercap.h is now obsolete and removed -- all is in surfdat.h * watercaptag initialized in surfini.F (up to 7 areas defined) instead of initracer.F - settings proposed by AS commented - experiments by TN decommented. use with caution. * water ice albedo and thermal inertia in callphys.def and inifis.F * water ice albedo in surfini.F * water ice albedo computation in albedocaps.F90 * alb_surfice is now obsolete in physiq.F, albedo_h2o_ice is used instead * frost_albedo_threshold defined in surfdat.h * water ice thermal inertia in soil.F TODO: * calibrate thermal inertia and ice albedo * have a look at subgrid-scale ice with dryness ? == 07/09/2011 == AC - Added new flag for the Richardson-based surface layer : callrichsl, can be changed in callphys.def One should always use the thermals model when using this surface layer model. Somes cases (weakly unstable with low winds), when not using thermals, won't be well represented by the Richardson surface layer. This stands for Mesoscale and Gcm but not for LES model. Correct configs : callrichsl = .true. calltherm = .true. callrichsl = .false. calltherm = .false. callrichsl = .false. calltherm = .true. Previously unstable config : callrichsl = .true. calltherm = .false. - To be able to run without thermals and with the new surface layer, a modification has been made to physiq.F to account for gustiness in GCM and MESOSCALE for negative Richardson, so that : callrichsl = .true. calltherm = .false. can now be used without problems, but is not recommended. - Consequently, callrichsl = .false. is now the default configuration for thermals. We recall the available options in callphys.def for thermals : outptherm = BOOLEAN (.false. by default) : outputs thermals related quantities (lots of diagfi) nsplit_thermals = INTEGER (50 by default in gcm, 2 in mesoscale) : subtimestep for thermals model. It is recommended to use at least 40 in the gcm, and at least 2 in the mesoscale. The user can lower these values but should check it's log for anomalies or errors regarding tracer transport in the thermals, or "granulosity" in the outputs for wmax, lmax and hfmax. == 08/09/11 == AS LMDZ.MARS + MESOSCALE. ---> Setting up a more realistic water ice source at the poles (notably outliers) [[ surfini.F ]] Main changes and bug fixes * reference to comcstfi.h was wrong. big problem because e.g. pi was not known. * commented about a problem to be fixed, due to surfini being called before initracer. * MESOSCALE: put the mesoscale north cap definition into a precompiling flag #MESOSCALE for the moment: if [alb_mean_TES > 0.26 and lat > 70] then outliers (previously done in meso_inc_caps.F) [[ inifis.F ]] Just changed a comment with wrong formatting --> below, only MESOSCALE [[ soil.F ]] if somewhere IT > IT_outliers, then makes it = IT_outliers [[ physiq.F ]] [[ meso_inc/meso_inc_caps.F ]] [[ meso_inc/meso_inc_ini.F ]] meso_inc_caps no longer called. keep for reference for the moment. [[ meso_inc/meso_inc_var.F ]] deleted lines with *_lim variables, now useless == 08/09/11 == EM -Fixed problem about undefined tracer names in 'surfini.F' by calling 'initracer' before 'surfini' in physiq.F == 09/09/11 == AS --> Fixed a problem with .eq. used with booleans in physiq.F [some good picky compilers complain about this] --> Added a warning in inifis about using callrichsl = .false. calltherm = .true. which is not recommended. The new surface layer has been built to go with the new mixing layer scheme (thermals). And anyway it is a much better scheme than the previous one. == 09/09/11 == AC This is a major update for thermals and richardson layer parametrization. This update stabilizes thermals (further studies might show that we can reduce the value of nsplit in gcm. To be tested.), and prevent downdrafts from descending into the surface layer, which was acting as a cold feedback on the thermals. The Richardson surface layer now features different gustiness coefficients for Richardson, heat and momentum so that u* and t* are correctly represented. Upcoming updates will change surflayer_interpol.F90 to implement those changes in the interpolation scheme as well. *************************** IMPORTANT : several parts of the vdifc code might want to use these new definitions for gustiness, u* and t*. exemple : dust devil routines that recompute u* ? lifting routines ? TODO ! ************************** M 289 libf/phymars/thermcell_main_mars.F90 ^-----------------> Added iterations to the velocity / buoyancy / entrainment / detrainment computation to ensure correct convergence. Iteration number is for now set to 4, which is probably too much. This will be changed once several cases are tested. The minimum is probably 2 iterations. M 289 libf/phymars/vdifc.F ^-----------------> Subgrid gustiness parametrization now uses different gustiness beta coefficients for heat and momentum. Comparisons with LES at Exomars landing site, matching u* and teta* suggests values of beta=0.7 for momentum and beta=1.2 for heat, where the formula for large scale horizontal winds in the first layer is : U0^2 = pu^2 + pv^2 + (beta*wmax_th)^2 M 289 libf/phymars/vdif_cd.F ^-----------------> LES data suggests that Richardson number distribution during daytime has a very large standart deviation due to highly negative Richardsons on several gridpoints. As a consequence, the mean Richardson in daytime in the LES is much more negative than in LES. In the gcm, parametrized subgrid turbulence prevents such cases, which can be dangerous in nearly unstable conditions. Hence, we use beta=0.5 instead of one, so that we keep the anti-crazy-hfx function of beta and we increase the norm of negative Richardsons in general for daytime conditions. This is set in conjonction with beta settings for heat and momentum in vdifc. M 289 libf/phymars/meso_inc/meso_inc_les.F ^-----------------> HFX and USTM computations now uses the different betas for heat and momentum. == 21/09/11 == AS --> Added MESOINI precompiler flag so that all fields needed to initialize the mesoscale model are being output in the diagfi.nc file. Previously it was made through a separate physiq.F which needed to be updated every time... --> This is completely transparent to the casual GCM users and only appears in the WRITEDIAGFI section of physiq.F == 21/09/11 == AC Revision on several settings for the thermals model. This version relies on fits done for an article to be published, and is more precise. M 299 libf/phymars/thermcell_main_mars.F90 ^-----------------> Changed downdraft to updraft mass flux ratio from -1.8 to -1.9 Changed first level for downdraft mixing from k=3 to k=2. Only level 1 is non-mixed now. Changed coefficients for downdraft to updraft thermal buoyancy ratios. M 299 libf/phymars/calltherm_mars.F90 ^-----------------> Changed r_aspect_thermals from 2. to 1.5 for the GCM version to better match buoyancy profiles. == 29/09/11 == AS --> To easily explore sensitivity to lifting thresholds: in dustlift.F, ustar_seuil=sqrt(stress/rho) and alpha_lift[dust_mass] can be prescribed through an external stress.def parameter file. --- alpha_lift[dust_number] is computed from alpha_lift[dust_mass] as in initracer.F --- ustar_seuil is more user-friendly than stress because direct comparison with ustar from model --> For the moment this is MESOSCALE only, but potentially useful to everyone == 30/09/11 == TN >> Bug correction in lect_start_archive.F ; in some cases layer(:) was not initialized. == 10/10/2011 == AC *********** This commit aims at increasing the thermals speed. Using these corrections, gcm performances in 64x48x32 using 1 tracer goes from 27.9% elapsed time in thermals to 18.76%. *********** Additional work needs to be done in tracer advection to gain speed in high tracer number configuration. (tracer advection (but not momentum nor temperature) could be decoupled from sub-timestep, as they do not act on the thermals scheme (water vapor is neglected as we use theta and not theta_v, and radiative effect of dust is not computed in the thermals.)) *********** => TOP 5 of routine contributions to gcm runtime : Each sample counts as 0.01 seconds. % cumulative self self total time seconds seconds calls s/call s/call name 18.76 6.33 6.33 960 0.01 0.01 thermcell_main_mars_ 17.19 12.13 5.80 __svml_powf4.A 13.72 16.76 4.63 10369 0.00 0.00 filtreg_ 3.94 18.09 1.33 __intel_new_memset 3.73 19.35 1.26 2880 0.00 0.00 thermcell_dqupdown_ note: thermcell_main_mars_ does call quite a lot power computations (see __svml_powf4.A), but this number will not increase with tracer numbers. *********** => LOG: M 312 libf/phymars/thermcell_main_mars.F90 ^------------------- removed (commented) computations on buoyancy which is purely diagnostic tuned internal convergence loop and added convergence criterion M 312 libf/phymars/thermcell_dqupdown.F90 ^------------------- removed (commented) downdraft-related if-loops (as we do not advect tracers and momentum in downdrafts for now) M 312 libf/phymars/calltherm_mars.F90 ^------------------- removed (commented) diagnostic-related computations changed default thermals spliting and aspect ratio corrected a bug where maximum height was not correctly computed and could result in convective adjustment used in place of thermals when using certains sets of nsplit and r_aspect (was not happening with the baseline version, so that this correction is transparent to users) ******************** == 14/10/2011 == AS + AC - Monin-Obukhov length is now output from surflay_interpol and written in diagfi if z_out not 0. - in calltherm_interface, defaut settings for qtransport_thermals and dtke_thermals == 20/10/2011 == EM - added FF's upgrade of writediagfi. Now, if at runtime there is a diagfi.def file, it should contain the list of variables (1 per line) than will be put in the diagfi.nc file. If there is no diagfi.def file, then all variables are put in the diagfi.nc file (as was the case before). == 21/10/2011 == EM - Corrected small bug in newstart: initracer was not always used and thus some tracer indexes (igm_co2, igcm_h2o_vap,...) were not set. This however means that we now also call inifis from newstart and that we read in flags set in 'callphys.def' (required for sanity checks in initracer). Also adapted 'inichim_readcallphys': removed some obsolescent tests on number of tracers for given combinations of options. == 21/10/2011 == AS - Added possibility for CH4 tracer in tracer.h and initracer.F == 27/10/11 == AS --> Added geticecover.F90 which computes seasonal ice cover given ls, lati(ngrid), long(ngrid) as proposed by T. Titus from TES observations [fitting functions for crocus line] ... output is icecover(ngrid) which value is 0 [no ice cover] or 1 [ice cover] ... no calculations are done for latitudes between -40 and +40 [ice cover is directly set to 0] --> In physiq.F, co2ice is set to a dummy high value to simulate a CO2 cap wherever icecover(ngrid) is 1. This is done at each timestep before newcondens is called. --> For the moment this is MESOSCALE only, but potentially useful to everyone. == 28/10/11 == FL + AS ADDED THE NEW FRAMEWORK FOR PHOTOCHEMISTRY This is not the last version. A new rewritten version of calchim.F [using LAPACK] is yet to be committed. --> A new version for photochemistry routines : *_B.F no longer exists, new routines in aeronomars D 333 libf/aeronomars/photochemist_B.F D 333 libf/aeronomars/init_chimie_B.F A 0 libf/aeronomars/read_phototable.F M 333 libf/aeronomars/calchim.F A 0 libf/aeronomars/photochemistry.F M 333 libf/aeronomars/chimiedata.h --> Changed calls to calchim and watercloud [surfdust and surfice needed] in physiq.F --> Left commented code for outputs in physiq.F [search for FL] --> Added settings which works for 35 levels in inidissip.F according to FL. Commented for the moment. --> Checked compilation and run, looks fine. Note that 'jmars.20101220' is needed. == 02/11/2011 == AC *********** This commit aims at increasing the thermals speed, especially for large tracer number configurations. The idea behind this commit is to advect non-active conserved variables outside of the sub-timestep of the thermals. Because these variables are not used in thermals computation, we can decouple them: momentum: can be decoupled because we assume a constant ratio between horizontal velocity in alimentation layer and maximum vertical velocity in the thermal s. tracer: can be decoupled because we do not take condensation of any tracer into account and hence do not liberate latent heat nor form clouds in the thermal s. temperature: cannot be decoupled (of course) *********** D 336 libf/phymars/thermcell_dqupdown.F90 ^---------------- Deleted and replaced by a simpler version. Notes about downdraft advection are still available from revision 336 of SVN in thermcell_dqu pdown. A 0 libf/phymars/thermcell_dqup.F90 ^---------------- New upward advection for tracers and momentum in thermals. Several changes are done compared to the old approach: - Updraft quantities are not longer computed by making hypothesis on the amount of advected air. - In general, the formalism for updraft computation is much simpler and clearer. - Tracer tendancies are no longer computed using the conservation equation. Instead, we use the divergence of an approximated turbulent flux of the concerned quantity, where downdraft are also neglected. M 336 libf/phymars/thermcell_main_mars.F90 ^---------------- The Main does not call anymore thermcell_dqupdown, which it was doing 2+tracer_number times per subtimestep (140 times per physical step for a 2 tracer config) M 336 libf/phymars/calltherm_mars.F90 ^---------------- Entrainment, detrainment and mass-fluxes are recomputed in the sub-timestep loop. Their final value after iterations is used by the new advection routine to compute tracer and momentum fluxes. *********** Results - Conservation of tracers has been assessed over 1 yr in 1D and found to be comparable to that obtained with the simple convective adjustment. (it actually seems to be better by a factor of 10%!) - GCM speed-up is of about 20% compared to the old thermal configuration, for a 2 tracer case. - Advection of sharp tracer profiles has been successfully observed, similar to the old method. == 07/11/11 == JBM >>> Changed watercloud.F to call two separate routines, simpleclouds.F or improvedclouds.F, which are a simplified and full microphysical scheme for cloud formation, respectively. Removed the tag called "improved" in watercloud.F, and added another tag called "microphys" which is defined in callphys.F instead. Changed routines: callkeys, inifis, physiq, watercloud. >>> Reimplemented the use of typical profiles of dust particle sizes and CCNs in simpleclouds.F. Corrected the previously used analytical CCN profile. Moved ccn_factor to simpleclouds.F, which won't be used in the new microphysical scheme. Changed routines: aeropacity, physiq, simpleclouds, watercloud. >>> Computed rdust at the end of callsedim instead of updatereffrad, except at firstcall. Renamed rfall into rsedcloud and moved it in simpleclouds. Moved nuice_sed to initracer.F and added it to tracer.h. Changed routines: callsedim, physiq, tracer.h, watercloud, initracer, simpleclouds, updatereffrad. >>> Added two tracers for the CCNs. Added the different tests in initracer.F to make sure that, for example, microphys is not used without doubleq, etc. Corrected an inconsistency in callsedim.F, and changed the way r0 is updated. Changes routines: callsedim, inifis, initracer, physiq, testphys1d, tracer.h. >>> Added the ability to have a spatially variable density in newsedim (same method as that used for the radius of sedimentation). Required the addition of one input to newsedim, which is the size of the density variable "rho". Changed routines: callsedim, newsedim. >>> Added an output to aeropacity called "tauscaling", which is a factor to convert dust_mass and dust_number to true quantities, based on tauvis. Removed ccn and qdust from aeropacity, which are obsolete now. >>> Wrote improvedcloud.F which includes all the microphysical scheme, and connected it to the sedimentation scheme. Added and changed routines: callsedim, physiq, growthrate, nuclea, improvedclouds, initracer, watercloud, microphys.h. == 07/11/11 == TN >> Added CCN & dust tracers updates in physiq.F Corrected a bug that can give negative CCN numbers, due to the use of single precision values (only 6 decimals) whereas up to 10E+10 CCN can be formed or disappears... >> Corrected physical bug that causes h2o_vap or h2o_ice to be negative in improvedclouds.F. >> Corrected physical bug that causes CCN & dust tracers to be negative in improvedclouds.F when all ice sublimates and releases dust >> Added parameter contact mteta in callphys.def Default value is still 0.95, see inifis.F >> Changed tendancy computation for dust_number in improvedclouds.F that was not the right one. Indeed, the scavenged dust_number tracer is computed from the dust_mass one, and its tendancy before scavenging must be taken into account to compute its scavenging's tendancy. == 09/11/11 == AS Added a more recent version of concentrations.F by FL == 22/11/11 == EM >> Added FL corrections to aeronomars/photochemistry.F : - Bug correction on the chemistry time step (now set to 1/3rd of the physics time step) - Updated kinetics coefficients (according to JPL 2011) == 22/11/11 == TN >>Changed scavenging in improvedclouds.F : CCNs are now virtual like dust to avoid tauscaling divergence. >>Water ice radius is now updated after sedimentation in callsedim.F. Added tau, tauscaling and zzlay as arguments >>Commented aggressive outputs at first call in suaer.F90, aeropacity.F == 23/11/11 == FGG + MALV >> New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified: -physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool. -inifis.F. Reading of nltemodel is added. -callkeys.h Declaration of nltemodel is added. The following lines should be added to callphys.def (ideally after setting callnlte): # NLTE 15um scheme to use. # 0-> Old scheme, static oxygen # 1-> Old scheme, dynamic oxygen # 2-> New scheme nltemodel = 2 A new directory, NLTEDAT, has to be included in datagcm. >> Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified: -nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file. -physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs -inifis: Reading new flag nircorr. -callkeys.h: Declaration of nircorr. The following lines have to be added to callphys.def (ideally after callnirco2): # NIR NLTE correction for variable SZA and O/CO2? # matters only if callnirco2=T # 0-> no correction # 1-> include correction nircorr=1 A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm. >> Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F == 24/11/11 == TN >> corrected minor bug in updatereffrad.F : reffdust was not saved >> ccn_factor as not correctly used in sedimentation. It is now initialized in inifis.F, declared in tracer.h and used in both simpleclouds.F & callsedim.F to update ice radius. >> Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users. == 06/12/11 == TN >> tiny change : nuice_sed initialisation is now done in inifis. Also changed initracer and improvedclouds. == 08/12/2011 == EM >> Include changes and updates for photochemistry by FL: - aeronomars/calchim.F : change in units of surface density. - aeronomars/surfacearea.F : new routine to compute ice and dust surface area (m2/m3) available for heterogeneous reactions. - phymars/initracer.F : bug correction: initialize igcm_ch4 and change loop bounds (when guessing tracer names/properties with old input files). - phymars/watercloud.F : cleanup. - phymars/physiq.F : add call to surfacearea; photochemistry is now called after sedimentation to take into acount updated rdust and rice. == 08/12/2011 == EM >> more updates for photochemistry: -aeronomars/inichim_newstart.F : add methane + bug correction in "oldnames" case -aeronomars/read_phototable.F : adapt routine so that photolysis table file name may be set in the def file with phototable = filename (default is "jmars.20111014") == 09/12/2011 == EM >> more updates for photochemistry from FL: improved aeronomars/surfacearea.F , with a change in arguments. == 14/12/2011 == EM >> fixed bug in callsedim.F: recomputation of rdust and rice should only be done if the appropriate tracers and flags are set. >> commented out call to uertst in ch.F (uertst crashes with gfortran and is moreover puerly diagnostic messages). >> In nltecool.F: Enforced using igcm_* indexes to identify tracers and not hard coded fixed numbers (which would lead to a disaster in cases where the order of tracers was changed). == 14/12/2011 == AC >> Improved computation of mixing in vdifc for temperatures less than the saturation temperature. The scheme works by computing mixing of T and then liberating latent heat when T> This scheme is based on initial work done by Melanie Vangvichith and has been adapted/tested/bug-fixed and improved to Mars. == 15/12/2011 == EM >> minor correction in aeronomars/calchim.F and aeronomars/surfacearea.F : calls to wstats should always be embeded in if (callstats) condition. == 19/12/2011 == TN >> new option in newstart: 'ini_h2osurf', water ice surface initialisation to keep seasonal frost but get rid of polar cap positive or negative values. == 20/12/2011 == AS >> new option possible in callphys.def: tituscap [logical] .true. will call functions geticecover to yield a nice evolving CO2 cap as measured by Titus (mostly useful for mesoscale) .false. is the default option, so this change is transparent to the casual user == 11/01/2012 == AC >> Corrections to the thermals scheme following latest revisions of the related paper. Replaced the surface layer interpolator by a more complete routine, that will ultimately become a post-processing utility (and disappear from libf) and a subroutine in the MCD. == 23/11/11 == FGG >> Cleanup of the NLTE routines which have been packed together to limit the number of files. Also enforced that file names are non-capitalized (needed by the create_make_gcm scripts to better evaluate dependencies when building the makefile). == 07/02/12 == JYC >> Adding the new molecular diffusion routines (old "moldiff" and "moldiffcoeff" are kept for now and if further tests are needed): moldiffcoeff_red.F, moldiff_red.F90, diffusion.h Also updated "phymars/comdiurn.h", "aeronomars/conc.h" and "aeronomars/chimiedata.h" to be fixed/free form compatible. == 09/02/12 == EM >> Updated the makegcm(s) so that default behaviour is to set LMDGCM env variable to be the directory in which the makegcm script is. Updated the makegcm_* to use "SOURCE" to identify main program to compile >> Adapted "create_make_gcm" script so that it can compile a main program that is either something.F or something.F90. Also added that object files are removed from archive before compiling a new version. == 10/02/12 == TN >> Major update on watercycle: a smaller integration timestep is now used in watercloud.F, sedimentation of clouds is done in watercloud instead of callsedim.F >> Temperature-dependant contact parameter in nuclea.F >> No dust lifting if CO2 ice in vdif.c >> Ice integrated column opacity is written in diagfi from physiq.F, instead of aeropacity.F. Mandatory if iradia is not 1. >> New definition of permanent ice in surfini.F and possibility to have an ice cap in it in 1d. >> Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging == 15/02/12 == AS >> Optimization of watercycle new capabilities, e.g.: ... in nuclea, removed 'if' from function fshape, modified **, decrease number of overall variables ... in growthrate and newsedim, modified ** ... in improvedclouds, modified **, improved recursive calls to erf (now loop w/ 2 erf + 1 log instead of 4 erf + 4 log) >> Checked consistency with previous formulations >> Code 35% faster with a Valles Marineris mesoscale run in springtime (set imicro=25) ... probably even better at cloudier seasons == 23/02/12 == AS >> In improvedclouds and waterclouds, added optimizations ... moved calculations to firstcall ... moved 'if' statement so that only useful calculations are performed >> Code about 10% faster [tests carried out with mesoscale summer Tharsis run] >> Check with a GCM run that results similar == 27/02/12 == AC >> Continuation of thermals setting, comparisons with mesoscale results on Case C >> Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES ... This is directly comparable to the variable tke_heat_flux in namelist.input ... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable) ... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results between the two models are comparable (might require a slitghly different tke_heat_flux between the two models due to difference in vertical diffusion schemes and subgrid effects) ... Syntax for use is to add "tke_heat_flux = ..." in callphys.def >> Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...) == 02/03/12 == EM >> 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 >> added missing call to surfacearea before call to calchim in physiq (mea culpa) and also moved photochemistry so it occurs after sedimentation. == 06/03/12 == EM >> Added option "q=profile" in newstart to initialize a given tracer with a profile specified in an ASCII file. >> Minor fix on lect_start_archive : when tracers were not found, a default value was asked to user, but twice! Once is enough. == 07/03/12 == AC ******************************************************** >> UPDATE ON PARAMETERS TO SET TO USE THE THERMALS MODEL ******************************************************** >> IN CALLPHYS.DEF << --------------------- - calltherm = .true. > activates the thermals - callrichsl = .true. > activate the Richardson Monin Obukhov surface layer with subgrid gustiness, and the yamada4 diffusion model from earth (corrected version of vdif_kc) - leave calladj = .true. > so that the convective adjustment may work ABOVE the PBL. >> IN RUN.DEF << ---------------- - 96 physicial timesteps per day (the internal sub-timestep of thermals will adapt so that it costs only 10% more CPU time to go from "48pdt+thermals" to "96pdt+thermals".) > this has to be respected because of a coupling between the radiative transfer and thermals > if you really wish to run with 48 pdt/day, you may want to deactive the thermals... >> USEFULL DIAGNOSTICS << ------------------------- - activating the thermals will automatically add 3 2D-outputs to diagfi.nc files: > zmax_th : interpolated boundary layer height (m) > hfmax_th : maximum vertical eddy heat flux reached in thermals (K.m/s) > wstar : free convection velocity reached in thermals (m/s) A rule of thumb for determining the mean velocities in subgrid-scale structures is: Mean updraft velocity : ~ wstar Maximum updraft velocity : ~ 2 or 3.*wstar Mean downdraft velocity : ~ 0.66*wstar Maximum absolute value of downdraft velocity : ~ 2.*wstar >> AT COMPILATION TIME << ------------------------- - +1 vertical level during the compilation, and then using the updated z2sig.def >> UPDATED Z2SIG.DEF << ----------------------- - one level has been added in the PBL, so that to reach 150km you have to compile with 33 vertical levels instead of 32. 10.00000 H: atmospheric scale height (km) (used as a reference only) 0.0020736 0.015735 0.0600151 0.163344 0.36358 0.708007 1.25334 2.06571 3.22069 4.80327 6.90787 9.63834 13.108 18.9666 25.0626 31.5527 38.4369 45.4369 52.4369 59.4369 66.4369 73.4369 80.4369 87.4369 94.4369 101.4369 108.437 115.437 122.437 129.437 136.437 143.437 150.437 157.437 164.437 171.437 178.437 185.437 192.437 199.437 206.437 213.437 220.437 227.437 234.437 241.437 248.437 255.437 262.437 269.437 276.437 283.437 290.437 297.437 304.437 311.437 318.437 325.437 332.437 339.437 346.437 353.437 360.437 367.437 374.437 381.437 388.437 395.437 ******************************************************** >> END OF UPDATE ON PARAMETERS FOR THE THERMALS MODEL ******************************************************** == 12/03/12 == EM >> correction in nirco2abs: io and ico2 should be declared as integers >> Update sponge: remove posibility of specifying 'hsponge', all modes apply to the last upper "nsponge" layers (default nsponge=3) and sponge time scale doubles with every level, starting from top. == 12/03/12 == FL + EM >> update a coefficient in photochemistry.F and some cleanup in watercloud.F and physiq.F. == 29/03/12 == JYC + EM >> Update vertical grid resolution for diffusion. == 30/03/12 == EM >> Updated model so that 1) reference pressure for opacity is now 610 Pa (and not 700 Pa anymore) and 2) the new MY24-MY30 dust scenarios (input files dust_MY24.nc, dust_MY25.nc, ..., dust_MY30.nc) can be used (when setting iaervar=24,25,...,30 in callphys.def): changed "readtesassim.F90" into "read_dust_scenario.F90" and apadpted aeropacity.F. one can still use the "old" MY24-MY25-MY26 (input files dust_tes_MY24.nc, dust_tes_MY25.nc and dust_tes_MY26.nc) scenarios by setting iaervar=124, 125 or 126 in callphys.def. == 31/03/12 == AS & FF >> Deleted obsolete lines in solarlong == 13/04/12 == FL >> Update of the chemistry package: - calchim.F90 : - upgraded from calchim.F - can run with or without CH4 - fixed the mass conservation scheme - photochemistry.F : -chemistry timestep is now independant from physics timestep -can run with or without a CH4 tracer - removed initial tests on species, since these are already done in calchim - removed inichim_readcallphys.F (not used any more) - concentrations.F: adaptated to better handle indexes and molar masses of tracers. "ncomp" (in chimiedata.h) no longer needed. - inichim_newstart.F90 : rewritten and cleaned (won't be compatible with old start files where tracer names are q01,...). Now handles hybrid levels and automaticaly adapts depending on which tracers are available. - newstart.F : adapted to followup changes in inchim_newstart.F90 and some cleanup around the initialization of tracer names and surface values. == 18/04/12 == TN >> New scheme in improvedclouds.F based on an analytical solution of the crystal radius growth equation >> Sedimentation of clouds done in callsedim.F, as it was done before version 520 >> Latent heat release of sublimating ground ice in vdifc.F >> Bug corrected in lwu.F, that was a source of NaN for very thick water ice clouds >> Bug corrected in updatereffrad.F, ice and dust radius are now updated each timestep before radiative transfer == 20/04/12 == EM >> Some cleanup on messages and comments in code about the reference pressure for dust opacity which is now 610 Pa. == 25/04/12 == TN >> The new scheme does not work for now. Instead, use of a subtimestep for nucleation and growth. == 26/04/12 == FGG+FL >> Update of the chemistry package, including: - 93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11) - computation of photoabsorption using raytracing - improved time stepping in the photochemistry - updated parameters (cross-sections); with this new version input files are in 'EUV/param_v5' of "datafile" directory. - transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90) - Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h - Upadted surfacearea.F - Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers). - Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense. == 03/05/12 == JYC+EM - Update of molecular diffusion routines (diffusion.h, moldiff_red.F90 and moldiffcoeff_red.F): optimized to run faster and diffuse only specific tracers (if they are present). - bugs fixed in simpleclouds.F (extra arguments) and watercloud.F (wrong arguments in call to simpleclouds and rdust() should only be recomputed if dust_mass & dust_number tracers are available). == 04/05/12 == JYC+EM - some syntax corrections in thermcall_main_mars, vdif_cd, pbl_parameters which cause problems when compiling with some strict compilers (g95, gfortran) - added an initialisation of 'varian' in initracer for cases when using conrath dust; because that value can be is used elsewhere (e.g. surfacearea) == 11/05/12 == FGG+FL - aeronomars/inichim_newstart.F : initialization of chemistry now handles nitrogen species and ions - aeronomars/conc.h : cleanup of obsolete variables - aeronomars/chimiedata.h : cleanup of obsolete variables == 12/05/12 == FGG+JYC+EM - updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F). - minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) . - minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed. - in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array). == 21/05/12 == JYC+FL - corrected bug in computation of HCO2+ in paramfoto_compact.F - avoid infinite H2O saturation ratio at very low temperatures in improvedclouds.F == 29/05/12 == EM - Added the possibility to use the "cold" (iaervar=6) or "warm" (iaervar=7) dust scenarios (derived from the 7 MY24 to MY30 scenarios) in aeropacity.F and read_dust_scenario.F == 04/06/12 == AC - Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration using the new physics. == 05/06/12 == AS - Change few function names in nlte_aux.F so that nesting is possible. This is completely transparent to GCM users. == 06/06/12 == EM - Minor cleanup in surfini.F - Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F) == 14/06/12 == EM - update z2sig.def file in deftank with improved vertical discretization for better thermals (this update should should have been done a while back...) == 14/06/12 == FGG - Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag. - Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir"). - Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars". == 15/06/12 == EM - Minor fix in "nuclea.F", max() and min() functions must have arguments of identical types. - Bug correction in "physiq.F": only water (and possibly dust and cnn if using microphysics) tracer tendencies should be updated after call to watercloud. == 22/06/12 == EM - Enforce that imposed TES ice albedo can never be greater than 0.9 == 28/06/12 == JYC+EM - Decreased time step for molecular diffusion (diffusion.h) for better stability - Added test in moldiff_red.F90 to enforce that the system can't be stuck in an infinite loop. - The creation and output of coefficients in file 'coeffs.dat' in moldiff_red.F is now controled by an internal flag (by default set to false). == 12/07/12 == TN - Use of stats.def to write variables in stats.nc, same as diagfi.def - some outputs in physiq.F