source: LMDZ5/trunk/libf/dyn3dpar/integrd_p.F @ 1680

Last change on this file since 1680 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
RevLine 
[630]1!
[1279]2! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z emillour $
[630]3!
4      SUBROUTINE integrd_p
5     $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
[1616]6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
[630]7      USE parallel
[1454]8      USE control_mod, only : planet_type
[630]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"
[1616]35#include "iniprint.h"
[630]36
37c   Arguments:
38c   ----------
39
[1616]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
[630]61
62c   Local:
63c   ------
64
65      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
[1616]66      REAL massescr( ip1jmp1,llm )
67!      REAL finvmasse(ip1jmp1,llm)
[764]68      REAL,SAVE :: p(ip1jmp1,llmp1)
[630]69      REAL tpn,tps,tppn(iim),tpps(iim)
70      REAL qpn,qps,qppn(iim),qpps(iim)
[985]71      REAL,SAVE :: deltap( ip1jmp1,llm )
[630]72
[1616]73      INTEGER  l,ij,iq,i,j
[630]74
75      REAL SSUM
76      EXTERNAL SSUM
77      INTEGER ijb,ije,jjb,jje
[764]78      REAL,SAVE :: ps(ip1jmp1)
[985]79      LOGICAL :: checksum
80      INTEGER :: stop_it
[630]81c-----------------------------------------------------------------------
[985]82c$OMP BARRIER     
[630]83      if (pole_nord) THEN
[764]84c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]85        DO  l = 1,llm
86          DO  ij = 1,iip1
87           ucov(    ij    , l) = 0.
88           uscr(     ij      ) = 0.
89           ENDDO
90        ENDDO
[764]91c$OMP END DO NOWAIT       
[630]92      ENDIF
93
94      if (pole_sud) THEN
[764]95c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]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
[764]102c$OMP END DO NOWAIT     
[630]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
[764]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
[985]117c$OMP DO SCHEDULE(STATIC)
[630]118      DO 2 ij = ijb,ije
[764]119       pscr (ij)    = ps0(ij)
[630]120       ps (ij)      = psm1(ij) + dt * dp(ij)
121   2  CONTINUE
[985]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)
[630]130      DO ij = ijb,ije
[985]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
[1616]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"
[630]147        ENDIF
[985]148
[630]149c
[985]150C$OMP MASTER
[630]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
[764]174c$OMP END MASTER
175c$OMP BARRIER
[630]176c
177c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
178c
[764]179
[630]180      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[764]181c$OMP BARRIER
[630]182      CALL massdair_p (     p  , masse         )
183
[1616]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  )
[630]199c
200
201c    ............   integration  de  ucov, vcov,  h     ..............
202
[764]203c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]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
[764]278c$OMP END DO NOWAIT
[630]279
280c
281c   .......  integration de   q   ......
282c
283      ijb=ij_begin
284      ije=ij_end
[1279]285
286         if (planet_type.eq."earth") then
287! Earth-specific treatment of first 2 tracers (water)
[985]288c$OMP BARRIER
289c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1279]290          DO l = 1, llm
291           DO ij = ijb, ije
292            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
293           ENDDO
[630]294          ENDDO
[985]295c$OMP END DO NOWAIT
296c$OMP BARRIER
[630]297
[1279]298          CALL qminimum_p( q, nq, deltap )
[630]299c
300c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
301c
[985]302c$OMP BARRIER
[630]303      IF (pole_nord) THEN
304     
305        DO iq = 1, nq
[985]306       
307c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]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
[985]320c$OMP END DO NOWAIT
321
[630]322        ENDDO
323     
324      ENDIF
325
326      IF (pole_sud) THEN
327     
328        DO iq = 1, nq
[985]329
330c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]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
[985]343c$OMP END DO NOWAIT
344
[630]345        ENDDO
346     
347      ENDIF
[764]348     
[1616]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
[764]357
[1454]358      endif ! of if (planet_type.eq."earth")
359
[630]360c
361c
362c     .....   FIN  de l'integration  de   q    .......
363
36415    continue
365
[985]366c$OMP DO SCHEDULE(STATIC)
367      DO ij=ijb,ije 
368        ps0(ij)=ps(ij)
369      ENDDO
370c$OMP END DO NOWAIT
371
[630]372c    .................................................................
373
374
375      IF( leapf )  THEN
376c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
377c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
[985]378c$OMP DO SCHEDULE(STATIC)
379      DO ij=ijb,ije 
380        psm1(ij)=pscr(ij)
381      ENDDO
382c$OMP END DO NOWAIT
[764]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       
[630]389      END IF
[985]390c$OMP BARRIER
[630]391      RETURN
392      END
Note: See TracBrowser for help on using the repository browser.