Index: trunk/libf/dyn3d/calfis.F =================================================================== --- trunk/libf/dyn3d/calfis.F (revision 101) +++ trunk/libf/dyn3d/calfis.F (revision 108) @@ -587,5 +587,4 @@ ! ADAPTATION GCM POUR CP(T) - ztfi=ztfi+zdtfi*dtphys call t2tpot(ngridmx*llm,ztfi,zteta,zpk) Index: trunk/libf/dyn3d/comconst.h =================================================================== --- trunk/libf/dyn3d/comconst.h (revision 101) +++ trunk/libf/dyn3d/comconst.h (revision 108) @@ -6,9 +6,9 @@ COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & - & iflag_top_bound + & iflag_top_bound,mode_top_bound COMMON/comconstr/dtvr,daysec, & & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & - & ,dissip_factz,dissip_deltaz,dissip_zref & - & ,tau_top_bound, & + & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta & + & ,dissip_pupstart ,tau_top_bound, & & daylen,year_day,molmass, ihf COMMON/cpdetvenus/nu_venus,t0_venus @@ -28,6 +28,7 @@ REAL g ! (m/s2) gravity REAL omeg ! (rad/s) rotation rate of the planet - REAL dissip_factz,dissip_deltaz,dissip_zref - INTEGER iflag_top_bound + REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta + REAL dissip_pupstart + INTEGER iflag_top_bound,mode_top_bound REAL tau_top_bound REAL daylen ! length of solar day, in 'standard' day length Index: trunk/libf/dyn3d/conf_gcm.F =================================================================== --- trunk/libf/dyn3d/conf_gcm.F (revision 101) +++ trunk/libf/dyn3d/conf_gcm.F (revision 108) @@ -303,17 +303,28 @@ ! Parametres controlant la variation sur la verticale des constantes de ! dissipation. -! Pour le moment actifs uniquement dans la version a 39 niveaux -! avec ok_strato=y - - dissip_factz=4. - dissip_deltaz=10. - dissip_zref=30. - CALL getin('dissip_factz',dissip_factz ) +! Actifs uniquement avec ok_strato=y + + dissip_fac_mid=2. + dissip_fac_up=10. + dissip_deltaz=10.! Intervalle (km) pour le changement mid / up + dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m) + dissip_pupstart=1.e3 ! pression (Pa) au bas la transition mid / up + CALL getin('dissip_fac_mid',dissip_fac_mid ) + CALL getin('dissip_fac_up',dissip_fac_up ) CALL getin('dissip_deltaz',dissip_deltaz ) - CALL getin('dissip_zref',dissip_zref ) - + CALL getin('dissip_hdelta',dissip_hdelta ) + CALL getin('dissip_pupstart',dissip_pupstart ) + +! Parametres controlant la couche eponge +! Actifs uniquement avec ok_strato=y +! mode = 0 : pas de sponge +! mode = 1 : u et v -> 0 +! mode = 2 : u et v -> moyenne zonale +! mode = 3 : u, v et h -> moyenne zonale iflag_top_bound=1 + mode_top_bound=3 tau_top_bound=1.e-5 CALL getin('iflag_top_bound',iflag_top_bound) + CALL getin('mode_top_bound',mode_top_bound) CALL getin('tau_top_bound',tau_top_bound) Index: trunk/libf/dyn3d/inidissip.F =================================================================== --- trunk/libf/dyn3d/inidissip.F (revision 101) +++ trunk/libf/dyn3d/inidissip.F (revision 108) @@ -31,5 +31,5 @@ INTEGER l,ij,idum,ii REAL tetamin - REAL pseudoz + REAL Pup REAL ran1 @@ -177,20 +177,38 @@ c -------------------------------------------------- - if (ok_strato .and. llm==39) then - do l=1,llm - pseudoz=8.*log(preff/presnivs(l)) - zvert(l)=1+ - s (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. - s *(dissip_factz-1.) - enddo - else - DO l=1,llm - zvert(l)=1. - ENDDO - fact=2. - DO l = 1, llm - zz = 1. - preff/presnivs(l) - zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) - ENDDO +c First step: going from 1 to dissip_fac_mid (in gcm.def) +c============ + DO l=1,llm + zz = 1. - preff/presnivs(l) + zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) + ENDDO + + write(*,*) 'Dissipation : ' + write(*,*) 'Multiplication de la dissipation en altitude :', + write(*,*) ' dissip_fac_mid =', dissip_fac_mid + +c Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) +c========================== +c Utilisation de la fonction tangente hyperbolique pour augmenter +c arbitrairement la dissipation et donc la stabilite du modele a +c partir d'une certaine altitude. + +c Le facteur multiplicatif de basse atmosphere etant deja pris +c en compte, il faut diviser le facteur multiplicatif de haute +c atmosphere par celui-ci. + + 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 =', + & dissip_pupstart,'Pa', dissip_deltaz , '(km)' + endif Index: trunk/libf/dyn3d/leapfrog.F =================================================================== --- trunk/libf/dyn3d/leapfrog.F (revision 101) +++ trunk/libf/dyn3d/leapfrog.F (revision 108) @@ -111,6 +111,4 @@ REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm) REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1) - - real :: duspg(ip1jmp1,llm) ! for bilan_dyn c variables pour le fichier histoire @@ -463,14 +461,11 @@ c ------------------- IF (ok_strato) THEN - CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi) + CALL top_bound( vcov,ucov,teta,masse,dutop,dvtop,dtetatop) +c dqtop=0, dptop=0 + CALL addfi( dtphys, leapf, forward , + $ ucov, vcov, teta , q ,ps , + $ dutop, dvtop, dtetatop , dqtop ,dptop ) ENDIF -c dqtop=0, dptop=0 - -c ajout des tendances physiques: -c ------------------------------ - CALL addfi( dtphys, leapf, forward , - $ ucov, vcov, teta , q ,ps , - $ dufi, dvfi, dtetafi , dqfi ,dpfi ) -c + c Diagnostique de conservation de l'énergie : difference IF (ip_ebil_dyn.ge.1 ) THEN @@ -512,10 +507,6 @@ ! Sponge layer (if any) IF (ok_strato) THEN - dutop(:,:)=0. - dvtop(:,:)=0. - dtetatop(:,:)=0. - dqtop(:,:,:)=0. - dptop(:)=0. CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop) +c dqtop=0, dptop=0 CALL addfi( dtvr, leapf, forward , $ ucov, vcov, teta , q ,ps , @@ -669,5 +660,5 @@ CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, - & du,dudis,duspg,dufi) + & du,dudis,dutop,dufi) #endif END IF @@ -801,5 +792,5 @@ CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, - & du,dudis,duspg,dufi) + & du,dudis,dutop,dufi) #endif ENDIF Index: trunk/libf/dyn3d/top_bound.F =================================================================== --- trunk/libf/dyn3d/top_bound.F (revision 101) +++ trunk/libf/dyn3d/top_bound.F (revision 108) @@ -52,4 +52,5 @@ LOGICAL,SAVE :: first=.true. + REAL zkm INTEGER j,l @@ -68,5 +69,5 @@ rdamp(llm-3)=tau_top_bound/8. else if (iflag_top_bound.eq.2) then -! couce eponge dans toutes les couches de pression plus faible que +! couche eponge dans toutes les couches de pression plus faible que ! 100 fois la pression de la derniere couche rdamp(:)=tau_top_bound @@ -74,10 +75,26 @@ endif first=.false. - print*,'TOP_BOUND rdamp=',rdamp + print*,'TOP_BOUND mode',mode_top_bound + print*,'Coeffs pour la couche eponge a l equateur' + print*,'p (Pa) z(km) tau (s) dt*rdamp' + do l=1,llm + if (rdamp(l).ne.0.) then + zkm = phi(iip1/2,jjp1/2,l)/(1000*g) + print*,presnivs(l),zkm, + . 1./rdamp(l), + . dt*rdamp(l) + endif + enddo endif CALL massbar(masse,massebx,masseby) - do l=1,llm +c mode = 0 : pas de sponge +c mode = 1 : u et v -> 0 +c mode = 2 : u et v -> moyenne zonale +c mode = 3 : u, v et h -> moyenne zonale + + if (mode_top_bound.ge.2) then + do l=1,llm do j=1,jjm vzon(j,l)=0. @@ -92,15 +109,7 @@ vzon(j,l)=vzon(j,l)/zm enddo - enddo + enddo - do l=1,llm - do i=1,iip1 - do j=1,jjm - dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l)) - enddo - enddo - enddo - - do l=1,llm + do l=1,llm do j=2,jjm uzon(j,l)=0. @@ -112,7 +121,12 @@ uzon(j,l)=uzon(j,l)/zm enddo - enddo + enddo + else + vzon(:,:)=0. + uzon(:,:)=0. + endif - do l=1,llm + if (mode_top_bound.ge.3) then + do l=1,llm do j=2,jjm zm=0. @@ -124,19 +138,37 @@ tzon(j,l)=tzon(j,l)/zm enddo - enddo + enddo + endif C AMORTISSEMENTS LINEAIRES: - do l=1,llm + if (mode_top_bound.ge.1) then + do l=1,llm + do i=1,iip1 + do j=1,jjm + dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l)) + enddo + enddo + enddo + + do l=1,llm do i=1,iip1 do j=2,jjm - du(i,j,l)=du(i,j,l) - s -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) - dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l)) + du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) enddo enddo - enddo + enddo + endif + + if (mode_top_bound.ge.3) then + do l=1,llm + do i=1,iip1 + do j=2,jjm + dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l)) + enddo + enddo + enddo + endif - RETURN END