source: LMDZ4/trunk/libf/dyn3d/disvert0.F @ 843

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