source: LMDZ6/branches/blowing_snow/libf/dyn3dmem/integrd_loc.F @ 5453

Last change on this file since 5453 was 4143, checked in by dcugnet, 3 years ago
  • Some variables are renamed or replaced by direct equivalents:
    • iso_indnum -> tracers(:)%iso_iName
    • niso_possibles -> niso
    • iqiso -> iqIsoPha ; index_trac -> itZonIso
    • ok_iso_verif -> isoCheck
    • ntraceurs_zone -> nzone ; ntraciso -> ntiso
    • qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
  • Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
  • Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
  • The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
  • 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.5 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 comconst_mod, ONLY: pi
14      USE logic_mod, ONLY: leapf
15      USE comvert_mod, ONLY: ap, bp
16      USE temps_mod, ONLY: dt
17     
18      IMPLICIT NONE
19
20
21c=======================================================================
22c
23c   Auteur:  P. Le Van
24c   -------
25c
26c   objet:
27c   ------
28c
29c   Incrementation des tendances dynamiques
30c
31c=======================================================================
32c-----------------------------------------------------------------------
33c   Declarations:
34c   -------------
35
36      include "dimensions.h"
37      include "paramet.h"
38      include "comgeom.h"
39      include "iniprint.h"
40
41c   Arguments:
42c   ----------
43
44      INTEGER,intent(in) :: nq ! number of tracers to handle in this routine
45
46      REAL,INTENT(INOUT) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
47      REAL,INTENT(INOUT) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
48      REAL,INTENT(INOUT) :: teta(ijb_u:ije_u,llm) ! potential temperature
49      REAL,INTENT(INOUT) :: q(ijb_u:ije_u,llm,nq) ! advected tracers
50      REAL,INTENT(INOUT) :: ps0(ijb_u:ije_u) ! surface pressure
51      REAL,INTENT(INOUT) :: masse(ijb_u:ije_u,llm) ! atmospheric mass
52      REAL,INTENT(INOUT) :: phis(ijb_u:ije_u) ! ground geopotential !!! unused
53      ! values at previous time step
54      REAL,INTENT(INOUT) :: vcovm1(ijb_v:ije_v,llm)
55      REAL,INTENT(INOUT) :: ucovm1(ijb_u:ije_u,llm)
56      REAL,INTENT(INOUT) :: tetam1(ijb_u:ije_u,llm)
57      REAL,INTENT(INOUT) :: psm1(ijb_u:ije_u)
58      REAL,INTENT(INOUT) :: massem1(ijb_u:ije_u,llm)
59      ! the tendencies to add
60      REAL,INTENT(INOUT) :: dv(ijb_v:ije_v,llm)
61      REAL,INTENT(INOUT) :: du(ijb_u:ije_u,llm)
62      REAL,INTENT(INOUT) :: dteta(ijb_u:ije_u,llm)
63      REAL,INTENT(INOUT) :: dp(ijb_u:ije_u)
64      REAL,INTENT(INOUT) :: dq(ijb_u:ije_u,llm,nq) !!! unused
65!      REAL,INTENT(INOUT) ::finvmaold(ijb_u:ije_u,llm) !!! unused
66
67c   Local:
68c   ------
69
70      REAL vscr( ijb_v:ije_v ),uscr( ijb_u:ije_u )
71      REAL hscr( ijb_u:ije_u ),pscr(ijb_u:ije_u)
72      REAL massescr( ijb_u:ije_u,llm )
73!      REAL finvmasse(ijb_u:ije_u,llm)
74      REAL tpn,tps,tppn(iim),tpps(iim)
75      REAL qpn,qps,qppn(iim),qpps(iim)
76
77      INTEGER  l,ij,iq,i,j
78
79      REAL SSUM
80      EXTERNAL SSUM
81      INTEGER ijb,ije,jjb,jje
82      LOGICAL :: checksum
83      LOGICAL,SAVE :: checksum_all=.TRUE.
84      INTEGER :: stop_it
85      INTEGER :: ierr
86
87      !write(*,*) 'integrd 88: entree, nq=',nq
88c-----------------------------------------------------------------------
89
90c$OMP BARRIER     
91      if (pole_nord) THEN
92c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
93        DO  l = 1,llm
94          DO  ij = 1,iip1
95           ucov(    ij    , l) = 0.
96           uscr(     ij      ) = 0.
97           ENDDO
98        ENDDO
99c$OMP END DO NOWAIT       
100      ENDIF
101
102      if (pole_sud) THEN
103c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
104        DO  l = 1,llm
105          DO  ij = 1,iip1
106           ucov( ij +ip1jm, l) = 0.
107           uscr( ij +ip1jm   ) = 0.
108          ENDDO
109        ENDDO
110c$OMP END DO NOWAIT     
111      ENDIF
112
113c    ............    integration  de       ps         ..............
114
115c      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
116
117      ijb=ij_begin
118      ije=ij_end
119c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
120      DO  l = 1,llm
121        massescr(ijb:ije,l)=masse(ijb:ije,l)
122      ENDDO
123c$OMP END DO NOWAIT
124
125c$OMP DO SCHEDULE(STATIC)
126      DO 2 ij = ijb,ije
127       pscr (ij)    = ps0(ij)
128       ps (ij)      = psm1(ij) + dt * dp(ij)     
129
130   2  CONTINUE
131
132c$OMP END DO 
133c$OMP BARRIER
134c --> ici synchro OPENMP pour ps
135       
136      checksum=.TRUE.
137      stop_it=0
138
139c$OMP MASTER
140!c$OMP DO SCHEDULE(STATIC)
141      DO ij = ijb,ije
142         IF( ps(ij).LT.0. ) THEN
143           IF (checksum) stop_it=ij
144           checksum=.FALSE.
145         ENDIF
146       ENDDO
147!c$OMP END DO NOWAIT
148       
149!      CALL MPI_ALLREDUCE(checksum,checksum_all,1,
150!     &                   MPI_LOGICAL,MPI_LOR,COMM_LMDZ,ierr)
151      IF( .NOT. checksum ) THEN
152         write(lunout,*) "integrd: ps = ", ps(stop_it)
153         write(lunout,*) " at node ij =", stop_it
154         ! since ij=j+(i-1)*jjp1 , we have
155          j=modulo(stop_it,jjp1)
156          i=1+(stop_it-j)/jjp1
157          write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
158     &                    " lat = ",rlatu(j)*180./pi, " deg"
159         call abort_gcm("integrd_loc", "negative surface pressure", 1)
160      ENDIF
161
162c$OMP END MASTER
163c$OMP BARRIER
164        !write(*,*) 'integrd 170'
165      IF (.NOT. Checksum_all) THEN
166        call WriteField_v('int_vcov',vcov)
167        call WriteField_u('int_ucov',ucov)
168        call WriteField_u('int_teta',teta)
169        call WriteField_u('int_ps0',ps0)
170        call WriteField_u('int_masse',masse)
171        call WriteField_u('int_phis',phis)
172        call WriteField_v('int_vcovm1',vcovm1)
173        call WriteField_u('int_ucovm1',ucovm1)
174        call WriteField_u('int_tetam1',tetam1)
175        call WriteField_u('int_psm1',psm1)
176        call WriteField_u('int_massem1',massem1)
177
178        call WriteField_v('int_dv',dv)
179        call WriteField_u('int_du',du)
180        call WriteField_u('int_dteta',dteta)
181        call WriteField_u('int_dp',dp)
182!        call WriteField_u('int_finvmaold',finvmaold)
183        do j=1,nq
184          call WriteField_u('int_q'//trim(int2str(j)),
185     .                q(:,:,j))
186          call WriteField_u('int_dq'//trim(int2str(j)),
187     .                dq(:,:,j))
188        enddo
189        call abort_gcm("integrd_loc", "", 1)
190      ENDIF
191   
192       
193c
194        !write(*,*) 'integrd 200'
195C$OMP MASTER
196      if (pole_nord) THEN
197     
198        DO  ij    = 1, iim
199         tppn(ij) = aire(   ij   ) * ps(  ij    )
200        ENDDO
201         tpn      = SSUM(iim,tppn,1)/apoln
202        DO ij   = 1, iip1
203         ps(   ij   )  = tpn
204        ENDDO
205     
206      ENDIF
207     
208      if (pole_sud) THEN
209     
210        DO  ij    = 1, iim
211         tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
212        ENDDO
213         tps      = SSUM(iim,tpps,1)/apols
214        DO ij   = 1, iip1
215         ps(ij+ip1jm)  = tps
216        ENDDO
217     
218      ENDIF
219c$OMP END MASTER
220c$OMP BARRIER
221      !write(*,*) 'integrd 217' 
222c
223c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
224c
225
226      CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
227
228c$OMP BARRIER
229      CALL massdair_loc (     p  , masse         )
230
231! Ehouarn : we don't use/need finvmaold and finvmasse,
232!           so might as well not compute them
233!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
234!      ijb=ij_begin
235!      ije=ij_end
236!     
237!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
238!      DO  l = 1,llm
239!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
240!      ENDDO
241!c$OMP END DO NOWAIT
242
243!      jjb=jj_begin
244!      jje=jj_end
245!      CALL filtreg_p( finvmasse,jjb_u,jje_u,jjb,jje, jjp1, llm,
246!     &                -2, 2, .TRUE., 1  )
247c
248
249c    ............   integration  de  ucov, vcov,  h     ..............
250
251c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
252      DO 10 l = 1,llm
253     
254      ijb=ij_begin
255      ije=ij_end
256      if (pole_nord) ijb=ij_begin+iip1
257      if (pole_sud)  ije=ij_end-iip1
258     
259      DO 4 ij = ijb,ije
260      uscr( ij )   =  ucov( ij,l )
261      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
262   4  CONTINUE
263
264      ijb=ij_begin
265      ije=ij_end
266      if (pole_sud)  ije=ij_end-iip1
267     
268      DO 5 ij = ijb,ije
269      vscr( ij )   =  vcov( ij,l )
270      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
271   5  CONTINUE
272     
273      ijb=ij_begin
274      ije=ij_end
275     
276      DO 6 ij = ijb,ije
277      hscr( ij )    =  teta(ij,l)
278      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
279     $                + dt * dteta(ij,l) / masse(ij,l)
280   6  CONTINUE
281
282c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
283c
284c
285        !write(*,*) 'integrd 291'
286      IF (pole_nord) THEN
287       
288        DO  ij   = 1, iim
289          tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
290        ENDDO
291          tpn      = SSUM(iim,tppn,1)/apoln
292
293        DO ij   = 1, iip1
294          teta(   ij   ,l)  = tpn
295        ENDDO
296     
297      ENDIF
298     
299      IF (pole_sud) THEN
300       
301        DO  ij   = 1, iim
302          tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
303        ENDDO
304          tps      = SSUM(iim,tpps,1)/apols
305
306        DO ij   = 1, iip1
307          teta(ij+ip1jm,l)  = tps
308        ENDDO
309     
310      ENDIF
311c
312
313      IF(leapf)  THEN
314c         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
315c         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
316c         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
317        ijb=ij_begin
318        ije=ij_end
319        ucovm1(ijb:ije,l)=uscr(ijb:ije)
320        tetam1(ijb:ije,l)=hscr(ijb:ije)
321        if (pole_sud) ije=ij_end-iip1
322        vcovm1(ijb:ije,l)=vscr(ijb:ije)
323     
324      END IF
325
326  10  CONTINUE
327c$OMP END DO NOWAIT
328
329c
330c   .......  integration de   q   ......
331c
332      ijb=ij_begin
333      ije=ij_end
334
335         if (planet_type.eq."earth") then
336! Earth-specific treatment of first 2 tracers (water)
337c$OMP BARRIER
338c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
339          DO l = 1, llm
340           DO ij = ijb, ije
341            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
342           ENDDO
343          ENDDO
344         
345c$OMP END DO NOWAIT
346c$OMP BARRIER
347
348        call check_isotopes(q,ijb,ije,'integrd 342')
349
350        !write(*,*) 'integrd 341'
351        CALL qminimum_loc( q, nq, deltap )
352        !write(*,*) 'integrd 343'
353
354        call check_isotopes(q,ijb,ije,'integrd 346')
355c
356c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
357c
358c$OMP BARRIER
359      IF (pole_nord) THEN
360     
361        DO iq = 1, nq
362       
363c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
364          DO l = 1, llm
365 
366             DO ij = 1, iim
367               qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
368             ENDDO
369               qpn  =  SSUM(iim,qppn,1)/apoln
370     
371             DO ij = 1, iip1
372               q(   ij   ,l,iq)  = qpn
373             ENDDO   
374 
375          ENDDO
376c$OMP END DO NOWAIT
377
378        ENDDO
379     
380      ENDIF
381
382      IF (pole_sud) THEN
383     
384        DO iq = 1, nq
385
386c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
387          DO l = 1, llm
388 
389             DO ij = 1, iim
390               qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
391             ENDDO
392               qps  =  SSUM(iim,qpps,1)/apols
393 
394             DO ij = 1, iip1
395               q(ij+ip1jm,l,iq)  = qps
396             ENDDO   
397 
398          ENDDO
399c$OMP END DO NOWAIT
400
401        ENDDO
402     
403      ENDIF
404
405      call check_isotopes(q,ijb,ije,'integrd 409')
406     
407! Ehouarn: forget about finvmaold
408!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
409
410!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
411!      DO l = 1, llm     
412!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
413!      ENDDO
414!c$OMP END DO NOWAIT
415
416      endif ! of if (planet_type.eq."earth")
417
418c
419c
420c     .....   FIN  de l'integration  de   q    .......
421
42215    continue
423          !write(*,*) 'integrd 410'
424
425c$OMP DO SCHEDULE(STATIC)
426      DO ij=ijb,ije 
427        ps0(ij)=ps(ij)
428      ENDDO
429c$OMP END DO NOWAIT
430
431c    .................................................................
432
433
434      IF( leapf )  THEN
435c       CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
436c       CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
437c$OMP DO SCHEDULE(STATIC)
438      DO ij=ijb,ije 
439        psm1(ij)=pscr(ij)
440      ENDDO
441c$OMP END DO NOWAIT
442
443c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
444          DO l = 1, llm
445            massem1(ijb:ije,l)=massescr(ijb:ije,l)
446          ENDDO
447c$OMP END DO NOWAIT         
448      END IF
449c$OMP BARRIER
450      RETURN
451      END
Note: See TracBrowser for help on using the repository browser.