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

Last change on this file was 2126, checked in by emillour, 6 years ago

Common dynamics:
Some work to enforce total tracer mass conservation in the dynamics.
Still to be further studied and validated.
For now these changes are triggered by setting a "force_conserv_tracer"
flag to ".true." in run.def (default is ".false." to not change anything
with respect to previous versions).
When force_conserv_tracer=.true. then:

  1. Rescale tracer mass in caladvtrac after tracer advection computations
  2. Recompute q ratios once atmospheric mass has been updated in integrd

These steps technically ensure total tracer mass conservation but it
might be the tracer advection scheme and/or time-stepping updating
sequence of fields that should be rethought or fixed.
EM
File size: 10.4 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     &                         omp_chunk
9      USE control_mod, only : planet_type,force_conserv_tracer
10      USE comvert_mod, ONLY: ap,bp
11      USE comconst_mod, ONLY: pi
12      USE logic_mod, ONLY: leapf
13      USE temps_mod, ONLY: dt
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 "comgeom.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 :: massratio(ip1jmp1,llm)
68!      REAL finvmasse(ip1jmp1,llm)
69      REAL,SAVE :: p(ip1jmp1,llmp1)
70      REAL tpn,tps,tppn(iim),tpps(iim)
71      REAL qpn,qps,qppn(iim),qpps(iim)
72      REAL,SAVE :: deltap( ip1jmp1,llm )
73
74      INTEGER  l,ij,iq,i,j
75
76      REAL SSUM
77      EXTERNAL SSUM
78      INTEGER ijb,ije,jjb,jje
79      REAL,SAVE :: ps(ip1jmp1)=0
80      LOGICAL :: checksum
81      INTEGER :: stop_it
82c-----------------------------------------------------------------------
83c$OMP BARRIER     
84      if (pole_nord) THEN
85c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
86        DO  l = 1,llm
87          DO  ij = 1,iip1
88           ucov(    ij    , l) = 0.
89           uscr(     ij      ) = 0.
90           ENDDO
91        ENDDO
92c$OMP END DO NOWAIT       
93      ENDIF
94
95      if (pole_sud) THEN
96c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
97        DO  l = 1,llm
98          DO  ij = 1,iip1
99           ucov( ij +ip1jm, l) = 0.
100           uscr( ij +ip1jm   ) = 0.
101          ENDDO
102        ENDDO
103c$OMP END DO NOWAIT     
104      ENDIF
105
106c    ............    integration  de       ps         ..............
107
108c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
109
110      ijb=ij_begin
111      ije=ij_end
112c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
113      DO  l = 1,llm
114        massescr(ijb:ije,l)=masse(ijb:ije,l)
115      ENDDO
116c$OMP END DO NOWAIT
117
118c$OMP DO SCHEDULE(STATIC)
119      DO 2 ij = ijb,ije
120       pscr (ij)    = ps0(ij)
121       ps (ij)      = psm1(ij) + dt * dp(ij)
122   2  CONTINUE
123c$OMP END DO 
124c$OMP BARRIER
125c --> ici synchro OPENMP pour ps
126       
127      checksum=.TRUE.
128      stop_it=0
129
130c$OMP DO SCHEDULE(STATIC)
131      DO ij = ijb,ije
132         IF( ps(ij).LT.0. ) THEN
133           IF (checksum) stop_it=ij
134           checksum=.FALSE.
135         ENDIF
136       ENDDO
137c$OMP END DO NOWAIT
138       
139        IF( .NOT. checksum ) THEN
140         write(lunout,*) "integrd: negative surface pressure ",
141     &                                                ps(stop_it)
142         write(lunout,*) " at node ij =", stop_it
143         ! since ij=j+(i-1)*jjp1 , we have
144         j=modulo(stop_it,jjp1)
145         i=1+(stop_it-j)/jjp1
146         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
147     &                   " lat = ",rlatu(j)*180./pi, " deg"
148         write(lunout,*) " psm1(ij)=",psm1(stop_it)," dt=",dt,
149     &                   " dp(ij)=",dp(stop_it)
150         call abort_gcm("integrd_p", "negative surface pressure", 1)
151        ENDIF
152
153c
154C$OMP MASTER
155      if (pole_nord) THEN
156     
157        DO  ij    = 1, iim
158         tppn(ij) = aire(   ij   ) * ps(  ij    )
159        ENDDO
160         tpn      = SSUM(iim,tppn,1)/apoln
161        DO ij   = 1, iip1
162         ps(   ij   )  = tpn
163        ENDDO
164     
165      ENDIF
166     
167      if (pole_sud) THEN
168     
169        DO  ij    = 1, iim
170         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
171        ENDDO
172         tps      = SSUM(iim,tpps,1)/apols
173        DO ij   = 1, iip1
174         ps(ij+ip1jm)  = tps
175        ENDDO
176     
177      ENDIF
178c$OMP END MASTER
179c$OMP BARRIER
180c
181c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
182c
183
184      CALL pression_p ( ip1jmp1, ap, bp, ps, p )
185c$OMP BARRIER
186      CALL massdair_p (     p  , masse         )
187
188! Ehouarn : we don't use/need finvmaold and finvmasse,
189!           so might as well not compute them
190!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
191!      ijb=ij_begin
192!      ije=ij_end
193!     
194!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
195!      DO  l = 1,llm
196!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
197!      ENDDO
198!c$OMP END DO NOWAIT
199!
200!      jjb=jj_begin
201!      jje=jj_end
202!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
203c
204
205c    ............   integration  de  ucov, vcov,  h     ..............
206
207c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
208      DO 10 l = 1,llm
209     
210      ijb=ij_begin
211      ije=ij_end
212      if (pole_nord) ijb=ij_begin+iip1
213      if (pole_sud)  ije=ij_end-iip1
214     
215      DO 4 ij = ijb,ije
216      uscr( ij )   =  ucov( ij,l )
217      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
218   4  CONTINUE
219
220      ijb=ij_begin
221      ije=ij_end
222      if (pole_sud)  ije=ij_end-iip1
223     
224      DO 5 ij = ijb,ije
225      vscr( ij )   =  vcov( ij,l )
226      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
227   5  CONTINUE
228     
229      ijb=ij_begin
230      ije=ij_end
231     
232      DO 6 ij = ijb,ije
233      hscr( ij )    =  teta(ij,l)
234      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
235     $                + dt * dteta(ij,l) / masse(ij,l)
236   6  CONTINUE
237
238c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
239c
240c
241      IF (pole_nord) THEN
242       
243        DO  ij   = 1, iim
244          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
245        ENDDO
246          tpn      = SSUM(iim,tppn,1)/apoln
247
248        DO ij   = 1, iip1
249          teta(   ij   ,l)  = tpn
250        ENDDO
251     
252      ENDIF
253     
254      IF (pole_sud) THEN
255       
256        DO  ij   = 1, iim
257          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
258        ENDDO
259          tps      = SSUM(iim,tpps,1)/apols
260
261        DO ij   = 1, iip1
262          teta(ij+ip1jm,l)  = tps
263        ENDDO
264     
265      ENDIF
266c
267
268      IF(leapf)  THEN
269c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
270c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
271c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
272        ijb=ij_begin
273        ije=ij_end
274        ucovm1(ijb:ije,l)=uscr(ijb:ije)
275        tetam1(ijb:ije,l)=hscr(ijb:ije)
276        if (pole_sud) ije=ij_end-iip1
277        vcovm1(ijb:ije,l)=vscr(ijb:ije)
278     
279      END IF
280
281  10  CONTINUE
282c$OMP END DO NOWAIT
283
284c
285c   .......  integration de   q   ......
286c
287      ijb=ij_begin
288      ije=ij_end
289
290         if (planet_type.eq."earth") then
291! Earth-specific treatment of first 2 tracers (water)
292c$OMP BARRIER
293c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
294          DO l = 1, llm
295           DO ij = ijb, ije
296            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
297           ENDDO
298          ENDDO
299c$OMP END DO NOWAIT
300c$OMP BARRIER
301
302          CALL qminimum_p( q, nq, deltap )
303c
304c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
305c
306c$OMP BARRIER
307      IF (pole_nord) THEN
308     
309        DO iq = 1, nq
310       
311c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312          DO l = 1, llm
313 
314             DO ij = 1, iim
315               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
316             ENDDO
317               qpn  =  SSUM(iim,qppn,1)/apoln
318     
319             DO ij = 1, iip1
320               q(   ij   ,l,iq)  = qpn
321             ENDDO   
322 
323          ENDDO
324c$OMP END DO NOWAIT
325
326        ENDDO
327     
328      ENDIF
329
330      IF (pole_sud) THEN
331     
332        DO iq = 1, nq
333
334c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
335          DO l = 1, llm
336 
337             DO ij = 1, iim
338               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
339             ENDDO
340               qps  =  SSUM(iim,qpps,1)/apols
341 
342             DO ij = 1, iip1
343               q(ij+ip1jm,l,iq)  = qps
344             ENDDO   
345 
346          ENDDO
347c$OMP END DO NOWAIT
348
349        ENDDO
350     
351      ENDIF
352     
353! Ehouarn: forget about finvmaold
354!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
355!
356!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
357!      DO l = 1, llm     
358!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
359!      ENDDO
360!c$OMP END DO NOWAIT
361
362      endif ! of if (planet_type.eq."earth")
363
364
365      if (force_conserv_tracer) then
366        ! Ehouarn: try to keep total amont of tracers fixed
367        ! by acounting for mass change in each cell
368        massratio(ijb:ije,1:llm)=massescr(ijb:ije,1:llm)
369     &                             /masse(ijb:ije,1:llm)
370        do iq=1,nq
371        q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)
372     &                        *massratio(ijb:ije,1:llm)
373        enddo
374      endif ! of if (force_conserv_tracer)
375c
376c
377c     .....   FIN  de l'integration  de   q    .......
378
37915    continue
380
381c$OMP DO SCHEDULE(STATIC)
382      DO ij=ijb,ije 
383        ps0(ij)=ps(ij)
384      ENDDO
385c$OMP END DO NOWAIT
386
387c    .................................................................
388
389
390      IF( leapf )  THEN
391c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
392c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
393c$OMP DO SCHEDULE(STATIC)
394      DO ij=ijb,ije 
395        psm1(ij)=pscr(ij)
396      ENDDO
397c$OMP END DO NOWAIT
398
399c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
400          DO l = 1, llm
401            massem1(ijb:ije,l)=massescr(ijb:ije,l)
402          ENDDO
403c$OMP END DO NOWAIT       
404      END IF
405c$OMP BARRIER
406      RETURN
407      END
Note: See TracBrowser for help on using the repository browser.