source: LMDZ6/trunk/libf/dyn3d_common/exner_hyb_m.F90 @ 3741

Last change on this file since 3741 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.3 KB
Line 
1module exner_hyb_m
2
3  IMPLICIT NONE
4
5contains
6
7  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
8
9    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
10    !    ..........
11    !
12    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
13    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
14    !
15    !   ************************************************************************
16    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
17    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
18    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
19    !   ************************************************************************
20    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
21    !    la pression et la fonction d'Exner  au  sol  .
22    !
23    !                                 -------- z
24    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
25    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
26    !    ( voir note de Fr.Hourdin )  ,
27    !
28    !    on determine successivement , du haut vers le bas des couches, les
29    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
30    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
31    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
32    !
33    !
34    !
35    USE comconst_mod, ONLY: jmp1, cpp, kappa, r
36    USE comvert_mod, ONLY: preff
37   
38    IMPLICIT NONE
39   
40    include "dimensions.h"
41    include "paramet.h"
42    include "comgeom.h"
43
44    INTEGER  ngrid
45    REAL p(ngrid,llmp1),pk(ngrid,llm)
46    real, optional:: pkf(ngrid,llm)
47    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
48
49    !    .... variables locales   ...
50
51    INTEGER l, ij
52    REAL unpl2k,dellta
53
54    logical,save :: firstcall=.true.
55    character(len=*),parameter :: modname="exner_hyb"
56
57    ! Sanity check
58    if (firstcall) then
59       ! sanity checks for Shallow Water case (1 vertical layer)
60       if (llm.eq.1) then
61          if (kappa.ne.1) then
62             call abort_gcm(modname, &
63                  "kappa!=1 , but running in Shallow Water mode!!",42)
64          endif
65          if (cpp.ne.r) then
66             call abort_gcm(modname, &
67                  "cpp!=r , but running in Shallow Water mode!!",42)
68          endif
69       endif ! of if (llm.eq.1)
70
71       firstcall=.false.
72    endif ! of if (firstcall)
73
74    ! Specific behaviour for Shallow Water (1 vertical layer) case:
75    if (llm.eq.1) then
76
77       ! Compute pks(:),pk(:),pkf(:)
78
79       DO   ij  = 1, ngrid
80          pks(ij) = (cpp/preff) * ps(ij)
81          pk(ij,1) = .5*pks(ij)
82       ENDDO
83
84       if (present(pkf)) then
85          pkf = pk
86          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
87       end if
88
89       ! our work is done, exit routine
90       return
91    endif ! of if (llm.eq.1)
92
93    ! General case:
94
95    unpl2k    = 1.+ 2.* kappa
96
97    !     -------------
98    !     Calcul de pks
99    !     -------------
100
101    DO   ij  = 1, ngrid
102       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
103    ENDDO
104
105    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
106    !
107    DO     ij      = 1, ngrid
108       alpha(ij,llm) = 0.
109       beta (ij,llm) = 1./ unpl2k
110    ENDDO
111    !
112    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
113    !
114    DO l = llm -1 , 2 , -1
115       !
116       DO ij = 1, ngrid
117          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
118          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
119          beta (ij,l)  =   p(ij,l  ) / dellta   
120       ENDDO
121    ENDDO
122
123    !  ***********************************************************************
124    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
125    !
126    DO   ij   = 1, ngrid
127       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
128            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
129    ENDDO
130    !
131    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
132    !
133    DO l = 2, llm
134       DO   ij   = 1, ngrid
135          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
136       ENDDO
137    ENDDO
138
139    if (present(pkf)) then
140       !    calcul de pkf
141       pkf = pk
142       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
143    end if
144
145  END SUBROUTINE exner_hyb
146
147end module exner_hyb_m
Note: See TracBrowser for help on using the repository browser.