source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/exner_milieu.F @ 3693

Last change on this file since 3693 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

File size: 4.3 KB
Line 
1!
2! $Id $
3!
4      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
5c
6c     Auteurs :  F. Forget , Y. Wanherdrick
7c P.Le Van  , Fr. Hourdin  .
8c    ..........
9c
10c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
11c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
12c
13c   ************************************************************************
14c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
15c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
16c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
17c   ************************************************************************
18c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
19c    la pression et la fonction d'Exner  au  sol  .
20c
21c     WARNING : CECI est une version speciale de exner_hyb originale
22c               Utilise dans la version martienne pour pouvoir
23c               tourner avec des coordonnees verticales complexe
24c              => Il ne verifie PAS la condition la proportionalite en
25c              energie totale/ interne / potentielle (F.Forget 2001)
26c    ( voir note de Fr.Hourdin )  ,
27c
28      IMPLICIT NONE
29c
30#include "dimensions.h"
31#include "paramet.h"
32#include "comconst.h"
33#include "comgeom.h"
34#include "comvert.h"
35#include "serre.h"
36
37      INTEGER  ngrid
38      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
39      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
40
41c    .... variables locales   ...
42
43      INTEGER l, ij
44      REAL dum1
45
46      REAL ppn(iim),pps(iim)
47      REAL xpn, xps
48      REAL SSUM
49      EXTERNAL SSUM
50      logical,save :: firstcall=.true.
51      character(len=*),parameter :: modname="exner_milieu"
52
53      ! Sanity check
54      if (firstcall) then
55        ! check that vertical discretization is compatible
56        ! with this routine
57        if (disvert_type.ne.2) then
58          call abort_gcm(modname,
59     &     "this routine should only be called if disvert_type==2",42)
60        endif
61       
62        ! sanity checks for Shallow Water case (1 vertical layer)
63        if (llm.eq.1) then
64          if (kappa.ne.1) then
65            call abort_gcm(modname,
66     &      "kappa!=1 , but running in Shallow Water mode!!",42)
67          endif
68          if (cpp.ne.r) then
69            call abort_gcm(modname,
70     &      "cpp!=r , but running in Shallow Water mode!!",42)
71          endif
72        endif ! of if (llm.eq.1)
73
74        firstcall=.false.
75      endif ! of if (firstcall)
76
77!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
78      if (llm.eq.1) then
79     
80        ! Compute pks(:),pk(:),pkf(:)
81       
82        DO   ij  = 1, ngrid
83          pks(ij) = (cpp/preff) * ps(ij)
84          pk(ij,1) = .5*pks(ij)
85        ENDDO
86       
87        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
88        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
89       
90        ! our work is done, exit routine
91        return
92
93      endif ! of if (llm.eq.1)
94
95!!!! General case:
96
97c     -------------
98c     Calcul de pks
99c     -------------
100   
101      DO   ij  = 1, ngrid
102        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
103      ENDDO
104
105      DO  ij   = 1, iim
106        ppn(ij) = aire(   ij   ) * pks(  ij     )
107        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
108      ENDDO
109      xpn      = SSUM(iim,ppn,1) /apoln
110      xps      = SSUM(iim,pps,1) /apols
111
112      DO ij   = 1, iip1
113        pks(   ij     )  =  xpn
114        pks( ij+ip1jm )  =  xps
115      ENDDO
116c
117c
118c    .... Calcul de pk  pour la couche l
119c    --------------------------------------------
120c
121      dum1 = cpp * (2*preff)**(-kappa)
122      DO l = 1, llm-1
123        DO   ij   = 1, ngrid
124         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
125        ENDDO
126      ENDDO
127
128c    .... Calcul de pk  pour la couche l = llm ..
129c    (on met la meme distance (en log pression)  entre Pk(llm)
130c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
131
132      DO   ij   = 1, ngrid
133         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
134      ENDDO
135
136
137c    calcul de pkf
138c    -------------
139      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
140      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
141     
142c    EST-CE UTILE ?? : calcul de beta
143c    --------------------------------
144      DO l = 2, llm
145        DO   ij   = 1, ngrid
146          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
147        ENDDO
148      ENDDO
149
150      RETURN
151      END
Note: See TracBrowser for help on using the repository browser.