source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/exner_hyb.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.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.6 KB
Line 
1!
2! $Id $
3!
4      SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
5c
6c     Auteurs :  P.Le Van  , Fr. Hourdin  .
7c    ..........
8c
9c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
10c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
11c
12c   ************************************************************************
13c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
14c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
15c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
16c   ************************************************************************
17c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
18c    la pression et la fonction d'Exner  au  sol  .
19c
20c                                 -------- z                                   
21c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
22c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
23c    ( voir note de Fr.Hourdin )  ,
24c
25c    on determine successivement , du haut vers le bas des couches, les
26c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
27c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
28c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
29c
30c
31      IMPLICIT NONE
32c
33#include "dimensions.h"
34#include "paramet.h"
35#include "comconst.h"
36#include "comgeom.h"
37#include "comvert.h"
38#include "serre.h"
39
40      INTEGER  ngrid
41      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
42      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
43
44c    .... variables locales   ...
45
46      INTEGER l, ij
47      REAL unpl2k,dellta
48
49      REAL ppn(iim),pps(iim)
50      REAL xpn, xps
51      REAL SSUM
52c
53      logical,save :: firstcall=.true.
54      character(len=*),parameter :: modname="exner_hyb"
55     
56      ! Sanity check
57      if (firstcall) then
58        ! check that vertical discretization is compatible
59        ! with this routine
60        if (disvert_type.ne.1) then
61          call abort_gcm(modname,
62     &     "this routine should only be called if disvert_type==1",42)
63        endif
64       
65        ! sanity checks for Shallow Water case (1 vertical layer)
66        if (llm.eq.1) then
67          if (kappa.ne.1) then
68            call abort_gcm(modname,
69     &      "kappa!=1 , but running in Shallow Water mode!!",42)
70          endif
71          if (cpp.ne.r) then
72            call abort_gcm(modname,
73     &      "cpp!=r , but running in Shallow Water mode!!",42)
74          endif
75        endif ! of if (llm.eq.1)
76
77        firstcall=.false.
78      endif ! of if (firstcall)
79
80      if (llm.eq.1) then
81       
82        ! Compute pks(:),pk(:),pkf(:)
83       
84        DO   ij  = 1, ngrid
85          pks(ij) = (cpp/preff) * ps(ij)
86          pk(ij,1) = .5*pks(ij)
87        ENDDO
88       
89        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
90        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
91       
92        ! our work is done, exit routine
93        return
94
95      endif ! of if (llm.eq.1)
96
97!!!! General case:
98     
99      unpl2k    = 1.+ 2.* kappa
100c
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 des coeff. alpha et beta  pour la couche l = llm ..
119c
120      DO     ij      = 1, ngrid
121       alpha(ij,llm) = 0.
122       beta (ij,llm) = 1./ unpl2k
123      ENDDO
124c
125c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
126c
127      DO l = llm -1 , 2 , -1
128c
129        DO ij = 1, ngrid
130        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
131        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
132        beta (ij,l)  =   p(ij,l  ) / dellta   
133        ENDDO
134c
135      ENDDO
136c
137c  ***********************************************************************
138c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
139c
140      DO   ij   = 1, ngrid
141       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
142     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
143      ENDDO
144c
145c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
146c
147      DO l = 2, llm
148        DO   ij   = 1, ngrid
149         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
150        ENDDO
151      ENDDO
152c
153c
154      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
155      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
156     
157
158      RETURN
159      END
Note: See TracBrowser for help on using the repository browser.