source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3d_common/disvert_noterre.F @ 5160

Last change on this file since 5160 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

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