source: LMDZ5/branches/LMDZ5_SPLA/libf/dyn3dpar/caldyn_p.F @ 4934

Last change on this file since 4934 was 1987, checked in by Ehouarn Millour, 11 years ago

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
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: 7.5 KB
Line 
1!
2! $Id: caldyn_p.F 1987 2014-02-24 15:05:47Z abarral $
3!
4#undef DEBUG_IO
5!#define DEBUG_IO
6
7      SUBROUTINE caldyn_p
8     $ (itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
9     $  phi,conser,du,dv,dteta,dp,w,pbaru,pbarv,time )
10      USE parallel_lmdz
11      USE Write_Field_p
12     
13      IMPLICIT NONE
14
15!=======================================================================
16!
17!  Auteur :  P. Le Van
18!
19!   Objet:
20!   ------
21!
22!   Calcul des tendances dynamiques.
23!
24! Modif 04/93 F.Forget
25!=======================================================================
26
27!-----------------------------------------------------------------------
28!   0. Declarations:
29!   ----------------
30
31#include "dimensions.h"
32#include "paramet.h"
33#include "comconst.h"
34#include "comvert.h"
35#include "comgeom.h"
36
37!   Arguments:
38!   ----------
39
40      LOGICAL,INTENT(IN) :: conser ! triggers printing some diagnostics
41      INTEGER,INTENT(IN) :: itau ! time step index
42      REAL,INTENT(IN) :: vcov(ip1jm,llm) ! covariant meridional wind
43      REAL,INTENT(IN) :: ucov(ip1jmp1,llm) ! covariant zonal wind
44      REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potential temperature
45      REAL,INTENT(IN) :: ps(ip1jmp1) ! surface pressure
46      REAL,INTENT(IN) :: phis(ip1jmp1) ! geopotential at the surface
47      REAL,INTENT(IN) :: pk(ip1jmp1,llm) ! Exner at mid-layer
48      REAL,INTENT(IN) :: pkf(ip1jmp1,llm) ! filtered Exner
49      REAL,INTENT(IN) :: phi(ip1jmp1,llm) ! geopotential
50      REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass
51      REAL,INTENT(OUT) :: dv(ip1jm,llm) ! tendency on vcov
52      REAL,INTENT(OUT) :: du(ip1jmp1,llm) ! tendency on ucov
53      REAL,INTENT(OUT) :: dteta(ip1jmp1,llm) ! tenddency on teta
54      REAL,INTENT(OUT) :: dp(ip1jmp1) ! tendency on ps
55      REAL,INTENT(OUT) :: w(ip1jmp1,llm) ! vertical velocity
56      REAL,INTENT(OUT) :: pbaru(ip1jmp1,llm) ! mass flux in the zonal direction
57      REAL,INTENT(OUT) :: pbarv(ip1jm,llm) ! mass flux in the meridional direction
58      REAL,INTENT(IN) :: time ! current time
59
60!   Local:
61!   ------
62
63      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
64      REAL,SAVE :: ang(ip1jmp1,llm)
65      REAL,SAVE :: p(ip1jmp1,llmp1)
66      REAL,SAVE :: massebx(ip1jmp1,llm),masseby(ip1jm,llm)
67      REAL,SAVE :: psexbarxy(ip1jm)
68      REAL,SAVE :: vorpot(ip1jm,llm)
69      REAL,SAVE :: ecin(ip1jmp1,llm)
70      REAL,SAVE :: bern(ip1jmp1,llm)
71      REAL,SAVE :: massebxy(ip1jm,llm)
72      REAL,SAVE :: convm(ip1jmp1,llm)
73      INTEGER   ij,l,ijb,ije,ierr
74
75!-----------------------------------------------------------------------
76!   Compute dynamical tendencies:
77!--------------------------------
78
79      ! compute contravariant winds ucont() and vcont
80      CALL covcont_p  ( llm    , ucov    , vcov , ucont, vcont        )
81      ! compute pressure p()
82      CALL pression_p ( ip1jmp1, ap      , bp   ,  ps  , p            )
83!ym      CALL psextbar (   ps   , psexbarxy                          )
84!$OMP BARRIER
85      ! compute mass in each atmospheric mesh: masse()
86      CALL massdair_p (    p   , masse                                )
87      ! compute X and Y-averages of mass, massebx() and masseby()
88      CALL massbar_p  (   masse, massebx , masseby                    )
89      ! compute XY-average of mass, massebxy()
90      call massbarxy_p(   masse, massebxy                             )
91      ! compute mass fluxes pbaru() and pbarv()
92      CALL flumass_p  ( massebx, masseby , vcont, ucont ,pbaru, pbarv )
93      ! compute dteta() , horizontal converging flux of theta
94      CALL dteta1_p   (   teta , pbaru   , pbarv, dteta               )
95      ! compute convm(), horizontal converging flux of mass
96      CALL convmas1_p  (   pbaru, pbarv   , convm                      )
97!$OMP BARRIER     
98      CALL convmas2_p  (   convm                      )
99!$OMP BARRIER
100#ifdef DEBUG_IO
101!$OMP BARRIER
102!$OMP MASTER
103      call WriteField_p('ucont',reshape(ucont,(/iip1,jmp1,llm/)))
104      call WriteField_p('vcont',reshape(vcont,(/iip1,jjm,llm/)))
105      call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
106      call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
107      call WriteField_p('massebx',reshape(massebx,(/iip1,jmp1,llm/)))
108      call WriteField_p('masseby',reshape(masseby,(/iip1,jjm,llm/)))
109      call WriteField_p('massebxy',reshape(massebxy,(/iip1,jjm,llm/)))
110      call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
111      call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
112      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
113      call WriteField_p('convm',reshape(convm,(/iip1,jmp1,llm/)))
114!$OMP END MASTER
115!$OMP BARRIER
116#endif     
117
118!$OMP BARRIER
119!$OMP MASTER
120      ijb=ij_begin
121      ije=ij_end
122      ! compute pressure variation due to mass convergence
123      DO ij =ijb, ije
124         dp( ij ) = convm( ij,1 ) / airesurg( ij )
125      ENDDO
126!$OMP END MASTER
127!$OMP BARRIER
128!$OMP FLUSH
129     
130      ! compute vertical velocity w()
131      CALL vitvert_p ( convm  , w                                  )
132      ! compute potential vorticity vorpot()
133      CALL tourpot_p ( vcov   , ucov  , massebxy  , vorpot         )
134      ! compute rotation induced du() and dv()
135      CALL dudv1_p   ( vorpot , pbaru , pbarv     , du     , dv    )
136
137#ifdef DEBUG_IO     
138!$OMP BARRIER
139!$OMP MASTER
140      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
141      call WriteField_p('vorpot',reshape(vorpot,(/iip1,jjm,llm/)))
142      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
143      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
144!$OMP END MASTER
145!$OMP BARRIER
146#endif     
147     
148      ! compute kinetic energy ecin()
149      CALL enercin_p ( vcov   , ucov  , vcont     , ucont  , ecin  )
150      ! compute Bernouilli function bern()
151      CALL bernoui_p ( ip1jmp1, llm   , phi       , ecin   , bern  )
152      ! compute and add du() and dv() contributions from Bernouilli and pressure
153      CALL dudv2_p   ( teta   , pkf   , bern      , du     , dv    )
154
155#ifdef DEBUG_IO
156!$OMP BARRIER
157!$OMP MASTER
158      call WriteField_p('ecin',reshape(ecin,(/iip1,jmp1,llm/)))
159      call WriteField_p('bern',reshape(bern,(/iip1,jmp1,llm/)))
160      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
161      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
162      call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
163!$OMP END MASTER
164!$OMP BARRIER
165#endif
166     
167      ijb=ij_begin-iip1
168      ije=ij_end+iip1
169     
170      if (pole_nord) ijb=ij_begin
171      if (pole_sud) ije=ij_end
172
173!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
174      DO l=1,llm
175         DO ij=ijb,ije
176            ang(ij,l) = ucov(ij,l) + constang(ij)
177        ENDDO
178      ENDDO
179!$OMP END DO
180
181      ! compute vertical advection contributions to du(), dv() and dteta()
182      CALL advect_new_p(ang,vcov,teta,w,massebx,masseby,du,dv,dteta)
183
184!  WARNING probleme de peridocite de dv sur les PC/linux. Pb d'arrondi
185!          probablement. Observe sur le code compile avec pgf90 3.0-1
186      ijb=ij_begin
187      ije=ij_end
188      if (pole_sud) ije=ij_end-iip1
189
190!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
191      DO l = 1, llm
192         DO ij = ijb, ije, iip1
193           IF( dv(ij,l).NE.dv(ij+iim,l) )  THEN
194!         PRINT *,'!!!ATTENTION!!! probleme de periodicite sur vcov', 
195!    ,   ' dans caldyn'
196!         PRINT *,' l,  ij = ', l, ij, ij+iim,dv(ij+iim,l),dv(ij,l)
197          dv(ij+iim,l) = dv(ij,l)
198          endif
199         enddo
200      enddo
201!$OMP END DO NOWAIT     
202!-----------------------------------------------------------------------
203!   Output some control variables:
204!---------------------------------
205
206      IF( conser )  THEN
207! ym ---> exige communication collective ( aussi dans advect)
208        CALL sortvarc
209     & ( itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time,vcov )
210
211      ENDIF
212
213      END
Note: See TracBrowser for help on using the repository browser.