source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dmem/integrd_loc.F @ 5452

Last change on this file since 5452 was 1705, checked in by Ehouarn Millour, 12 years ago

Added arch files for ADA (IDRIS IBMx3750) and made the following code modifications:
phylmd/printflag.F : removed "print" of unset variable (radpas0)
dyn3dmem/integrd_loc.F : removed unecessary "include mpif.h"
dyn3dmem/leapfrog_loc.F : removed unecessary "include mpif.h" and allocate saved variables at first call
dyn3dmem/mod_filtreg_p.F : added matmul() alternatives to call to BLAS routine SGEMM (which was incorectly set as DGEMM; which would fail if running with -r4)
filtrez/filtreg.F: changed calls to DGEMM into calls to SGEMM, so that code works with either -r4 or -r8 (the later being used in conjunction with "BLAS SGEMV=DGEMV SGEMM=DGEMM" preprocessing statements)
EM
File size: 11.0 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
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"
39#include "iniprint.h"
40!      include 'mpif.h'
41
42c   Arguments:
43c   ----------
44
45      INTEGER,intent(in) :: nq ! number of tracers to handle in this routine
46
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
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)
73      REAL massescr( ijb_u:ije_u,llm )
74!      REAL finvmasse(ijb_u:ije_u,llm)
75      REAL tpn,tps,tppn(iim),tpps(iim)
76      REAL qpn,qps,qppn(iim),qpps(iim)
77
78      INTEGER  l,ij,iq,i,j
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
86      INTEGER :: ierr
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
149          write(lunout,*) "integrd: negative surface pressure ",
150     &                                                ps(stop_it)
151         write(lunout,*) " at node ij =", stop_it
152         ! since ij=j+(i-1)*jjp1 , we have
153!         j=modulo(stop_it,jjp1)
154!         i=1+(stop_it-j)/jjp1
155!         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
156!     &                   " lat = ",rlatu(j)*180./pi, " deg"
157      ENDIF
158
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)
178!        call WriteField_u('int_finvmaold',finvmaold)
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
185      STOP
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
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
235
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  )
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     
390! Ehouarn: forget about finvmaold
391!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
392
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
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.