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

Last change on this file since 5158 was 5158, checked in by abarral, 11 months ago

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables

  • 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  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
43  !   Arguments:
44  !   ----------
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
69  !   Local:
70  !   ------
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  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
88  !-----------------------------------------------------------------------
89
90!$OMP BARRIER
91  IF (pole_nord) THEN
92!$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
99!$OMP END DO NOWAIT
100  ENDIF
101
102  IF (pole_sud) THEN
103!$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
110!$OMP END DO NOWAIT
111  ENDIF
112
113  !    ............    integration  de       ps         ..............
114
115   ! CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
116
117  ijb=ij_begin
118  ije=ij_end
119!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
120  DO  l = 1,llm
121    massescr(ijb:ije,l)=masse(ijb:ije,l)
122  ENDDO
123!$OMP END DO NOWAIT
124
125!$OMP DO SCHEDULE(STATIC)
126  DO ij = ijb,ije
127   pscr (ij)    = ps0(ij)
128   ps (ij)      = psm1(ij) + dt * dp(ij)
129
130  END DO
131
132!$OMP END DO
133!$OMP BARRIER
134  ! --> ici synchro OPENMP pour ps
135
136  checksum=.TRUE.
137  stop_it=0
138
139!$OMP MASTER
140  !c$OMP DO SCHEDULE(STATIC)
141  DO ij = ijb,ije
142     IF( ps(ij)<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
162!$OMP END MASTER
163!$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
193  !
194  !   !WRITE(*,*) 'integrd 200'
195!$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
219!$OMP END MASTER
220!$OMP BARRIER
221  !WRITE(*,*) 'integrd 217'
222  !
223  !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
224  !
225
226  CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
227
228!$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  )
247  !
248
249  !    ............   integration  de  ucov, vcov,  h     ..............
250
251!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
252  DO 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 ij = ijb,ije
260  uscr( ij )   =  ucov( ij,l )
261  ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
262  END DO
263
264  ijb=ij_begin
265  ije=ij_end
266  IF (pole_sud)  ije=ij_end-iip1
267
268  DO ij = ijb,ije
269  vscr( ij )   =  vcov( ij,l )
270  vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
271  END DO
272
273  ijb=ij_begin
274  ije=ij_end
275
276  DO 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  END DO
281
282  !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
283  !
284  !
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
311  !
312
313  IF(leapf)  THEN
314      ! CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
315      ! CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
316      ! 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  END DO
327!$OMP END DO NOWAIT
328
329  !
330  !   .......  integration de   q   ......
331  !
332  ijb=ij_begin
333  ije=ij_end
334
335     IF (planet_type=="earth") THEN
336  ! Earth-specific treatment of first 2 tracers (water)
337!$OMP BARRIER
338!$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
345!$OMP END DO NOWAIT
346!$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')
355  !
356  !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
357  !
358!$OMP BARRIER
359  IF (pole_nord) THEN
360
361    DO iq = 1, nq
362
363!$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
376!$OMP END DO NOWAIT
377
378    ENDDO
379
380  ENDIF
381
382  IF (pole_sud) THEN
383
384    DO iq = 1, nq
385
386!$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
399!$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
418  !
419  !
420  ! .....   FIN  de l'integration  de   q    .......
421
422      !WRITE(*,*) 'integrd 410'
423
424!$OMP DO SCHEDULE(STATIC)
425  DO ij=ijb,ije
426    ps0(ij)=ps(ij)
427  ENDDO
428!$OMP END DO NOWAIT
429
430  !    .................................................................
431
432
433  IF( leapf )  THEN
434    ! CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
435    ! CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
436!$OMP DO SCHEDULE(STATIC)
437  DO ij=ijb,ije
438    psm1(ij)=pscr(ij)
439  ENDDO
440!$OMP END DO NOWAIT
441
442!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
443      DO l = 1, llm
444        massem1(ijb:ije,l)=massescr(ijb:ije,l)
445      ENDDO
446!$OMP END DO NOWAIT
447  END IF
448!$OMP BARRIER
449
450END SUBROUTINE integrd_loc
Note: See TracBrowser for help on using the repository browser.