source: LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p_m.F90 @ 3667

Last change on this file since 3667 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: 5.2 KB
RevLine 
[2021]1module exner_hyb_p_m
[630]2
[2021]3  IMPLICIT NONE
[630]4
[2021]5contains
[630]6
[2021]7  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
[630]8
[2021]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    USE parallel_lmdz
[2597]35    USE comconst_mod, ONLY: cpp, kappa, r, jmp1
[2600]36    USE comvert_mod, ONLY: preff
[2021]37    !
38    include "dimensions.h"
39    include "paramet.h"
40    include "comgeom.h"
[1520]41
[2021]42    INTEGER  ngrid
43    REAL p(ngrid,llmp1),pk(ngrid,llm)
44    REAL, optional:: pkf(ngrid,llm)
45    REAL ps(ngrid),pks(ngrid)
46    REAL alpha(ngrid,llm),beta(ngrid,llm)
47
48    !    .... variables locales   ...
49
50    INTEGER l, ij
51    REAL unpl2k,dellta
52
53    INTEGER ije,ijb,jje,jjb
54    logical,save :: firstcall=.true.
55    !$OMP THREADPRIVATE(firstcall)
56    character(len=*),parameter :: modname="exner_hyb_p"
57
58    ! Sanity check
59    if (firstcall) then
60       ! sanity checks for Shallow Water case (1 vertical layer)
61       if (llm.eq.1) then
[1520]62          if (kappa.ne.1) then
[2021]63             call abort_gcm(modname, &
64                  "kappa!=1 , but running in Shallow Water mode!!",42)
[1520]65          endif
66          if (cpp.ne.r) then
[2021]67             call abort_gcm(modname, &
68                  "cpp!=r , but running in Shallow Water mode!!",42)
[1520]69          endif
[2021]70       endif ! of if (llm.eq.1)
[1520]71
[2021]72       firstcall=.false.
73    endif ! of if (firstcall)
[1520]74
[2021]75    !$OMP BARRIER
[1403]76
[2021]77    ! Specific behaviour for Shallow Water (1 vertical layer) case:
78    if (llm.eq.1) then
[1403]79
[2021]80       ! Compute pks(:),pk(:),pkf(:)
81       ijb=ij_begin
82       ije=ij_end
83       !$OMP DO SCHEDULE(STATIC)
84       DO ij=ijb, ije
85          pks(ij) = (cpp/preff) * ps(ij)
[1403]86          pk(ij,1) = .5*pks(ij)
[2021]87          if (present(pkf)) pkf(ij,1)=pk(ij,1)
88       ENDDO
89       !$OMP ENDDO
[1403]90
[2021]91       !$OMP BARRIER
92       if (present(pkf)) then
93          jjb=jj_begin
94          jje=jj_end
95          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
96       end if
[1403]97
[2021]98       ! our work is done, exit routine
99       return
100    endif ! of if (llm.eq.1)
[1403]101
[2021]102    ! General case:
[774]103
[2021]104    unpl2k    = 1.+ 2.* kappa
[630]105
[2021]106    !     -------------
107    !     Calcul de pks
108    !     -------------
109
110    ijb=ij_begin
111    ije=ij_end
112
113    !$OMP DO SCHEDULE(STATIC)
114    DO   ij  = ijb, ije
115       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
116    ENDDO
117    !$OMP ENDDO
118    ! Synchro OPENMP ici
119
120    !$OMP BARRIER
121    !
122    !
123    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
124    !
125    !$OMP DO SCHEDULE(STATIC)
126    DO     ij      = ijb,ije
[630]127       alpha(ij,llm) = 0.
128       beta (ij,llm) = 1./ unpl2k
[2021]129    ENDDO
130    !$OMP ENDDO NOWAIT
131    !
132    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
133    !
134    DO l = llm -1 , 2 , -1
135       !
136       !$OMP DO SCHEDULE(STATIC)
137       DO ij = ijb, ije
138          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
139          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
140          beta (ij,l)  =   p(ij,l  ) / dellta   
141       ENDDO
142       !$OMP ENDDO NOWAIT
143    ENDDO
[774]144
[2021]145    !  ***********************************************************************
146    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
147    !
148    !$OMP DO SCHEDULE(STATIC)
149    DO   ij   = ijb, ije
150       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
151            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
152    ENDDO
153    !$OMP ENDDO NOWAIT
154    !
155    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
156    !
157    DO l = 2, llm
158       !$OMP DO SCHEDULE(STATIC)
159       DO   ij   = ijb, ije
160          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
161       ENDDO
162       !$OMP ENDDO NOWAIT       
163    ENDDO
[985]164
[2021]165    if (present(pkf)) then
166       !    calcul de pkf
[630]167
[2021]168       DO l = 1, llm
169          !$OMP DO SCHEDULE(STATIC)
170          DO   ij   = ijb, ije
171             pkf(ij,l)=pk(ij,l)
172          ENDDO
173          !$OMP ENDDO NOWAIT             
174       ENDDO
175
176       !$OMP BARRIER
177
178       jjb=jj_begin
179       jje=jj_end
180       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
181    end if
182
183  END SUBROUTINE exner_hyb_p
184
185end module exner_hyb_p_m
Note: See TracBrowser for help on using the repository browser.