source: LMDZ6/trunk/libf/dyn3dmem/integrd_loc.f90 @ 5281

Last change on this file since 5281 was 5281, checked in by abarral, 8 hours ago

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