source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3dpar/integrd_p.F @ 4493

Last change on this file since 4493 was 1550, checked in by lguez, 13 years ago

Bug fix in "bilan_dyn_p". The index was out of bounds in the removed
assignment . Also, the removed assignment was useless.

Bug fix in "coefkzmin". The size of a dummy array cannot exceed the
size of the associated actual array. ("coefkzmin" is called by
"coef_diff_turb".) "km(:, klev+1)" and "kn(:, klev+1)" were not
defined in "coefkzmin" so this was maybe an innocuous bug.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.5 KB
Line 
1!
2! $Id: integrd_p.F 1550 2011-07-05 09:44:55Z crisi $
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
36c   Arguments:
37c   ----------
38
39      INTEGER nq
40
41      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
42      REAL q(ip1jmp1,llm,nq)
43      REAL ps0(ip1jmp1),masse(ip1jmp1,llm),phis(ip1jmp1)
44
45      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
46      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1),massem1(ip1jmp1,llm)
47
48      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
49      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
50      REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm)
51
52c   Local:
53c   ------
54
55      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
56      REAL massescr( ip1jmp1,llm ), finvmasse(ip1jmp1,llm)
57      REAL,SAVE :: p(ip1jmp1,llmp1)
58      REAL tpn,tps,tppn(iim),tpps(iim)
59      REAL qpn,qps,qppn(iim),qpps(iim)
60      REAL,SAVE :: deltap( ip1jmp1,llm )
61
62      INTEGER  l,ij,iq
63
64      REAL SSUM
65      EXTERNAL SSUM
66      INTEGER ijb,ije,jjb,jje
67      REAL,SAVE :: ps(ip1jmp1)
68      LOGICAL :: checksum
69      INTEGER :: stop_it
70c-----------------------------------------------------------------------
71c$OMP BARRIER     
72      if (pole_nord) THEN
73c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
74        DO  l = 1,llm
75          DO  ij = 1,iip1
76           ucov(    ij    , l) = 0.
77           uscr(     ij      ) = 0.
78           ENDDO
79        ENDDO
80c$OMP END DO NOWAIT       
81      ENDIF
82
83      if (pole_sud) THEN
84c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
85        DO  l = 1,llm
86          DO  ij = 1,iip1
87           ucov( ij +ip1jm, l) = 0.
88           uscr( ij +ip1jm   ) = 0.
89          ENDDO
90        ENDDO
91c$OMP END DO NOWAIT     
92      ENDIF
93
94c    ............    integration  de       ps         ..............
95
96c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
97
98      ijb=ij_begin
99      ije=ij_end
100c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
101      DO  l = 1,llm
102        massescr(ijb:ije,l)=masse(ijb:ije,l)
103      ENDDO
104c$OMP END DO NOWAIT
105
106c$OMP DO SCHEDULE(STATIC)
107      DO 2 ij = ijb,ije
108       pscr (ij)    = ps0(ij)
109       ps (ij)      = psm1(ij) + dt * dp(ij)
110   2  CONTINUE
111c$OMP END DO 
112c$OMP BARRIER
113c --> ici synchro OPENMP pour ps
114       
115      checksum=.TRUE.
116      stop_it=0
117
118c$OMP DO SCHEDULE(STATIC)
119      DO ij = ijb,ije
120         IF( ps(ij).LT.0. ) THEN
121           IF (checksum) stop_it=ij
122           checksum=.FALSE.
123         ENDIF
124       ENDDO
125c$OMP END DO NOWAIT
126       
127        IF( .NOT. checksum ) THEN
128         PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. '
129     &         , ps(stop_it)
130         print *, ' dans integrd'
131         stop 1
132        ENDIF
133
134c
135C$OMP MASTER
136      if (pole_nord) THEN
137     
138        DO  ij    = 1, iim
139         tppn(ij) = aire(   ij   ) * ps(  ij    )
140        ENDDO
141         tpn      = SSUM(iim,tppn,1)/apoln
142        DO ij   = 1, iip1
143         ps(   ij   )  = tpn
144        ENDDO
145     
146      ENDIF
147     
148      if (pole_sud) THEN
149     
150        DO  ij    = 1, iim
151         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
152        ENDDO
153         tps      = SSUM(iim,tpps,1)/apols
154        DO ij   = 1, iip1
155         ps(ij+ip1jm)  = tps
156        ENDDO
157     
158      ENDIF
159c$OMP END MASTER
160c$OMP BARRIER
161c
162c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
163c
164
165      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
166c$OMP BARRIER
167      CALL massdair_p (     p  , masse         )
168
169c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
170      ijb=ij_begin
171      ije=ij_end
172     
173c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
174      DO  l = 1,llm
175        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
176      ENDDO
177c$OMP END DO NOWAIT
178
179      jjb=jj_begin
180      jje=jj_end
181      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
182c
183
184c    ............   integration  de  ucov, vcov,  h     ..............
185
186c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
187      DO 10 l = 1,llm
188     
189      ijb=ij_begin
190      ije=ij_end
191      if (pole_nord) ijb=ij_begin+iip1
192      if (pole_sud)  ije=ij_end-iip1
193     
194      DO 4 ij = ijb,ije
195      uscr( ij )   =  ucov( ij,l )
196      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
197   4  CONTINUE
198
199      ijb=ij_begin
200      ije=ij_end
201      if (pole_sud)  ije=ij_end-iip1
202     
203      DO 5 ij = ijb,ije
204      vscr( ij )   =  vcov( ij,l )
205      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
206   5  CONTINUE
207     
208      ijb=ij_begin
209      ije=ij_end
210     
211      DO 6 ij = ijb,ije
212      hscr( ij )    =  teta(ij,l)
213      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
214     $                + dt * dteta(ij,l) / masse(ij,l)
215   6  CONTINUE
216
217c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
218c
219c
220      IF (pole_nord) THEN
221       
222        DO  ij   = 1, iim
223          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
224        ENDDO
225          tpn      = SSUM(iim,tppn,1)/apoln
226
227        DO ij   = 1, iip1
228          teta(   ij   ,l)  = tpn
229        ENDDO
230     
231      ENDIF
232     
233      IF (pole_sud) THEN
234       
235        DO  ij   = 1, iim
236          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
237        ENDDO
238          tps      = SSUM(iim,tpps,1)/apols
239
240        DO ij   = 1, iip1
241          teta(ij+ip1jm,l)  = tps
242        ENDDO
243     
244      ENDIF
245c
246
247      IF(leapf)  THEN
248c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
249c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
250c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
251        ijb=ij_begin
252        ije=ij_end
253        ucovm1(ijb:ije,l)=uscr(ijb:ije)
254        tetam1(ijb:ije,l)=hscr(ijb:ije)
255        if (pole_sud) ije=ij_end-iip1
256        vcovm1(ijb:ije,l)=vscr(ijb:ije)
257     
258      END IF
259
260  10  CONTINUE
261c$OMP END DO NOWAIT
262
263c
264c   .......  integration de   q   ......
265c
266      ijb=ij_begin
267      ije=ij_end
268
269         if (planet_type.eq."earth") then
270! Earth-specific treatment of first 2 tracers (water)
271c$OMP BARRIER
272c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
273          DO l = 1, llm
274           DO ij = ijb, ije
275            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
276           ENDDO
277          ENDDO
278c$OMP END DO NOWAIT
279c$OMP BARRIER
280
281          CALL qminimum_p( q, nq, deltap )
282c
283c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
284c
285c$OMP BARRIER
286      IF (pole_nord) THEN
287     
288        DO iq = 1, nq
289       
290c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
291          DO l = 1, llm
292 
293             DO ij = 1, iim
294               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
295             ENDDO
296               qpn  =  SSUM(iim,qppn,1)/apoln
297     
298             DO ij = 1, iip1
299               q(   ij   ,l,iq)  = qpn
300             ENDDO   
301 
302          ENDDO
303c$OMP END DO NOWAIT
304
305        ENDDO
306     
307      ENDIF
308
309      IF (pole_sud) THEN
310     
311        DO iq = 1, nq
312
313c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
314          DO l = 1, llm
315 
316             DO ij = 1, iim
317               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
318             ENDDO
319               qps  =  SSUM(iim,qpps,1)/apols
320 
321             DO ij = 1, iip1
322               q(ij+ip1jm,l,iq)  = qps
323             ENDDO   
324 
325          ENDDO
326c$OMP END DO NOWAIT
327
328        ENDDO
329     
330      ENDIF
331     
332c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
333
334c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
335      DO l = 1, llm     
336        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
337      ENDDO
338c$OMP END DO NOWAIT
339
340      endif ! of if (planet_type.eq."earth")
341
342c
343c
344c     .....   FIN  de l'integration  de   q    .......
345
34615    continue
347
348c$OMP DO SCHEDULE(STATIC)
349      DO ij=ijb,ije 
350        ps0(ij)=ps(ij)
351      ENDDO
352c$OMP END DO NOWAIT
353
354c    .................................................................
355
356
357      IF( leapf )  THEN
358c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
359c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
360c$OMP DO SCHEDULE(STATIC)
361      DO ij=ijb,ije 
362        psm1(ij)=pscr(ij)
363      ENDDO
364c$OMP END DO NOWAIT
365
366c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
367          DO l = 1, llm
368            massem1(ijb:ije,l)=massescr(ijb:ije,l)
369          ENDDO
370c$OMP END DO NOWAIT       
371      END IF
372c$OMP BARRIER
373      RETURN
374      END
Note: See TracBrowser for help on using the repository browser.