Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py (revision 310) @@ -26,7 +26,17 @@ #title = "Atmospheric temperature" +name = "wrfout_d01_2024-03-22_01:00:00_z" +itime = 12 +ndiv = 10 +zey = 120 +var = "VMR_ICE" +title = "Volume mixing ratio of water ice [ppm]" +vmin = -0.5 +vmax = 400. + nc = Dataset(name) what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d2=zey ) + y = nc.variables["vert"][:] @@ -34,4 +44,10 @@ horinp = len(what_I_plot[0,:]) x = np.linspace(0.,horinp*500.,horinp) / 1000. +xlabel("Horizontal coordinate (km)") + +horinp = len(what_I_plot[0,:]) +x = np.linspace(0.,horinp,horinp) +xlabel("x grid points") + zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) @@ -46,5 +62,4 @@ extend='neither',spacing='proportional') ptitle(title) -xlabel("Horizontal coordinate (km)") ylabel("Altitude (m)") makeplotres(var+str(itime),res=200.,disp=False) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py (revision 310) @@ -23,5 +23,7 @@ from scipy.ndimage.measurements import minimum_position from netCDF4 import Dataset - + import matplotlib.pyplot as plt + import myplot as myp + ### LOAD NETCDF DATA filename = "psfc.nc" @@ -34,24 +36,51 @@ allsizesx = [] allsizesy = [] - for i in range(0,shape[0]): - #for i in range(0,2): - - print i, ' on ', shape[0] + stride = 1 #5 + for i in range(0,shape[0],stride): + psfc2d = np.array ( psfc [ i, : , : ] ) ############### CRITERION - ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy - where = np.where(psfc2d - ave < -0.3) ## comme le papier Phoenix + ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy + + #limdp = -0.2 ## on en loupe pas mal + #where = np.where(psfc2d - ave < limdp) ## comme le papier Phoenix - std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy - fac = 2.5 ## not too low, otherwise vortices are not caught. 4 good choice. + std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy + fac = 4. ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice. + ## 2.5 clearly too low, 3.5 not too bad, 4 probably good + fac = 3.5 + #fac = 2.5 lim = ave - fac*std - print lim, ave, std where = np.where(psfc2d < lim) ############### END CRITERION - + lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine lab[where] = 1. ## do not treat points close to 'mean' (background) pressure - + + draw = False + if draw: + ################################################################################## + vmin = -0.3 + vmax = 0.0 + ndiv = 3 + palette = plt.get_cmap(name="YlGnBu") + what_I_plot = psfc2d-ave + zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) + what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax) + zelevels = np.linspace(zevmin,zevmax) + fig = plt.figure() + sub = myp.definesubplot(2,fig) + plt.subplot(sub) + plt.contourf(what_I_plot,zelevels,cmap=palette) + plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\ + ticks=np.linspace(zevmin,zevmax,ndiv+1),\ + extend='both',spacing='proportional') + plt.subplot(sub+1) + palette = plt.get_cmap(name="binary") + plt.contourf(lab,2,cmap=palette) + plt.show() + ################################################################################## + xx = [] yy = [] @@ -69,9 +98,18 @@ sizex = detsize( xx, res = 10, loga=False, thres=2 ) sizey = detsize( yy, res = 10, loga=False, thres=2 ) + #sizex = detsize( xx, res = 10, loga=False, thres=3 ) + #sizey = detsize( yy, res = 10, loga=False, thres=3 ) ### + #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex ### un peu limite dans certains cas + if (len(sizex) > len(sizey)) : sizey = sizex ### plus fidele mais petit souci lorsque PBC + elif (len(sizex) == len(sizey)) : + if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex + else : sizex = sizey + else : sizex = sizey allsizesx = np.append(allsizesx,sizex) allsizesy = np.append(allsizesy,sizey) + print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex ######## - + allsizesx.sort() allsizesy.sort() @@ -79,11 +117,6 @@ return allsizesx, allsizesy -def writeascii ( tab, filename ): - mydata = tab - myfile = open(filename, 'w') - for line in mydata: - myfile.write(str(line) + '\n') - myfile.close() - return +######################################################################### +######################################################################### import matplotlib.pyplot as plt @@ -91,21 +124,23 @@ import numpy as np import matplotlib.mlab as mlab +import mymath as mym save = True -#save = False +save = False + if save: allsizesx, allsizesy = getsize("psfc.nc") - writeascii(allsizesx,'allsizex.txt') - writeascii(allsizesy,'allsizey.txt') - + ### sauvegarde texte pour inspection + mym.writeascii(allsizesx,'allsizex.txt') + mym.writeascii(allsizesy,'allsizey.txt') + ### sauvegarde binaire pour utilisation python myfile = open('allsizex.bin', 'wb') pickle.dump(allsizesx, myfile) myfile.close() - myfile = open('allsizey.bin', 'wb') pickle.dump(allsizesy, myfile) myfile.close() - +### load files myfile = open('allsizex.bin', 'r') allsizesx = pickle.load(myfile) @@ -113,16 +148,33 @@ allsizesy = pickle.load(myfile) +### append sizes plothist = np.append(allsizesx,allsizesy) +plothist = allsizesx +#plothist = allsizesy plothist.sort() -nbins=100 +### make bins +nbins = 100 zebins = [2.0] +#nbins = 8 +#zebins = [19.0] +#nbins = 15 +#zebins = [11.] ##12 non mais donne un peu la meme chose +#zebins = [20.] ##20 non car trop pres du premier +nbins = 100 +zebins = [2./np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde +nbins = 200 +zebins = [0.2/np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde + for i in range(0,nbins): zebins.append(zebins[i]*np.sqrt(2)) zebins = np.array(zebins) -print zebins - +### select reasonable bins for DD zebins = zebins [ zebins > 10. ] zebins = zebins [ zebins < 1000. ] -print zebins +#print 'bins ',zebins +print 'corrected bins ',zebins +print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy) + +#plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.] ## pour calibrer #### HISTOGRAM @@ -130,8 +182,7 @@ log=True,\ bins=zebins,\ - #cumulative=-1,\ +# cumulative=-1,\ ) plt.xscale('log') - plt.show() Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 310) @@ -14,2 +14,10 @@ return u'\u00b0' +def writeascii ( tab, filename ): + mydata = tab + myfile = open(filename, 'w') + for line in mydata: + myfile.write(str(line) + '\n') + myfile.close() + return + Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 310) @@ -450,9 +450,12 @@ from mymath import max,min,mean ### might be convenient to add the missing value in arguments - what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) + #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7) + if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7) + else: what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7) print "new min ", min(what_I_plot) what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 - what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) + what_I_plot[ what_I_plot > zevmax ] = zevmax print "new max ", max(what_I_plot) + return what_I_plot @@ -508,5 +511,5 @@ "USTM": "YlOrRd",\ "HFX": "RdYlBu",\ - "ICETOT": "YlGnBu",\ + "ICETOT": "YlGnBu_r",\ "MTOT": "PuBu",\ "TAU_ICE": "Blues",\ Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F (revision 310) @@ -497,6 +497,6 @@ ztoploc = 10 ! target pseudo-altitude of local storm (km) radloc = 0.5 ! radius of dust storm (degree) - lonloc = 80 ! center longitude of storm (deg) - latloc = -25 ! center latitude of storm (deg) + lonloc = 25. ! center longitude of storm (deg) + latloc = -3. ! center latitude of storm (deg) DO ig=1,ngrid Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F (revision 310) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F (revision 310) @@ -0,0 +1,749 @@ + SUBROUTINE initracer(qsurf,co2ice) + + IMPLICIT NONE +c======================================================================= +c subject: +c -------- +c Initialization related to tracer +c (transported dust, water, chemical species, ice...) +c +c Name of the tracer +c +c Test of dimension : +c Initialize COMMON tracer in tracer.h, using tracer names provided +c by the dynamics in "advtrac.h" +c +c Old conventions: (not used any more) +c +c If water=T : q(iq=nqmx) is the water mass mixing ratio +c and q(iq=nqmx-1) is the ice mass mixing ratio + +c If there is transported dust, it uses iq=1 to iq=dustbin +c If there is no transported dust : dustbin=0 +c If doubleq=T : q(iq=1) is the dust mass mixing ratio +c q(iq=2) is the dust number mixing ratio + +c If (photochem.or.thermochem) there is "ncomp" chemical species (ncomp +c is set in aeronomars/chimiedata.h) using the ncomp iq values starting at +c iq=nqchem_min = dustbin+1 (nqchem_min is defined in inifis.F) +c +c +c author: F.Forget +c ------ +c Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003) +c Ehouarn Millour (oct. 2008) identify tracers by their names +c======================================================================= + + +#include "dimensions.h" +#include "dimphys.h" +#include "comcstfi.h" +#include "callkeys.h" +#include "tracer.h" +#include "advtrac.h" +#include "comgeomfi.h" +#include "chimiedata.h" + +#include "surfdat.h" + + real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) + real co2ice(ngridmx) ! co2 ice mass on surface (e.g. kg.m-2) + integer iq,ig,count + real r0_lift , reff_lift, nueff_lift +c Ratio of small over large dust particles (used when both +c doubleq and the submicron mode are active); In Montmessin +c et al. (2002), a value of 25 has been deduced; + real, parameter :: popratio = 25. + logical :: oldnames ! =.true. if old tracer naming convention (q01,...) + character(len=20) :: txt ! to store some text + +c----------------------------------------------------------------------- +c radius(nqmx) ! aerosol particle radius (m) +c rho_q(nqmx) ! tracer densities (kg.m-3) +c alpha_lift(nqmx) ! saltation vertical flux/horiz flux ratio (m-1) +c alpha_devil(nqmx) ! lifting coeeficient by dust devil +c rho_dust ! Mars dust density +c rho_ice ! Water ice density +c nuice_ref ! Effective variance nueff of the +c ! water-ice size distributions +c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio +c varian ! Characteristic variance of log-normal distribution +c----------------------------------------------------------------------- + + integer :: nqchem_start + +! Initialization: get tracer names from the dynamics and check if we are +! using 'old' tracer convention ('q01',q02',...) +! or new convention (full tracer names) + ! check if tracers have 'old' names + oldnames=.false. + count=0 + do iq=1,nqmx + txt=" " + write(txt,'(a1,i2.2)') 'q',iq + if (txt.eq.tnom(iq)) then + count=count+1 + endif + enddo ! of do iq=1,nqmx + + if (count.eq.nqmx) then + write(*,*) "initracer: tracers seem to follow old naming ", + & "convention (q01,q02,...)" + write(*,*) " => will work for now ... " + write(*,*) " but you should run newstart to rename them" + oldnames=.true. + endif + + ! copy/set tracer names + if (oldnames) then ! old names (derived from iq & option values) + ! 1. dust: + if (dustbin.ne.0) then ! transported dust + do iq=1,dustbin + txt=" " + write(txt,'(a4,i2.2)') 'dust',iq + noms(iq)=txt + mmol(iq)=100. + enddo + ! special case if doubleq + if (doubleq) then + if (dustbin.ne.2) then + write(*,*) 'initracer: error expected dustbin=2' + else +! noms(1)='dust_mass' ! dust mass mixing ratio +! noms(2)='dust_number' ! dust number mixing ratio +! same as above, but avoid explict possible out-of-bounds on array + noms(1)='dust_mass' ! dust mass mixing ratio + do iq=2,2 + noms(iq)='dust_number' ! dust number mixing ratio + enddo + endif + endif + endif + ! 2. water & ice + if (water) then + noms(nqmx)='h2o_vap' + mmol(nqmx)=18. +! noms(nqmx-1)='h2o_ice' +! mmol(nqmx-1)=18. + !"loop" to avoid potential out-of-bounds in array + do iq=nqmx-1,nqmx-1 + noms(iq)='h2o_ice' + mmol(iq)=18. + enddo + endif + ! 3. Chemistry + if (photochem .or. callthermos) then + if (doubleq) then + nqchem_start=3 + else + nqchem_start=dustbin+1 + end if + do iq=nqchem_start,nqchem_start+ncomp-1 + noms(iq)=nomchem(iq-nqchem_start+1) + mmol(iq)=mmolchem(iq-nqchem_start+1) + enddo + endif ! of if (photochem .or. callthermos) + ! 4. Other tracers ???? + if ((dustbin.eq.0).and.(.not.water)) then + noms(1)='co2' + mmol(1)=44 + if (nqmx.eq.2) then + noms(nqmx)='Ar_N2' + mmol(nqmx)=30 + endif + endif + else + ! copy tracer names from dynamics + do iq=1,nqmx + noms(iq)=tnom(iq) + enddo + ! mmol(:) array is initialized later (see below) + endif ! of if (oldnames) + + ! special modification when using 'old' tracers: + ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour + ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice + if (oldnames.and.water) then + write(*,*)'initracer: moving surface water ice to index ',nqmx-1 + ! "loop" to avoid potential out-of-bounds on arrays + do iq=nqmx-1,nqmx-1 + qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1) + enddo + qsurf(1:ngridmx,nqmx)=0 + endif + +c------------------------------------------------------------ +c Test Dimensions tracers +c------------------------------------------------------------ + +! Ehouarn: testing number of tracers is obsolete... +! if(photochem.or.thermochem) then +! if (water) then +! if ((dustbin+ncomp+2).ne.nqmx) then +! print*,'initracer: tracer dimension problem:' +! print*,"(dustbin+ncomp+2).ne.nqmx" +! print*,"ncomp: ",ncomp +! print*,"dustbin: ",dustbin +! print*,"nqmx: ",nqmx +! print*,'Change ncomp in chimiedata.h' +! endif +! else +! if ((dustbin+ncomp+1).ne.nqmx) then +! print*,'initracer: tracer dimension problem:' +! print*,"(dustbin+ncomp+1).ne.nqmx" +! print*,"ncomp: ",ncomp +! print*,"dustbin: ",dustbin +! print*,"nqmx: ",nqmx +! print*,'Change ncomp in chimiedata.h' +! STOP +! endif +! endif +! endif + +c------------------------------------------------------------ +c NAME and molar mass of the tracer +c------------------------------------------------------------ + + +! Identify tracers by their names: (and set corresponding values of mmol) + ! 0. initialize tracer indexes to zero: + do iq=1,nqmx + igcm_dustbin(iq)=0 + enddo + igcm_dust_mass=0 + igcm_dust_number=0 + igcm_dust_submicron=0 + igcm_h2o_vap=0 + igcm_h2o_ice=0 + igcm_co2=0 + igcm_co=0 + igcm_o=0 + igcm_o1d=0 + igcm_o2=0 + igcm_o3=0 + igcm_h=0 + igcm_h2=0 + igcm_oh=0 + igcm_ho2=0 + igcm_h2o2=0 + igcm_n2=0 + igcm_ar=0 + igcm_ar_n2=0 + igcm_n=0 + igcm_no=0 + igcm_no2=0 + igcm_n2d=0 + igcm_co2plus=0 + igcm_oplus=0 + igcm_o2plus=0 + igcm_coplus=0 + igcm_cplus=0 + igcm_nplus=0 + igcm_noplus=0 + igcm_n2plus=0 + igcm_hplus=0 + igcm_elec=0 + + + ! 1. find dust tracers + count=0 + if (dustbin.gt.0) then + do iq=1,nqmx + txt=" " + write(txt,'(a4,i2.2)')'dust',count+1 + if (noms(iq).eq.txt) then + count=count+1 + igcm_dustbin(count)=iq + mmol(iq)=100. + endif + enddo !do iq=1,nqmx + endif ! of if (dustbin.gt.0) + if (doubleq) then + do iq=1,nqmx + if (noms(iq).eq."dust_mass") then + igcm_dust_mass=iq + count=count+1 + endif + if (noms(iq).eq."dust_number") then + igcm_dust_number=iq + count=count+1 + endif + enddo + endif ! of if (doubleq) + if (submicron) then + do iq=1,nqmx + if (noms(iq).eq."dust_submicron") then + igcm_dust_submicron=iq + mmol(iq)=100. + count=count+1 + endif + enddo + endif ! of if (submicron) + ! 2. find chemistry and water tracers + do iq=1,nqmx + if (noms(iq).eq."co2") then + igcm_co2=iq + mmol(igcm_co2)=44. + count=count+1 + endif + if (noms(iq).eq."co") then + igcm_co=iq + mmol(igcm_co)=28. + count=count+1 + endif + if (noms(iq).eq."o") then + igcm_o=iq + mmol(igcm_o)=16. + count=count+1 + endif + if (noms(iq).eq."o1d") then + igcm_o1d=iq + mmol(igcm_o1d)=16. + count=count+1 + endif + if (noms(iq).eq."o2") then + igcm_o2=iq + mmol(igcm_o2)=32. + count=count+1 + endif + if (noms(iq).eq."o3") then + igcm_o3=iq + mmol(igcm_o3)=48. + count=count+1 + endif + if (noms(iq).eq."h") then + igcm_h=iq + mmol(igcm_h)=1. + count=count+1 + endif + if (noms(iq).eq."h2") then + igcm_h2=iq + mmol(igcm_h2)=2. + count=count+1 + endif + if (noms(iq).eq."oh") then + igcm_oh=iq + mmol(igcm_oh)=17. + count=count+1 + endif + if (noms(iq).eq."ho2") then + igcm_ho2=iq + mmol(igcm_ho2)=33. + count=count+1 + endif + if (noms(iq).eq."h2o2") then + igcm_h2o2=iq + mmol(igcm_h2o2)=34. + count=count+1 + endif + if (noms(iq).eq."n2") then + igcm_n2=iq + mmol(igcm_n2)=28. + count=count+1 + endif + if (noms(iq).eq."ar") then + igcm_ar=iq + mmol(igcm_ar)=40. + count=count+1 + endif + if (noms(iq).eq."n") then + igcm_n=iq + mmol(igcm_n)=14. + count=count+1 + endif + if (noms(iq).eq."no") then + igcm_no=iq + mmol(igcm_no)=30. + count=count+1 + endif + if (noms(iq).eq."no2") then + igcm_no2=iq + mmol(igcm_no2)=46. + count=count+1 + endif + if (noms(iq).eq."n2d") then + igcm_n2d=iq + mmol(igcm_n2d)=28. + count=count+1 + endif + if (noms(iq).eq."co2plus") then + igcm_co2plus=iq + mmol(igcm_co2plus)=44. + count=count+1 + endif + if (noms(iq).eq."oplus") then + igcm_oplus=iq + mmol(igcm_oplus)=16. + count=count+1 + endif + if (noms(iq).eq."o2plus") then + igcm_o2plus=iq + mmol(igcm_o2plus)=32. + count=count+1 + endif + if (noms(iq).eq."coplus") then + igcm_coplus=iq + mmol(igcm_coplus)=28. + count=count+1 + endif + if (noms(iq).eq."cplus") then + igcm_cplus=iq + mmol(igcm_cplus)=12. + count=count+1 + endif + if (noms(iq).eq."nplus") then + igcm_nplus=iq + mmol(igcm_nplus)=14. + count=count+1 + endif + if (noms(iq).eq."noplus") then + igcm_noplus=iq + mmol(igcm_noplus)=30. + count=count+1 + endif + if (noms(iq).eq."n2plus") then + igcm_n2plus=iq + mmol(igcm_n2plus)=28. + count=count+1 + endif + if (noms(iq).eq."hplus") then + igcm_hplus=iq + mmol(igcm_hplus)=1. + count=count+1 + endif + if (noms(iq).eq."elec") then + igcm_elec=iq + mmol(igcm_elec)=1./1822.89 + count=count+1 + endif + if (noms(iq).eq."h2o_vap") then + igcm_h2o_vap=iq + mmol(igcm_h2o_vap)=18. + count=count+1 + endif + if (noms(iq).eq."h2o_ice") then + igcm_h2o_ice=iq + mmol(igcm_h2o_ice)=18. + count=count+1 + endif + ! Other stuff: e.g. for simulations using co2 + neutral gaz + if (noms(iq).eq."Ar_N2") then + igcm_ar_n2=iq + mmol(igcm_ar_n2)=30. + count=count+1 + endif + + + enddo ! of do iq=1,nqmx +! count=count+nbqchem + + ! check that we identified all tracers: + if (count.ne.nqmx) then + write(*,*) "initracer: found only ",count," tracers" + write(*,*) " expected ",nqmx + do iq=1,count + write(*,*)' ',iq,' ',trim(noms(iq)) + enddo + stop + else + write(*,*) "initracer: found all expected tracers, namely:" + do iq=1,nqmx + write(*,*)' ',iq,' ',trim(noms(iq)) + enddo + endif + + ! if water cycle but iceparty=.false., there will nevertheless be + ! water ice at the surface (iceparty is not used anymore, but this + ! part is still relevant, as we want to stay compatible with the + ! older versions). + if (water.and.(igcm_h2o_ice.eq.0)) then + igcm_h2o_ice=igcm_h2o_vap ! so that qsurf(i_h2o_ice) is identified + ! even though there is no q(i_h2o_ice) + else + ! surface ice qsurf(i_h2o_ice) was loaded twice by phyetat0, + ! as qsurf(i_h2o_vap) & as qsurf(i_h2o_ice), so to be clean: + if (igcm_h2o_vap.ne.0) then + qsurf(1:ngridmx,igcm_h2o_vap)=0 + endif + endif + +c------------------------------------------------------------ +c Initialisation tracers .... +c------------------------------------------------------------ + call zerophys(nqmx,rho_q) + + rho_dust=2500. ! Mars dust density (kg.m-3) + rho_ice=920. ! Water ice density (kg.m-3) + nuice_ref=0.1 ! Effective variance nueff of the + ! water-ice size distributions + + if (doubleq) then +c "doubleq" technique +c ------------------- +c (transport of mass and number mixing ratio) +c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio + + if( (nqmx.lt.2).or.(water.and.(nqmx.lt.4)) ) then + write(*,*)'initracer: nqmx is too low : nqmx=', nqmx + write(*,*)'water= ',water,' doubleq= ',doubleq + end if + + nueff_lift = 0.5 + varian=sqrt(log(1.+nueff_lift)) + + rho_q(igcm_dust_mass)=rho_dust + rho_q(igcm_dust_number)=rho_dust + +c Intermediate calcul for computing geometric mean radius r0 +c as a function of mass and number mixing ratio Q and N +c (r0 = (r3n_q * Q/ N)^(1/3)) + r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust) + +c Intermediate calcul for computing effective radius reff +c from geometric mean radius r0 +c (reff = ref_r0 * r0) + ref_r0 = exp(2.5*varian**2) + +c lifted dust : +c ''''''''''' + reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m) + alpha_devil(igcm_dust_mass)=9.e-9 ! dust devil lift mass coeff +c alpha_lift(igcm_dust_mass)=3.0e-15 ! Lifted mass coeff + alpha_lift(igcm_dust_mass)= 3.3e-10!1.e-6 !Lifted mass coeff + + r0_lift = reff_lift/ref_r0 + alpha_devil(igcm_dust_number)=r3n_q* + & alpha_devil(igcm_dust_mass)/r0_lift**3 + alpha_lift(igcm_dust_number)=r3n_q* + & alpha_lift(igcm_dust_mass)/r0_lift**3 + + radius(igcm_dust_mass) = reff_lift + radius(igcm_dust_number) = reff_lift + + write(*,*) "initracer: doubleq_param reff_lift:", reff_lift + write(*,*) "initracer: doubleq_param nueff_lift:", nueff_lift + write(*,*) "initracer: doubleq_param alpha_lift:", + & alpha_lift(igcm_dust_mass) + + else + + if (dustbin.gt.1) then + print*,'initracer: STOP!', + $ ' properties of dust need to be set in initracer !!!' + stop + + else if (dustbin.eq.1) then + +c This will be used for 1 dust particle size: +c ------------------------------------------ + radius(igcm_dustbin(1))=3.e-6 + alpha_lift(igcm_dustbin(1))=0.0e-6 + alpha_devil(igcm_dustbin(1))=7.65e-9 + rho_q(igcm_dustbin(1))=rho_dust + + endif + end if ! (doubleq) + +c Submicron dust mode: +c -------------------- + + if (submicron) then + radius(igcm_dust_submicron)=0.1e-6 + rho_q(igcm_dust_submicron)=rho_dust + if (doubleq) then +c If doubleq is also active, we use the population ratio: + alpha_lift(igcm_dust_submicron) = + & alpha_lift(igcm_dust_number)*popratio* + & rho_q(igcm_dust_submicron)*4./3.*pi* + & radius(igcm_dust_submicron)**3. + alpha_devil(igcm_dust_submicron)=1.e-30 + else + alpha_lift(igcm_dust_submicron)=1e-6 + alpha_devil(igcm_dust_submicron)=1.e-30 + endif ! (doubleq) + end if ! (submicron) + +c Initialization for photochemistry: +c --------------------------------- + if (photochem) then + ! initialize chemistry+water (water will be correctly initialized below) + ! by initializing everything which is not dust ... + do iq=1,nqmx + txt=noms(iq) + if (txt(1:4).ne."dust") then + radius(iq)=0. + alpha_lift(iq) =0. + alpha_devil(iq)=0. + endif + enddo ! do iq=1,nqmx + endif + +c Initialization for water vapor +c ------------------------------ + if(water) then + radius(igcm_h2o_vap)=0. + alpha_lift(igcm_h2o_vap) =0. + alpha_devil(igcm_h2o_vap)=0. + if(water.and.(nqmx.ge.2)) then + radius(igcm_h2o_ice)=3.e-6 + rho_q(igcm_h2o_ice)=rho_ice + alpha_lift(igcm_h2o_ice) =0. + alpha_devil(igcm_h2o_ice)=0. + elseif(water.and.(nqmx.lt.2)) then + write(*,*) 'nqmx is too low : nqmx=', nqmx + write(*,*) 'water= ',water + endif + + end if ! (water) + +c Output for records: +c ~~~~~~~~~~~~~~~~~~ + write(*,*) + Write(*,*) '******** initracer : dust transport parameters :' + write(*,*) 'alpha_lift = ', alpha_lift + write(*,*) 'alpha_devil = ', alpha_devil + write(*,*) 'radius = ', radius + if(doubleq) then + write(*,*) 'reff_lift (um) = ', reff_lift + write(*,*) 'size distribution variance = ', varian + write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0 + end if + +! +! some extra (possibly redundant) sanity checks for tracers: +! --------------------------------------------------------- + + if (doubleq) then + ! verify that we indeed have dust_mass and dust_number tracers + if (igcm_dust_mass.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use doubleq option without ", + & "a dust_mass tracer !" + stop + endif + if (igcm_dust_number.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use doubleq option without ", + & "a dust_number tracer !" + stop + endif + endif + + if ((.not.doubleq).and.(dustbin.gt.0)) then + ! verify that we indeed have 'dustbin' dust tracers + count=0 + do iq=1,dustbin + if (igcm_dustbin(iq).ne.0) then + count=count+1 + endif + enddo + if (count.ne.dustbin) then + write(*,*) "initracer: error !!" + write(*,*) " dusbin is set to ",dustbin, + & " but we only have the following dust tracers:" + do iq=1,count + write(*,*)" ",trim(noms(igcm_dustbin(iq))) + enddo + stop + endif + endif + + if (water) then + ! verify that we indeed have h2o_vap and h2o_ice tracers + if (igcm_h2o_vap.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use water option without ", + & "an h2o_vap tracer !" + stop + endif + if (igcm_h2o_ice.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use water option without ", + & "an h2o_ice tracer !" + stop + endif + endif + + if (photochem .or. callthermos) then + ! verify that we indeed have the chemistry tracers + if (igcm_co2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a co2 tracer !" + stop + endif + if (igcm_co.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a co tracer !" + stop + endif + if (igcm_o.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a o tracer !" + stop + endif + if (igcm_o1d.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a o1d tracer !" + stop + endif + if (igcm_o2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "an o2 tracer !" + stop + endif + if (igcm_o3.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "an o3 tracer !" + stop + endif + if (igcm_h.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a h tracer !" + stop + endif + if (igcm_h2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a h2 tracer !" + stop + endif + if (igcm_oh.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "an oh tracer !" + stop + endif + if (igcm_ho2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a ho2 tracer !" + stop + endif + if (igcm_h2o2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a h2o2 tracer !" + stop + endif + if (igcm_n2.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "a n2 tracer !" + stop + endif + if (igcm_ar.eq.0) then + write(*,*) "initracer: error !!" + write(*,*) " cannot use chemistry option without ", + & "an ar tracer !" + stop + endif + endif ! of if (photochem .or. callthermos) + + end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F (revision 309) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F (revision 310) @@ -1,1 +1,1 @@ -link ../../../mars_lmd_new/libf/phymars/initracer.F +link STORM_JULIEN_LAST/initracer.F