source: LMDZ4/trunk/libf/dyn3dpar/disvert.F @ 1353

Last change on this file since 1353 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.8 KB
Line 
1!
2! $Id: disvert.F 1279 2009-12-10 09:02:56Z musat $
3!
4      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
5
6c    Auteur :  P. Le Van .
7c
8      IMPLICIT NONE
9
10#include "dimensions.h"
11#include "paramet.h"
12#include "iniprint.h"
13#include "logic.h"
14c
15c=======================================================================
16c
17c
18c    s = sigma ** kappa   :  coordonnee  verticale
19c    dsig(l)            : epaisseur de la couche l ds la coord.  s
20c    sig(l)             : sigma a l'interface des couches l et l-1
21c    ds(l)              : distance entre les couches l et l-1 en coord.s
22c
23c=======================================================================
24c
25      REAL pa,preff
26      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
27      REAL presnivs(llm)
28c
29c   declarations:
30c   -------------
31c
32      REAL sig(llm+1),dsig(llm)
33       real zzz(1:llm+1)
34       real dzz(1:llm)
35      real zk,zkm1,dzk1,dzk2,k0,k1
36c
37      INTEGER l
38      REAL snorm,dsigmin
39      REAL alpha,beta,gama,delta,deltaz,h
40      INTEGER np,ierr
41      REAL pi,x
42
43      REAL SSUM
44c
45c-----------------------------------------------------------------------
46c
47      pi=2.*ASIN(1.)
48
49      OPEN(99,file='sigma.def',status='old',form='formatted',
50     s   iostat=ierr)
51
52c-----------------------------------------------------------------------
53c   cas 1 on lit les options dans sigma.def:
54c   ----------------------------------------
55
56      IF (ierr.eq.0) THEN
57
58      READ(99,*) h           ! hauteur d'echelle 8.
59      READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
60      READ(99,*) beta        ! facteur d'acroissement en haut 1.3
61      READ(99,*) k0          ! nombre de couches dans la transition surf
62      READ(99,*) k1          ! nombre de couches dans la transition haute
63      CLOSE(99)
64      alpha=deltaz/(llm*h)
65      write(lunout,*)'h,alpha,k0,k1,beta'
66
67c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
68
69      alpha=deltaz/tanh(1./k0)*2.
70      zkm1=0.
71      sig(1)=1.
72      do l=1,llm
73        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
74     + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
75        zk=-h*log(sig(l+1))
76
77        dzk1=alpha*tanh(l/k0)
78        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
79        write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
80        zkm1=zk
81      enddo
82
83      sig(llm+1)=0.
84
85c
86       DO 2  l = 1, llm
87       dsig(l) = sig(l)-sig(l+1)
88   2   CONTINUE
89c
90
91      ELSE
92c-----------------------------------------------------------------------
93c   cas 2 ancienne discretisation (LMD5...):
94c   ----------------------------------------
95
96      WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
97
98      if (ok_strato) then
99         if (llm==39) then
100            dsigmin=0.3
101         else if (llm==50) then
102            dsigmin=1.
103         else
104            WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
105            dsigmin=1.
106         endif
107         WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
108      endif
109
110      h=7.
111      snorm  = 0.
112      DO l = 1, llm
113         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
114
115         IF (ok_strato) THEN
116           dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
117     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
118         ELSE
119           dsig(l) = 1.0 + 7.0 * SIN(x)**2
120         ENDIF
121
122         snorm = snorm + dsig(l)
123      ENDDO
124      snorm = 1./snorm
125      DO l = 1, llm
126         dsig(l) = dsig(l)*snorm
127      ENDDO
128      sig(llm+1) = 0.
129      DO l = llm, 1, -1
130         sig(l) = sig(l+1) + dsig(l)
131      ENDDO
132
133      ENDIF
134
135
136      DO l=1,llm
137        nivsigs(l) = FLOAT(l)
138      ENDDO
139
140      DO l=1,llmp1
141        nivsig(l)= FLOAT(l)
142      ENDDO
143
144c
145c    ....  Calculs  de ap(l) et de bp(l)  ....
146c    .........................................
147c
148c
149c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
150c
151
152      bp(llmp1) =   0.
153
154      DO l = 1, llm
155cc
156ccc    ap(l) = 0.
157ccc    bp(l) = sig(l)
158
159      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
160      ap(l) = pa * ( sig(l) - bp(l) )
161c
162      ENDDO
163
164      bp(1)=1.
165      ap(1)=0.
166
167      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
168
169      write(lunout,*)' BP '
170      write(lunout,*)  bp
171      write(lunout,*)' AP '
172      write(lunout,*)  ap
173
174      write(lunout,*)
175     .'Niveaux de pressions approximatifs aux centres des'
176      write(lunout,*)'couches calcules pour une pression de surface =',
177     .                 preff
178      write(lunout,*)
179     .     'et altitudes equivalentes pour une hauteur d echelle de'
180      write(lunout,*)'8km'
181      DO l = 1, llm
182       dpres(l) = bp(l) - bp(l+1)
183       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
184       write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
185     .        log(preff/presnivs(l))*8.
186     .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
187     .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
188      ENDDO
189
190      write(lunout,*)' PRESNIVS '
191      write(lunout,*)presnivs
192
193      RETURN
194      END
Note: See TracBrowser for help on using the repository browser.