source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/integrd_p.F @ 3069

Last change on this file since 3069 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

  • 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 1446 2010-10-22 09:27:25Z dcugnet $
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         STOP' dans integrd'
131        ENDIF
132
133c
134C$OMP MASTER
135      if (pole_nord) THEN
136     
137        DO  ij    = 1, iim
138         tppn(ij) = aire(   ij   ) * ps(  ij    )
139        ENDDO
140         tpn      = SSUM(iim,tppn,1)/apoln
141        DO ij   = 1, iip1
142         ps(   ij   )  = tpn
143        ENDDO
144     
145      ENDIF
146     
147      if (pole_sud) THEN
148     
149        DO  ij    = 1, iim
150         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
151        ENDDO
152         tps      = SSUM(iim,tpps,1)/apols
153        DO ij   = 1, iip1
154         ps(ij+ip1jm)  = tps
155        ENDDO
156     
157      ENDIF
158c$OMP END MASTER
159c$OMP BARRIER
160c
161c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
162c
163
164      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
165c$OMP BARRIER
166      CALL massdair_p (     p  , masse         )
167
168c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
169      ijb=ij_begin
170      ije=ij_end
171     
172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173      DO  l = 1,llm
174        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
175      ENDDO
176c$OMP END DO NOWAIT
177
178      jjb=jj_begin
179      jje=jj_end
180      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
181c
182
183c    ............   integration  de  ucov, vcov,  h     ..............
184
185c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
186      DO 10 l = 1,llm
187     
188      ijb=ij_begin
189      ije=ij_end
190      if (pole_nord) ijb=ij_begin+iip1
191      if (pole_sud)  ije=ij_end-iip1
192     
193      DO 4 ij = ijb,ije
194      uscr( ij )   =  ucov( ij,l )
195      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
196   4  CONTINUE
197
198      ijb=ij_begin
199      ije=ij_end
200      if (pole_sud)  ije=ij_end-iip1
201     
202      DO 5 ij = ijb,ije
203      vscr( ij )   =  vcov( ij,l )
204      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
205   5  CONTINUE
206     
207      ijb=ij_begin
208      ije=ij_end
209     
210      DO 6 ij = ijb,ije
211      hscr( ij )    =  teta(ij,l)
212      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
213     $                + dt * dteta(ij,l) / masse(ij,l)
214   6  CONTINUE
215
216c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
217c
218c
219      IF (pole_nord) THEN
220       
221        DO  ij   = 1, iim
222          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
223        ENDDO
224          tpn      = SSUM(iim,tppn,1)/apoln
225
226        DO ij   = 1, iip1
227          teta(   ij   ,l)  = tpn
228        ENDDO
229     
230      ENDIF
231     
232      IF (pole_sud) THEN
233       
234        DO  ij   = 1, iim
235          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
236        ENDDO
237          tps      = SSUM(iim,tpps,1)/apols
238
239        DO ij   = 1, iip1
240          teta(ij+ip1jm,l)  = tps
241        ENDDO
242     
243      ENDIF
244c
245
246      IF(leapf)  THEN
247c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
248c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
249c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
250        ijb=ij_begin
251        ije=ij_end
252        ucovm1(ijb:ije,l)=uscr(ijb:ije)
253        tetam1(ijb:ije,l)=hscr(ijb:ije)
254        if (pole_sud) ije=ij_end-iip1
255        vcovm1(ijb:ije,l)=vscr(ijb:ije)
256     
257      END IF
258
259  10  CONTINUE
260c$OMP END DO NOWAIT
261
262c
263c   .......  integration de   q   ......
264c
265      ijb=ij_begin
266      ije=ij_end
267
268         if (planet_type.eq."earth") then
269! Earth-specific treatment of first 2 tracers (water)
270c$OMP BARRIER
271c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
272          DO l = 1, llm
273           DO ij = ijb, ije
274            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
275           ENDDO
276          ENDDO
277c$OMP END DO NOWAIT
278c$OMP BARRIER
279
280          CALL qminimum_p( q, nq, deltap )
281c
282c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
283c
284c$OMP BARRIER
285      IF (pole_nord) THEN
286     
287        DO iq = 1, nq
288       
289c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
290          DO l = 1, llm
291 
292             DO ij = 1, iim
293               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
294             ENDDO
295               qpn  =  SSUM(iim,qppn,1)/apoln
296     
297             DO ij = 1, iip1
298               q(   ij   ,l,iq)  = qpn
299             ENDDO   
300 
301          ENDDO
302c$OMP END DO NOWAIT
303
304        ENDDO
305     
306      ENDIF
307
308      IF (pole_sud) THEN
309     
310        DO iq = 1, nq
311
312c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
313          DO l = 1, llm
314 
315             DO ij = 1, iim
316               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
317             ENDDO
318               qps  =  SSUM(iim,qpps,1)/apols
319 
320             DO ij = 1, iip1
321               q(ij+ip1jm,l,iq)  = qps
322             ENDDO   
323 
324          ENDDO
325c$OMP END DO NOWAIT
326
327        ENDDO
328     
329      ENDIF
330     
331c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
332
333c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
334      DO l = 1, llm     
335        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
336      ENDDO
337c$OMP END DO NOWAIT
338
339      endif ! of if (planet_type.eq."earth")
340
341c
342c
343c     .....   FIN  de l'integration  de   q    .......
344
34515    continue
346
347c$OMP DO SCHEDULE(STATIC)
348      DO ij=ijb,ije 
349        ps0(ij)=ps(ij)
350      ENDDO
351c$OMP END DO NOWAIT
352
353c    .................................................................
354
355
356      IF( leapf )  THEN
357c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
358c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
359c$OMP DO SCHEDULE(STATIC)
360      DO ij=ijb,ije 
361        psm1(ij)=pscr(ij)
362      ENDDO
363c$OMP END DO NOWAIT
364
365c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
366          DO l = 1, llm
367            massem1(ijb:ije,l)=massescr(ijb:ije,l)
368          ENDDO
369c$OMP END DO NOWAIT       
370      END IF
371c$OMP BARRIER
372      RETURN
373      END
Note: See TracBrowser for help on using the repository browser.