source: LMDZ5/branches/LMDZ5-dev032011/libf/dyn3d/disvert.F90 @ 2069

Last change on this file since 2069 was 1492, checked in by Laurent Fairhead, 14 years ago

Merge of development branch LMDZ5V2.0-dev r1455:r1491 into the trunk.
Validation made locally: restart files are strictly equal between the HEAD of the trunk
and r1491 of LMDZ5V2.0-dev


Synchro de la branche de développement LMDZ5V2.0-dev r1455:r1491 et de la trunk
Validation faite en local: les fichiers restart sont équivalents entre la HEAD de la trunk
et la révision r1491 de LMDZ5V2.0-dev

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.6 KB
Line 
1! $Id: disvert.F90 1492 2011-03-08 08:10:25Z fairhead $
2
3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
4
5  ! Auteur : P. Le Van
6
7  IMPLICIT NONE
8
9  include "dimensions.h"
10  include "paramet.h"
11  include "iniprint.h"
12  include "logic.h"
13
14  ! s = sigma ** kappa : coordonnee verticale
15  ! dsig(l) : epaisseur de la couche l ds la coord. s
16  ! sig(l) : sigma a l'interface des couches l et l-1
17  ! ds(l) : distance entre les couches l et l-1 en coord.s
18
19  REAL pa, preff
20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
21  REAL presnivs(llm)
22
23  REAL sig(llm+1), dsig(llm)
24  real zk, zkm1, dzk1, dzk2, k0, k1
25
26  INTEGER l
27  REAL dsigmin
28  REAL alpha, beta, deltaz, h
29  INTEGER iostat
30  REAL pi, x
31
32  !-----------------------------------------------------------------------
33
34  pi = 2 * ASIN(1.)
35
36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
37
38  IF (iostat == 0) THEN
39     ! cas 1 on lit les options dans sigma.def:
40     READ(99, *) h ! hauteur d'echelle 8.
41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
43     READ(99, *) k0 ! nombre de couches dans la transition surf
44     READ(99, *) k1 ! nombre de couches dans la transition haute
45     CLOSE(99)
46     alpha=deltaz/(llm*h)
47     write(lunout, *)'h, alpha, k0, k1, beta'
48
49     alpha=deltaz/tanh(1./k0)*2.
50     zkm1=0.
51     sig(1)=1.
52     do l=1, llm
53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
55        zk=-h*log(sig(l+1))
56
57        dzk1=alpha*tanh(l/k0)
58        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
60        zkm1=zk
61     enddo
62
63     sig(llm+1)=0.
64
65     DO l = 1, llm
66        dsig(l) = sig(l)-sig(l+1)
67     end DO
68  ELSE
69     if (ok_strato) then
70        if (llm==39) then
71           dsigmin=0.3
72        else if (llm==50) then
73           dsigmin=1.
74        else
75           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
76           dsigmin=1.
77        endif
78        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
79     endif
80
81     DO l = 1, llm
82        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
83
84        IF (ok_strato) THEN
85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
86                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
87        ELSE
88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
89        ENDIF
90     ENDDO
91     dsig = dsig / sum(dsig)
92     sig(llm+1) = 0.
93     DO l = llm, 1, -1
94        sig(l) = sig(l+1) + dsig(l)
95     ENDDO
96  ENDIF
97
98  DO l=1, llm
99     nivsigs(l) = REAL(l)
100  ENDDO
101
102  DO l=1, llmp1
103     nivsig(l)= REAL(l)
104  ENDDO
105
106  ! .... Calculs de ap(l) et de bp(l) ....
107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
108
109  bp(llmp1) = 0.
110
111  DO l = 1, llm
112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
113     ap(l) = pa * ( sig(l) - bp(l) )
114  ENDDO
115
116  bp(1)=1.
117  ap(1)=0.
118
119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
120
121  write(lunout, *)' BP '
122  write(lunout, *) bp
123  write(lunout, *)' AP '
124  write(lunout, *) ap
125
126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
127  write(lunout, *)'couches calcules pour une pression de surface =', preff
128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
129  write(lunout, *)'8km'
130  DO l = 1, llm
131     dpres(l) = bp(l) - bp(l+1)
132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
134          log(preff/presnivs(l))*8. &
135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
137  ENDDO
138
139  write(lunout, *) 'PRESNIVS '
140  write(lunout, *) presnivs
141
142END SUBROUTINE disvert
Note: See TracBrowser for help on using the repository browser.