source: LMDZ4/trunk/libf/dyn3d/disvert.F @ 802

Last change on this file since 802 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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