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

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod.
EM

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