source: LMDZ.3.3/trunk/libf/dyn3d/disvert.F @ 507

Last change on this file since 507 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 3.2 KB
Line 
1      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
2
3c    Auteur :  P. Le Van .
4c
5      IMPLICIT NONE
6
7#include "dimensions.h"
8#include "paramet.h"
9c
10c=======================================================================
11c
12c
13c    s = sigma ** kappa   :  coordonnee  verticale
14c    dsig(l)            : epaisseur de la couche l ds la coord.  s
15c    sig(l)             : sigma a l'interface des couches l et l-1
16c    ds(l)              : distance entre les couches l et l-1 en coord.s
17c
18c=======================================================================
19c
20      REAL pa,preff
21      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
22      REAL presnivs(llm)
23c
24c   declarations:
25c   -------------
26c
27      REAL sig(llm+1),dsig(llm)
28c
29      INTEGER l
30      REAL snorm
31      REAL alpha,beta,gama,delta,deltaz,h
32      INTEGER np,ierr
33      REAL pi,x
34
35      REAL SSUM
36      EXTERNAL SSUM
37c
38c-----------------------------------------------------------------------
39c
40      pi=2.*ASIN(1.)
41
42      OPEN(99,file='sigma.def',status='old',form='formatted',
43     s   iostat=ierr)
44
45c-----------------------------------------------------------------------
46c   cas 1 on lit les options dans sigma.def:
47c   ----------------------------------------
48
49      IF (ierr.eq.0) THEN
50
51      READ(99,*) deltaz
52      READ(99,*) h
53      READ(99,*) beta
54      READ(99,*) gama
55      READ(99,*) delta
56      READ(99,*) np
57      CLOSE(99)
58      alpha=deltaz/(llm*h)
59c
60
61       DO 1  l = 1, llm
62       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
63     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
64     $            (1.-l/FLOAT(llm))*delta )
65   1   CONTINUE
66
67       sig(1)=1.
68       DO 101 l=1,llm-1
69          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
70101    CONTINUE
71       sig(llm+1)=0.
72
73       DO 2  l = 1, llm
74       dsig(l) = sig(l)-sig(l+1)
75   2   CONTINUE
76c
77
78      ELSE
79c-----------------------------------------------------------------------
80c   cas 2 ancienne discretisation (LMD5...):
81c   ----------------------------------------
82
83      PRINT*,'WARNING!!! Ancienne discretisation verticale'
84
85      h=7.
86      snorm  = 0.
87      DO l = 1, llm
88         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
89         dsig(l) = 1.0 + 7.0 * SIN(x)**2
90         snorm = snorm + dsig(l)
91      ENDDO
92      snorm = 1./snorm
93      DO l = 1, llm
94         dsig(l) = dsig(l)*snorm
95      ENDDO
96      sig(llm+1) = 0.
97      DO l = llm, 1, -1
98         sig(l) = sig(l+1) + dsig(l)
99      ENDDO
100
101      ENDIF
102
103
104      DO l=1,llm
105        nivsigs(l) = FLOAT(l)
106      ENDDO
107
108      DO l=1,llmp1
109        nivsig(l)= FLOAT(l)
110      ENDDO
111
112c
113c    ....  Calculs  de ap(l) et de bp(l)  ....
114c    .........................................
115c
116c
117c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
118c
119
120      bp(llmp1) =   0.
121
122      DO l = 1, llm
123cc
124ccc    ap(l) = 0.
125ccc    bp(l) = sig(l)
126
127      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
128      ap(l) = pa * ( sig(l) - bp(l) )
129c
130      ENDDO
131      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
132
133      PRINT *,' BP '
134      PRINT *,  bp
135      PRINT *,' AP '
136      PRINT *,  ap
137
138      DO l = 1, llm
139       dpres(l) = bp(l) - bp(l+1)
140       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
141      ENDDO
142
143      PRINT *,' PRESNIVS '
144      PRINT *,presnivs
145
146      RETURN
147      END
Note: See TracBrowser for help on using the repository browser.