source: trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F @ 1959

Last change on this file since 1959 was 1959, checked in by emillour, 6 years ago

Common (parallel) dynamics:
Add some initializations of global arrays. This is not a bug fix, strictly
speaking, as results are unchanged. But in some cases computations extend
into halos and these cases are identified as errors with recent versions
of ifort (2018) with the -init=arrays,snan debugging option.
EM

File size: 10.0 KB
Line 
1!
2! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z emillour $
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, ONLY: ij_begin, ij_end, pole_nord, pole_sud,
8     &                         omp_chunk
9      USE control_mod, only : planet_type
10      USE comvert_mod, ONLY: ap,bp
11      USE comconst_mod, ONLY: pi
12      USE logic_mod, ONLY: leapf
13      USE temps_mod, ONLY: dt
14      IMPLICIT NONE
15
16
17c=======================================================================
18c
19c   Auteur:  P. Le Van
20c   -------
21c
22c   objet:
23c   ------
24c
25c   Incrementation des tendances dynamiques
26c
27c=======================================================================
28c-----------------------------------------------------------------------
29c   Declarations:
30c   -------------
31
32#include "dimensions.h"
33#include "paramet.h"
34#include "comgeom.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)=0
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         write(lunout,*) " psm1(ij)=",psm1(stop_it)," dt=",dt,
148     &                   " dp(ij)=",dp(stop_it)
149         call abort_gcm("integrd_p", "negative surface pressure", 1)
150        ENDIF
151
152c
153C$OMP MASTER
154      if (pole_nord) THEN
155     
156        DO  ij    = 1, iim
157         tppn(ij) = aire(   ij   ) * ps(  ij    )
158        ENDDO
159         tpn      = SSUM(iim,tppn,1)/apoln
160        DO ij   = 1, iip1
161         ps(   ij   )  = tpn
162        ENDDO
163     
164      ENDIF
165     
166      if (pole_sud) THEN
167     
168        DO  ij    = 1, iim
169         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
170        ENDDO
171         tps      = SSUM(iim,tpps,1)/apols
172        DO ij   = 1, iip1
173         ps(ij+ip1jm)  = tps
174        ENDDO
175     
176      ENDIF
177c$OMP END MASTER
178c$OMP BARRIER
179c
180c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
181c
182
183      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
184c$OMP BARRIER
185      CALL massdair_p (     p  , masse         )
186
187! Ehouarn : we don't use/need finvmaold and finvmasse,
188!           so might as well not compute them
189!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
190!      ijb=ij_begin
191!      ije=ij_end
192!     
193!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
194!      DO  l = 1,llm
195!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
196!      ENDDO
197!c$OMP END DO NOWAIT
198!
199!      jjb=jj_begin
200!      jje=jj_end
201!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
202c
203
204c    ............   integration  de  ucov, vcov,  h     ..............
205
206c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
207      DO 10 l = 1,llm
208     
209      ijb=ij_begin
210      ije=ij_end
211      if (pole_nord) ijb=ij_begin+iip1
212      if (pole_sud)  ije=ij_end-iip1
213     
214      DO 4 ij = ijb,ije
215      uscr( ij )   =  ucov( ij,l )
216      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
217   4  CONTINUE
218
219      ijb=ij_begin
220      ije=ij_end
221      if (pole_sud)  ije=ij_end-iip1
222     
223      DO 5 ij = ijb,ije
224      vscr( ij )   =  vcov( ij,l )
225      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
226   5  CONTINUE
227     
228      ijb=ij_begin
229      ije=ij_end
230     
231      DO 6 ij = ijb,ije
232      hscr( ij )    =  teta(ij,l)
233      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
234     $                + dt * dteta(ij,l) / masse(ij,l)
235   6  CONTINUE
236
237c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
238c
239c
240      IF (pole_nord) THEN
241       
242        DO  ij   = 1, iim
243          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
244        ENDDO
245          tpn      = SSUM(iim,tppn,1)/apoln
246
247        DO ij   = 1, iip1
248          teta(   ij   ,l)  = tpn
249        ENDDO
250     
251      ENDIF
252     
253      IF (pole_sud) THEN
254       
255        DO  ij   = 1, iim
256          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
257        ENDDO
258          tps      = SSUM(iim,tpps,1)/apols
259
260        DO ij   = 1, iip1
261          teta(ij+ip1jm,l)  = tps
262        ENDDO
263     
264      ENDIF
265c
266
267      IF(leapf)  THEN
268c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
269c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
270c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
271        ijb=ij_begin
272        ije=ij_end
273        ucovm1(ijb:ije,l)=uscr(ijb:ije)
274        tetam1(ijb:ije,l)=hscr(ijb:ije)
275        if (pole_sud) ije=ij_end-iip1
276        vcovm1(ijb:ije,l)=vscr(ijb:ije)
277     
278      END IF
279
280  10  CONTINUE
281c$OMP END DO NOWAIT
282
283c
284c   .......  integration de   q   ......
285c
286      ijb=ij_begin
287      ije=ij_end
288
289         if (planet_type.eq."earth") then
290! Earth-specific treatment of first 2 tracers (water)
291c$OMP BARRIER
292c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
293          DO l = 1, llm
294           DO ij = ijb, ije
295            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
296           ENDDO
297          ENDDO
298c$OMP END DO NOWAIT
299c$OMP BARRIER
300
301          CALL qminimum_p( q, nq, deltap )
302c
303c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
304c
305c$OMP BARRIER
306      IF (pole_nord) THEN
307     
308        DO iq = 1, nq
309       
310c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
311          DO l = 1, llm
312 
313             DO ij = 1, iim
314               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
315             ENDDO
316               qpn  =  SSUM(iim,qppn,1)/apoln
317     
318             DO ij = 1, iip1
319               q(   ij   ,l,iq)  = qpn
320             ENDDO   
321 
322          ENDDO
323c$OMP END DO NOWAIT
324
325        ENDDO
326     
327      ENDIF
328
329      IF (pole_sud) THEN
330     
331        DO iq = 1, nq
332
333c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
334          DO l = 1, llm
335 
336             DO ij = 1, iim
337               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
338             ENDDO
339               qps  =  SSUM(iim,qpps,1)/apols
340 
341             DO ij = 1, iip1
342               q(ij+ip1jm,l,iq)  = qps
343             ENDDO   
344 
345          ENDDO
346c$OMP END DO NOWAIT
347
348        ENDDO
349     
350      ENDIF
351     
352! Ehouarn: forget about finvmaold
353!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
354!
355!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
356!      DO l = 1, llm     
357!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
358!      ENDDO
359!c$OMP END DO NOWAIT
360
361      endif ! of if (planet_type.eq."earth")
362
363c
364c
365c     .....   FIN  de l'integration  de   q    .......
366
36715    continue
368
369c$OMP DO SCHEDULE(STATIC)
370      DO ij=ijb,ije 
371        ps0(ij)=ps(ij)
372      ENDDO
373c$OMP END DO NOWAIT
374
375c    .................................................................
376
377
378      IF( leapf )  THEN
379c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
380c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
381c$OMP DO SCHEDULE(STATIC)
382      DO ij=ijb,ije 
383        psm1(ij)=pscr(ij)
384      ENDDO
385c$OMP END DO NOWAIT
386
387c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
388          DO l = 1, llm
389            massem1(ijb:ije,l)=massescr(ijb:ije,l)
390          ENDDO
391c$OMP END DO NOWAIT       
392      END IF
393c$OMP BARRIER
394      RETURN
395      END
Note: See TracBrowser for help on using the repository browser.