Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F =================================================================== --- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F (revision 1019) +++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F (revision 1021) @@ -277,4 +277,6 @@ !Config dissip_period>0 => on prend cette valeur dissip_period = 0 + call getin('idissip',dissip_period) ! old Mars/Genreic model parameter + ! if there is a "dissip_period" in run.def, it overrides "idissip" CALL getin('dissip_period',dissip_period) @@ -601,5 +603,8 @@ $ "bad value for vert_prof_dissip") else - vert_prof_dissip=0 + vert_prof_dissip=0 ! default for planets ! + if (planet_type=="mars") then + vert_prof_dissip=1 ! use fac_mid & fac_up & startalt & delta + endif endif Index: trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90 (revision 1019) +++ trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90 (revision 1021) @@ -5,9 +5,22 @@ tetagdiv,tetagrot,tetatemp, vert_prof_dissip) !======================================================================= - ! initialisation de la dissipation horizontale + ! Initialization for horizontal (lateral) dissipation + ! - in all cases, there is a multiplicative coefficient which increases + ! the dissipation in the higher levels of the atmosphere, but there + ! are different ways of seting the vertical profile of this coefficient + ! (see code below). + ! - the call to dissipation, every 'dissip_period' dynamical time step, + ! can be imposed via 'dissip_period=...' in run.def (or via the + ! 'idissip=...' flag, but this flag should become obsolete, and is + ! overridden by the 'dissip_period' flag). Note that setting + ! dissip_period=0 has the special meaning of requesting an "optimal" + ! value for "dissip_period" (then taken as the largest possible value) + ! - the three characteristic time scales (relative to the smallest + ! longitudinal grid mesh size), which are privided in run.def, which + ! are used for the dissipations steps are: + ! tetagdiv : time scale for the gradient of the divergence of velocities + ! tetagrot : time scale for the curl of the curl of velocities + ! tetatemp : time scale for the laplacian of the potential temperature !======================================================================= - !----------------------------------------------------------------------- - ! declarations: - ! ------------- USE control_mod, only : dissip_period,iperiod,planet_type @@ -26,10 +39,12 @@ REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp - integer, INTENT(in):: vert_prof_dissip ! for the Earth model !! - ! Vertical profile of horizontal dissipation + integer, INTENT(in):: vert_prof_dissip ! Vertical profile of horizontal dissipation + ! For the Earth model: ! Allowed values: ! 0: rational fraction, function of pressure ! 1: tanh of altitude - + ! For planets: + ! 0: use fac_mid (read from run.def) + ! 1: use fac_mid, fac_up, startalt, delta (hard coded in inidissip.F90) ! Local variables: REAL fact,zvert(llm),zz @@ -46,5 +61,23 @@ REAL ran1 - + logical,save :: firstcall=.true. + real,save :: fac_mid,fac_up,startalt,delta,middle + + if (firstcall) then + firstcall=.false. + if ((planet_type.ne."earth").and.(vert_prof_dissip.eq.1)) then + ! initialize values for dissipation variation along the vertical (Mars) + fac_mid=3 ! coefficient for lower/middle atmosphere + fac_up=30 ! coefficient for upper atmosphere + startalt=70. ! altitude (in km) for the transition from middle to upper atm. + delta=30.! Size (in km) of the transition region between middle/upper atm. + + middle=startalt+delta/2 + write(lunout,*)"inidissip: multiplicative factors in altitude:", & + " fac_mid=",fac_mid," fac_up=",fac_up + write(lunout,*)" transition mid/up : startalt (km) =",startalt, & + " delta (km) =",delta + endif + endif !of if firstcall !----------------------------------------------------------------------- @@ -196,14 +229,15 @@ else ! other planets + if (vert_prof_dissip == 0) then ! First step: going from 1 to dissip_fac_mid (in gcm.def) !============ - DO l=1,llm + DO l=1,llm zz = 1. - preff/presnivs(l) zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) - ENDDO - - write(lunout,*) 'Dissipation : ' - write(lunout,*) 'Multiplication de la dissipation en altitude :' - write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid + ENDDO + + write(lunout,*) 'Dissipation : ' + write(lunout,*) 'Multiplication de la dissipation en altitude :' + write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) @@ -217,23 +251,38 @@ ! atmosphere par celui-ci. - if (ok_strato) then - - Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) - do l=1,llm + if (ok_strato) then + + Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) + do l=1,llm zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & *(1-(0.5*(1+tanh(-6./dissip_deltaz* & (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) - enddo - - write(*,*) ' dissip_fac_up =', dissip_fac_up - write(*,*) 'Transition mid /up: Pupstart,delta =', & + enddo + + write(*,*) ' dissip_fac_up =', dissip_fac_up + write(*,*) 'Transition mid /up: Pupstart,delta =', & dissip_pupstart,'Pa', dissip_deltaz , '(km)' - endif ! of if (ok_strato) - + endif ! of if (ok_strato) + elseif (vert_prof_dissip==1) then + DO l=1,llm + zz = 1. - preff/presnivs(l) + zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) + + zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)* & + (1-(0.5*(1+tanh(-6./ & + delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) & + )) + ENDDO + + else + write(lunout,*) 'wrong value for vert_prof_dissip:',vert_prof_dissip + abort_message='wrong value for vert_prof_dissip' + call abort_gcm('inidissip',abort_message,1) + endif ! of (vert_prof_dissip == 0) endif ! of if (planet_type.eq."earth") - write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' + write(lunout,*)'inidissip: Time constants for lateral dissipation' tetamin = 1.e+6 @@ -259,8 +308,28 @@ write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr + write(lunout,*)'pseudoZ(km) zvert dt(tetagdiv) dt(tetagrot) dt(divgrad)' DO l = 1,llm - write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & - dtdiss*tetah(l) - ENDDO + write(lunout,'(f6.1,x,4(1pe14.7,x))') & + pseudoalt(l),zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),dtdiss*tetah(l) + ! test if disipation is not too strong (for Explicit Euler time marching) + if (dtdiss*tetaudiv(l).gt.1.9) then + write(lunout,*)"STOP : lateral dissipation is too intense and will" + write(lunout,*)" generate instabilities in the model !" + write(lunout,*)" You must increase tetagdiv (or increase dissip_period" + write(lunout,*)" or increase day_stap)" + endif + if (dtdiss*tetaurot(l).gt.1.9) then + write(lunout,*)"STOP : lateral dissipation is too intense and will" + write(lunout,*)" generate instabilities in the model !" + write(lunout,*)" You must increase tetaurot (or increase dissip_period" + write(lunout,*)" or increase day_stap)" + endif + if (dtdiss*tetah(l).gt.1.9) then + write(lunout,*)"STOP : lateral dissipation is too intense and will" + write(lunout,*)" generate instabilities in the model !" + write(lunout,*)" You must increase tetah (or increase dissip_period" + write(lunout,*)" or increase day_stap)" + endif + ENDDO ! of DO l=1,llm END SUBROUTINE inidissip