source: LMDZ5/branches/LF-private/libf/dyn3d/disvert_noterre.F @ 4585

Last change on this file since 4585 was 1654, checked in by aslmd, 12 years ago

  • Changements dans calfis pour appeler la physique generique planetaire NB: dans le futur il faudrait un appel commun mais actuellement differences dans les arguments...
  • Changements dans dynetat0 pour s'accommoder des tableaux "controle" terrestre et planeto
  • Passage de la variable pseudoalt dans comvert.h (compatibilite planeto)
  • arch : ajout des fichiers pour faire fonctionner sur ciclad avec ifort

verifie : compatibilite bit-par-bit des resultats et initialisations, ainsi que 1+1=2 avec unpun.sh


  • Changes in calfis so that generic physics for planets could be interfaced NB: ideally a common call should be implemented but for now there are differences in arguments...
  • Changes in dynetat0 to cope with difference in "controle" between earth and planeto
  • Transfer of variable pseudoalt into comvert.h (this is only used by planeto and done for compatibility)
  • arch: added files to allow use of ifort compiler on ciclad cluster

checked : bit-by-bit compatibility of results and starts, as well as 1+1=2 with unpun.sh


Author : AS equipe planeto

File size: 9.9 KB
Line 
1! $Id: $
2      SUBROUTINE disvert_noterre
3
4c    Auteur :  F. Forget Y. Wanherdrick, P. Levan
5c    Nouvelle version 100% Mars !!
6c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
7
8#ifdef CPP_IOIPSL
9      use IOIPSL
10#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      use ioipsl_getincom
13#endif
14
15      IMPLICIT NONE
16
17#include "dimensions.h"
18#include "paramet.h"
19#include "comvert.h"
20#include "comconst.h"
21#include "logic.h"
22#include "iniprint.h"
23c
24c=======================================================================
25c    Discretisation verticale en coordonnée hybride (ou sigma)
26c
27c=======================================================================
28c
29c   declarations:
30c   -------------
31c
32c
33      INTEGER l,ll
34      REAL snorm
35      REAL alpha,beta,gama,delta,deltaz
36      real quoi,quand
37      REAL zsig(llm),sig(llm+1)
38      INTEGER np,ierr
39      integer :: ierr1,ierr2,ierr3,ierr4
40      REAL x
41
42      REAL SSUM
43      EXTERNAL SSUM
44      real newsig
45      REAL dz0,dz1,nhaut,sig1,esig,csig,zz
46      real tt,rr,gg, prevz
47      real s(llm),dsig(llm)
48
49      integer iz
50      real z, ps,p
51      character(len=*),parameter :: modname="disvert_noterre"
52
53c
54c-----------------------------------------------------------------------
55c
56! Initializations:
57!      pi=2.*ASIN(1.) ! already done in iniconst
58     
59      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
60      CALL getin('hybrid',hybrid)
61      write(lunout,*) trim(modname),': hybrid=',hybrid
62
63! Ouverture possible de fichiers typiquement E.T.
64
65         open(99,file="esasig.def",status='old',form='formatted',
66     s   iostat=ierr2)
67         if(ierr2.ne.0) then
68              close(99)
69              open(99,file="z2sig.def",status='old',form='formatted',
70     s        iostat=ierr4)
71         endif
72
73c-----------------------------------------------------------------------
74c   cas 1 on lit les options dans esasig.def:
75c   ----------------------------------------
76
77      IF(ierr2.eq.0) then
78
79c        Lecture de esasig.def :
80c        Systeme peu souple, mais qui respecte en theorie
81c        La conservation de l'energie (conversion Energie potentielle
82c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
83
84         write(lunout,*)'*****************************'
85         write(lunout,*)'WARNING reading esasig.def'
86         write(lunout,*)'*****************************'
87         READ(99,*) scaleheight
88         READ(99,*) dz0
89         READ(99,*) dz1
90         READ(99,*) nhaut
91         CLOSE(99)
92
93         dz0=dz0/scaleheight
94         dz1=dz1/scaleheight
95
96         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
97
98         esig=1.
99
100         do l=1,20
101            esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
102         enddo
103         csig=(1./sig1-1.)/(exp(esig)-1.)
104
105         DO L = 2, llm
106            zz=csig*(exp(esig*(l-1.))-1.)
107            sig(l) =1./(1.+zz)
108     &      * tanh(.5*(llm+1-l)/nhaut)
109         ENDDO
110         sig(1)=1.
111         sig(llm+1)=0.
112         quoi      = 1. + 2.* kappa
113         s( llm )  = 1.
114         s(llm-1) = quoi
115         IF( llm.gt.2 )  THEN
116            DO  ll = 2, llm-1
117               l         = llm+1 - ll
118               quand     = sig(l+1)/ sig(l)
119               s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
120            ENDDO
121         END IF
122c
123         snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
124         DO l = 1, llm
125            s(l)    = s(l)/ snorm
126         ENDDO
127
128c-----------------------------------------------------------------------
129c   cas 2 on lit les options dans z2sig.def:
130c   ----------------------------------------
131
132      ELSE IF(ierr4.eq.0) then
133         write(lunout,*)'****************************'
134         write(lunout,*)'Reading z2sig.def'
135         write(lunout,*)'****************************'
136
137         READ(99,*) scaleheight
138         do l=1,llm
139            read(99,*) zsig(l)
140         end do
141         CLOSE(99)
142
143         sig(1) =1
144         do l=2,llm
145           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) +
146     &                      exp(-zsig(l-1)/scaleheight) )
147         end do
148         sig(llm+1) =0
149
150c-----------------------------------------------------------------------
151      ELSE
152         write(lunout,*) 'didn t you forget something ??? '
153         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
154         stop
155      ENDIF
156c-----------------------------------------------------------------------
157
158      DO l=1,llm
159        nivsigs(l) = REAL(l)
160      ENDDO
161
162      DO l=1,llmp1
163        nivsig(l)= REAL(l)
164      ENDDO
165
166 
167c-----------------------------------------------------------------------
168c    ....  Calculs  de ap(l) et de bp(l)  ....
169c    .........................................
170c
171c   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
172c-----------------------------------------------------------------------
173c
174
175      if (hybrid) then  ! use hybrid coordinates
176         write(lunout,*) "*********************************"
177         write(lunout,*) "Using hybrid vertical coordinates"
178         write(lunout,*)
179c        Coordonnees hybrides avec mod
180         DO l = 1, llm
181
182         call sig_hybrid(sig(l),pa,preff,newsig)
183            bp(l) = EXP( 1. - 1./(newsig**2)  )
184            ap(l) = pa * (newsig - bp(l) )
185         enddo
186         bp(llmp1) = 0.
187         ap(llmp1) = 0.
188      else ! use sigma coordinates
189         write(lunout,*) "********************************"
190         write(lunout,*) "Using sigma vertical coordinates"
191         write(lunout,*)
192c        Pour ne pas passer en coordonnees hybrides
193         DO l = 1, llm
194            ap(l) = 0.
195            bp(l) = sig(l)
196         ENDDO
197         ap(llmp1) = 0.
198      endif
199
200      bp(llmp1) =   0.
201
202      write(lunout,*) trim(modname),': BP '
203      write(lunout,*)  bp
204      write(lunout,*) trim(modname),': AP '
205      write(lunout,*)  ap
206
207c     Calcul au milieu des couches :
208c     WARNING : le choix de placer le milieu des couches au niveau de
209c     pression intermédiaire est arbitraire et pourrait etre modifié.
210c     Le calcul du niveau pour la derniere couche
211c     (on met la meme distance (en log pression)  entre P(llm)
212c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
213c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
214
215      DO l = 1, llm-1
216       aps(l) =  0.5 *( ap(l) +ap(l+1))
217       bps(l) =  0.5 *( bp(l) +bp(l+1))
218      ENDDO
219     
220      if (hybrid) then
221         aps(llm) = aps(llm-1)**2 / aps(llm-2)
222         bps(llm) = 0.5*(bp(llm) + bp(llm+1))
223      else
224         bps(llm) = bps(llm-1)**2 / bps(llm-2)
225         aps(llm) = 0.
226      end if
227
228      write(lunout,*) trim(modname),': BPs '
229      write(lunout,*)  bps
230      write(lunout,*) trim(modname),': APs'
231      write(lunout,*)  aps
232
233      DO l = 1, llm
234       presnivs(l) = aps(l)+bps(l)*preff
235       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
236      ENDDO
237
238      write(lunout,*)trim(modname),' : PRESNIVS'
239      write(lunout,*)presnivs
240      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
241     &                'height of ',scaleheight,' km)'
242      write(lunout,*)pseudoalt
243
244c     --------------------------------------------------
245c     This can be used to plot the vertical discretization
246c     (> xmgrace -nxy testhybrid.tab )
247c     --------------------------------------------------
248c     open (53,file='testhybrid.tab')
249c     scaleheight=15.5
250c     do iz=0,34
251c       z = -5 + min(iz,34-iz)
252c     approximation of scale height for Venus
253c       scaleheight = 15.5 - z/55.*10.
254c       ps = preff*exp(-z/scaleheight)
255c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
256c       do l=2,llm
257c     approximation of scale height for Venus
258c          if (zsig(l-1).le.55.) then
259c             scaleheight = 15.5 - zsig(l-1)/55.*10.
260c          else
261c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
262c          endif
263c          zsig(l)= zsig(l-1)-scaleheight*
264c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
265c       end do
266c       write(53,'(I3,50F10.5)') iz, zsig
267c      end do
268c      close(53)
269c     --------------------------------------------------
270
271
272      RETURN
273      END
274
275c ************************************************************
276      subroutine sig_hybrid(sig,pa,preff,newsig)
277c     ----------------------------------------------
278c     Subroutine utilisee pour calculer des valeurs de sigma modifie
279c     pour conserver les coordonnees verticales decrites dans
280c     esasig.def/z2sig.def lors du passage en coordonnees hybrides
281c     F. Forget 2002
282c     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
283c     L'objectif est de calculer newsig telle que
284c       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
285c     Cela ne se résoud pas analytiquement:
286c     => on résoud par iterration bourrine
287c     ----------------------------------------------
288c     Information  : where exp(1-1./x**2) become << x
289c           x      exp(1-1./x**2) /x
290c           1           1
291c           0.68       0.5
292c           0.5        1.E-1
293c           0.391      1.E-2
294c           0.333      1.E-3
295c           0.295      1.E-4
296c           0.269      1.E-5
297c           0.248      1.E-6
298c        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
299
300
301      implicit none
302      real x1, x2, sig,pa,preff, newsig, F
303      integer j
304
305      newsig = sig
306      x1=0
307      x2=1
308      if (sig.ge.1) then
309            newsig= sig
310      else if (sig*preff/pa.ge.0.25) then
311        DO J=1,9999  ! nombre d''iteration max
312          F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
313c         write(0,*) J, ' newsig =', newsig, ' F= ', F
314          if (F.gt.1) then
315              X2 = newsig
316              newsig=(X1+newsig)*0.5
317          else
318              X1 = newsig
319              newsig=(X2+newsig)*0.5
320          end if
321c         Test : on arete lorsque on approxime sig à moins de 0.01 m près
322c         (en pseudo altitude) :
323          IF(abs(10.*log(F)).LT.1.E-5) goto 999
324        END DO
325       else   !    if (sig*preff/pa.le.0.25) then
326             newsig= sig*preff/pa
327       end if
328 999   continue
329       Return
330      END
Note: See TracBrowser for help on using the repository browser.