! $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 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 lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet 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 :: 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/=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==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>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==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>=1) THEN newsig= sig ELSE IF (sig*preff/pa>=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>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))<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