source: trunk/LMDZ.MARS/libf/dyn3d/disvert.F @ 1594

Last change on this file since 1594 was 1422, checked in by milmd, 10 years ago

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