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

Last change on this file since 5225 was 5159, checked in by abarral, 3 months ago

Put dimensions.h and paramet.h into modules

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