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 598) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F (revision 609) @@ -1,5 +1,8 @@ SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, & pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, - & QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d, + & QVISsQREF3d,omegaVIS3d,gVIS3d, + & QIRsQREF3d,omegaIR3d,gIR3d, + & QREFvis3d,QREFir3d, + & omegaREFvis3d,omegaREFir3d, & zdqnorm,dsodust,pt) @@ -80,10 +83,25 @@ REAL reffrad(ngrid,nlayer,naerkind) REAL nueffrad(ngrid,nlayer,naerkind) - REAL QREFvis3d(ngridmx,nlayermx,naerkind) - REAL QREFir3d(ngridmx,nlayermx,naerkind) - REAL omegaREFvis3d(ngridmx,nlayermx,naerkind) - REAL omegaREFir3d(ngridmx,nlayermx,naerkind) +! REAL QREFvis3d(ngridmx,nlayermx,naerkind) +! REAL QREFir3d(ngridmx,nlayermx,naerkind) +! REAL omegaREFvis3d(ngridmx,nlayermx,naerkind) +! REAL omegaREFir3d(ngridmx,nlayermx,naerkind) REAL zdqnorm(ngridmx,nlayermx,2) !mass mixing ratio perturbation due to the dust storm. Output for meso_physiq.F REAL pt(ngrid,nlayer) !input: temperature. usefull to compute precisely the l_top parameter + +c Aerosol optical properties + REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind) + REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind) + REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind) + + REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind) + REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind) + REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind) + + REAL :: QREFvis3d(ngridmx,nlayermx,naerkind) + REAL :: QREFir3d(ngridmx,nlayermx,naerkind) + + REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind) + REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) @@ -142,4 +160,5 @@ logical localstorm ! =true to create a local dust storm real taulocref,ztoploc,radloc,lonloc,latloc !local dust storm parameters + real reffstorm, yeah REAL ray(ngridmx) !distance from dust storm center REAL tauuser(ngridmx) !opacity perturbation due to dust storm @@ -488,5 +507,4 @@ IF (firstcall) THEN - WRITE(*,*) " RENORMALISATION !!! " c-------------------------------------------------- c Parameters of the opacity perturbation @@ -494,11 +512,59 @@ iaer=1 !!!! PROVISOIRE !!!! - taulocref = 4.25 !10 ! ref optical depth of the local dust storm - ztoploc = 10 ! target pseudo-altitude of local storm (km) - radloc = 0.5 ! radius of dust storm (degree) - lonloc = 25. ! center longitude of storm (deg) - latloc = -3. ! center latitude of storm (deg) + + write(*,*) "Add a local storm ?" + localstorm=.true. ! default value + call getin("localstorm",localstorm) + write(*,*) " localstorm = ",localstorm + + IF (localstorm) THEN + WRITE(*,*) "********************" + WRITE(*,*) "ADDING A LOCAL STORM" + WRITE(*,*) "********************" + + !!! DEFAULT CASE + !4.25 ! ref optical depth of the local dust storm + !10 ! target pseudo-altitude of local storm (km) + !0.5 ! radius of dust storm (degree) (actual radius is twice this) + !25. ! center longitude of storm (deg) + !-2.5 ! center latitude of storm (deg) + !0.0 ! reff in storm (microns) 0. for background + + write(*,*) "ref opacity of local dust storm" + taulocref = 4.25 ! default value + call getin("taulocref",taulocref) + write(*,*) " taulocref = ",taulocref + + write(*,*) "target altitude of local storm (km)" + ztoploc = 10.0 ! default value + call getin("ztoploc",ztoploc) + write(*,*) " ztoploc = ",ztoploc + + write(*,*) "radius of dust storm (degree)" + radloc = 0.5 ! default value + call getin("radloc",radloc) + write(*,*) " radloc = ",radloc + + write(*,*) "center longitude of storm (deg)" + lonloc = 25.0 ! default value + call getin("lonloc",lonloc) + write(*,*) " lonloc = ",lonloc + + write(*,*) "center latitude of storm (deg)" + latloc = -2.5 ! default value + call getin("latloc",latloc) + write(*,*) " latloc = ",latloc + + write(*,*) "reff storm (mic) 0. for background" + reffstorm = 0.0 ! default value + call getin("reffstorm",reffstorm) + write(*,*) " reffstorm = ",reffstorm + ELSE + write(*,*) 'I need localstorm = T' + STOP + ENDIF DO ig=1,ngrid + c--------------------------------------- @@ -509,10 +575,14 @@ & (long(ig)*180./pi -lonloc)**2) + !! transition factor for storm + !! increase factor ray diff for steepness + yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2. + c------------------------------------------------- c Tau's new map: c------------------------------------------ - tauuser(ig)=max(tauref(ig) * pplev(ig,1) / 700.E0 , taulocref - & * (TANH(2.+(radloc-ray(ig))*2)+1.)/2.) + tauuser(ig)=max(tauref(ig) * pplev(ig,1) / 700.E0 , + & taulocref * yeah) c--------------------------------------------------------- @@ -529,10 +599,41 @@ ENDDO - ENDDO + +c--------------------------------------------------------- +c change reffrad if ever needed +c---------------------------------------------------------- + + IF (reffstorm .gt. 0.) THEN + DO l=1,nlayer + IF (l .lt. l_top+1) THEN + + WRITE(*,*) "OLD REFFRAD", reffrad(ig,l,iaer) + + reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm + & * 1.e-6 * yeah ) + + WRITE(*,*) "NEW REFFRAD", reffrad(ig,l,iaer), reffstorm + & * 1.e-6 * yeah + + ENDIF + ENDDO + ENDIF + +c--------------------------------------------------------------------------- + ENDDO !! boucle sur ngridmx +c--------------------------------------------------------------------------- + + !!! ReComputing 3D scattering parameters: if needed. + IF (reffstorm .gt. 0.) THEN + CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad, + & QVISsQREF3d,omegaVIS3d,gVIS3d, + & QIRsQREF3d,omegaIR3d,gIR3d, + & QREFvis3d,QREFir3d, + & omegaREFvis3d,omegaREFir3d) + ENDIF c---------------------------------------------------------------------------- c compute int_factor c--------------------------------------------------------------------------- - DO ig=1,ngrid Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F (revision 598) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F (revision 609) @@ -391,5 +391,8 @@ CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, & pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, - & QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d, + & QVISsQREF3d,omegaVIS3d,gVIS3d, + & QIRsQREF3d,omegaIR3d,gIR3d, + & QREFvis3d,QREFir3d, + & omegaREFvis3d,omegaREFir3d, & zdqnorm,dsodust,pt) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F (revision 598) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F (revision 609) @@ -1,1 +1,1 @@ -link ../../../mars_lmd_new/libf/phymars/callsedim.F +link callsedim.F.modified Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified (revision 609) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified (revision 609) @@ -0,0 +1,370 @@ + SUBROUTINE callsedim(ngrid,nlay, ptimestep, + $ pplev,zlev, pt, rdust, rice, + & pq, pdqfi, pdqsed,pdqs_sed,nq) + +! to use 'getin' + USE ioipsl_getincom + + IMPLICIT NONE + +c======================================================================= +c Sedimentation of the Martian aerosols +c depending on their density and radius +c +c F.Forget 1999 +c +c Modified by J.-B. Madeleine 2010: Now includes the doubleq +c technique in order to have only one call to callsedim in +c physiq.F. +c +c======================================================================= + +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "dimensions.h" +#include "dimphys.h" +#include "comcstfi.h" +#include "tracer.h" +#include "callkeys.h" + +c +c arguments: +c ---------- + + INTEGER ngrid ! number of horizontal grid points + INTEGER nlay ! number of atmospheric layers + REAL ptimestep ! physics time step (s) + REAL pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa) + REAL pt(ngrid,nlay) ! temperature at mid-layer (K) + REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries +c Aerosol radius provided by the water ice microphysical scheme: + REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m) + REAL rice(ngrid,nlay) ! Ice geometric mean radius (m) + +c Traceurs : + integer nq ! number of tracers + real pq(ngrid,nlay,nq) ! tracers (kg/kg) + real pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1) + real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1) + real pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1) + +c local: +c ------ + + REAL CBRT + EXTERNAL CBRT + + INTEGER l,ig, iq + + real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers + real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2) + real epaisseur (ngridmx,nlayermx) ! Layer thickness (m) + real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2) + real r0(ngridmx,nlayermx) ! geometric mean doubleq radius (m) +c Sedimentation radius of water ice + real rfall(ngridmx,nlayermx) +c Sedimentation effective variance of water ice + REAL, PARAMETER :: nuice_sed = 0.45 !! TESTS_JB !! 0.1 avant + +c Discrete size distributions (doubleq) +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c 1) Parameters used to represent the changes in fall +c velocity as a function of particle size; + integer nr,ir + parameter (nr=12) !(nr=7) ! number of bins + real rd(nr),qr(ngridmx,nlayermx,nr) + real rdi(nr+1) ! extreme and intermediate radii + real Sq(ngridmx,nlayermx) + real rdmin,rdmax,rdimin,rdimax + data rdmin/1.e-8/ !/1.e-7/ + data rdmax/30.e-6/ + data rdimin/1.e-10/ + data rdimax/1e-4/ + save rd, rdi + +c 2) Second size distribution for the log-normal integration +c (the mass mixing ratio is computed for each radius) + + integer ninter, iint + parameter (ninter=4) ! nombre de point entre chaque rayon rdi + real rr(ninter,nr) + save rr + integer radpower + +c 3) Other local variables used in doubleq + + real reff(ngridmx,nlayermx,2) ! for diagnostic only + INTEGER idust_mass ! index of tracer containing dust mass + ! mix. ratio + INTEGER idust_number ! index of tracer containing dust number + ! mix. ratio + SAVE idust_mass,idust_number + +c Firstcall: + + LOGICAL firstcall + SAVE firstcall + DATA firstcall/.true./ + + REAL betaset + SAVE betaset + + +c ** un petit test de coherence +c -------------------------- + + IF (firstcall) THEN + IF(ngrid.NE.ngridmx) THEN + PRINT*,'STOP dans callsedim' + PRINT*,'probleme de dimensions :' + PRINT*,'ngrid =',ngrid + PRINT*,'ngridmx =',ngridmx + STOP + ENDIF + + write(*,*) "beta coefficient for sedimentation" + betaset = 0.5 ! default value + call getin("betaset",betaset) + write(*,*) " betaset = ",betaset + +c Doubleq: initialization + IF (doubleq) THEN + do ir=1,nr + rd(ir)= rdmin*(rdmax/rdmin)**(float(ir-1)/float(nr-1)) + end do + rdi(1)=rdimin + do ir=2,nr + rdi(ir)= sqrt(rd(ir-1)*rd(ir)) + end do + rdi(nr+1)=rdimax + + do ir=1,nr + do iint=1,ninter + rr(iint,ir)= + & rdi(ir)* + & (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1)) +c write(*,*) rr(iint,ir) + end do + end do + + ! identify tracers corresponding to mass mixing ratio and + ! number mixing ratio + idust_mass=0 ! dummy initialization + idust_number=0 ! dummy initialization + + do iq=1,nq + if (noms(iq).eq."dust_mass") then + idust_mass=iq + endif + if (noms(iq).eq."dust_number") then + idust_number=iq + endif + enddo + + ! check that we did find the tracers + if ((idust_mass.eq.0).or.(idust_number.eq.0)) then + write(*,*) 'callsedim: error! could not identify' + write(*,*) ' tracers for dust mass and number mixing' + write(*,*) ' ratio and doubleq is activated!' + stop + endif + ENDIF !of if (doubleq) + + IF (water) THEN + write(*,*) "water_param nueff Sedimentation:", nuice_sed + IF (activice) THEN + write(*,*) "water_param nueff Radiative:", nuice_ref + ENDIF + ENDIF + + firstcall=.false. + ENDIF ! of IF (firstcall) + +c----------------------------------------------------------------------- +c 1. Initialization +c ----------------- + + zqi(1:ngrid,1:nlay,1:nqmx) = 0. +c Updating the mass mixing ratio with the tendencies coming +c from other parameterizations: +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + do iq=1,nq + do l=1,nlay + do ig=1,ngrid + zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep + enddo + enddo + enddo + +c Computing the different layer properties +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c Mass (kg.m-2), thickness(m), crossing time (s) etc. + + do l=1,nlay + do ig=1, ngrid + masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g + epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l) + end do + end do + +c ================================================================= + do iq=1,nq + if(radius(iq).gt.1.e-9) then ! no sedim for gaz + +c ----------------------------------------------------------------- +c DOUBLEQ CASE +c ----------------------------------------------------------------- + if (doubleq.and. + & ((iq.eq.idust_mass).or. + & (iq.eq.idust_number))) then + +c Computing size distribution: +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + do l=1,nlay + do ig=1, ngrid + r0(ig,l)= + & CBRT(r3n_q*zqi(ig,l,idust_mass)/ + & max(zqi(ig,l,idust_number),0.01)) + r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6) + end do + end do + +c Computing mass mixing ratio for each particle size +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (iq.EQ.idust_mass) then + radpower = 2 + ELSE + radpower = -1 + ENDIF + Sq(1:ngrid,1:nlay) = 0. + do ir=1,nr + do l=1,nlay + do ig=1,ngrid +c **************** +c Size distribution integration +c (Trapezoid Integration Method) + qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))* + & (rr(1,ir)**radpower)* + & exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*varian**2)) + do iint=2,ninter-1 + qr(ig,l,ir)=qr(ig,l,ir) + + & 0.5*(rr(iint+1,ir)-rr(iint-1,ir))* + & (rr(iint,ir)**radpower)* + & exp(-(log(rr(iint,ir)/r0(ig,l)))**2/ + & (2*varian**2)) + end do + qr(ig,l,ir)=qr(ig,l,ir) + + & 0.5*(rr(ninter,ir)-rr(ninter-1,ir))* + & (rr(ninter,ir)**radpower)* + & exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/ + & (2*varian**2)) + +c **************** old method (not recommended!) +c qr(ig,l,ir)=(rd(ir)**(5-3*iq))* +c & exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*varian**2) ) +c ****************************** + + Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir) + enddo + enddo + enddo + + do ir=1,nr + do l=1,nlay + do ig=1,ngrid + qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l) + enddo + enddo + enddo + +c Computing sedimentation for each tracer +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +c call zerophys(ngridmx*nlayermx,zqi(1,1,iq)) + zqi(1:ngrid,1:nlay,iq) = 0. +c call zerophys(ngridmx,pdqs_sed(1,iq)) + pdqs_sed(1:ngrid,iq) = 0. + + do ir=1,nr + call newsedim(ngrid,nlay,1,ptimestep, + & pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir), + & wq,betaset) + +c Tendencies +c ~~~~~~~~~~ + do ig=1,ngrid + pdqs_sed(ig,iq) = pdqs_sed(ig,iq) + & + wq(ig,1)/ptimestep + end do + DO l = 1, nlay + DO ig=1,ngrid + zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir) + ENDDO + ENDDO + enddo ! of do ir=1,nr +c ----------------------------------------------------------------- +c WATER CYCLE CASE +c ----------------------------------------------------------------- + else if (water.and.(iq.eq.igcm_h2o_ice)) then +c if (doubleq) then +c do l=1,nlay +c do ig=1,ngrid +c rfall(ig,l)=max( rice(ig,l),rdust(ig,l) ) +c rfall(ig,l)=min(rfall(ig,l),1.e-4) +c enddo +c enddo ! of do l=1,nlay +c else + do l=1,nlay + do ig=1,ngrid +c For water cycle, a typical dust size is assumed: +c r(z)=r0*exp(-z/H) with r0=0.8 micron and H=18 km. +c rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) ) + rfall(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3., + & rdust(ig,l) ) +c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns +c mars commente pour l'instant rfall(ig,l)=20.e-6 + rfall(ig,l)=min(rfall(ig,l),1.e-4) + enddo + enddo ! of do l=1,nlay +c endif + call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, + & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi(1,1,iq), + & wq,1.0) +c Tendencies +c ~~~~~~~~~~ + do ig=1,ngrid + pdqs_sed(ig,iq)=wq(ig,1)/ptimestep + end do +c ----------------------------------------------------------------- +c GENERAL CASE +c ----------------------------------------------------------------- + else + call newsedim(ngrid,nlay,1,ptimestep, + & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq), + & zqi(1,1,iq),wq,1.0) +c Tendencies +c ~~~~~~~~~~ + do ig=1,ngrid + pdqs_sed(ig,iq)=wq(ig,1)/ptimestep + end do + endif ! of if doubleq and if water +c ----------------------------------------------------------------- + + DO l = 1, nlay + DO ig=1,ngrid + pdqsed(ig,l,iq)=(zqi(ig,l,iq)- + $ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep + ENDDO + ENDDO + + endif ! of if(radius(iq).gt.1.e-9) +c ================================================================= + enddo ! of do iq=1,nq + + RETURN + END + Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def (revision 609) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def (revision 609) @@ -0,0 +1,39 @@ +### "Add a local storm ?" +### default is T +######################### +localstorm = .true. + +### "ref opacity of local dust storm" +### default is 4.25 +##################################### +taulocref = 4.25 + +### "target altitude of local storm (km)" +### default is 10.0 +######################################### +ztoploc = 10.0 + +### "radius of dust storm (degree)" +### default is 0.5 +################################### +radloc = 0.5 + +### "center longitude of storm (deg)" +### default is 25.0 +##################################### +lonloc = 25.0 + +### "center latitude of storm (deg)" +### default is -2.5 +#################################### +latloc = -2.5 + +### "reff storm (mic) 0. for background" +### default is 0.0 (use background reff) +############################################# +reffstorm = 0.0 + +### "beta coefficient for sedimentation" +### default is 0.5 +######################################## +betaset = 0.5 Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F (revision 598) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F (revision 609) @@ -296,4 +296,5 @@ SAVE ccn !! in case iradia != 1 real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m) + real reffdust(ngridmx,nlayermx) ! effective radius (m) real qtot1,qtot2 ! total aerosol mass integer igmin, lmin @@ -1046,4 +1047,9 @@ ENDDO ENDDO + + DO ig=1,ngrid + zdqssed_diag(ig)= zdqssed(ig,igcm_dust_mass) + ENDDO + END IF ! of IF (sedimentation) @@ -1456,6 +1462,8 @@ . * mugaz / mmol(igcm_h2o_ice) ENDIF + !! Dust quantity integration along the vertical axe dustot(:)=0 + reffdust(:,:)=0 do ig=1,ngrid do l=1,nlayermx @@ -1463,4 +1471,5 @@ & zq(ig,l,igcm_dust_mass) & * (pplev(ig,l) - pplev(ig,l+1)) / g + reffdust(ig,l) = rdust(ig,l)*ref_r0 enddo enddo Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def (revision 609) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def (revision 609) @@ -0,0 +1,4 @@ +# some definitions for the physics, in file 'callphys.def' +INCLUDEDEF=callphys.def + +INCLUDEDEF=localstorm.def Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def (revision 609) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def (revision 609) @@ -0,0 +1,29 @@ +10. +0.000001 + +0.5 !! ulim +0.002 !! alpha + +10. +0.000001 + +0.5 +0.002 + +0.5 +0.005 + +0.5 +0.001 + +0.5 +0.05 + +0.6 +0.1 + +0.6 +1. + +1.061 +1.