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

Last change on this file since 2397 was 2286, checked in by Ehouarn Millour, 10 years ago
  • fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
  • comment out many invasive debug writes.

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