Index: /trunk/chantiers/commit_importants.log =================================================================== --- /trunk/chantiers/commit_importants.log (revision 108) +++ /trunk/chantiers/commit_importants.log (revision 109) @@ -802,4 +802,38 @@ Ajout dissip_horiz.tex (et .pdf) dans la documentation. -Reste la discretisation verticale a voir. - +********************* +**** commit_v109 **** +********************* + +** Discretisation verticale +**----------------------- + +* Modif de dyn3d[par]/iniconst.F + Appel de disvert_[no]terre mis sous flag planet_type + +* Modif de dyn3d[par]/disvert.F90 + Renomme disvert_terre.F90 + +* Ajout de dyn3d[par]/disvert_noterre.F + Correspond au disvert pour Titan et Venus (et Mars ? A verifier) + Declaration aps et bps enlevee (=> comvert.h) + +* Modif de dyn3d[par]/comvert.h + Ajout de aps,bps + => entraine la modif de limy.F(+par), vlsplt[_p].F, vlspltqs[_p].F, + +* Modif de dyn3d[par]/limy.F, vlsplt[_p].F, vlspltqs[_p].F, + Changement de aps en apps (et du coup, apn en appn) + +* Ajout dyn3d[par]/exner_milieu.F (vient directement de mars/libf/dyn3d/) + pour remplacer exner_hyb dans cas noterre + +* Ajout dyn3dpar/exner_milieu_p.F + adapte de exner_milieu.F + pour remplacer exner_hyb_p dans cas noterre + +* Modif de dyn3d[par]/leapfrog[_p].F + Appels à exner_hyb[_p]/exner_milieu[_p] mis sous flag planet_type. + +* Ajout de disvert.tex[/.pdf] dans la doc + Index: /trunk/chantiers/meschantiers-Seb.txt =================================================================== --- /trunk/chantiers/meschantiers-Seb.txt (revision 108) +++ /trunk/chantiers/meschantiers-Seb.txt (revision 109) @@ -10,14 +10,11 @@ Reflechir au controle des frequences de sorties -- dynamique: introduire les modifs Titan/Venus dans dyn3d[par] - Quid de la discretisation verticale ? Comment rendre souple ? +- dynamique: + Il reste a voir tout ce qui est lie a l'etat initial + et au newstart - I/O: evolution des I/O de Venus et Titan vers Terre ? souplesse pour I/O Mars ? -- newstart... - - testphys1d... -- couche eponge... - Index: /trunk/documentation/disvert.tex =================================================================== --- /trunk/documentation/disvert.tex (revision 109) +++ /trunk/documentation/disvert.tex (revision 109) @@ -0,0 +1,98 @@ +\documentclass[a4paper,10pt]{article} +%\usepackage{graphicx} +\usepackage{natbib} % si appel à bibtex +%\usepackage[francais]{babel} +%\usepackage[latin1]{inputenc} % accents directs (é...), avec babel +%\usepackage{rotating} + +\setlength{\hoffset}{-1.in} +\setlength{\oddsidemargin}{3.cm} +\setlength{\textwidth}{15.cm} +\setlength{\marginparsep}{0.mm} +\setlength{\marginparwidth}{0.mm} + +\setlength{\voffset}{-1.in} +\setlength{\topmargin}{0.mm} +\setlength{\headheight}{0.mm} +\setlength{\headsep}{30.mm} +\setlength{\textheight}{24.cm} +\setlength{\footskip}{1.cm} + +\setlength{\parindent}{0.mm} +\setlength{\parskip}{1 em} +\newcommand{\ten}[1]{$\times 10^{#1}$~} +\renewcommand{\baselinestretch}{1.} + +\begin{document} +\pagestyle{plain} + +\begin{center} +{\bf \LARGE +Documentation for LMDZ, Planets version + +\vspace{1cm} +\Large +The vertical discretization +} + +\vspace{1cm} +S\'ebastien Lebonnois + +\vspace{1cm} +Latest version: \today +\end{center} + + +\section{Theoretical aspects} + +The position of the layers: +\begin{itemize} +\item pressure limit between two layers, +\item pressure within the layers +\end{itemize} + +The Exner function: definition. +It corresponds to the pressure levels within the layers. +Used for the computation of the potential temperature. +For the Earth, we use a specific scheme that computes these positions so that +it maintains a condition of proportionality between total, +internal and potential energy (cf. a note from F. Hourdin). + +\section{Pratical aspects in the code} + +\begin{itemize} +\item \textsf{disvert\_[no]terre.F[90]}: +position of the interface pressure levels from an input file +(several possibilities). +Definition of ap, bp and presnivs. +In the planetary version, definition of aps and bps. + +This is done only once, called at the beginning from \textsf{iniconst.F}. + +\item Interface pressures: +computed in \textsf{caldyn0.F, caldyn.F, integrd.F, leapfrog.F} +through the \textsf{pression.F} routine. + +\item Exner function (and therefore pressure within the layers): +computed at three different places in \textsf{leapfrog.F} through the +\textsf{exner\_[hyb/milieu].F} routine. +For the Earth, we use \textsf{exner\_hyb.F}, that computes the positions in a +specific way to maintain a condition of proportionality between total, +internal and potential energy (cf. a note from F. Hourdin). +For other planets, we use \textsf{exner\_milieu.F}, that computes the positions +of these pressure levels exactly in the middle of each layer. +Though this fails to maintain the previous condition, there is no evidence of +any significant influence on the results, and it makes it a lot easier to +define correctly the level positions with the input file. +\end{itemize} + +%\begin{thebibliography}{2} +%\providecommand{\natexlab}[1]{#1} +%\expandafter\ifx\csname urlstyle\endcsname\relax +% \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else +% \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup +% \urlstyle{rm}\Url}\fi + +%\end{thebibliography} + +\end{document} Index: /trunk/libf/dyn3d/comvert.h =================================================================== --- /trunk/libf/dyn3d/comvert.h (revision 108) +++ /trunk/libf/dyn3d/comvert.h (revision 109) @@ -6,7 +6,8 @@ COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & - & pa,preff,nivsigs(llm),nivsig(llm+1) + & pa,preff,nivsigs(llm),nivsig(llm+1), & + & aps(llm),bps(llm) - REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig + REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps !----------------------------------------------------------------------- Index: unk/libf/dyn3d/disvert.F90 =================================================================== --- /trunk/libf/dyn3d/disvert.F90 (revision 108) +++ (revision ) @@ -1,142 +1,0 @@ -! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ - -SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) - - ! Auteur : P. Le Van - - IMPLICIT NONE - - include "dimensions.h" - include "paramet.h" - include "iniprint.h" - include "logic.h" - - ! s = sigma ** kappa : coordonnee verticale - ! dsig(l) : epaisseur de la couche l ds la coord. s - ! sig(l) : sigma a l'interface des couches l et l-1 - ! ds(l) : distance entre les couches l et l-1 en coord.s - - REAL pa, preff - REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) - REAL presnivs(llm) - - REAL sig(llm+1), dsig(llm) - real zk, zkm1, dzk1, dzk2, k0, k1 - - INTEGER l - REAL dsigmin - REAL alpha, beta, deltaz, h - INTEGER iostat - REAL pi, x - - !----------------------------------------------------------------------- - - pi = 2 * ASIN(1.) - - OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) - - IF (iostat == 0) THEN - ! cas 1 on lit les options dans sigma.def: - READ(99, *) h ! hauteur d'echelle 8. - READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 - READ(99, *) beta ! facteur d'acroissement en haut 1.3 - READ(99, *) k0 ! nombre de couches dans la transition surf - READ(99, *) k1 ! nombre de couches dans la transition haute - CLOSE(99) - alpha=deltaz/(llm*h) - write(lunout, *)'disvert: h, alpha, k0, k1, beta' - - alpha=deltaz/tanh(1./k0)*2. - zkm1=0. - sig(1)=1. - do l=1, llm - sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & - *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) - zk=-h*log(sig(l+1)) - - dzk1=alpha*tanh(l/k0) - dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) - write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 - zkm1=zk - enddo - - sig(llm+1)=0. - - DO l = 1, llm - dsig(l) = sig(l)-sig(l+1) - end DO - ELSE - if (ok_strato) then - if (llm==39) then - dsigmin=0.3 - else if (llm==50) then - dsigmin=1. - else - WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster' - dsigmin=1. - endif - WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin - endif - - DO l = 1, llm - x = 2*asin(1.) * (l - 0.5) / (llm + 1) - - IF (ok_strato) THEN - dsig(l) =(dsigmin + 7. * SIN(x)**2) & - *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 - ELSE - dsig(l) = 1.0 + 7.0 * SIN(x)**2 - ENDIF - ENDDO - dsig = dsig / sum(dsig) - sig(llm+1) = 0. - DO l = llm, 1, -1 - sig(l) = sig(l+1) + dsig(l) - ENDDO - 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 lectba ..... - - bp(llmp1) = 0. - - DO l = 1, llm - bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) - ap(l) = pa * ( sig(l) - bp(l) ) - ENDDO - - bp(1)=1. - ap(1)=0. - - ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) - - write(lunout, *)'disvert: BP ' - write(lunout, *) bp - write(lunout, *)'disvert: AP ' - write(lunout, *) ap - - write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' - write(lunout, *)'couches calcules pour une pression de surface =', preff - write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' - write(lunout, *)'8km' - DO l = 1, llm - dpres(l) = bp(l) - bp(l+1) - presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) - write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & - log(preff/presnivs(l))*8. & - , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & - max(ap(l+1)+bp(l+1)*preff, 1.e-10)) - ENDDO - - write(lunout, *) 'disvert: PRESNIVS ' - write(lunout, *) presnivs - -END SUBROUTINE disvert Index: /trunk/libf/dyn3d/disvert_noterre.F =================================================================== --- /trunk/libf/dyn3d/disvert_noterre.F (revision 109) +++ /trunk/libf/dyn3d/disvert_noterre.F (revision 109) @@ -0,0 +1,315 @@ + SUBROUTINE disvert_noterre + +c Auteur : F. Forget Y. Wanherdrick, P. Levan +c Nouvelle version 100% Mars !! +c On l'utilise aussi pour Venus et Titan, legerment modifiee. + + IMPLICIT NONE + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "comconst.h" +#include "logic.h" +c +c======================================================================= +c Discretisation verticale en coordonnée hybride +c +c======================================================================= +c +c declarations: +c ------------- +c +c + INTEGER l,ll + REAL snorm + REAL alpha,beta,gama,delta,deltaz,h,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) + real pseudoalt(llm) + + integer iz + real z, ps,p + +c +c----------------------------------------------------------------------- +c + pi=2.*ASIN(1.) + +! Ouverture possible de fichiers typiquement E.T. + + open(99,file="esasig.def",status='old',form='formatted', + s iostat=ierr2) + if(ierr2.ne.0) then + close(99) + open(99,file="z2sig.def",status='old',form='formatted', + s iostat=ierr4) + endif + +c----------------------------------------------------------------------- +c cas 1 on lit les options dans esasig.def: +c ---------------------------------------- + + IF(ierr2.eq.0) then + +c Lecture de esasig.def : +c Systeme peu souple, mais qui respecte en theorie +c La conservation de l'energie (conversion Energie potentielle +c <-> energie cinetique, d'apres la note de Frederic Hourdin... + + PRINT*,'*****************************' + PRINT*,'WARNING Lecture de esasig.def' + PRINT*,'*****************************' + READ(99,*) h + READ(99,*) dz0 + READ(99,*) dz1 + READ(99,*) nhaut + CLOSE(99) + + dz0=dz0/h + dz1=dz1/h + + 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 +c + 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 + +c----------------------------------------------------------------------- +c cas 2 on lit les options dans z2sig.def: +c ---------------------------------------- + + ELSE IF(ierr4.eq.0) then + PRINT*,'****************************' + PRINT*,'Lecture de z2sig.def' + PRINT*,'****************************' + + READ(99,*) h + 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)/h) + exp(-zsig(l-1)/h) ) + end do + sig(llm+1) =0 + +c----------------------------------------------------------------------- + ELSE + write(*,*) 'didn t you forget something ??? ' + write(*,*) 'We need the file z2sig.def ! (OR esasig.def) ' + stop + ENDIF +c----------------------------------------------------------------------- + + DO l=1,llm + nivsigs(l) = FLOAT(l) + ENDDO + + DO l=1,llmp1 + nivsig(l)= FLOAT(l) + ENDDO + + +c----------------------------------------------------------------------- +c .... Calculs de ap(l) et de bp(l) .... +c ......................................... +c +c ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... +c----------------------------------------------------------------------- +c +c se trouvait dans Titan... Verifier... +c h = 40. + + if (1.eq.1) then ! toujours hybrides + write(*,*) "*******************************" + write(*,*) "Systeme en coordonnees hybrides" + write(*,*) +c 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 + write(*,*) "****************************" + write(*,*) "Systeme en coordonnees sigma" + write(*,*) +c 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. + + PRINT *,' BP ' + PRINT *, bp + PRINT *,' AP ' + PRINT *, ap + +c Calcul au milieu des couches : +c WARNING : le choix de placer le milieu des couches au niveau de +c pression intermédiaire est arbitraire et pourrait etre modifié. +c Le calcul du niveau pour la derniere couche +c (on met la meme distance (en log pression) entre P(llm) +c et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est +c Specifique. Ce choix est spécifié ici ET dans exner_hyb.F + + DO l = 1, llm + 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 + + PRINT *,' BPs ' + PRINT *, bps + PRINT *,' APs' + PRINT *, aps + + DO l = 1, llm + presnivs(l) = aps(l)+bps(l)*preff + pseudoalt(l) = -h*log(presnivs(l)/preff) + ENDDO + + PRINT *,' PRESNIVS' + PRINT *,presnivs + PRINT *,'Pseudo altitude des Presnivs : (avec H = ',h,' km)' + PRINT *,pseudoalt + +c -------------------------------------------------- +c This can be used to plot the vertical discretization +c (> xmgrace -nxy testhybrid.tab ) +c -------------------------------------------------- +c open (53,file='testhybrid.tab') +c h=15.5 +c do iz=0,34 +c z = -5 + min(iz,34-iz) +c approximation of scale height for Venus +c h = 15.5 - z/55.*10. +c ps = preff*exp(-z/h) +c zsig(1)= -h*log((aps(1) + bps(1)*ps)/preff) +c do l=2,llm +c approximation of scale height for Venus +c if (zsig(l-1).le.55.) then +c h = 15.5 - zsig(l-1)/55.*10. +c else +c h = 5.5 - (zsig(l-1)-55.)/35.*2. +c endif +c zsig(l)= zsig(l-1)-h* +c . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) +c end do +c write(53,'(I3,50F10.5)') iz, zsig +c end do +c close(53) +c -------------------------------------------------- + + + RETURN + END + +c ************************************************************ + subroutine sig_hybrid(sig,pa,preff,newsig) +c ---------------------------------------------- +c Subroutine utilisee pour calculer des valeurs de sigma modifie +c pour conserver les coordonnees verticales decrites dans +c esasig.def/z2sig.def lors du passage en coordonnees hybrides +c F. Forget 2002 +c Connaissant sig (niveaux "sigma" ou on veut mettre les couches) +c L'objectif est de calculer newsig telle que +c (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig +c Cela ne se résoud pas analytiquement: +c => on résoud par iterration bourrine +c ---------------------------------------------- +c Information : where exp(1-1./x**2) become << x +c x exp(1-1./x**2) /x +c 1 1 +c 0.68 0.5 +c 0.5 1.E-1 +c 0.391 1.E-2 +c 0.333 1.E-3 +c 0.295 1.E-4 +c 0.269 1.E-5 +c 0.248 1.E-6 +c => 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 +c 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 +c Test : on arete lorsque on approxime sig à moins de 0.01 m près +c (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 Index: /trunk/libf/dyn3d/disvert_terre.F90 =================================================================== --- /trunk/libf/dyn3d/disvert_terre.F90 (revision 109) +++ /trunk/libf/dyn3d/disvert_terre.F90 (revision 109) @@ -0,0 +1,142 @@ +! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ + +SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) + + ! Auteur : P. Le Van + + IMPLICIT NONE + + include "dimensions.h" + include "paramet.h" + include "iniprint.h" + include "logic.h" + + ! s = sigma ** kappa : coordonnee verticale + ! dsig(l) : epaisseur de la couche l ds la coord. s + ! sig(l) : sigma a l'interface des couches l et l-1 + ! ds(l) : distance entre les couches l et l-1 en coord.s + + REAL pa, preff + REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) + REAL presnivs(llm) + + REAL sig(llm+1), dsig(llm) + real zk, zkm1, dzk1, dzk2, k0, k1 + + INTEGER l + REAL dsigmin + REAL alpha, beta, deltaz, h + INTEGER iostat + REAL pi, x + + !----------------------------------------------------------------------- + + pi = 2 * ASIN(1.) + + OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) + + IF (iostat == 0) THEN + ! cas 1 on lit les options dans sigma.def: + READ(99, *) h ! hauteur d'echelle 8. + READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 + READ(99, *) beta ! facteur d'acroissement en haut 1.3 + READ(99, *) k0 ! nombre de couches dans la transition surf + READ(99, *) k1 ! nombre de couches dans la transition haute + CLOSE(99) + alpha=deltaz/(llm*h) + write(lunout, *)'disvert: h, alpha, k0, k1, beta' + + alpha=deltaz/tanh(1./k0)*2. + zkm1=0. + sig(1)=1. + do l=1, llm + sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) + zk=-h*log(sig(l+1)) + + dzk1=alpha*tanh(l/k0) + dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) + write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 + zkm1=zk + enddo + + sig(llm+1)=0. + + DO l = 1, llm + dsig(l) = sig(l)-sig(l+1) + end DO + ELSE + if (ok_strato) then + if (llm==39) then + dsigmin=0.3 + else if (llm==50) then + dsigmin=1. + else + WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster' + dsigmin=1. + endif + WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin + endif + + DO l = 1, llm + x = 2*asin(1.) * (l - 0.5) / (llm + 1) + + IF (ok_strato) THEN + dsig(l) =(dsigmin + 7. * SIN(x)**2) & + *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 + ELSE + dsig(l) = 1.0 + 7.0 * SIN(x)**2 + ENDIF + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + 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 lectba ..... + + bp(llmp1) = 0. + + DO l = 1, llm + bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) + ap(l) = pa * ( sig(l) - bp(l) ) + ENDDO + + bp(1)=1. + ap(1)=0. + + ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) + + write(lunout, *)'disvert: BP ' + write(lunout, *) bp + write(lunout, *)'disvert: AP ' + write(lunout, *) ap + + write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' + write(lunout, *)'couches calcules pour une pression de surface =', preff + write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' + write(lunout, *)'8km' + DO l = 1, llm + dpres(l) = bp(l) - bp(l+1) + presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) + write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & + log(preff/presnivs(l))*8. & + , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & + max(ap(l+1)+bp(l+1)*preff, 1.e-10)) + ENDDO + + write(lunout, *) 'disvert: PRESNIVS ' + write(lunout, *) presnivs + +END SUBROUTINE disvert Index: /trunk/libf/dyn3d/exner_milieu.F =================================================================== --- /trunk/libf/dyn3d/exner_milieu.F (revision 109) +++ /trunk/libf/dyn3d/exner_milieu.F (revision 109) @@ -0,0 +1,102 @@ + SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) +c +c Auteurs : F. Forget , Y. Wanherdrick +c P.Le Van , Fr. Hourdin . +c .......... +c +c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... +c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... +c +c ************************************************************************ +c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des +c couches . Pk(l) sera calcule aux milieux des couches l ,entre les +c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . +c ************************************************************************ +c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont +c la pression et la fonction d'Exner au sol . +c +c WARNING : CECI est une version speciale de exner_hyb originale +c Utilise dans la version martienne pour pouvoir +c tourner avec des coordonnees verticales complexe +c => Il ne verifie PAS la condition la proportionalite en +c energie totale/ interne / potentielle (F.Forget 2001) +c ( voir note de Fr.Hourdin ) , +c + IMPLICIT NONE +c +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comgeom.h" +#include "comvert.h" +#include "serre.h" + + INTEGER ngrid + REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm) + REAL ps(ngrid),pks(ngrid), beta(ngrid,llm) + +c .... variables locales ... + + INTEGER l, ij + REAL dum1 + + REAL ppn(iim),pps(iim) + REAL xpn, xps + REAL SSUM + EXTERNAL filtreg, SSUM + +c ------------- +c Calcul de pks +c ------------- + + DO ij = 1, ngrid + pks(ij) = cpp * ( ps(ij)/preff ) ** kappa + ENDDO + + DO ij = 1, iim + ppn(ij) = aire( ij ) * pks( ij ) + pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) + ENDDO + xpn = SSUM(iim,ppn,1) /apoln + xps = SSUM(iim,pps,1) /apols + + DO ij = 1, iip1 + pks( ij ) = xpn + pks( ij+ip1jm ) = xps + ENDDO +c +c +c .... Calcul de pk pour la couche l +c -------------------------------------------- +c + dum1 = cpp * (2*preff)**(-kappa) + DO l = 1, llm-1 + DO ij = 1, ngrid + pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa + ENDDO + ENDDO + +c .... Calcul de pk pour la couche l = llm .. +c (on met la meme distance (en log pression) entre Pk(llm) +c et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) + + DO ij = 1, ngrid + pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2) + ENDDO + + +c calcul de pkf +c ------------- + CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) + CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) + +c EST-CE UTILE ?? : calcul de beta +c -------------------------------- + DO l = 2, llm + DO ij = 1, ngrid + beta(ij,l) = pk(ij,l) / pk(ij,l-1) + ENDDO + ENDDO + + RETURN + END Index: /trunk/libf/dyn3d/iniconst.F =================================================================== --- /trunk/libf/dyn3d/iniconst.F (revision 108) +++ /trunk/libf/dyn3d/iniconst.F (revision 109) @@ -53,6 +53,9 @@ c----------------------------------------------------------------------- - CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) -c + if (planet_type.eq."earth") then + CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) + else + CALL disvert_noterre + endif c RETURN Index: /trunk/libf/dyn3d/leapfrog.F =================================================================== --- /trunk/libf/dyn3d/leapfrog.F (revision 108) +++ /trunk/libf/dyn3d/leapfrog.F (revision 109) @@ -238,5 +238,10 @@ dq(:,:,:)=0. CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif + c------------------ c TEST PK MONOTONE @@ -404,5 +409,9 @@ CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif ! rdaym_ini = itau * dtvr / daysec @@ -519,5 +528,9 @@ CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif Index: /trunk/libf/dyn3d/limy.F =================================================================== --- /trunk/libf/dyn3d/limy.F (revision 108) +++ /trunk/libf/dyn3d/limy.F (revision 109) @@ -40,5 +40,5 @@ REAL qbyv(ip1jm,llm) - REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2 + REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2 Logical extremum,first save first @@ -117,14 +117,14 @@ c print*,dyqv(iip1+1) -c apn=abs(dyq(1)/dyqv(iip1+1)) +c appn=abs(dyq(1)/dyqv(iip1+1)) c print*,dyq(ip1jm+1) c print*,dyqv(ip1jm-iip1+1) -c aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +c apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) c do ij=2,iim -c apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -c aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +c appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +c apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) c enddo -c apn=min(pente_max/apn,1.) -c aps=min(pente_max/aps,1.) +c appn=min(pente_max/appn,1.) +c apps=min(pente_max/apps,1.) @@ -132,13 +132,13 @@ c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -c & apn=0. +c & appn=0. c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -c & aps=0. +c & apps=0. c limitation des pentes aux poles c do ij=1,iip1 -c dyq(ij)=apn*dyq(ij) -c dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +c dyq(ij)=appn*dyq(ij) +c dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) c enddo Index: /trunk/libf/dyn3d/vlsplt.F =================================================================== --- /trunk/libf/dyn3d/vlsplt.F (revision 108) +++ /trunk/libf/dyn3d/vlsplt.F (revision 109) @@ -478,5 +478,5 @@ REAL qbyv(ip1jm,llm) - REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs + REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs c REAL newq,oldmasse Logical extremum,first,testcpu @@ -602,14 +602,14 @@ C PRINT*,dyq(1) C PRINT*,dyqv(iip1+1) -C apn=abs(dyq(1)/dyqv(iip1+1)) +C appn=abs(dyq(1)/dyqv(iip1+1)) C PRINT*,dyq(ip1jm+1) C PRINT*,dyqv(ip1jm-iip1+1) -C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) C DO ij=2,iim -C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) C ENDDO -C apn=min(pente_max/apn,1.) -C aps=min(pente_max/aps,1.) +C appn=min(pente_max/appn,1.) +C apps=min(pente_max/apps,1.) C C @@ -617,13 +617,13 @@ C C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -C & apn=0. +C & appn=0. C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -C & aps=0. +C & apps=0. C C limitation des pentes aux poles C DO ij=1,iip1 -C dyq(ij)=apn*dyq(ij) -C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +C dyq(ij)=appn*dyq(ij) +C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) C ENDDO C Index: /trunk/libf/dyn3d/vlspltqs.F =================================================================== --- /trunk/libf/dyn3d/vlspltqs.F (revision 108) +++ /trunk/libf/dyn3d/vlspltqs.F (revision 109) @@ -635,14 +635,14 @@ C PRINT*,dyq(1) C PRINT*,dyqv(iip1+1) -C apn=abs(dyq(1)/dyqv(iip1+1)) +C appn=abs(dyq(1)/dyqv(iip1+1)) C PRINT*,dyq(ip1jm+1) C PRINT*,dyqv(ip1jm-iip1+1) -C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) C DO ij=2,iim -C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) C ENDDO -C apn=min(pente_max/apn,1.) -C aps=min(pente_max/aps,1.) +C appn=min(pente_max/appn,1.) +C apps=min(pente_max/apps,1.) C C @@ -650,13 +650,13 @@ C C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -C & apn=0. +C & appn=0. C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -C & aps=0. +C & apps=0. C C limitation des pentes aux poles C DO ij=1,iip1 -C dyq(ij)=apn*dyq(ij) -C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +C dyq(ij)=appn*dyq(ij) +C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) C ENDDO C Index: /trunk/libf/dyn3dpar/comvert.h =================================================================== --- /trunk/libf/dyn3dpar/comvert.h (revision 108) +++ /trunk/libf/dyn3dpar/comvert.h (revision 109) @@ -1,12 +1,13 @@ ! -! $Header$ +! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $ ! !----------------------------------------------------------------------- ! INCLUDE 'comvert.h' - COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & - & pa,preff,nivsigs(llm),nivsig(llm+1) + COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & + & pa,preff,nivsigs(llm),nivsig(llm+1), & + & aps(llm),bps(llm) - REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig + REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps -!----------------------------------------------------------------------- + !----------------------------------------------------------------------- Index: unk/libf/dyn3dpar/disvert.F90 =================================================================== --- /trunk/libf/dyn3dpar/disvert.F90 (revision 108) +++ (revision ) @@ -1,142 +1,0 @@ -! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ - -SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) - - ! Auteur : P. Le Van - - IMPLICIT NONE - - include "dimensions.h" - include "paramet.h" - include "iniprint.h" - include "logic.h" - - ! s = sigma ** kappa : coordonnee verticale - ! dsig(l) : epaisseur de la couche l ds la coord. s - ! sig(l) : sigma a l'interface des couches l et l-1 - ! ds(l) : distance entre les couches l et l-1 en coord.s - - REAL pa, preff - REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) - REAL presnivs(llm) - - REAL sig(llm+1), dsig(llm) - real zk, zkm1, dzk1, dzk2, k0, k1 - - INTEGER l - REAL dsigmin - REAL alpha, beta, deltaz, h - INTEGER iostat - REAL pi, x - - !----------------------------------------------------------------------- - - pi = 2 * ASIN(1.) - - OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) - - IF (iostat == 0) THEN - ! cas 1 on lit les options dans sigma.def: - READ(99, *) h ! hauteur d'echelle 8. - READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 - READ(99, *) beta ! facteur d'acroissement en haut 1.3 - READ(99, *) k0 ! nombre de couches dans la transition surf - READ(99, *) k1 ! nombre de couches dans la transition haute - CLOSE(99) - alpha=deltaz/(llm*h) - write(lunout, *)'disvert: h, alpha, k0, k1, beta' - - alpha=deltaz/tanh(1./k0)*2. - zkm1=0. - sig(1)=1. - do l=1, llm - sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & - *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) - zk=-h*log(sig(l+1)) - - dzk1=alpha*tanh(l/k0) - dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) - write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 - zkm1=zk - enddo - - sig(llm+1)=0. - - DO l = 1, llm - dsig(l) = sig(l)-sig(l+1) - end DO - ELSE - if (ok_strato) then - if (llm==39) then - dsigmin=0.3 - else if (llm==50) then - dsigmin=1. - else - WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster' - dsigmin=1. - endif - WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin - endif - - DO l = 1, llm - x = 2*asin(1.) * (l - 0.5) / (llm + 1) - - IF (ok_strato) THEN - dsig(l) =(dsigmin + 7. * SIN(x)**2) & - *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 - ELSE - dsig(l) = 1.0 + 7.0 * SIN(x)**2 - ENDIF - ENDDO - dsig = dsig / sum(dsig) - sig(llm+1) = 0. - DO l = llm, 1, -1 - sig(l) = sig(l+1) + dsig(l) - ENDDO - 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 lectba ..... - - bp(llmp1) = 0. - - DO l = 1, llm - bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) - ap(l) = pa * ( sig(l) - bp(l) ) - ENDDO - - bp(1)=1. - ap(1)=0. - - ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) - - write(lunout, *)'disvert: BP ' - write(lunout, *) bp - write(lunout, *)'disvert: AP ' - write(lunout, *) ap - - write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' - write(lunout, *)'couches calcules pour une pression de surface =', preff - write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' - write(lunout, *)'8km' - DO l = 1, llm - dpres(l) = bp(l) - bp(l+1) - presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) - write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & - log(preff/presnivs(l))*8. & - , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & - max(ap(l+1)+bp(l+1)*preff, 1.e-10)) - ENDDO - - write(lunout, *) 'disvert: PRESNIVS ' - write(lunout, *) presnivs - -END SUBROUTINE disvert Index: /trunk/libf/dyn3dpar/disvert_noterre.F =================================================================== --- /trunk/libf/dyn3dpar/disvert_noterre.F (revision 109) +++ /trunk/libf/dyn3dpar/disvert_noterre.F (revision 109) @@ -0,0 +1,315 @@ + SUBROUTINE disvert_noterre + +c Auteur : F. Forget Y. Wanherdrick, P. Levan +c Nouvelle version 100% Mars !! +c On l'utilise aussi pour Venus et Titan, legerment modifiee. + + IMPLICIT NONE + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "comconst.h" +#include "logic.h" +c +c======================================================================= +c Discretisation verticale en coordonnée hybride +c +c======================================================================= +c +c declarations: +c ------------- +c +c + INTEGER l,ll + REAL snorm + REAL alpha,beta,gama,delta,deltaz,h,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) + real pseudoalt(llm) + + integer iz + real z, ps,p + +c +c----------------------------------------------------------------------- +c + pi=2.*ASIN(1.) + +! Ouverture possible de fichiers typiquement E.T. + + open(99,file="esasig.def",status='old',form='formatted', + s iostat=ierr2) + if(ierr2.ne.0) then + close(99) + open(99,file="z2sig.def",status='old',form='formatted', + s iostat=ierr4) + endif + +c----------------------------------------------------------------------- +c cas 1 on lit les options dans esasig.def: +c ---------------------------------------- + + IF(ierr2.eq.0) then + +c Lecture de esasig.def : +c Systeme peu souple, mais qui respecte en theorie +c La conservation de l'energie (conversion Energie potentielle +c <-> energie cinetique, d'apres la note de Frederic Hourdin... + + PRINT*,'*****************************' + PRINT*,'WARNING Lecture de esasig.def' + PRINT*,'*****************************' + READ(99,*) h + READ(99,*) dz0 + READ(99,*) dz1 + READ(99,*) nhaut + CLOSE(99) + + dz0=dz0/h + dz1=dz1/h + + 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 +c + 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 + +c----------------------------------------------------------------------- +c cas 2 on lit les options dans z2sig.def: +c ---------------------------------------- + + ELSE IF(ierr4.eq.0) then + PRINT*,'****************************' + PRINT*,'Lecture de z2sig.def' + PRINT*,'****************************' + + READ(99,*) h + 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)/h) + exp(-zsig(l-1)/h) ) + end do + sig(llm+1) =0 + +c----------------------------------------------------------------------- + ELSE + write(*,*) 'didn t you forget something ??? ' + write(*,*) 'We need the file z2sig.def ! (OR esasig.def) ' + stop + ENDIF +c----------------------------------------------------------------------- + + DO l=1,llm + nivsigs(l) = FLOAT(l) + ENDDO + + DO l=1,llmp1 + nivsig(l)= FLOAT(l) + ENDDO + + +c----------------------------------------------------------------------- +c .... Calculs de ap(l) et de bp(l) .... +c ......................................... +c +c ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... +c----------------------------------------------------------------------- +c +c se trouvait dans Titan... Verifier... +c h = 40. + + if (1.eq.1) then ! toujours hybrides + write(*,*) "*******************************" + write(*,*) "Systeme en coordonnees hybrides" + write(*,*) +c 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 + write(*,*) "****************************" + write(*,*) "Systeme en coordonnees sigma" + write(*,*) +c 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. + + PRINT *,' BP ' + PRINT *, bp + PRINT *,' AP ' + PRINT *, ap + +c Calcul au milieu des couches : +c WARNING : le choix de placer le milieu des couches au niveau de +c pression intermédiaire est arbitraire et pourrait etre modifié. +c Le calcul du niveau pour la derniere couche +c (on met la meme distance (en log pression) entre P(llm) +c et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est +c Specifique. Ce choix est spécifié ici ET dans exner_hyb.F + + DO l = 1, llm + 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 + + PRINT *,' BPs ' + PRINT *, bps + PRINT *,' APs' + PRINT *, aps + + DO l = 1, llm + presnivs(l) = aps(l)+bps(l)*preff + pseudoalt(l) = -h*log(presnivs(l)/preff) + ENDDO + + PRINT *,' PRESNIVS' + PRINT *,presnivs + PRINT *,'Pseudo altitude des Presnivs : (avec H = ',h,' km)' + PRINT *,pseudoalt + +c -------------------------------------------------- +c This can be used to plot the vertical discretization +c (> xmgrace -nxy testhybrid.tab ) +c -------------------------------------------------- +c open (53,file='testhybrid.tab') +c h=15.5 +c do iz=0,34 +c z = -5 + min(iz,34-iz) +c approximation of scale height for Venus +c h = 15.5 - z/55.*10. +c ps = preff*exp(-z/h) +c zsig(1)= -h*log((aps(1) + bps(1)*ps)/preff) +c do l=2,llm +c approximation of scale height for Venus +c if (zsig(l-1).le.55.) then +c h = 15.5 - zsig(l-1)/55.*10. +c else +c h = 5.5 - (zsig(l-1)-55.)/35.*2. +c endif +c zsig(l)= zsig(l-1)-h* +c . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) +c end do +c write(53,'(I3,50F10.5)') iz, zsig +c end do +c close(53) +c -------------------------------------------------- + + + RETURN + END + +c ************************************************************ + subroutine sig_hybrid(sig,pa,preff,newsig) +c ---------------------------------------------- +c Subroutine utilisee pour calculer des valeurs de sigma modifie +c pour conserver les coordonnees verticales decrites dans +c esasig.def/z2sig.def lors du passage en coordonnees hybrides +c F. Forget 2002 +c Connaissant sig (niveaux "sigma" ou on veut mettre les couches) +c L'objectif est de calculer newsig telle que +c (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig +c Cela ne se résoud pas analytiquement: +c => on résoud par iterration bourrine +c ---------------------------------------------- +c Information : where exp(1-1./x**2) become << x +c x exp(1-1./x**2) /x +c 1 1 +c 0.68 0.5 +c 0.5 1.E-1 +c 0.391 1.E-2 +c 0.333 1.E-3 +c 0.295 1.E-4 +c 0.269 1.E-5 +c 0.248 1.E-6 +c => 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 +c 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 +c Test : on arete lorsque on approxime sig à moins de 0.01 m près +c (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 Index: /trunk/libf/dyn3dpar/disvert_terre.F90 =================================================================== --- /trunk/libf/dyn3dpar/disvert_terre.F90 (revision 109) +++ /trunk/libf/dyn3dpar/disvert_terre.F90 (revision 109) @@ -0,0 +1,142 @@ +! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ + +SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) + + ! Auteur : P. Le Van + + IMPLICIT NONE + + include "dimensions.h" + include "paramet.h" + include "iniprint.h" + include "logic.h" + + ! s = sigma ** kappa : coordonnee verticale + ! dsig(l) : epaisseur de la couche l ds la coord. s + ! sig(l) : sigma a l'interface des couches l et l-1 + ! ds(l) : distance entre les couches l et l-1 en coord.s + + REAL pa, preff + REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) + REAL presnivs(llm) + + REAL sig(llm+1), dsig(llm) + real zk, zkm1, dzk1, dzk2, k0, k1 + + INTEGER l + REAL dsigmin + REAL alpha, beta, deltaz, h + INTEGER iostat + REAL pi, x + + !----------------------------------------------------------------------- + + pi = 2 * ASIN(1.) + + OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) + + IF (iostat == 0) THEN + ! cas 1 on lit les options dans sigma.def: + READ(99, *) h ! hauteur d'echelle 8. + READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 + READ(99, *) beta ! facteur d'acroissement en haut 1.3 + READ(99, *) k0 ! nombre de couches dans la transition surf + READ(99, *) k1 ! nombre de couches dans la transition haute + CLOSE(99) + alpha=deltaz/(llm*h) + write(lunout, *)'disvert: h, alpha, k0, k1, beta' + + alpha=deltaz/tanh(1./k0)*2. + zkm1=0. + sig(1)=1. + do l=1, llm + sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) + zk=-h*log(sig(l+1)) + + dzk1=alpha*tanh(l/k0) + dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) + write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 + zkm1=zk + enddo + + sig(llm+1)=0. + + DO l = 1, llm + dsig(l) = sig(l)-sig(l+1) + end DO + ELSE + if (ok_strato) then + if (llm==39) then + dsigmin=0.3 + else if (llm==50) then + dsigmin=1. + else + WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster' + dsigmin=1. + endif + WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin + endif + + DO l = 1, llm + x = 2*asin(1.) * (l - 0.5) / (llm + 1) + + IF (ok_strato) THEN + dsig(l) =(dsigmin + 7. * SIN(x)**2) & + *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 + ELSE + dsig(l) = 1.0 + 7.0 * SIN(x)**2 + ENDIF + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + 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 lectba ..... + + bp(llmp1) = 0. + + DO l = 1, llm + bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) + ap(l) = pa * ( sig(l) - bp(l) ) + ENDDO + + bp(1)=1. + ap(1)=0. + + ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) + + write(lunout, *)'disvert: BP ' + write(lunout, *) bp + write(lunout, *)'disvert: AP ' + write(lunout, *) ap + + write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' + write(lunout, *)'couches calcules pour une pression de surface =', preff + write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' + write(lunout, *)'8km' + DO l = 1, llm + dpres(l) = bp(l) - bp(l+1) + presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) + write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & + log(preff/presnivs(l))*8. & + , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & + max(ap(l+1)+bp(l+1)*preff, 1.e-10)) + ENDDO + + write(lunout, *) 'disvert: PRESNIVS ' + write(lunout, *) presnivs + +END SUBROUTINE disvert Index: /trunk/libf/dyn3dpar/exner_milieu.F =================================================================== --- /trunk/libf/dyn3dpar/exner_milieu.F (revision 109) +++ /trunk/libf/dyn3dpar/exner_milieu.F (revision 109) @@ -0,0 +1,102 @@ + SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) +c +c Auteurs : F. Forget , Y. Wanherdrick +c P.Le Van , Fr. Hourdin . +c .......... +c +c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... +c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... +c +c ************************************************************************ +c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des +c couches . Pk(l) sera calcule aux milieux des couches l ,entre les +c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . +c ************************************************************************ +c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont +c la pression et la fonction d'Exner au sol . +c +c WARNING : CECI est une version speciale de exner_hyb originale +c Utilise dans la version martienne pour pouvoir +c tourner avec des coordonnees verticales complexe +c => Il ne verifie PAS la condition la proportionalite en +c energie totale/ interne / potentielle (F.Forget 2001) +c ( voir note de Fr.Hourdin ) , +c + IMPLICIT NONE +c +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comgeom.h" +#include "comvert.h" +#include "serre.h" + + INTEGER ngrid + REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm) + REAL ps(ngrid),pks(ngrid), beta(ngrid,llm) + +c .... variables locales ... + + INTEGER l, ij + REAL dum1 + + REAL ppn(iim),pps(iim) + REAL xpn, xps + REAL SSUM + EXTERNAL filtreg, SSUM + +c ------------- +c Calcul de pks +c ------------- + + DO ij = 1, ngrid + pks(ij) = cpp * ( ps(ij)/preff ) ** kappa + ENDDO + + DO ij = 1, iim + ppn(ij) = aire( ij ) * pks( ij ) + pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) + ENDDO + xpn = SSUM(iim,ppn,1) /apoln + xps = SSUM(iim,pps,1) /apols + + DO ij = 1, iip1 + pks( ij ) = xpn + pks( ij+ip1jm ) = xps + ENDDO +c +c +c .... Calcul de pk pour la couche l +c -------------------------------------------- +c + dum1 = cpp * (2*preff)**(-kappa) + DO l = 1, llm-1 + DO ij = 1, ngrid + pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa + ENDDO + ENDDO + +c .... Calcul de pk pour la couche l = llm .. +c (on met la meme distance (en log pression) entre Pk(llm) +c et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) + + DO ij = 1, ngrid + pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2) + ENDDO + + +c calcul de pkf +c ------------- + CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) + CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) + +c EST-CE UTILE ?? : calcul de beta +c -------------------------------- + DO l = 2, llm + DO ij = 1, ngrid + beta(ij,l) = pk(ij,l) / pk(ij,l-1) + ENDDO + ENDDO + + RETURN + END Index: /trunk/libf/dyn3dpar/exner_milieu_p.F =================================================================== --- /trunk/libf/dyn3dpar/exner_milieu_p.F (revision 109) +++ /trunk/libf/dyn3dpar/exner_milieu_p.F (revision 109) @@ -0,0 +1,142 @@ + SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) +c +c Auteurs : F. Forget , Y. Wanherdrick +c P.Le Van , Fr. Hourdin . +c .......... +c +c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... +c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... +c +c ************************************************************************ +c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des +c couches . Pk(l) sera calcule aux milieux des couches l ,entre les +c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . +c ************************************************************************ +c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont +c la pression et la fonction d'Exner au sol . +c +c WARNING : CECI est une version speciale de exner_hyb originale +c Utilise dans la version martienne pour pouvoir +c tourner avec des coordonnees verticales complexe +c => Il ne verifie PAS la condition la proportionalite en +c energie totale/ interne / potentielle (F.Forget 2001) +c ( voir note de Fr.Hourdin ) , +c + USE parallel + IMPLICIT NONE +c +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comgeom.h" +#include "comvert.h" +#include "serre.h" + + INTEGER ngrid + REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm) + REAL ps(ngrid),pks(ngrid), beta(ngrid,llm) + +c .... variables locales ... + + INTEGER l, ij + REAL dum1 + + REAL ppn(iim),pps(iim) + REAL xpn, xps + REAL SSUM + EXTERNAL SSUM + INTEGER ije,ijb,jje,jjb + +c$OMP BARRIER + +c ------------- +c Calcul de pks +c ------------- + + ijb=ij_begin + ije=ij_end + +c$OMP DO SCHEDULE(STATIC) + DO ij = ijb, ije + pks(ij) = cpp * ( ps(ij)/preff ) ** kappa + ENDDO +c$OMP ENDDO +c Synchro OPENMP ici + +c$OMP MASTER + if (pole_nord) then + DO ij = 1, iim + ppn(ij) = aire( ij ) * pks( ij ) + ENDDO + xpn = SSUM(iim,ppn,1) /apoln + + DO ij = 1, iip1 + pks( ij ) = xpn + ENDDO + endif + + if (pole_sud) then + DO ij = 1, iim + pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) + ENDDO + xps = SSUM(iim,pps,1) /apols + + DO ij = 1, iip1 + pks( ij+ip1jm ) = xps + ENDDO + endif +c$OMP END MASTER +c +c +c .... Calcul de pk pour la couche l +c -------------------------------------------- +c + dum1 = cpp * (2*preff)**(-kappa) + DO l = 1, llm-1 +c$OMP DO SCHEDULE(STATIC) + DO ij = ijb, ije + pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa + ENDDO +c$OMP ENDDO NOWAIT + ENDDO + +c .... Calcul de pk pour la couche l = llm .. +c (on met la meme distance (en log pression) entre Pk(llm) +c et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) + +c$OMP DO SCHEDULE(STATIC) + DO ij = ijb, ije + pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2) + ENDDO +c$OMP ENDDO NOWAIT + + +c calcul de pkf +c ------------- +c CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) + DO l = 1, llm +c$OMP DO SCHEDULE(STATIC) + DO ij = ijb, ije + pkf(ij,l)=pk(ij,l) + ENDDO +c$OMP ENDDO NOWAIT + ENDDO + +c$OMP BARRIER + + jjb=jj_begin + jje=jj_end + CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 ) + +c EST-CE UTILE ?? : calcul de beta +c -------------------------------- + DO l = 2, llm +c$OMP DO SCHEDULE(STATIC) + DO ij = ijb, ije + beta(ij,l) = pk(ij,l) / pk(ij,l-1) + ENDDO +c$OMP ENDDO NOWAIT + ENDDO + + RETURN + END Index: /trunk/libf/dyn3dpar/iniconst.F =================================================================== --- /trunk/libf/dyn3dpar/iniconst.F (revision 108) +++ /trunk/libf/dyn3dpar/iniconst.F (revision 109) @@ -53,6 +53,9 @@ c----------------------------------------------------------------------- - CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) -c + if (planet_type.eq."earth") then + CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) + else + CALL disvert_noterre + endif c RETURN Index: /trunk/libf/dyn3dpar/inidissip.F =================================================================== --- /trunk/libf/dyn3dpar/inidissip.F (revision 108) +++ /trunk/libf/dyn3dpar/inidissip.F (revision 109) @@ -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/dyn3dpar/leapfrog_p.F =================================================================== --- /trunk/libf/dyn3dpar/leapfrog_p.F (revision 108) +++ /trunk/libf/dyn3dpar/leapfrog_p.F (revision 109) @@ -271,5 +271,9 @@ CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif c$OMP END MASTER c----------------------------------------------------------------------- @@ -746,5 +750,9 @@ c$OMP BARRIER - CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif c$OMP BARRIER jD_cur = jD_ref + day_ini - day_ref @@ -1094,5 +1102,9 @@ CALL pression_p ( ip1jmp1, ap, bp, ps, p ) c$OMP BARRIER - CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + if (planet_type.eq."earth") then + CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + else + CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf ) + endif c$OMP BARRIER Index: /trunk/libf/dyn3dpar/limy.F =================================================================== --- /trunk/libf/dyn3dpar/limy.F (revision 108) +++ /trunk/libf/dyn3dpar/limy.F (revision 109) @@ -40,5 +40,5 @@ REAL qbyv(ip1jm,llm) - REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2 + REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2 Logical extremum,first save first @@ -116,14 +116,14 @@ c print*,dyqv(iip1+1) -c apn=abs(dyq(1)/dyqv(iip1+1)) +c appn=abs(dyq(1)/dyqv(iip1+1)) c print*,dyq(ip1jm+1) c print*,dyqv(ip1jm-iip1+1) -c aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +c apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) c do ij=2,iim -c apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -c aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +c appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +c apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) c enddo -c apn=min(pente_max/apn,1.) -c aps=min(pente_max/aps,1.) +c appn=min(pente_max/appn,1.) +c apps=min(pente_max/apps,1.) @@ -131,13 +131,13 @@ c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -c & apn=0. +c & appn=0. c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -c & aps=0. +c & apps=0. c limitation des pentes aux poles c do ij=1,iip1 -c dyq(ij)=apn*dyq(ij) -c dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +c dyq(ij)=appn*dyq(ij) +c dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) c enddo Index: /trunk/libf/dyn3dpar/vlsplt_p.F =================================================================== --- /trunk/libf/dyn3dpar/vlsplt_p.F (revision 108) +++ /trunk/libf/dyn3dpar/vlsplt_p.F (revision 109) @@ -561,5 +561,5 @@ REAL qbyv(ip1jm,llm) - REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs + REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs c REAL newq,oldmasse Logical extremum,first,testcpu @@ -732,14 +732,14 @@ C PRINT*,dyq(1) C PRINT*,dyqv(iip1+1) -C apn=abs(dyq(1)/dyqv(iip1+1)) +C appn=abs(dyq(1)/dyqv(iip1+1)) C PRINT*,dyq(ip1jm+1) C PRINT*,dyqv(ip1jm-iip1+1) -C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) C DO ij=2,iim -C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) C ENDDO -C apn=min(pente_max/apn,1.) -C aps=min(pente_max/aps,1.) +C appn=min(pente_max/appn,1.) +C apps=min(pente_max/apps,1.) C C @@ -747,13 +747,13 @@ C C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -C & apn=0. +C & appn=0. C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -C & aps=0. +C & apps=0. C C limitation des pentes aux poles C DO ij=1,iip1 -C dyq(ij)=apn*dyq(ij) -C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +C dyq(ij)=appn*dyq(ij) +C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) C ENDDO C Index: /trunk/libf/dyn3dpar/vlspltqs_p.F =================================================================== --- /trunk/libf/dyn3dpar/vlspltqs_p.F (revision 108) +++ /trunk/libf/dyn3dpar/vlspltqs_p.F (revision 109) @@ -780,14 +780,14 @@ C PRINT*,dyq(1) C PRINT*,dyqv(iip1+1) -C apn=abs(dyq(1)/dyqv(iip1+1)) +C appn=abs(dyq(1)/dyqv(iip1+1)) C PRINT*,dyq(ip1jm+1) C PRINT*,dyqv(ip1jm-iip1+1) -C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) +C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) C DO ij=2,iim -C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn) -C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps) +C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) +C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) C ENDDO -C apn=min(pente_max/apn,1.) -C aps=min(pente_max/aps,1.) +C appn=min(pente_max/appn,1.) +C apps=min(pente_max/apps,1.) C C @@ -795,13 +795,13 @@ C C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) -C & apn=0. +C & appn=0. C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) -C & aps=0. +C & apps=0. C C limitation des pentes aux poles C DO ij=1,iip1 -C dyq(ij)=apn*dyq(ij) -C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij) +C dyq(ij)=appn*dyq(ij) +C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) C ENDDO C