source: LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90 @ 2965

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