source: trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F @ 1243

Last change on this file since 1243 was 1019, checked in by emillour, 11 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

  • General stuff:
  • makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

  • bld.cfg: adaptations to enable compiling using multiple threads
  • build_gcm: adaptations to enable compiling using multiple threads
  • makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change
  • bibio:
  • wxios.F90 : Added for possible future use of XIOS library
  • filtrez:
  • mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
  • filtreg_mod : limit use of FFT to parallel mode
  • mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz
  • dyn3d:
  • abort_gcm.F : add things for xios
  • advtrac.F90 : minor change in CFL outputs
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • inigeom.F : move two computations outside loop
  • dyn3dpar:
  • parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • leapfrog_p.F : add xios + correction for times in Newtonian case
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • inigeom.F : move two computations outside loop

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