source: LMDZ5/trunk/libf/dyn3dmem/integrd_loc.F @ 2261

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