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

Last change on this file since 2146 was 2110, checked in by lguez, 10 years ago

Abort if surface pressure becomes negative. The call to abort_gcm
already was in the sequential version but not in dyn3dpar nor in
dyn3dmem. In dyn3dmem, there is a variable checksum_all, which is
always true. The call to MPI_ALLREDUCE which should update
checksum_all is commented out. I do not know why (performance?).

Non-ASCII characters in comments are not always rendered properly and
they risk being lost. See revision [1740].

Bug fix in procedure convect2: the dimension len must be declared
before the array idcum which has this dimension. Bug fix in procedure
icefrac: the dimensions nl and len must be declared before the arrays
which have these dimensions.
  • 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: 9.6 KB
Line 
1!
2! $Id: integrd_p.F 2110 2014-08-27 15:54:44Z 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_lmdz
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: ps = ", ps(stop_it)
140         write(lunout,*) " at node ij =", stop_it
141         ! since ij=j+(i-1)*jjp1 , we have
142         j=modulo(stop_it,jjp1)
143         i=1+(stop_it-j)/jjp1
144         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
145     &                   " lat = ",rlatu(j)*180./pi, " deg"
146         call abort_gcm("integrd_p", "negative surface pressure", 1)
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.