source: LMDZ6/trunk/libf/dyn3d_common/disvert_noterre.f90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 23 hours ago

Turn paramet.h into a module
  • 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: 9.4 KB
Line 
1! $Id: $
2SUBROUTINE disvert_noterre
3
4  !    Auteur :  F. Forget Y. Wanherdrick, P. Levan
5  !    Nouvelle version 100% Mars !!
6  !    On l'utilise aussi pour Venus et Titan, legerment modifiee.
7
8  use IOIPSL
9
10  USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt, &
11        nivsig,nivsigs,pa,preff,scaleheight
12  USE comconst_mod, ONLY: kappa
13  USE logic_mod, ONLY: hybrid
14
15  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
16USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
17          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
18IMPLICIT NONE
19
20
21
22  include "iniprint.h"
23  !
24  !=======================================================================
25  !    Discretisation verticale en coordonn�e hybride (ou sigma)
26  !
27  !=======================================================================
28  !
29  !   declarations:
30  !   -------------
31  !
32  !
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
53  !
54  !-----------------------------------------------------------------------
55  !
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           iostat=ierr2)
67     if(ierr2.ne.0) then
68          close(99)
69          open(99,file="z2sig.def",status='old',form='formatted', &
70                iostat=ierr4)
71     endif
72
73  !-----------------------------------------------------------------------
74  !   cas 1 on lit les options dans esasig.def:
75  !   ----------------------------------------
76
77  IF(ierr2.eq.0) then
78
79     ! Lecture de esasig.def :
80     ! Systeme peu souple, mais qui respecte en theorie
81     ! La conservation de l'energie (conversion Energie potentielle
82     ! <-> 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
122  !
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
128  !-----------------------------------------------------------------------
129  !   cas 2 on lit les options dans z2sig.def:
130  !   ----------------------------------------
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
150  !-----------------------------------------------------------------------
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
156  !-----------------------------------------------------------------------
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
167  !-----------------------------------------------------------------------
168  !    ....  Calculs  de ap(l) et de bp(l)  ....
169  !    .........................................
170  !
171  !   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
172  !-----------------------------------------------------------------------
173  !
174
175  if (hybrid) then  ! use hybrid coordinates
176     write(lunout,*) "*********************************"
177     write(lunout,*) "Using hybrid vertical coordinates"
178     write(lunout,*)
179     ! 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,*)
192     ! 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
207  ! Calcul au milieu des couches :
208  ! WARNING : le choix de placer le milieu des couches au niveau de
209  ! pression interm�diaire est arbitraire et pourrait etre modifi�.
210  ! Le calcul du niveau pour la derniere couche
211  ! (on met la meme distance (en log pression)  entre P(llm)
212  ! et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
213  ! 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
244  ! --------------------------------------------------
245  ! This can be used to plot the vertical discretization
246  ! (> xmgrace -nxy testhybrid.tab )
247  ! --------------------------------------------------
248  ! open (53,file='testhybrid.tab')
249  ! scaleheight=15.5
250  ! do iz=0,34
251  !   z = -5 + min(iz,34-iz)
252  ! approximation of scale height for Venus
253  !   scaleheight = 15.5 - z/55.*10.
254  !   ps = preff*exp(-z/scaleheight)
255  !   zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
256  !   do l=2,llm
257  ! approximation of scale height for Venus
258  !      if (zsig(l-1).le.55.) then
259  !         scaleheight = 15.5 - zsig(l-1)/55.*10.
260  !      else
261  !         scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
262  !      endif
263  !      zsig(l)= zsig(l-1)-scaleheight*
264  !    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
265  !   end do
266  !   write(53,'(I3,50F10.5)') iz, zsig
267  !  end do
268  !  close(53)
269  ! --------------------------------------------------
270
271
272  RETURN
273END SUBROUTINE disvert_noterre
274
275! ************************************************************
276subroutine sig_hybrid(sig,pa,preff,newsig)
277  ! ----------------------------------------------
278  ! Subroutine utilisee pour calculer des valeurs de sigma modifie
279  ! pour conserver les coordonnees verticales decrites dans
280  ! esasig.def/z2sig.def lors du passage en coordonnees hybrides
281  ! F. Forget 2002
282  ! Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
283  ! L'objectif est de calculer newsig telle que
284  !   (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
285  ! Cela ne se r�soud pas analytiquement:
286  ! => on r�soud par iterration bourrine
287  ! ----------------------------------------------
288  ! Information  : where exp(1-1./x**2) become << x
289  !       x      exp(1-1./x**2) /x
290  !       1           1
291  !       0.68       0.5
292  !       0.5        1.E-1
293  !       0.391      1.E-2
294  !       0.333      1.E-3
295  !       0.295      1.E-4
296  !       0.269      1.E-5
297  !       0.248      1.E-6
298  !    => 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
313      ! 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
321      ! Test : on arete lorsque on approxime sig � moins de 0.01 m pr�s
322      ! (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
330END SUBROUTINE sig_hybrid
Note: See TracBrowser for help on using the repository browser.