Index: /trunk/LMDZ.MARS/README =================================================================== --- /trunk/LMDZ.MARS/README (revision 2627) +++ /trunk/LMDZ.MARS/README (revision 2628) @@ -3556,12 +3556,12 @@ Also cleaned up and commented comsaison_h in the process. -== 04/01/2021 == CM +== 04/01/2022 == CM - Clean co2condens_mod.F - remove dqsurf duplication after call co2condens -== 04/01/2021 == AB +== 04/01/2022 == AB Some further cleaning of co2condens following r2599 and r2600 that fixed the bug appearing when the scavenging by CO2 ice was activated. -== 04/01/2021 == CM +== 04/01/2022 == CM following r2600, remove use co2condens_mod4micro in physiq_mod.F @@ -3608,2 +3608,19 @@ Open_MP :Small OpenMP fixes in conf_phys for reading radia.def with ifort +== 28/02/2022 == AB +Big changes on mountain top dust flows for GCM6: +- the scheme now activates only in grid meshes that contain a mountain among + a hard-written list, instead of every meshes. This is done to prevent strong + artificial reinjections of dust in places that don't present huge converging + slopes that concentrate dust (ex: Valles Marineris, Hellas). + Topdust is now also detrained as soon as it leaves the column it originated from. +- the list of the 19 allowed mountains is used by the subroutine topmons_setup + in module topmons_mod, to compute a logical array contains_mons(ngrid). + alpha_hmons and hsummit are also set up once and for all by this subroutine, + which is called in physiq_mod's firstcall. +- contains_mons, alpha_hmons and hsummit are now saved variables of the module + surfdat_h, and are called as such and not as arguments in the subroutines + using them. +- the logical flag "slpwind" and the comments in the code have also been updated + to the new terminology "mountain top dust flows", accordingly to ticket #71. + The new flag read in callphys.def is "topflows". Index: /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F (revision 2628) @@ -14,5 +14,5 @@ & QREFvis3d,QREFir3d,omegaREFir3d, & totstormfract,clearatm,dsords,dsotop, - & alpha_hmons,nohmons, + & nohmons, & clearsky,totcloudfrac) @@ -38,4 +38,6 @@ use dust_scaling_mod, only: compute_dustscaling use density_co2_ice_mod, only: density_co2_ice + use surfdat_h,only: alpha_hmons,contains_mons + IMPLICIT NONE c======================================================================= @@ -104,7 +106,6 @@ REAL, INTENT(IN) :: totstormfract(ngrid) ! mesh fraction with a rocket ! dust storm - LOGICAL, INTENT(IN) :: nohmons ! true to compute RT without slope wind - ! topdust, false to compute RT in the topdust - REAL, INTENT(IN) :: alpha_hmons(ngrid) + LOGICAL, INTENT(IN) :: nohmons ! true to compute RT without topdust, + ! false to compute RT in the topdust REAL,INTENT(OUT) :: tauscaling(ngrid) ! Scaling factor for qdust and Ndust REAL,INTENT(OUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment @@ -596,7 +597,7 @@ c and once in the part of the mesh without the sub-grid mountain (nohmons=true) aerosol(1:ngrid,1:nlayer,iaer) = 0. - IF (nohmons) THEN ! considering part of the mesh without storm + IF (nohmons) THEN ! considering part of the mesh without top flows aerosol(1:ngrid,1:nlayer,iaer)=1.e-25 - ELSE ! part of the mesh with concentred dust storm + ELSE ! part of the mesh with concentrated dust flows DO l=1,nlayer IF (l.LE.cstdustlevel) THEN @@ -883,9 +884,9 @@ c aerosol/X for topdust to prepare calculation of radiative transfer c ----------------------------------------------------------------- - IF (slpwind) THEN + IF (topflows) THEN DO ig=1,ngrid - IF (alpha_hmons(ig) .gt. 0.) THEN + IF (contains_mons(ig)) THEN ! contains_mons=True ensures that alpha_hmons>0 DO l=1,nlayer - ! topdust: opacity relative to the storm fraction (topdust/x) + ! topdust: opacity relative to the mons fraction (topdust/x) aerosol(ig,l,iaer_topdust_doubleq) = & aerosol(ig,l,iaer_topdust_doubleq)/alpha_hmons(ig) Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 2628) @@ -15,5 +15,5 @@ & ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds & & ,co2useh2o,meteo_flux,activeco2ice,CLFvaryingCO2,spantCO2 & - & ,CLFvarying,satindexco2,rdstorm,slpwind,calllott_nonoro & + & ,CLFvarying,satindexco2,rdstorm,topflows,calllott_nonoro & & ,latentheat_surfwater,gwd_convective_source,startphy_file & & ,hdo,hdofrac,cst_cap_albedo,temp_dependant_m,refill_watercap @@ -72,5 +72,5 @@ logical scavenging logical rdstorm ! rocket dust storm parametrization - logical slpwind ! entrainment by slope wind parametrization + logical topflows ! entrainment by mountain top dust flows parametrization logical latentheat_surfwater ! latent heat release from ground water ice sublimation/condensation logical cst_cap_albedo ! polar cap albedo remains unchanged by water frost deposition Index: /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F (revision 2628) @@ -12,5 +12,5 @@ $ taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice, $ rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop, - $ alpha_hmons,nohmons,clearsky,totcloudfrac) + $ nohmons,clearsky,totcloudfrac) use aeropacity_mod, only: aeropacity @@ -213,7 +213,6 @@ REAL,INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust -c entrainment by slope wind - LOGICAL, INTENT(IN) :: nohmons ! true for background dust - REAL, INTENT(IN) :: alpha_hmons(ngrid) ! sub-grid scale topography mesh fraction +c entrainment by mountain top dust flows + LOGICAL, INTENT(IN) :: nohmons ! true for background dust REAL,INTENT(OUT) :: rtopdust(ngrid,nlayer) ! Topdust geometric mean radius (m) REAL,INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust @@ -376,5 +375,5 @@ enddo end if - if (slpwind.AND.active) then + if (topflows.AND.active) then do iaer=1,naerkind if (name_iaer(iaer).eq."topdust_doubleq") then @@ -448,6 +447,6 @@ & QREFvis3d,QREFir3d,omegaREFir3d, & totstormfract,clearatm,dsords,dsotop, - & alpha_hmons,nohmons, - & clearsky,totcloudfrac) + & nohmons,clearsky,totcloudfrac) + c Starting loop on sub-domain c ---------------------------- Index: /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F (revision 2628) @@ -372,5 +372,5 @@ ENDIF !of if (rdstorm) - IF (slpwind) THEN ! identifying topdust tracers for sedimentation + IF (topflows) THEN ! identifying topdust tracers for sedimentation itopdust_mass=0 ! dummy initialization itopdust_number=0 ! dummy initialization @@ -392,8 +392,8 @@ write(*,*) 'callsedim: error! could not identify' write(*,*) ' tracers for topdust mass and number mixing' - write(*,*) ' ratio and slpwind is activated!' + write(*,*) ' ratio and topflows is activated!' call abort_physic(modname,"missing topdust tracers",1) endif - ENDIF !of if (slpwind) + ENDIF !of if (topflows) firstcall=.false. @@ -452,6 +452,6 @@ end do endif -c entrainment by slope wind - if (slpwind) then +c entrainment by mountain top dust flows + if (topflows) then do l=1,nlay do ig=1, ngrid @@ -710,5 +710,5 @@ endif ! of if (rdstorm) - if (slpwind) then + if (topflows) then DO l = 1, nlay DO ig=1,ngrid @@ -718,5 +718,5 @@ ENDDO ENDDO - endif ! of if (slpwind) + endif ! of if (topflows) c Update the ice particle size "rice" Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 2628) @@ -314,9 +314,9 @@ write(*,*)" coeff_detrainment = ",coeff_detrainment -! entrainment by slope wind scheme - write(*,*)"call slope wind lifting parametrization" - slpwind=.false. ! default value - call getin_p("slpwind",slpwind) - write(*,*)" slpwind = ",slpwind +! entrainment by mountain top dust flows scheme + write(*,*)"call mountain top dust flows parametrization" + topflows=.false. ! default value + call getin_p("topflows",topflows) + write(*,*)" topflows = ",topflows ! latent heat release from ground water ice sublimation/condensation @@ -471,18 +471,18 @@ endif endif -! rocket dust storm and entrainment by slope wind +! rocket dust storm and entrainment by top flows ! Test of incompatibility: -! if rdstorm or slpwind is used, then doubleq should be true - if ((rdstorm.or.slpwind).and..not.doubleq) then - print*,'if rdstorm or slpwind is used, then doubleq +! if rdstorm or topflows is used, then doubleq should be true + if ((rdstorm.or.topflows).and..not.doubleq) then + print*,'if rdstorm or topflows is used, then doubleq & should be used !' call abort_physic(modname, - & "rdstorm or slpwind requires doubleq",1) - endif - if ((rdstorm.or.slpwind).and..not.active) then - print*,'if rdstorm or slpwind is used, then active + & "rdstorm or topflows requires doubleq",1) + endif + if ((rdstorm.or.topflows).and..not.active) then + print*,'if rdstorm or topflows is used, then active & should be used !' call abort_physic(modname, - & "rdstorm or slpwind requires activ",1) + & "rdstorm or topflows requires activ",1) endif if (rdstorm.and..not.lifting) then @@ -492,9 +492,9 @@ & "rdstorm requires lifting",1) endif - if ((rdstorm.or.slpwind).and..not.freedust) then - print*,'if rdstorm or slpwind is used, then freedust + if ((rdstorm.or.topflows).and..not.freedust) then + print*,'if rdstorm or topflows is used, then freedust & should be used !' call abort_physic(modname, - & "rdstorm or slpwind requires freedust",1) + & "rdstorm or topflows requires freedust",1) endif if (rdstorm.and.(dustinjection.eq.0)) then @@ -900,25 +900,25 @@ ! and picky compilers who know name_iaer(2) is out of bounds j=2 - IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) then + IF (rdstorm.AND..NOT.activice.AND..NOT.topflows) then name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme j = j+1 END IF - IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind) then + IF (rdstorm.AND.water.AND.activice.AND..NOT.topflows) then name_iaer(j) = "stormdust_doubleq" j = j+1 END IF - IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) then + IF (topflows.AND..NOT.activice.AND..NOT.rdstorm) then name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme j = j+1 END IF - IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm) then + IF (topflows.AND.water.AND.activice.AND..NOT.rdstorm) then name_iaer(j) = "topdust_doubleq" j = j+1 END IF - IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN + IF (rdstorm.AND.topflows.AND..NOT.activice) THEN name_iaer(j) = "stormdust_doubleq" name_iaer(j+1) = "topdust_doubleq" @@ -926,5 +926,5 @@ ENDIF - IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN + IF (rdstorm.AND.topflows.AND.water.AND.activice) THEN name_iaer(j) = "stormdust_doubleq" name_iaer(j+1) = "topdust_doubleq" Index: /trunk/LMDZ.MARS/libf/phymars/initracer.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/initracer.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/initracer.F (revision 2628) @@ -181,5 +181,5 @@ enddo endif ! of if (rdstorm) - if (slpwind) then + if (topflows) then do iq=1,nq if (noms(iq).eq."topdust_mass") then @@ -192,5 +192,5 @@ endif enddo - endif ! of if (slpwind) + endif ! of if (topflows) ! 2. find chemistry and water tracers do iq=1,nq @@ -613,13 +613,13 @@ end if !(rdstorm) !c ---------------------------------------------------------------------- -!c slope wind scheme +!c mountain top dust flows scheme !c you need a radius value for topdust to active its sedimentation !c we take the same value as for the normal dust - if (slpwind) then + if (topflows) then rho_q(igcm_topdust_mass)=rho_dust rho_q(igcm_topdust_number)=rho_dust radius(igcm_topdust_mass) = 3.e-6 radius(igcm_topdust_number) = 3.e-6 - end if !(slpwind) + end if !(topflows) !c ---------------------------------------------------------------------- @@ -859,17 +859,17 @@ endif - if (slpwind) then + if (topflows) then ! verify that we indeed have topdust_mass and topdust_number tracers if (igcm_topdust_mass.eq.0) then write(*,*) "initracer: error !!" - write(*,*) " cannot use slpwind option without ", + write(*,*) " cannot use topflows option without ", & "a topdust_mass tracer !" - call abort_physic("initracer","slpwind issue",1) + call abort_physic("initracer","topflows issue",1) endif if (igcm_topdust_number.eq.0) then write(*,*) "initracer: error !!" - write(*,*) " cannot use slpwind option without ", + write(*,*) " cannot use topflows option without ", & "a topdust_number tracer !" - call abort_physic("initracer","slpwind issue",1) + call abort_physic("initracer","topflows issue",1) endif endif Index: /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90 (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90 (revision 2628) @@ -52,6 +52,4 @@ use rocketduststorm_mod, only: ini_rocketduststorm_mod, & end_rocketduststorm_mod - use topmons_mod, only: ini_topmons_mod, & - end_topmons_mod use calchim_mod, only: ini_calchim_mod,end_calchim_mod use watercloud_mod, only: ini_watercloud_mod, & @@ -137,8 +135,4 @@ call ini_rocketduststorm_mod(ngrid) - ! allocate arrays in "topmons_mod": - call end_topmons_mod - call ini_topmons_mod(ngrid,nlayer) - ! allocate arrays in "calchim_mod" (aeronomars) call end_calchim_mod Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2628) @@ -23,5 +23,5 @@ use rocketduststorm_mod, only: rocketduststorm, dustliftday use calcstormfract_mod, only: calcstormfract - use topmons_mod, only: topmons,alpha_hmons + use topmons_mod, only: topmons,topmons_setup use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice, & igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice, @@ -308,7 +308,7 @@ ! - in a mesh with stormdust and background dust (false) ! - in a mesh with background dust only (true) -c entrainment by slope winds +c entrainment by mountain top dust flows logical nohmons ! nohmons used to calculate twice the radiative - ! transfer when slpwind is active : + ! transfer when topflows is active : ! - in a mesh with topdust and background dust (false) ! - in a mesh with background dust only (true) @@ -507,8 +507,6 @@ REAL co2totB -c entrainment by slope winds above sb-grid scale topography +c entrainment by mountain top dust flows above sub-grid scale topography REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons - REAL hmax,hmin - REAL hsummit(ngrid) c when no startfi file is asked for init @@ -718,23 +716,7 @@ endif -c Initialize mountain mesh fraction for the entrainment by slope wind param. +c Initialize mountain mesh fraction for the entrainment by top flows param. c ~~~~~~~~~~~~~~~ - if (slpwind) then - !! alpha_hmons calculation - if (ngrid.gt.1) then - call planetwide_maxval(hmons,hmax ) - call planetwide_minval(hmons,hmin ) - do ig=1,ngrid - alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin) - enddo - else - hmin=0. - hmax=23162.1 !set here the height of the sub-grid scaled topography - do ig=1,ngrid - alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin) - print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig) - enddo - endif ! (ngrid.gt.1) - endif ! if (slpwind) + if (topflows) call topmons_setup(ngrid) #endif @@ -972,7 +954,7 @@ c Transfer through CO2 (except NIR CO2 absorption) c and aerosols (dust and water ice) - ! callradite for background dust + ! callradite for background dust (out of the rdstorm fraction) clearatm=.true. - !! callradite for background dust in the case of slope wind entrainment + !! callradite for background dust (out of the topflows fraction) nohmons=.true. @@ -988,5 +970,5 @@ & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice, & rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop, - & alpha_hmons,nohmons,clearsky,totcloudfrac) + & nohmons,clearsky,totcloudfrac) ! case of sub-grid water ice clouds: callradite for the clear case @@ -1006,6 +988,6 @@ & rice,nuice,riceco2, nuiceco2,co2ice,rstormdust, & rtopdust,totstormfract, - & clearatm,dsords,dsotop,alpha_hmons,nohmons, - & clearsky,totcloudfrac) + & clearatm,dsords,dsotop, + & nohmons,clearsky,totcloudfrac) clearsky = .false. ! just in case. ! Sum the fluxes and heating rates from cloudy/clear @@ -1200,5 +1182,5 @@ & clearsky,totcloudfrac, c input sub-grid scale topography - & nohmons,alpha_hmons, + & nohmons, c output & pdqrds,wspeed,dsodust,dsords,dsotop, @@ -1250,10 +1232,5 @@ c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography c ------------------------------------------- - IF (slpwind) THEN - if (ngrid.gt.1) then - hsummit(:)=summit(:)-phisfi(:)/g - else - hsummit(:)=14000. - endif + IF (topflows) THEN clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains nohmons=.false. @@ -1266,5 +1243,5 @@ & totstormfract,clearatm, & clearsky,totcloudfrac, - & nohmons,hsummit, + & nohmons, & pdqtop,wtop,dsodust,dsords,dsotop, & tau_pref_scenario,tau_pref_gcm) @@ -1287,5 +1264,5 @@ ENDDO - ENDIF ! end of if (slpwind) + ENDIF ! end of if (topflows) c 3.3 Dust injection from the surface @@ -3387,5 +3364,5 @@ endif ! (rdstorm) - if (slpwind) then + if (topflows) then call WRITEDIAGFI(ngrid,'refftopdust','refftopdust', & 'm',3,rtopdust*ref_r0) @@ -3404,5 +3381,5 @@ & 'm2.kg-1',3,dsotop) end select - endif ! (slpwind) + endif ! (topflows) if (dustscaling_mode==2) then Index: /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 (revision 2628) @@ -30,5 +30,5 @@ clearsky,totcloudfrac, & ! input sub-grid scale topography - nohmons,alpha_hmons, & + nohmons, & ! output pdqrds,wrad,dsodust,dsords,dsotop, & @@ -81,11 +81,10 @@ REAL,INTENT(OUT) :: dust_rad_adjust(ngrid) -! sbgrid scale water ice clouds +! subgrid scale water ice clouds logical, intent(in) :: clearsky real, intent(in) :: totcloudfrac(ngrid) -! sbgrid scale topography +! subgrid scale topography LOGICAL, INTENT(IN) :: nohmons - REAL, INTENT(IN) :: alpha_hmons(ngrid) !-------------------------------------------------------- @@ -257,5 +256,5 @@ tau,aerosol,dsodust,tauscaling,dust_rad_adjust, & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, & - totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,& + totstormfract,clearatm,dsords,dsotop,nohmons,& clearsky,totcloudfrac) Index: /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90 (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90 (revision 2628) @@ -35,4 +35,11 @@ !$OMP TESice_Scoef,iceradius,dtemisice, & !$OMP zmea,zstd,zsig,zgam,zthe,hmons,summit,base,z0,z0_default ) + + !! mountain top dust flows + REAL,SAVE,ALLOCATABLE :: alpha_hmons(:) ! sub-grid scale mountain mesh fraction + REAL,SAVE,ALLOCATABLE :: hsummit(:) ! mountain height above the GCM surface + LOGICAL,SAVE,ALLOCATABLE :: contains_mons(:) ! is there a mountain in the grid mesh ? + +!$OMP THREADPRIVATE(alpha_hmons,hsummit,contains_mons) !! variables @@ -75,5 +82,8 @@ allocate(summit(ngrid)) allocate(base(ngrid)) - + allocate(alpha_hmons(ngrid)) + allocate(hsummit(ngrid)) + allocate(contains_mons(ngrid)) + end subroutine ini_surfdat_h @@ -83,24 +93,27 @@ implicit none - if (allocated(albedodat)) deallocate(albedodat) - if (allocated(phisfi)) deallocate(phisfi) - if (allocated(watercaptag)) deallocate(watercaptag) - if (allocated(dryness)) deallocate(dryness) - if (allocated(zmea)) deallocate(zmea) - if (allocated(zstd)) deallocate(zstd) - if (allocated(zsig)) deallocate(zsig) - if (allocated(zgam)) deallocate(zgam) - if (allocated(zthe)) deallocate(zthe) - if (allocated(z0)) deallocate(z0) - if (allocated(qsurf)) deallocate(qsurf) - if (allocated(tsurf)) deallocate(tsurf) - if (allocated(co2ice)) deallocate(co2ice) - if (allocated(watercap)) deallocate(watercap) - if (allocated(emis)) deallocate(emis) - if (allocated(capcal)) deallocate(capcal) - if (allocated(fluxgrd)) deallocate(fluxgrd) - if (allocated(hmons)) deallocate(hmons) - if (allocated(summit)) deallocate(summit) - if (allocated(base)) deallocate(base) + if (allocated(albedodat)) deallocate(albedodat) + if (allocated(phisfi)) deallocate(phisfi) + if (allocated(watercaptag)) deallocate(watercaptag) + if (allocated(dryness)) deallocate(dryness) + if (allocated(zmea)) deallocate(zmea) + if (allocated(zstd)) deallocate(zstd) + if (allocated(zsig)) deallocate(zsig) + if (allocated(zgam)) deallocate(zgam) + if (allocated(zthe)) deallocate(zthe) + if (allocated(z0)) deallocate(z0) + if (allocated(qsurf)) deallocate(qsurf) + if (allocated(tsurf)) deallocate(tsurf) + if (allocated(co2ice)) deallocate(co2ice) + if (allocated(watercap)) deallocate(watercap) + if (allocated(emis)) deallocate(emis) + if (allocated(capcal)) deallocate(capcal) + if (allocated(fluxgrd)) deallocate(fluxgrd) + if (allocated(hmons)) deallocate(hmons) + if (allocated(summit)) deallocate(summit) + if (allocated(base)) deallocate(base) + if (allocated(alpha_hmons)) deallocate(alpha_hmons) + if (allocated(hsummit)) deallocate(hsummit) + if (allocated(contains_mons)) deallocate(contains_mons) end subroutine end_surfdat_h Index: /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 (revision 2628) @@ -2,9 +2,4 @@ IMPLICIT NONE - -! sub-grid scale mountain mesh fraction - REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) - -!$OMP THREADPRIVATE(alpha_hmons) CONTAINS @@ -32,5 +27,5 @@ clearsky,totcloudfrac, & ! input sub-grid scale mountain - nohmons,hsummit, & + nohmons, & ! output pdqtop,wfin,dsodust,dsords,dsotop, & @@ -43,5 +38,6 @@ USE dimradmars_mod, only: albedo,naerkind USE comsaison_h, only: dist_sol,mu0,fract - USE surfdat_h, only: emis,co2ice,hmons,summit + USE surfdat_h, only: emis,co2ice,hmons,summit,alpha_hmons, & + hsummit,contains_mons USE callradite_mod, only: callradite @@ -87,5 +83,4 @@ ! input sub-grid scale mountain LOGICAL, INTENT(IN) :: nohmons - REAL, INTENT(IN) :: hsummit(ngrid) !-------------------------------------------------------- @@ -196,5 +191,5 @@ ! Detrainment - REAL coefdetrain(ngrid,nlayer) ! coefficient for detrainment : % of stormdust detrained + REAL coefdetrain(ngrid,nlayer) ! coefficient for detrainment : % of topdust detrained REAL dqdet_topdust_mass(ngrid,nlayer) ! tendancy pdq topdust mass after detrainment only REAL dqdet_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after detrainment only @@ -206,6 +201,6 @@ ! ********************************************************************** ! ********************************************************************** - ! Parametrization of the entrainment by slope wind above the sub-grid - ! scale topography + ! Parametrization of the entrainment of dust by slope winds above the + ! converging sub-grid scale topography ("mountain top dust flows") ! ********************************************************************** ! ********************************************************************** @@ -279,4 +274,10 @@ ! 1.1. Call the second radiative transfer for topdust, obtain the extra heating ! ********************************************************************* + + ! NB: since the only grid meshes that matter for the radiative transfer + ! are the ones that contain a mount (contains_mons(ig)=.true.), + ! it could be relevant to optimize the code by calling the RT + ! only for those grid meshes. + ! See Ticket #92 on trac.lmd.jussieu.fr/Planeto CALL callradite(icount,ngrid,nlayer,nq,zday,zls,zq,albedo, & emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & @@ -285,5 +286,5 @@ tau,aerosol,dsodust,tauscaling,dust_rad_adjust, & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, & - totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,& + totstormfract,clearatm,dsords,dsotop,nohmons,& clearsky,totcloudfrac) ! ********************************************************************** @@ -291,5 +292,8 @@ ! ********************************************************************** DO ig=1,ngrid - IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN + IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN + !! mu0(ig)>mu0lim ensures that it is daytime + !! contains_mons=True ensures that there is a mount in the mesh and alpha_hmons>0 + !! ********************************************************************** !! Temperature profile above the mountain and in the close environment @@ -353,5 +357,5 @@ endif ENDDO - ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. alpha_hmons(ig) .gt. 0.) + ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. contains_mons(ig)) ENDDO ! DO ig=1,ngrid @@ -364,5 +368,5 @@ ! ********************************************************************** DO ig=1,ngrid - IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN + IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN !! Positive buoyancy: negative vertical velocity entrains UP IF (dt_top(ig) .gt. 0.) THEN @@ -477,5 +481,5 @@ DO ig=1,ngrid - IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN + IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN !! Total air mass within the PBL before entrainment (=> by PBL we mean between the surface and the layer where the vertical wind is maximum) masse_pbl(ig)=0. @@ -489,5 +493,5 @@ ! ********************************************************************** DO ig=1,ngrid - IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) .and. (masse_pbl(ig) .gt. 0.) ) THEN + IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) .and. (masse_pbl(ig) .gt. 0.) ) THEN !! Transport of background dust + concentrated topdust above lwmax DO l=lwmax(ig),nlayer @@ -583,5 +587,5 @@ DO l=1,nlayer!-1 !! Detrainment during the day - IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01)) THEN + IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01) .and. (contains_mons(ig)) ) THEN coefdetrain(ig,l)=1.*( rhobarz(ig,l+1)*abs(wfin(ig,l+1)) - rhobarz(ig,l)*abs(wfin(ig,l)) ) / masse(ig,l) !! Detrainment when abs(w(l)) > abs(w(l+1)), i.e. coefdetrain < 0 @@ -601,5 +605,5 @@ ! dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_dust_number(ig,l)/ptimestep endif - !! Full detrainment during the night imposed + !! Full detrainment imposed during the night or when topdust leaves its origin column (contains_mons=False) ELSE dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep @@ -681,4 +685,8 @@ CALL WRITEDIAGFI(ngrid,'wfin_top', & 'wfin_top','',3,wfin(:,:)) + CALL WRITEDIAGFI(ngrid,'hmons', & + 'hmons','',2,hmons) + CALL WRITEDIAGFI(ngrid,'hsummit', & + 'hsummit','',2,hsummit) CALL WRITEDIAGFI(ngrid,'alpha_hmons', & 'alpha_hmons','',2,alpha_hmons) @@ -961,26 +969,116 @@ - end subroutine van_leer - + end subroutine van_leer + !******************************************************************************** - ! initialization module variables - subroutine ini_topmons_mod(ngrid,nlayer) - - implicit none - - integer, intent(in) :: ngrid - integer, intent(in) :: nlayer - - allocate(alpha_hmons(ngrid)) - - end subroutine ini_topmons_mod - - subroutine end_topmons_mod - - implicit none - - if (allocated(alpha_hmons)) deallocate(alpha_hmons) - - end subroutine end_topmons_mod + subroutine topmons_setup(ngrid) + ! Purpose: + ! 1) Fill the logical array contains_mons(:), + ! with contains_mons(ig)=True if there is + ! a mountain in the mesh ig + ! 2) Compute alpha_hmons(:) and hsummit(:) + use surfdat_h,only: phisfi,hmons,summit,base,& + alpha_hmons,hsummit,contains_mons + use comcstfi_h,only: pi,g + use planetwide_mod,only: planetwide_maxval, planetwide_minval + use geometry_mod,only: longitude_deg,latitude_deg,& + boundslon,boundslat !boundslon/lat(ngrid,4) : + ! |------------------------------| + ! |north_west=2 north_east=1| + ! | | + ! | (ig) | + ! | | + ! |south_west=3 south_east=4| + ! |------------------------------| + + implicit none + integer,intent(in) :: ngrid + + ! Local variables + integer,parameter :: ntop_max = 19 ! total number of mounts written in the hmons list + real lon_top(ntop_max),lat_top(ntop_max) ! coordinates of the mounts (in deg) + ! Mountains list : + ! Olympus Mons,Ascraeus Mons,Elysium Mons,Arsia Mons,Pavonis Mons, + ! Hecates Tholus,Tharsis Tholus,Ceraunius Tholus,Alba Mons,Apollinaris Mons, + ! Albor Tholus,Biblis Tholus,Anseris Mons,Ulysses Tholus,Aeolis Mons, + ! Euripus Mons,Hadriacus Mons,Tyrrhenus Mons,Uranius Mons + ! + ! NB: in 64x48 horiz. resolution, Biblis Tholus & Ulysses Tholus fall in the + ! same mesh, hence only Biblis Tholus is kept by the alpha_hmons computation + data lon_top/-134.,-104.5,146.9,-121.1,-113.4,& + 150.2,-90.8,-97.4,-109.6,174.4,& + 150.4,-124.6,86.6,-121.6,137.8,& + 105.,91.8,106.5,-92.2/ + data lat_top/18.4,11.8,24.8,-8.4,-0.8,& + 31.8,13.4,24.,40.4,-9.3,& + 18.8,2.6,-29.8,2.9,-5.4,& + -44.8,-32.1,-21.1,26.8/ + integer,parameter :: ntop = 19 ! the topmons scheme is limited to the first ntop mounts + real :: boundslon_deg(4),boundslat_deg(4) + real :: hmax,hmin + integer :: ig,itop + + + + IF (ngrid.gt.1) THEN + ! Sanity check + if (ntop.gt.ntop_max) then + call abort_physic("topmons_setup","Number of mountains ntop greater than ntop_max",1) + endif + + ! Determine contains_mons + contains_mons(:)=.false. + + do ig=1,ngrid + boundslon_deg(:)=boundslon(ig,:)/pi*180. + boundslat_deg(:)=boundslat(ig,:)/pi*180. + + do itop=1,ntop + if ( (lon_top(itop).gt.boundslon_deg(2)).and.(lon_top(itop).lt.boundslon_deg(1)) & + .and.(lat_top(itop).gt.boundslat_deg(3)).and.(lat_top(itop).lt.boundslat_deg(2)) ) then + contains_mons(ig)=.true. + write(*,*) "topmons_setup: Found a mount at:" + write(*,*) "(",boundslon_deg(2),",",boundslat_deg(2),") (",boundslon_deg(1),",",boundslat_deg(1),")" + write(*,*) " ((",lon_top(itop),",",lat_top(itop),"))" + write(*,*) "(",boundslon_deg(3),",",boundslat_deg(3),") (",boundslon_deg(4),",",boundslat_deg(4),")" + endif + enddo + enddo + + ! Compute alpha_hmons + call planetwide_maxval(hmons,hmax) + call planetwide_minval(hmons,hmin) + do ig=1,ngrid + if (contains_mons(ig)) then + ! the mesh ig contains a mountain + alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin) + ! Sanity check + if (alpha_hmons(ig).le.0) then + call abort_physic("topmons_setup","ERROR: alpha_hmons cannot be <0 "//& + "if the mesh contains a mountain. Please check your "//& + "formula or your start files.",1) + endif + else + ! the mesh ig doesn't contain a mountain + alpha_hmons(ig)= 0 + endif + enddo + + ! Compute hsummit + hsummit(:)=summit(:)-phisfi(:)/g + + ELSE ! 1D case + hmin=0. + hmax=23162.1 !set here the height of the sub-grid scale topography + do ig=1,ngrid + alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin) + print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig) + enddo + + hsummit(:)=14000. + ENDIF ! (ngrid.gt.1) + + end subroutine topmons_setup +!******************************************************************************** END MODULE topmons_mod Index: /trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F (revision 2627) +++ /trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F (revision 2628) @@ -142,5 +142,5 @@ ! updating radius of topdust particles - IF (slpwind.AND.active) THEN + IF (topflows.AND.active) THEN DO l=1,nlayer DO ig=1, ngrid