source: LMDZ5/trunk/libf/dyn3dpar/caldyn_p.F @ 3024

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