source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/integrd_loc.f90 @ 5133

Last change on this file since 5133 was 5128, checked in by abarral, 5 months ago

Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments.
Merge the three different versions of abort_gcm into one
Fix seq, para 3D compilation broken from r5107 onwards
(lint) usual + Remove uneeded fixed-form continuations

  • 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: 10.4 KB
Line 
1
2! $Id: integrd_p.F 1299 2010-01-20 14:27:21Z fairhead $
3
4SUBROUTINE 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 lmdz_filtreg_p
10  USE write_field_loc
11  USE lmdz_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  USE lmdz_strings, ONLY: int2str
18  USE lmdz_iniprint, ONLY: lunout, prt_level
19  USE lmdz_ssum_scopy, ONLY: ssum
20
21
22  IMPLICIT NONE
23
24
25  !=======================================================================
26  !
27  !   Auteur:  P. Le Van
28  !   -------
29  !
30  !   objet:
31  !   ------
32  !
33  !   Incrementation des tendances dynamiques
34  !
35  !=======================================================================
36  !-----------------------------------------------------------------------
37  !   Declarations:
38  !   -------------
39
40  include "dimensions.h"
41  include "paramet.h"
42  include "comgeom.h"
43
44  !   Arguments:
45  !   ----------
46
47  INTEGER,INTENT(IN) :: nq ! number of tracers to handle in this routine
48
49  REAL,INTENT(INOUT) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
50  REAL,INTENT(INOUT) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
51  REAL,INTENT(INOUT) :: teta(ijb_u:ije_u,llm) ! potential temperature
52  REAL,INTENT(INOUT) :: q(ijb_u:ije_u,llm,nq) ! advected tracers
53  REAL,INTENT(INOUT) :: ps0(ijb_u:ije_u) ! surface pressure
54  REAL,INTENT(INOUT) :: masse(ijb_u:ije_u,llm) ! atmospheric mass
55  REAL,INTENT(INOUT) :: phis(ijb_u:ije_u) ! ground geopotential !!! unused
56  ! values at previous time step
57  REAL,INTENT(INOUT) :: vcovm1(ijb_v:ije_v,llm)
58  REAL,INTENT(INOUT) :: ucovm1(ijb_u:ije_u,llm)
59  REAL,INTENT(INOUT) :: tetam1(ijb_u:ije_u,llm)
60  REAL,INTENT(INOUT) :: psm1(ijb_u:ije_u)
61  REAL,INTENT(INOUT) :: massem1(ijb_u:ije_u,llm)
62  ! the tendencies to add
63  REAL,INTENT(INOUT) :: dv(ijb_v:ije_v,llm)
64  REAL,INTENT(INOUT) :: du(ijb_u:ije_u,llm)
65  REAL,INTENT(INOUT) :: dteta(ijb_u:ije_u,llm)
66  REAL,INTENT(INOUT) :: dp(ijb_u:ije_u)
67  REAL,INTENT(INOUT) :: dq(ijb_u:ije_u,llm,nq) !!! unused
68   ! REAL,INTENT(INOUT) ::finvmaold(ijb_u:ije_u,llm) !!! unused
69
70  !   Local:
71  !   ------
72
73  REAL :: vscr( ijb_v:ije_v ),uscr( ijb_u:ije_u )
74  REAL :: hscr( ijb_u:ije_u ),pscr(ijb_u:ije_u)
75  REAL :: massescr( ijb_u:ije_u,llm )
76   ! REAL finvmasse(ijb_u:ije_u,llm)
77  REAL :: tpn,tps,tppn(iim),tpps(iim)
78  REAL :: qpn,qps,qppn(iim),qpps(iim)
79
80  INTEGER :: l,ij,iq,i,j
81
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
89  !-----------------------------------------------------------------------
90
91!$OMP BARRIER
92  IF (pole_nord) THEN
93!$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
100!$OMP END DO NOWAIT
101  ENDIF
102
103  IF (pole_sud) THEN
104!$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
111!$OMP END DO NOWAIT
112  ENDIF
113
114  !    ............    integration  de       ps         ..............
115
116   ! CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
117
118  ijb=ij_begin
119  ije=ij_end
120!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
121  DO  l = 1,llm
122    massescr(ijb:ije,l)=masse(ijb:ije,l)
123  ENDDO
124!$OMP END DO NOWAIT
125
126!$OMP DO SCHEDULE(STATIC)
127  DO ij = ijb,ije
128   pscr (ij)    = ps0(ij)
129   ps (ij)      = psm1(ij) + dt * dp(ij)
130
131  END DO
132
133!$OMP END DO
134!$OMP BARRIER
135  ! --> ici synchro OPENMP pour ps
136
137  checksum=.TRUE.
138  stop_it=0
139
140!$OMP MASTER
141  !c$OMP DO SCHEDULE(STATIC)
142  DO ij = ijb,ije
143     IF( ps(ij)<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
163!$OMP END MASTER
164!$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
194  !
195  !   !WRITE(*,*) 'integrd 200'
196!$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
220!$OMP END MASTER
221!$OMP BARRIER
222  !WRITE(*,*) 'integrd 217'
223  !
224  !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
225  !
226
227  CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
228
229!$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  )
248  !
249
250  !    ............   integration  de  ucov, vcov,  h     ..............
251
252!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
253  DO 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 ij = ijb,ije
261  uscr( ij )   =  ucov( ij,l )
262  ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
263  END DO
264
265  ijb=ij_begin
266  ije=ij_end
267  IF (pole_sud)  ije=ij_end-iip1
268
269  DO ij = ijb,ije
270  vscr( ij )   =  vcov( ij,l )
271  vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
272  END DO
273
274  ijb=ij_begin
275  ije=ij_end
276
277  DO 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  END DO
282
283  !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
284  !
285  !
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
312  !
313
314  IF(leapf)  THEN
315      ! CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
316      ! CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
317      ! 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  END DO
328!$OMP END DO NOWAIT
329
330  !
331  !   .......  integration de   q   ......
332  !
333  ijb=ij_begin
334  ije=ij_end
335
336     IF (planet_type=="earth") THEN
337  ! Earth-specific treatment of first 2 tracers (water)
338!$OMP BARRIER
339!$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
346!$OMP END DO NOWAIT
347!$OMP BARRIER
348
349    CALL check_isotopes(q,ijb,ije,'integrd 342')
350
351    !WRITE(*,*) 'integrd 341'
352    CALL qminimum_loc( q, nq, deltap )
353    !WRITE(*,*) 'integrd 343'
354
355    CALL check_isotopes(q,ijb,ije,'integrd 346')
356  !
357  !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
358  !
359!$OMP BARRIER
360  IF (pole_nord) THEN
361
362    DO iq = 1, nq
363
364!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
365      DO l = 1, llm
366
367         DO ij = 1, iim
368           qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
369         ENDDO
370           qpn  =  SSUM(iim,qppn,1)/apoln
371
372         DO ij = 1, iip1
373           q(   ij   ,l,iq)  = qpn
374         ENDDO
375
376      ENDDO
377!$OMP END DO NOWAIT
378
379    ENDDO
380
381  ENDIF
382
383  IF (pole_sud) THEN
384
385    DO iq = 1, nq
386
387!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
388      DO l = 1, llm
389
390         DO ij = 1, iim
391           qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
392         ENDDO
393           qps  =  SSUM(iim,qpps,1)/apols
394
395         DO ij = 1, iip1
396           q(ij+ip1jm,l,iq)  = qps
397         ENDDO
398
399      ENDDO
400!$OMP END DO NOWAIT
401
402    ENDDO
403
404  ENDIF
405
406  CALL check_isotopes(q,ijb,ije,'integrd 409')
407
408  ! Ehouarn: forget about finvmaold
409  !c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
410
411  !c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
412   ! DO l = 1, llm
413   !   finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)
414   ! ENDDO
415  !c$OMP END DO NOWAIT
416
417  ENDIF ! of if (planet_type.EQ."earth")
418
419  !
420  !
421  ! .....   FIN  de l'integration  de   q    .......
422
423      !WRITE(*,*) 'integrd 410'
424
425!$OMP DO SCHEDULE(STATIC)
426  DO ij=ijb,ije
427    ps0(ij)=ps(ij)
428  ENDDO
429!$OMP END DO NOWAIT
430
431  !    .................................................................
432
433
434  IF( leapf )  THEN
435    ! CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
436    ! CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
437!$OMP DO SCHEDULE(STATIC)
438  DO ij=ijb,ije
439    psm1(ij)=pscr(ij)
440  ENDDO
441!$OMP END DO NOWAIT
442
443!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
444      DO l = 1, llm
445        massem1(ijb:ije,l)=massescr(ijb:ije,l)
446      ENDDO
447!$OMP END DO NOWAIT
448  END IF
449!$OMP BARRIER
450
451END SUBROUTINE integrd_loc
Note: See TracBrowser for help on using the repository browser.