source: LMDZ6/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90

Last change on this file 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
Line 
1module exner_hyb_loc_m
2
3  IMPLICIT NONE
4
5contains
6
7  SUBROUTINE  exner_hyb_loc(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    USE parallel_lmdz
35    USE mod_filtreg_p
36    USE write_field_loc
37    USE comconst_mod, ONLY: cpp, kappa, r, jmp1
38    USE comvert_mod, ONLY: preff
39   
40    IMPLICIT NONE
41    !
42    include "dimensions.h"
43    include "paramet.h"
44    include "comgeom.h"
45
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
68          if (kappa.ne.1) then
69             call abort_gcm(modname, &
70                  "kappa!=1 , but running in Shallow Water mode!!",42)
71          endif
72          if (cpp.ne.r) then
73             call abort_gcm(modname, &
74                  "cpp!=r , but running in Shallow Water mode!!",42)
75          endif
76       endif ! of if (llm.eq.1)
77
78       firstcall=.false.
79    endif ! of if (firstcall)
80
81    !$OMP BARRIER
82
83    ! Specific behaviour for Shallow Water (1 vertical layer) case:
84    if (llm.eq.1) then
85
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)
92          pk(ij,1) = .5*pks(ij)
93          if (present(pkf)) pkf(ij,1)=pk(ij,1)
94       ENDDO
95       !$OMP ENDDO
96
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
104
105       ! our work is done, exit routine
106       return
107    endif ! of if (llm.eq.1)
108
109    ! General case:
110
111    unpl2k    = 1.+ 2.* kappa
112
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
134       alpha(ij,llm) = 0.
135       beta (ij,llm) = 1./ unpl2k
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
151
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
171
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
187#ifdef DEBUG_IO   
188       call WriteField_u('pkf',pkf)
189#endif
190       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
191            2, 1, .TRUE., 1 )
192#ifdef DEBUG_IO   
193       call WriteField_u('pkf',pkf)
194#endif     
195    end if
196
197  END SUBROUTINE exner_hyb_loc
198
199end module exner_hyb_loc_m
Note: See TracBrowser for help on using the repository browser.