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

Last change on this file since 1033 was 999, checked in by Laurent Fairhead, 16 years ago
  • Modifs sur le parallelisme: masquage dans la physique
  • Inclusion strato
  • mise en coherence etat0
  • le mode offline fonctionne maintenant en parallele,
  • les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
  • la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
  • L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.5 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"
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
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      h=7.
99      snorm  = 0.
100      DO l = 1, llm
101         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
102
103         IF (ok_strato) THEN
104           dsig(l) =(1.0 + 7.0 * SIN(x)**2)
105     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
106         ELSE
107           dsig(l) = 1.0 + 7.0 * SIN(x)**2
108         ENDIF
109
110         snorm = snorm + dsig(l)
111      ENDDO
112      snorm = 1./snorm
113      DO l = 1, llm
114         dsig(l) = dsig(l)*snorm
115      ENDDO
116      sig(llm+1) = 0.
117      DO l = llm, 1, -1
118         sig(l) = sig(l+1) + dsig(l)
119      ENDDO
120
121      ENDIF
122
123
124      DO l=1,llm
125        nivsigs(l) = FLOAT(l)
126      ENDDO
127
128      DO l=1,llmp1
129        nivsig(l)= FLOAT(l)
130      ENDDO
131
132c
133c    ....  Calculs  de ap(l) et de bp(l)  ....
134c    .........................................
135c
136c
137c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
138c
139
140      bp(llmp1) =   0.
141
142      DO l = 1, llm
143cc
144ccc    ap(l) = 0.
145ccc    bp(l) = sig(l)
146
147      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
148      ap(l) = pa * ( sig(l) - bp(l) )
149c
150      ENDDO
151      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
152
153      write(lunout,*)' BP '
154      write(lunout,*)  bp
155      write(lunout,*)' AP '
156      write(lunout,*)  ap
157
158      write(lunout,*)
159     .'Niveaux de pressions approximatifs aux centres des'
160      write(lunout,*)'couches calcules pour une pression de surface =',
161     .                 preff
162      write(lunout,*)
163     .     'et altitudes equivalentes pour une hauteur d echelle de'
164      write(lunout,*)'8km'
165      DO l = 1, llm
166       dpres(l) = bp(l) - bp(l+1)
167       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
168       write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
169     .        log(preff/presnivs(l))*8.
170     .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
171     .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
172      ENDDO
173
174      write(lunout,*)' PRESNIVS '
175      write(lunout,*)presnivs
176
177      RETURN
178      END
Note: See TracBrowser for help on using the repository browser.