source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dpar/integrd_p.F @ 4785

Last change on this file since 4785 was 1616, checked in by Ehouarn Millour, 13 years ago

Some cleanup around what is done during the integration step of dynamical tendencies and namely removed computation of (unused) finvmaold, thereby saving us the expense of a call to the (costly) filter at every dynamical time step.
Checked (on Vargas, in seq, omp, mpi and mixed mode) that this doesn't change the GCM results, as expected.
EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
1!
2! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z idelkadi $
3!
4      SUBROUTINE integrd_p
5     $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
7      USE parallel
8      USE control_mod, only : planet_type
9      IMPLICIT NONE
10
11
12c=======================================================================
13c
14c   Auteur:  P. Le Van
15c   -------
16c
17c   objet:
18c   ------
19c
20c   Incrementation des tendances dynamiques
21c
22c=======================================================================
23c-----------------------------------------------------------------------
24c   Declarations:
25c   -------------
26
27#include "dimensions.h"
28#include "paramet.h"
29#include "comconst.h"
30#include "comgeom.h"
31#include "comvert.h"
32#include "logic.h"
33#include "temps.h"
34#include "serre.h"
35#include "iniprint.h"
36
37c   Arguments:
38c   ----------
39
40      integer,intent(in) :: nq ! number of tracers to handle in this routine
41      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
42      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
43      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
44      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
45      real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
46      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
47      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
48      ! values at previous time step
49      real,intent(inout) :: vcovm1(ip1jm,llm)
50      real,intent(inout) :: ucovm1(ip1jmp1,llm)
51      real,intent(inout) :: tetam1(ip1jmp1,llm)
52      real,intent(inout) :: psm1(ip1jmp1)
53      real,intent(inout) :: massem1(ip1jmp1,llm)
54      ! the tendencies to add
55      real,intent(in) :: dv(ip1jm,llm)
56      real,intent(in) :: du(ip1jmp1,llm)
57      real,intent(in) :: dteta(ip1jmp1,llm)
58      real,intent(in) :: dp(ip1jmp1)
59      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
60!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
61
62c   Local:
63c   ------
64
65      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
66      REAL massescr( ip1jmp1,llm )
67!      REAL finvmasse(ip1jmp1,llm)
68      REAL,SAVE :: p(ip1jmp1,llmp1)
69      REAL tpn,tps,tppn(iim),tpps(iim)
70      REAL qpn,qps,qppn(iim),qpps(iim)
71      REAL,SAVE :: deltap( ip1jmp1,llm )
72
73      INTEGER  l,ij,iq,i,j
74
75      REAL SSUM
76      EXTERNAL SSUM
77      INTEGER ijb,ije,jjb,jje
78      REAL,SAVE :: ps(ip1jmp1)
79      LOGICAL :: checksum
80      INTEGER :: stop_it
81c-----------------------------------------------------------------------
82c$OMP BARRIER     
83      if (pole_nord) THEN
84c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
85        DO  l = 1,llm
86          DO  ij = 1,iip1
87           ucov(    ij    , l) = 0.
88           uscr(     ij      ) = 0.
89           ENDDO
90        ENDDO
91c$OMP END DO NOWAIT       
92      ENDIF
93
94      if (pole_sud) THEN
95c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
96        DO  l = 1,llm
97          DO  ij = 1,iip1
98           ucov( ij +ip1jm, l) = 0.
99           uscr( ij +ip1jm   ) = 0.
100          ENDDO
101        ENDDO
102c$OMP END DO NOWAIT     
103      ENDIF
104
105c    ............    integration  de       ps         ..............
106
107c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
108
109      ijb=ij_begin
110      ije=ij_end
111c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
112      DO  l = 1,llm
113        massescr(ijb:ije,l)=masse(ijb:ije,l)
114      ENDDO
115c$OMP END DO NOWAIT
116
117c$OMP DO SCHEDULE(STATIC)
118      DO 2 ij = ijb,ije
119       pscr (ij)    = ps0(ij)
120       ps (ij)      = psm1(ij) + dt * dp(ij)
121   2  CONTINUE
122c$OMP END DO 
123c$OMP BARRIER
124c --> ici synchro OPENMP pour ps
125       
126      checksum=.TRUE.
127      stop_it=0
128
129c$OMP DO SCHEDULE(STATIC)
130      DO ij = ijb,ije
131         IF( ps(ij).LT.0. ) THEN
132           IF (checksum) stop_it=ij
133           checksum=.FALSE.
134         ENDIF
135       ENDDO
136c$OMP END DO NOWAIT
137       
138        IF( .NOT. checksum ) THEN
139         write(lunout,*) "integrd: negative surface pressure ",
140     &                                                ps(stop_it)
141         write(lunout,*) " at node ij =", stop_it
142         ! since ij=j+(i-1)*jjp1 , we have
143         j=modulo(stop_it,jjp1)
144         i=1+(stop_it-j)/jjp1
145         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
146     &                   " lat = ",rlatu(j)*180./pi, " deg"
147        ENDIF
148
149c
150C$OMP MASTER
151      if (pole_nord) THEN
152     
153        DO  ij    = 1, iim
154         tppn(ij) = aire(   ij   ) * ps(  ij    )
155        ENDDO
156         tpn      = SSUM(iim,tppn,1)/apoln
157        DO ij   = 1, iip1
158         ps(   ij   )  = tpn
159        ENDDO
160     
161      ENDIF
162     
163      if (pole_sud) THEN
164     
165        DO  ij    = 1, iim
166         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
167        ENDDO
168         tps      = SSUM(iim,tpps,1)/apols
169        DO ij   = 1, iip1
170         ps(ij+ip1jm)  = tps
171        ENDDO
172     
173      ENDIF
174c$OMP END MASTER
175c$OMP BARRIER
176c
177c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
178c
179
180      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
181c$OMP BARRIER
182      CALL massdair_p (     p  , masse         )
183
184! Ehouarn : we don't use/need finvmaold and finvmasse,
185!           so might as well not compute them
186!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
187!      ijb=ij_begin
188!      ije=ij_end
189!     
190!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
191!      DO  l = 1,llm
192!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
193!      ENDDO
194!c$OMP END DO NOWAIT
195!
196!      jjb=jj_begin
197!      jje=jj_end
198!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
199c
200
201c    ............   integration  de  ucov, vcov,  h     ..............
202
203c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
204      DO 10 l = 1,llm
205     
206      ijb=ij_begin
207      ije=ij_end
208      if (pole_nord) ijb=ij_begin+iip1
209      if (pole_sud)  ije=ij_end-iip1
210     
211      DO 4 ij = ijb,ije
212      uscr( ij )   =  ucov( ij,l )
213      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
214   4  CONTINUE
215
216      ijb=ij_begin
217      ije=ij_end
218      if (pole_sud)  ije=ij_end-iip1
219     
220      DO 5 ij = ijb,ije
221      vscr( ij )   =  vcov( ij,l )
222      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
223   5  CONTINUE
224     
225      ijb=ij_begin
226      ije=ij_end
227     
228      DO 6 ij = ijb,ije
229      hscr( ij )    =  teta(ij,l)
230      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
231     $                + dt * dteta(ij,l) / masse(ij,l)
232   6  CONTINUE
233
234c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
235c
236c
237      IF (pole_nord) THEN
238       
239        DO  ij   = 1, iim
240          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
241        ENDDO
242          tpn      = SSUM(iim,tppn,1)/apoln
243
244        DO ij   = 1, iip1
245          teta(   ij   ,l)  = tpn
246        ENDDO
247     
248      ENDIF
249     
250      IF (pole_sud) THEN
251       
252        DO  ij   = 1, iim
253          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
254        ENDDO
255          tps      = SSUM(iim,tpps,1)/apols
256
257        DO ij   = 1, iip1
258          teta(ij+ip1jm,l)  = tps
259        ENDDO
260     
261      ENDIF
262c
263
264      IF(leapf)  THEN
265c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
266c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
267c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
268        ijb=ij_begin
269        ije=ij_end
270        ucovm1(ijb:ije,l)=uscr(ijb:ije)
271        tetam1(ijb:ije,l)=hscr(ijb:ije)
272        if (pole_sud) ije=ij_end-iip1
273        vcovm1(ijb:ije,l)=vscr(ijb:ije)
274     
275      END IF
276
277  10  CONTINUE
278c$OMP END DO NOWAIT
279
280c
281c   .......  integration de   q   ......
282c
283      ijb=ij_begin
284      ije=ij_end
285
286         if (planet_type.eq."earth") then
287! Earth-specific treatment of first 2 tracers (water)
288c$OMP BARRIER
289c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
290          DO l = 1, llm
291           DO ij = ijb, ije
292            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
293           ENDDO
294          ENDDO
295c$OMP END DO NOWAIT
296c$OMP BARRIER
297
298          CALL qminimum_p( q, nq, deltap )
299c
300c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
301c
302c$OMP BARRIER
303      IF (pole_nord) THEN
304     
305        DO iq = 1, nq
306       
307c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
308          DO l = 1, llm
309 
310             DO ij = 1, iim
311               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
312             ENDDO
313               qpn  =  SSUM(iim,qppn,1)/apoln
314     
315             DO ij = 1, iip1
316               q(   ij   ,l,iq)  = qpn
317             ENDDO   
318 
319          ENDDO
320c$OMP END DO NOWAIT
321
322        ENDDO
323     
324      ENDIF
325
326      IF (pole_sud) THEN
327     
328        DO iq = 1, nq
329
330c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
331          DO l = 1, llm
332 
333             DO ij = 1, iim
334               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
335             ENDDO
336               qps  =  SSUM(iim,qpps,1)/apols
337 
338             DO ij = 1, iip1
339               q(ij+ip1jm,l,iq)  = qps
340             ENDDO   
341 
342          ENDDO
343c$OMP END DO NOWAIT
344
345        ENDDO
346     
347      ENDIF
348     
349! Ehouarn: forget about finvmaold
350!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
351!
352!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
353!      DO l = 1, llm     
354!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
355!      ENDDO
356!c$OMP END DO NOWAIT
357
358      endif ! of if (planet_type.eq."earth")
359
360c
361c
362c     .....   FIN  de l'integration  de   q    .......
363
36415    continue
365
366c$OMP DO SCHEDULE(STATIC)
367      DO ij=ijb,ije 
368        ps0(ij)=ps(ij)
369      ENDDO
370c$OMP END DO NOWAIT
371
372c    .................................................................
373
374
375      IF( leapf )  THEN
376c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
377c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
378c$OMP DO SCHEDULE(STATIC)
379      DO ij=ijb,ije 
380        psm1(ij)=pscr(ij)
381      ENDDO
382c$OMP END DO NOWAIT
383
384c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
385          DO l = 1, llm
386            massem1(ijb:ije,l)=massescr(ijb:ije,l)
387          ENDDO
388c$OMP END DO NOWAIT       
389      END IF
390c$OMP BARRIER
391      RETURN
392      END
Note: See TracBrowser for help on using the repository browser.