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

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM
  • 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 2600 2016-07-23 05:45:38Z 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
8      USE control_mod, only : planet_type
9      USE comconst_mod, ONLY: pi
10      USE comvert_mod, ONLY: ap, bp
11      IMPLICIT NONE
12
13
14c=======================================================================
15c
16c   Auteur:  P. Le Van
17c   -------
18c
19c   objet:
20c   ------
21c
22c   Incrementation des tendances dynamiques
23c
24c=======================================================================
25c-----------------------------------------------------------------------
26c   Declarations:
27c   -------------
28
29#include "dimensions.h"
30#include "paramet.h"
31#include "comgeom.h"
32#include "logic.h"
33#include "temps.h"
34#include "iniprint.h"
35
36c   Arguments:
37c   ----------
38
39      integer,intent(in) :: nq ! number of tracers to handle in this routine
40      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
41      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
42      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
43      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
44      real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
45      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
46      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
47      ! values at previous time step
48      real,intent(inout) :: vcovm1(ip1jm,llm)
49      real,intent(inout) :: ucovm1(ip1jmp1,llm)
50      real,intent(inout) :: tetam1(ip1jmp1,llm)
51      real,intent(inout) :: psm1(ip1jmp1)
52      real,intent(inout) :: massem1(ip1jmp1,llm)
53      ! the tendencies to add
54      real,intent(in) :: dv(ip1jm,llm)
55      real,intent(in) :: du(ip1jmp1,llm)
56      real,intent(in) :: dteta(ip1jmp1,llm)
57      real,intent(in) :: dp(ip1jmp1)
58      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
59!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
60
61c   Local:
62c   ------
63
64      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
65      REAL massescr( ip1jmp1,llm )
66!      REAL finvmasse(ip1jmp1,llm)
67      REAL,SAVE :: p(ip1jmp1,llmp1)
68      REAL tpn,tps,tppn(iim),tpps(iim)
69      REAL qpn,qps,qppn(iim),qpps(iim)
70      REAL,SAVE :: deltap( ip1jmp1,llm )
71
72      INTEGER  l,ij,iq,i,j
73
74      REAL SSUM
75      EXTERNAL SSUM
76      INTEGER ijb,ije,jjb,jje
77      REAL,SAVE :: ps(ip1jmp1)
78      LOGICAL :: checksum
79      INTEGER :: stop_it
80c-----------------------------------------------------------------------
81c$OMP BARRIER     
82      if (pole_nord) THEN
83c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
84        DO  l = 1,llm
85          DO  ij = 1,iip1
86           ucov(    ij    , l) = 0.
87           uscr(     ij      ) = 0.
88           ENDDO
89        ENDDO
90c$OMP END DO NOWAIT       
91      ENDIF
92
93      if (pole_sud) THEN
94c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
95        DO  l = 1,llm
96          DO  ij = 1,iip1
97           ucov( ij +ip1jm, l) = 0.
98           uscr( ij +ip1jm   ) = 0.
99          ENDDO
100        ENDDO
101c$OMP END DO NOWAIT     
102      ENDIF
103
104c    ............    integration  de       ps         ..............
105
106c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
107
108      ijb=ij_begin
109      ije=ij_end
110c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
111      DO  l = 1,llm
112        massescr(ijb:ije,l)=masse(ijb:ije,l)
113      ENDDO
114c$OMP END DO NOWAIT
115
116c$OMP DO SCHEDULE(STATIC)
117      DO 2 ij = ijb,ije
118       pscr (ij)    = ps0(ij)
119       ps (ij)      = psm1(ij) + dt * dp(ij)
120   2  CONTINUE
121c$OMP END DO 
122c$OMP BARRIER
123c --> ici synchro OPENMP pour ps
124       
125      checksum=.TRUE.
126      stop_it=0
127
128c$OMP DO SCHEDULE(STATIC)
129      DO ij = ijb,ije
130         IF( ps(ij).LT.0. ) THEN
131           IF (checksum) stop_it=ij
132           checksum=.FALSE.
133         ENDIF
134       ENDDO
135c$OMP END DO NOWAIT
136       
137        IF( .NOT. checksum ) THEN
138           write(lunout,*) "integrd: ps = ", ps(stop_it)
139         write(lunout,*) " at node ij =", stop_it
140         ! since ij=j+(i-1)*jjp1 , we have
141         j=modulo(stop_it,jjp1)
142         i=1+(stop_it-j)/jjp1
143         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
144     &                   " lat = ",rlatu(j)*180./pi, " deg"
145         call abort_gcm("integrd_p", "negative surface pressure", 1)
146        ENDIF
147
148c
149C$OMP MASTER
150      if (pole_nord) THEN
151     
152        DO  ij    = 1, iim
153         tppn(ij) = aire(   ij   ) * ps(  ij    )
154        ENDDO
155         tpn      = SSUM(iim,tppn,1)/apoln
156        DO ij   = 1, iip1
157         ps(   ij   )  = tpn
158        ENDDO
159     
160      ENDIF
161     
162      if (pole_sud) THEN
163     
164        DO  ij    = 1, iim
165         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
166        ENDDO
167         tps      = SSUM(iim,tpps,1)/apols
168        DO ij   = 1, iip1
169         ps(ij+ip1jm)  = tps
170        ENDDO
171     
172      ENDIF
173c$OMP END MASTER
174c$OMP BARRIER
175c
176c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
177c
178
179      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
180c$OMP BARRIER
181      CALL massdair_p (     p  , masse         )
182
183! Ehouarn : we don't use/need finvmaold and finvmasse,
184!           so might as well not compute them
185!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
186!      ijb=ij_begin
187!      ije=ij_end
188!     
189!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
190!      DO  l = 1,llm
191!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
192!      ENDDO
193!c$OMP END DO NOWAIT
194!
195!      jjb=jj_begin
196!      jje=jj_end
197!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
198c
199
200c    ............   integration  de  ucov, vcov,  h     ..............
201
202c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
203      DO 10 l = 1,llm
204     
205      ijb=ij_begin
206      ije=ij_end
207      if (pole_nord) ijb=ij_begin+iip1
208      if (pole_sud)  ije=ij_end-iip1
209     
210      DO 4 ij = ijb,ije
211      uscr( ij )   =  ucov( ij,l )
212      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
213   4  CONTINUE
214
215      ijb=ij_begin
216      ije=ij_end
217      if (pole_sud)  ije=ij_end-iip1
218     
219      DO 5 ij = ijb,ije
220      vscr( ij )   =  vcov( ij,l )
221      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
222   5  CONTINUE
223     
224      ijb=ij_begin
225      ije=ij_end
226     
227      DO 6 ij = ijb,ije
228      hscr( ij )    =  teta(ij,l)
229      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
230     $                + dt * dteta(ij,l) / masse(ij,l)
231   6  CONTINUE
232
233c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
234c
235c
236      IF (pole_nord) THEN
237       
238        DO  ij   = 1, iim
239          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
240        ENDDO
241          tpn      = SSUM(iim,tppn,1)/apoln
242
243        DO ij   = 1, iip1
244          teta(   ij   ,l)  = tpn
245        ENDDO
246     
247      ENDIF
248     
249      IF (pole_sud) THEN
250       
251        DO  ij   = 1, iim
252          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
253        ENDDO
254          tps      = SSUM(iim,tpps,1)/apols
255
256        DO ij   = 1, iip1
257          teta(ij+ip1jm,l)  = tps
258        ENDDO
259     
260      ENDIF
261c
262
263      IF(leapf)  THEN
264c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
265c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
266c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
267        ijb=ij_begin
268        ije=ij_end
269        ucovm1(ijb:ije,l)=uscr(ijb:ije)
270        tetam1(ijb:ije,l)=hscr(ijb:ije)
271        if (pole_sud) ije=ij_end-iip1
272        vcovm1(ijb:ije,l)=vscr(ijb:ije)
273     
274      END IF
275
276  10  CONTINUE
277c$OMP END DO NOWAIT
278
279c
280c   .......  integration de   q   ......
281c
282      ijb=ij_begin
283      ije=ij_end
284
285         if (planet_type.eq."earth") then
286! Earth-specific treatment of first 2 tracers (water)
287c$OMP BARRIER
288c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
289          DO l = 1, llm
290           DO ij = ijb, ije
291            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
292           ENDDO
293          ENDDO
294c$OMP END DO NOWAIT
295c$OMP BARRIER
296
297          CALL qminimum_p( q, nq, deltap )
298c
299c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
300c
301c$OMP BARRIER
302      IF (pole_nord) THEN
303     
304        DO iq = 1, nq
305       
306c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
307          DO l = 1, llm
308 
309             DO ij = 1, iim
310               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
311             ENDDO
312               qpn  =  SSUM(iim,qppn,1)/apoln
313     
314             DO ij = 1, iip1
315               q(   ij   ,l,iq)  = qpn
316             ENDDO   
317 
318          ENDDO
319c$OMP END DO NOWAIT
320
321        ENDDO
322     
323      ENDIF
324
325      IF (pole_sud) THEN
326     
327        DO iq = 1, nq
328
329c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
330          DO l = 1, llm
331 
332             DO ij = 1, iim
333               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
334             ENDDO
335               qps  =  SSUM(iim,qpps,1)/apols
336 
337             DO ij = 1, iip1
338               q(ij+ip1jm,l,iq)  = qps
339             ENDDO   
340 
341          ENDDO
342c$OMP END DO NOWAIT
343
344        ENDDO
345     
346      ENDIF
347     
348! Ehouarn: forget about finvmaold
349!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
350!
351!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
352!      DO l = 1, llm     
353!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
354!      ENDDO
355!c$OMP END DO NOWAIT
356
357      endif ! of if (planet_type.eq."earth")
358
359c
360c
361c     .....   FIN  de l'integration  de   q    .......
362
36315    continue
364
365c$OMP DO SCHEDULE(STATIC)
366      DO ij=ijb,ije 
367        ps0(ij)=ps(ij)
368      ENDDO
369c$OMP END DO NOWAIT
370
371c    .................................................................
372
373
374      IF( leapf )  THEN
375c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
376c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
377c$OMP DO SCHEDULE(STATIC)
378      DO ij=ijb,ije 
379        psm1(ij)=pscr(ij)
380      ENDDO
381c$OMP END DO NOWAIT
382
383c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
384          DO l = 1, llm
385            massem1(ijb:ije,l)=massescr(ijb:ije,l)
386          ENDDO
387c$OMP END DO NOWAIT       
388      END IF
389c$OMP BARRIER
390      RETURN
391      END
Note: See TracBrowser for help on using the repository browser.