source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/disvert.F @ 1523

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

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper


General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

  • 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 1299 2010-01-20 14:27:21Z idelkadi $
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.) * (REAL(l)-0.5) / REAL(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) = REAL(l)
138      ENDDO
139
140      DO l=1,llmp1
141        nivsig(l)= REAL(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.