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

Last change on this file since 2461 was 2461, checked in by fhourdin, 10 years ago

Impression de la latitude et de la longitude en cas de plantage dans
dyn3dmeme/intergrd_loc

  • 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
Line 
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,
6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
7      USE parallel_lmdz
8      USE control_mod
9      USE mod_filtreg_p
10      USE write_field_loc
11      USE write_field
12      USE integrd_mod
13      USE infotrac, ONLY: ok_iso_verif ! ajout CRisi
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"
40#include "iniprint.h"
41!      include 'mpif.h'
42
43c   Arguments:
44c   ----------
45
46      INTEGER,intent(in) :: nq ! number of tracers to handle in this routine
47
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
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)
74      REAL massescr( ijb_u:ije_u,llm )
75!      REAL finvmasse(ijb_u:ije_u,llm)
76      REAL tpn,tps,tppn(iim),tpps(iim)
77      REAL qpn,qps,qppn(iim),qpps(iim)
78
79      INTEGER  l,ij,iq,i,j
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
87      INTEGER :: ierr
88
89      !write(*,*) 'integrd 88: entree, nq=',nq
90c-----------------------------------------------------------------------
91
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)
130       ps (ij)      = psm1(ij) + dt * dp(ij)     
131
132   2  CONTINUE
133
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
154         write(lunout,*) "integrd: ps = ", ps(stop_it)
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"
161         call abort_gcm("integrd_loc", "negative surface pressure", 1)
162      ENDIF
163
164c$OMP END MASTER
165c$OMP BARRIER
166        !write(*,*) 'integrd 170'
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)
184!        call WriteField_u('int_finvmaold',finvmaold)
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
191        call abort_gcm("integrd_loc", "", 1)
192      ENDIF
193   
194       
195c
196        !write(*,*) 'integrd 200'
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
223      !write(*,*) 'integrd 217' 
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 )
229
230c$OMP BARRIER
231      CALL massdair_loc (     p  , masse         )
232
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
244
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  )
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
287        !write(*,*) 'integrd 291'
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
346         
347c$OMP END DO NOWAIT
348c$OMP BARRIER
349
350        if (ok_iso_verif) then
351           call check_isotopes(q,ijb,ije,'integrd 342')
352        endif !if (ok_iso_verif) then
353
354        !write(*,*) 'integrd 341'
355        CALL qminimum_loc( q, nq, deltap )
356        !write(*,*) 'integrd 343'
357
358        if (ok_iso_verif) then
359           call check_isotopes(q,ijb,ije,'integrd 346')
360        endif !if (ok_iso_verif) then
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
410
411        if (ok_iso_verif) then
412           call check_isotopes(q,ijb,ije,'integrd 409')
413        endif !if (ok_iso_verif) then
414     
415! Ehouarn: forget about finvmaold
416!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
417
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
426c
427c
428c     .....   FIN  de l'integration  de   q    .......
429
43015    continue
431          !write(*,*) 'integrd 410'
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.