! $Id: $ SUBROUTINE disvert_noterre ! Auteur : F. Forget Y. Wanherdrick, P. Levan ! Nouvelle version 100% Mars !! ! On l'utilise aussi pour Venus et Titan, legerment modifiee. USE iniprint_mod_h use IOIPSL USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt, & nivsig,nivsigs,pa,preff,scaleheight USE comconst_mod, ONLY: kappa USE logic_mod, ONLY: hybrid USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h IMPLICIT NONE ! !======================================================================= ! Discretisation verticale en coordonn�e hybride (ou sigma) ! !======================================================================= ! ! declarations: ! ------------- ! ! INTEGER :: l,ll REAL :: snorm REAL :: alpha,beta,gama,delta,deltaz real :: quoi,quand REAL :: zsig(llm),sig(llm+1) INTEGER :: np,ierr integer :: ierr1,ierr2,ierr3,ierr4 REAL :: x REAL :: SSUM EXTERNAL SSUM real :: newsig REAL :: dz0,dz1,nhaut,sig1,esig,csig,zz real :: tt,rr,gg, prevz real :: s(llm),dsig(llm) integer :: iz real :: z, ps,p character(len=*),parameter :: modname="disvert_noterre" ! !----------------------------------------------------------------------- ! ! Initializations: ! pi=2.*ASIN(1.) ! already done in iniconst hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates) CALL getin('hybrid',hybrid) write(lunout,*) trim(modname),': hybrid=',hybrid ! Ouverture possible de fichiers typiquement E.T. open(99,file="esasig.def",status='old',form='formatted', & iostat=ierr2) if(ierr2.ne.0) then close(99) open(99,file="z2sig.def",status='old',form='formatted', & iostat=ierr4) endif !----------------------------------------------------------------------- ! cas 1 on lit les options dans esasig.def: ! ---------------------------------------- IF(ierr2.eq.0) then ! Lecture de esasig.def : ! Systeme peu souple, mais qui respecte en theorie ! La conservation de l'energie (conversion Energie potentielle ! <-> energie cinetique, d'apres la note de Frederic Hourdin... write(lunout,*)'*****************************' write(lunout,*)'WARNING reading esasig.def' write(lunout,*)'*****************************' READ(99,*) scaleheight READ(99,*) dz0 READ(99,*) dz1 READ(99,*) nhaut CLOSE(99) dz0=dz0/scaleheight dz1=dz1/scaleheight sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut) esig=1. do l=1,20 esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.) enddo csig=(1./sig1-1.)/(exp(esig)-1.) DO L = 2, llm zz=csig*(exp(esig*(l-1.))-1.) sig(l) =1./(1.+zz) & * tanh(.5*(llm+1-l)/nhaut) ENDDO sig(1)=1. sig(llm+1)=0. quoi = 1. + 2.* kappa s( llm ) = 1. s(llm-1) = quoi IF( llm.gt.2 ) THEN DO ll = 2, llm-1 l = llm+1 - ll quand = sig(l+1)/ sig(l) s(l-1) = quoi * (1.-quand) * s(l) + quand * s(l+1) ENDDO END IF ! snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2) DO l = 1, llm s(l) = s(l)/ snorm ENDDO !----------------------------------------------------------------------- ! cas 2 on lit les options dans z2sig.def: ! ---------------------------------------- ELSE IF(ierr4.eq.0) then write(lunout,*)'****************************' write(lunout,*)'Reading z2sig.def' write(lunout,*)'****************************' READ(99,*) scaleheight do l=1,llm read(99,*) zsig(l) end do CLOSE(99) sig(1) =1 do l=2,llm sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + & exp(-zsig(l-1)/scaleheight) ) end do sig(llm+1) =0 !----------------------------------------------------------------------- ELSE write(lunout,*) 'didn t you forget something ??? ' write(lunout,*) 'We need file z2sig.def ! (OR esasig.def)' stop ENDIF !----------------------------------------------------------------------- DO l=1,llm nivsigs(l) = REAL(l) ENDDO DO l=1,llmp1 nivsig(l)= REAL(l) ENDDO !----------------------------------------------------------------------- ! .... Calculs de ap(l) et de bp(l) .... ! ......................................... ! ! ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... !----------------------------------------------------------------------- ! if (hybrid) then ! use hybrid coordinates write(lunout,*) "*********************************" write(lunout,*) "Using hybrid vertical coordinates" write(lunout,*) ! Coordonnees hybrides avec mod DO l = 1, llm call sig_hybrid(sig(l),pa,preff,newsig) bp(l) = EXP( 1. - 1./(newsig**2) ) ap(l) = pa * (newsig - bp(l) ) enddo bp(llmp1) = 0. ap(llmp1) = 0. else ! use sigma coordinates write(lunout,*) "********************************" write(lunout,*) "Using sigma vertical coordinates" write(lunout,*) ! Pour ne pas passer en coordonnees hybrides DO l = 1, llm ap(l) = 0. bp(l) = sig(l) ENDDO ap(llmp1) = 0. endif bp(llmp1) = 0. write(lunout,*) trim(modname),': BP ' write(lunout,*) bp write(lunout,*) trim(modname),': AP ' write(lunout,*) ap ! Calcul au milieu des couches : ! WARNING : le choix de placer le milieu des couches au niveau de ! pression interm�diaire est arbitraire et pourrait etre modifi�. ! Le calcul du niveau pour la derniere couche ! (on met la meme distance (en log pression) entre P(llm) ! et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est ! Specifique. Ce choix est sp�cifi� ici ET dans exner_milieu.F DO l = 1, llm-1 aps(l) = 0.5 *( ap(l) +ap(l+1)) bps(l) = 0.5 *( bp(l) +bp(l+1)) ENDDO if (hybrid) then aps(llm) = aps(llm-1)**2 / aps(llm-2) bps(llm) = 0.5*(bp(llm) + bp(llm+1)) else bps(llm) = bps(llm-1)**2 / bps(llm-2) aps(llm) = 0. end if write(lunout,*) trim(modname),': BPs ' write(lunout,*) bps write(lunout,*) trim(modname),': APs' write(lunout,*) aps DO l = 1, llm presnivs(l) = aps(l)+bps(l)*preff pseudoalt(l) = -scaleheight*log(presnivs(l)/preff) ENDDO write(lunout,*)trim(modname),' : PRESNIVS' write(lunout,*)presnivs write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', & 'height of ',scaleheight,' km)' write(lunout,*)pseudoalt ! -------------------------------------------------- ! This can be used to plot the vertical discretization ! (> xmgrace -nxy testhybrid.tab ) ! -------------------------------------------------- ! open (53,file='testhybrid.tab') ! scaleheight=15.5 ! do iz=0,34 ! z = -5 + min(iz,34-iz) ! approximation of scale height for Venus ! scaleheight = 15.5 - z/55.*10. ! ps = preff*exp(-z/scaleheight) ! zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff) ! do l=2,llm ! approximation of scale height for Venus ! if (zsig(l-1).le.55.) then ! scaleheight = 15.5 - zsig(l-1)/55.*10. ! else ! scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2. ! endif ! zsig(l)= zsig(l-1)-scaleheight* ! . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) ! end do ! write(53,'(I3,50F10.5)') iz, zsig ! end do ! close(53) ! -------------------------------------------------- RETURN END SUBROUTINE disvert_noterre ! ************************************************************ subroutine sig_hybrid(sig,pa,preff,newsig) ! ---------------------------------------------- ! Subroutine utilisee pour calculer des valeurs de sigma modifie ! pour conserver les coordonnees verticales decrites dans ! esasig.def/z2sig.def lors du passage en coordonnees hybrides ! F. Forget 2002 ! Connaissant sig (niveaux "sigma" ou on veut mettre les couches) ! L'objectif est de calculer newsig telle que ! (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig ! Cela ne se r�soud pas analytiquement: ! => on r�soud par iterration bourrine ! ---------------------------------------------- ! Information : where exp(1-1./x**2) become << x ! x exp(1-1./x**2) /x ! 1 1 ! 0.68 0.5 ! 0.5 1.E-1 ! 0.391 1.E-2 ! 0.333 1.E-3 ! 0.295 1.E-4 ! 0.269 1.E-5 ! 0.248 1.E-6 ! => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25 implicit none real :: x1, x2, sig,pa,preff, newsig, F integer :: j newsig = sig x1=0 x2=1 if (sig.ge.1) then newsig= sig else if (sig*preff/pa.ge.0.25) then DO J=1,9999 ! nombre d''iteration max F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig ! write(0,*) J, ' newsig =', newsig, ' F= ', F if (F.gt.1) then X2 = newsig newsig=(X1+newsig)*0.5 else X1 = newsig newsig=(X2+newsig)*0.5 end if ! Test : on arete lorsque on approxime sig � moins de 0.01 m pr�s ! (en pseudo altitude) : IF(abs(10.*log(F)).LT.1.E-5) goto 999 END DO else ! if (sig*preff/pa.le.0.25) then newsig= sig*preff/pa end if 999 continue Return END SUBROUTINE sig_hybrid