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

Last change on this file since 1508 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

  • in dynlonlat_phylonlat:
  • add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)
  • in misc:
  • follow up updates on wxios.F (add missing_val module variable)
  • in dyn3d_common:
  • pression.F => pression.F90
  • misc_mod.F90: moved from misc to dyn3d_common
  • added new iso_verif_dyn.F
  • covcont.F => covcont.F90
  • infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
  • dynetat0.F => dynetat0.F90 with some code factorization
  • dynredem.F => dynredem.F90 with some code factorization
  • added dynredem_mod.F90: routines used by dynredem
  • iniacademic.F90 : added isotopes-related initialization for Earth case
  • in dyn3d:
  • added check_isotopes.F
  • modified (isotopes) advtrac.F90, caladvtrac.F
  • guide_mod.F90: ported updates
  • leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
  • qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
  • vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
  • vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)
  • in dyn3dpar:
  • leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
  • dynredem_p.F => dynredem_p.F90 and some code factorization
  • and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

EM
File size: 9.9 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, ONLY: ij_begin, ij_end, pole_nord, pole_sud
8      USE control_mod, only : planet_type
9      USE comvert_mod, ONLY: ap,bp
10      USE comconst_mod, ONLY: pi
11      USE logic_mod, ONLY: leapf
12      USE temps_mod, ONLY: dt
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 "comgeom.h"
34#include "iniprint.h"
35
36c   Arguments:
37c   ----------
38
39      integer,intent(in) :: nq ! number of tracers to handle in this routine
40      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
41      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
42      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
43      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
44      real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
45      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
46      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
47      ! values at previous time step
48      real,intent(inout) :: vcovm1(ip1jm,llm)
49      real,intent(inout) :: ucovm1(ip1jmp1,llm)
50      real,intent(inout) :: tetam1(ip1jmp1,llm)
51      real,intent(inout) :: psm1(ip1jmp1)
52      real,intent(inout) :: massem1(ip1jmp1,llm)
53      ! the tendencies to add
54      real,intent(in) :: dv(ip1jm,llm)
55      real,intent(in) :: du(ip1jmp1,llm)
56      real,intent(in) :: dteta(ip1jmp1,llm)
57      real,intent(in) :: dp(ip1jmp1)
58      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
59!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
60
61c   Local:
62c   ------
63
64      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
65      REAL massescr( ip1jmp1,llm )
66!      REAL finvmasse(ip1jmp1,llm)
67      REAL,SAVE :: p(ip1jmp1,llmp1)
68      REAL tpn,tps,tppn(iim),tpps(iim)
69      REAL qpn,qps,qppn(iim),qpps(iim)
70      REAL,SAVE :: deltap( ip1jmp1,llm )
71
72      INTEGER  l,ij,iq,i,j
73
74      REAL SSUM
75      EXTERNAL SSUM
76      INTEGER ijb,ije,jjb,jje
77      REAL,SAVE :: ps(ip1jmp1)
78      LOGICAL :: checksum
79      INTEGER :: stop_it
80c-----------------------------------------------------------------------
81c$OMP BARRIER     
82      if (pole_nord) THEN
83c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
84        DO  l = 1,llm
85          DO  ij = 1,iip1
86           ucov(    ij    , l) = 0.
87           uscr(     ij      ) = 0.
88           ENDDO
89        ENDDO
90c$OMP END DO NOWAIT       
91      ENDIF
92
93      if (pole_sud) THEN
94c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
95        DO  l = 1,llm
96          DO  ij = 1,iip1
97           ucov( ij +ip1jm, l) = 0.
98           uscr( ij +ip1jm   ) = 0.
99          ENDDO
100        ENDDO
101c$OMP END DO NOWAIT     
102      ENDIF
103
104c    ............    integration  de       ps         ..............
105
106c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
107
108      ijb=ij_begin
109      ije=ij_end
110c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
111      DO  l = 1,llm
112        massescr(ijb:ije,l)=masse(ijb:ije,l)
113      ENDDO
114c$OMP END DO NOWAIT
115
116c$OMP DO SCHEDULE(STATIC)
117      DO 2 ij = ijb,ije
118       pscr (ij)    = ps0(ij)
119       ps (ij)      = psm1(ij) + dt * dp(ij)
120   2  CONTINUE
121c$OMP END DO 
122c$OMP BARRIER
123c --> ici synchro OPENMP pour ps
124       
125      checksum=.TRUE.
126      stop_it=0
127
128c$OMP DO SCHEDULE(STATIC)
129      DO ij = ijb,ije
130         IF( ps(ij).LT.0. ) THEN
131           IF (checksum) stop_it=ij
132           checksum=.FALSE.
133         ENDIF
134       ENDDO
135c$OMP END DO NOWAIT
136       
137        IF( .NOT. checksum ) THEN
138         write(lunout,*) "integrd: negative surface pressure ",
139     &                                                ps(stop_it)
140         write(lunout,*) " at node ij =", stop_it
141         ! since ij=j+(i-1)*jjp1 , we have
142         j=modulo(stop_it,jjp1)
143         i=1+(stop_it-j)/jjp1
144         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
145     &                   " lat = ",rlatu(j)*180./pi, " deg"
146         write(lunout,*) " psm1(ij)=",psm1(stop_it)," dt=",dt,
147     &                   " dp(ij)=",dp(stop_it)
148         call abort_gcm("integrd_p", "negative surface pressure", 1)
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.