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

Last change on this file since 2278 was 2270, checked in by crisi, 9 years ago

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

  • 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
92        if (ok_iso_verif) then
93           call check_isotopes(q,ijb,ije,'integrd 92')
94        endif !if (ok_iso_verif) then
95
96c$OMP BARRIER     
97      if (pole_nord) THEN
98c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
99        DO  l = 1,llm
100          DO  ij = 1,iip1
101           ucov(    ij    , l) = 0.
102           uscr(     ij      ) = 0.
103           ENDDO
104        ENDDO
105c$OMP END DO NOWAIT       
106      ENDIF
107
108      if (pole_sud) THEN
109c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
110        DO  l = 1,llm
111          DO  ij = 1,iip1
112           ucov( ij +ip1jm, l) = 0.
113           uscr( ij +ip1jm   ) = 0.
114          ENDDO
115        ENDDO
116c$OMP END DO NOWAIT     
117      ENDIF
118
119c    ............    integration  de       ps         ..............
120
121c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
122
123      ijb=ij_begin
124      ije=ij_end
125c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
126      DO  l = 1,llm
127        massescr(ijb:ije,l)=masse(ijb:ije,l)
128      ENDDO
129c$OMP END DO NOWAIT
130
131c$OMP DO SCHEDULE(STATIC)
132      DO 2 ij = ijb,ije
133       pscr (ij)    = ps0(ij)
134       ps (ij)      = psm1(ij) + dt * dp(ij)     
135
136   2  CONTINUE
137
138c$OMP END DO 
139c$OMP BARRIER
140c --> ici synchro OPENMP pour ps
141       
142      checksum=.TRUE.
143      stop_it=0
144
145c$OMP MASTER
146!c$OMP DO SCHEDULE(STATIC)
147      DO ij = ijb,ije
148         IF( ps(ij).LT.0. ) THEN
149           IF (checksum) stop_it=ij
150           checksum=.FALSE.
151         ENDIF
152       ENDDO
153!c$OMP END DO NOWAIT
154       
155!      CALL MPI_ALLREDUCE(checksum,checksum_all,1,
156!     &                   MPI_LOGICAL,MPI_LOR,COMM_LMDZ,ierr)
157      IF( .NOT. checksum ) THEN
158         write(lunout,*) "integrd: ps = ", ps(stop_it)
159         write(lunout,*) " at node ij =", stop_it
160         ! since ij=j+(i-1)*jjp1 , we have
161!         j=modulo(stop_it,jjp1)
162!         i=1+(stop_it-j)/jjp1
163!         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
164!     &                   " lat = ",rlatu(j)*180./pi, " deg"
165         call abort_gcm("integrd_loc", "negative surface pressure", 1)
166      ENDIF
167
168c$OMP END MASTER
169c$OMP BARRIER
170      IF (.NOT. Checksum_all) THEN
171        call WriteField_v('int_vcov',vcov)
172        call WriteField_u('int_ucov',ucov)
173        call WriteField_u('int_teta',teta)
174        call WriteField_u('int_ps0',ps0)
175        call WriteField_u('int_masse',masse)
176        call WriteField_u('int_phis',phis)
177        call WriteField_v('int_vcovm1',vcovm1)
178        call WriteField_u('int_ucovm1',ucovm1)
179        call WriteField_u('int_tetam1',tetam1)
180        call WriteField_u('int_psm1',psm1)
181        call WriteField_u('int_massem1',massem1)
182
183        call WriteField_v('int_dv',dv)
184        call WriteField_u('int_du',du)
185        call WriteField_u('int_dteta',dteta)
186        call WriteField_u('int_dp',dp)
187!        call WriteField_u('int_finvmaold',finvmaold)
188        do j=1,nq
189          call WriteField_u('int_q'//trim(int2str(j)),
190     .                q(:,:,j))
191          call WriteField_u('int_dq'//trim(int2str(j)),
192     .                dq(:,:,j))
193        enddo
194        call abort_gcm("integrd_loc", "", 1)
195      ENDIF
196   
197       
198c
199C$OMP MASTER
200      if (pole_nord) THEN
201     
202        DO  ij    = 1, iim
203         tppn(ij) = aire(   ij   ) * ps(  ij    )
204        ENDDO
205         tpn      = SSUM(iim,tppn,1)/apoln
206        DO ij   = 1, iip1
207         ps(   ij   )  = tpn
208        ENDDO
209     
210      ENDIF
211     
212      if (pole_sud) THEN
213     
214        DO  ij    = 1, iim
215         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
216        ENDDO
217         tps      = SSUM(iim,tpps,1)/apols
218        DO ij   = 1, iip1
219         ps(ij+ip1jm)  = tps
220        ENDDO
221     
222      ENDIF
223c$OMP END MASTER
224c$OMP BARRIER
225      !write(*,*) 'integrd 217' 
226c
227c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
228c
229
230      CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
231
232c$OMP BARRIER
233      CALL massdair_loc (     p  , masse         )
234
235! Ehouarn : we don't use/need finvmaold and finvmasse,
236!           so might as well not compute them
237!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
238!      ijb=ij_begin
239!      ije=ij_end
240!     
241!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
242!      DO  l = 1,llm
243!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
244!      ENDDO
245!c$OMP END DO NOWAIT
246
247!      jjb=jj_begin
248!      jje=jj_end
249!      CALL filtreg_p( finvmasse,jjb_u,jje_u,jjb,jje, jjp1, llm,
250!     &                -2, 2, .TRUE., 1  )
251c
252
253c    ............   integration  de  ucov, vcov,  h     ..............
254
255c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
256      DO 10 l = 1,llm
257     
258      ijb=ij_begin
259      ije=ij_end
260      if (pole_nord) ijb=ij_begin+iip1
261      if (pole_sud)  ije=ij_end-iip1
262     
263      DO 4 ij = ijb,ije
264      uscr( ij )   =  ucov( ij,l )
265      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
266   4  CONTINUE
267
268      ijb=ij_begin
269      ije=ij_end
270      if (pole_sud)  ije=ij_end-iip1
271     
272      DO 5 ij = ijb,ije
273      vscr( ij )   =  vcov( ij,l )
274      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
275   5  CONTINUE
276     
277      ijb=ij_begin
278      ije=ij_end
279     
280      DO 6 ij = ijb,ije
281      hscr( ij )    =  teta(ij,l)
282      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
283     $                + dt * dteta(ij,l) / masse(ij,l)
284   6  CONTINUE
285
286c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
287c
288c
289      IF (pole_nord) THEN
290       
291        DO  ij   = 1, iim
292          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
293        ENDDO
294          tpn      = SSUM(iim,tppn,1)/apoln
295
296        DO ij   = 1, iip1
297          teta(   ij   ,l)  = tpn
298        ENDDO
299     
300      ENDIF
301     
302      IF (pole_sud) THEN
303       
304        DO  ij   = 1, iim
305          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
306        ENDDO
307          tps      = SSUM(iim,tpps,1)/apols
308
309        DO ij   = 1, iip1
310          teta(ij+ip1jm,l)  = tps
311        ENDDO
312     
313      ENDIF
314c
315
316      IF(leapf)  THEN
317c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
318c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
319c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
320        ijb=ij_begin
321        ije=ij_end
322        ucovm1(ijb:ije,l)=uscr(ijb:ije)
323        tetam1(ijb:ije,l)=hscr(ijb:ije)
324        if (pole_sud) ije=ij_end-iip1
325        vcovm1(ijb:ije,l)=vscr(ijb:ije)
326     
327      END IF
328
329  10  CONTINUE
330c$OMP END DO NOWAIT
331
332c
333c   .......  integration de   q   ......
334c
335      ijb=ij_begin
336      ije=ij_end
337
338         if (planet_type.eq."earth") then
339! Earth-specific treatment of first 2 tracers (water)
340c$OMP BARRIER
341c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
342          DO l = 1, llm
343           DO ij = ijb, ije
344            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
345           ENDDO
346          ENDDO
347         
348c$OMP END DO NOWAIT
349c$OMP BARRIER
350
351        if (ok_iso_verif) then
352           call check_isotopes(q,ijb,ije,'integrd 342')
353        endif !if (ok_iso_verif) then
354
355          !write(*,*) 'integrd 341'
356          CALL qminimum_loc( q, nq, deltap )
357          !write(*,*) 'integrd 343'
358
359        if (ok_iso_verif) then
360           call check_isotopes(q,ijb,ije,'integrd 346')
361        endif !if (ok_iso_verif) then
362c
363c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
364c
365c$OMP BARRIER
366      IF (pole_nord) THEN
367     
368        DO iq = 1, nq
369       
370c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
371          DO l = 1, llm
372 
373             DO ij = 1, iim
374               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
375             ENDDO
376               qpn  =  SSUM(iim,qppn,1)/apoln
377     
378             DO ij = 1, iip1
379               q(   ij   ,l,iq)  = qpn
380             ENDDO   
381 
382          ENDDO
383c$OMP END DO NOWAIT
384
385        ENDDO
386     
387      ENDIF
388
389      IF (pole_sud) THEN
390     
391        DO iq = 1, nq
392
393c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
394          DO l = 1, llm
395 
396             DO ij = 1, iim
397               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
398             ENDDO
399               qps  =  SSUM(iim,qpps,1)/apols
400 
401             DO ij = 1, iip1
402               q(ij+ip1jm,l,iq)  = qps
403             ENDDO   
404 
405          ENDDO
406c$OMP END DO NOWAIT
407
408        ENDDO
409     
410      ENDIF
411
412        if (ok_iso_verif) then
413           call check_isotopes(q,ijb,ije,'integrd 409')
414        endif !if (ok_iso_verif) then
415     
416! Ehouarn: forget about finvmaold
417!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
418
419!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
420!      DO l = 1, llm     
421!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
422!      ENDDO
423!c$OMP END DO NOWAIT
424
425      endif ! of if (planet_type.eq."earth")
426
427c
428c
429c     .....   FIN  de l'integration  de   q    .......
430
43115    continue
432          !write(*,*) 'integrd 410'
433
434c$OMP DO SCHEDULE(STATIC)
435      DO ij=ijb,ije 
436        ps0(ij)=ps(ij)
437      ENDDO
438c$OMP END DO NOWAIT
439
440c    .................................................................
441
442
443      IF( leapf )  THEN
444c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
445c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
446c$OMP DO SCHEDULE(STATIC)
447      DO ij=ijb,ije 
448        psm1(ij)=pscr(ij)
449      ENDDO
450c$OMP END DO NOWAIT
451
452c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
453          DO l = 1, llm
454            massem1(ijb:ije,l)=massescr(ijb:ije,l)
455          ENDDO
456c$OMP END DO NOWAIT       
457      END IF
458c$OMP BARRIER
459      RETURN
460      END
Note: See TracBrowser for help on using the repository browser.