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
RevLine 
[5099]1
[1632]2! $Id: integrd_p.F 1299 2010-01-20 14:27:21Z fairhead $
[5099]3
[5105]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
[5106]9  USE lmdz_filtreg_p
[5105]10  USE write_field_loc
[5117]11  USE lmdz_write_field
[5105]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
[5117]17  USE lmdz_strings, ONLY: int2str
[5118]18  USE lmdz_iniprint, ONLY: lunout, prt_level
[5123]19  USE lmdz_ssum_scopy, ONLY: ssum
[5136]20  USE lmdz_comgeom
[1632]21
[5159]22  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
23  USE lmdz_paramet
[5105]24  IMPLICIT NONE
[1632]25
26
[5105]27  !=======================================================================
[5159]28
[5105]29  !   Auteur:  P. Le Van
30  !   -------
[5159]31
[5105]32  !   objet:
33  !   ------
[5159]34
[5105]35  !   Incrementation des tendances dynamiques
[5159]36
[5105]37  !=======================================================================
38  !-----------------------------------------------------------------------
39  !   Declarations:
40  !   -------------
[1632]41
42
[5159]43
44
[5105]45  !   Arguments:
46  !   ----------
[1632]47
[5117]48  INTEGER,INTENT(IN) :: nq ! number of tracers to handle in this routine
[1632]49
[5105]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
[5113]57  ! values at previous time step
[5105]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)
[5113]63  ! the tendencies to add
[5105]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
[1632]70
[5105]71  !   Local:
72  !   ------
[1632]73
[5105]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)
[1632]80
[5105]81  INTEGER :: l,ij,iq,i,j
[1632]82
[5105]83  INTEGER :: ijb,ije,jjb,jje
84  LOGICAL :: checksum
85  LOGICAL,SAVE :: checksum_all=.TRUE.
86  INTEGER :: stop_it
87  INTEGER :: ierr
[2270]88
[5116]89  !WRITE(*,*) 'integrd 88: entree, nq=',nq
[5105]90  !-----------------------------------------------------------------------
[1632]91
[5105]92!$OMP BARRIER
[5117]93  IF (pole_nord) THEN
[5105]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
[1632]103
[5117]104  IF (pole_sud) THEN
[5105]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
[1632]114
[5105]115  !    ............    integration  de       ps         ..............
[1632]116
[5105]117   ! CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
[1632]118
[5105]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
[2270]126
[5105]127!$OMP DO SCHEDULE(STATIC)
128  DO ij = ijb,ije
129   pscr (ij)    = ps0(ij)
130   ps (ij)      = psm1(ij) + dt * dp(ij)
[2270]131
[5105]132  END DO
[1632]133
[5105]134!$OMP END DO
135!$OMP BARRIER
136  ! --> ici synchro OPENMP pour ps
[1673]137
[5105]138  checksum=.TRUE.
139  stop_it=0
[1632]140
[5105]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
[1632]150
[5105]151   ! CALL MPI_ALLREDUCE(checksum,checksum_all,1,
152  ! &                   MPI_LOGICAL,MPI_LOR,COMM_LMDZ,ierr)
153  IF( .NOT. checksum ) THEN
[5116]154     WRITE(lunout,*) "integrd: ps = ", ps(stop_it)
155     WRITE(lunout,*) " at node ij =", stop_it
[5113]156     ! since ij=j+(i-1)*jjp1 , we have
[5105]157      j=modulo(stop_it,jjp1)
158      i=1+(stop_it-j)/jjp1
[5116]159      WRITE(lunout,*) " lon = ",rlonv(i)*180./pi, " deg", &
[5105]160            " lat = ",rlatu(j)*180./pi, " deg"
161     CALL abort_gcm("integrd_loc", "negative surface pressure", 1)
162  ENDIF
[2270]163
[5105]164!$OMP END MASTER
165!$OMP BARRIER
[5116]166    !WRITE(*,*) 'integrd 170'
[5105]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)
[1632]179
[5105]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)
[5158]185    DO j=1,nq
[5105]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
[5099]193
[1632]194
[5159]195
[5116]196  !   !WRITE(*,*) 'integrd 200'
[5105]197!$OMP MASTER
[5117]198  IF (pole_nord) THEN
[1632]199
[5105]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
[1632]207
[5105]208  ENDIF
[1632]209
[5117]210  IF (pole_sud) THEN
[1632]211
[5105]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
[1632]219
[5105]220  ENDIF
221!$OMP END MASTER
222!$OMP BARRIER
[5116]223  !WRITE(*,*) 'integrd 217'
[5159]224
[5105]225  !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
226  !
[1632]227
[5105]228  CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
[1632]229
[5105]230!$OMP BARRIER
231  CALL massdair_loc (     p  , masse         )
[1632]232
[5105]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
[1632]238
[5105]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
[1632]244
[5105]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  !
[1632]250
[5105]251  !    ............   integration  de  ucov, vcov,  h     ..............
[2270]252
[5105]253!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
254  DO l = 1,llm
[2270]255
[5105]256  ijb=ij_begin
257  ije=ij_end
[5117]258  IF (pole_nord) ijb=ij_begin+iip1
259  IF (pole_sud)  ije=ij_end-iip1
[1632]260
[5105]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
[1632]265
[5105]266  ijb=ij_begin
267  ije=ij_end
[5117]268  IF (pole_sud)  ije=ij_end-iip1
[1632]269
[5105]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
[1632]274
[5105]275  ijb=ij_begin
276  ije=ij_end
[2270]277
[5105]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
[1632]283
[5105]284  !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
[5159]285
286
[5116]287  !   !WRITE(*,*) 'integrd 291'
[5105]288  IF (pole_nord) THEN
[1673]289
[5105]290    DO  ij   = 1, iim
291      tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
292    ENDDO
293      tpn      = SSUM(iim,tppn,1)/apoln
[1673]294
[5105]295    DO ij   = 1, iip1
296      teta(   ij   ,l)  = tpn
297    ENDDO
[1632]298
[5105]299  ENDIF
[1632]300
[5105]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)
[5117]323    IF (pole_sud) ije=ij_end-iip1
[5105]324    vcovm1(ijb:ije,l)=vscr(ijb:ije)
325
326  END IF
327
328  END DO
329!$OMP END DO NOWAIT
330
[5159]331
[5105]332  !   .......  integration de   q   ......
[5159]333
[5105]334  ijb=ij_begin
335  ije=ij_end
336
[5117]337     IF (planet_type=="earth") THEN
[5105]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
[1632]345      ENDDO
346
[5105]347!$OMP END DO NOWAIT
348!$OMP BARRIER
[1632]349
[5105]350    CALL check_isotopes(q,ijb,ije,'integrd 342')
[1632]351
[5116]352    !WRITE(*,*) 'integrd 341'
[5105]353    CALL qminimum_loc( q, nq, deltap )
[5116]354    !WRITE(*,*) 'integrd 343'
[5105]355
356    CALL check_isotopes(q,ijb,ije,'integrd 346')
[5159]357
[5105]358  !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
[5159]359
[5105]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
[1632]377      ENDDO
[5105]378!$OMP END DO NOWAIT
[1632]379
[5105]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
[5117]418  ENDIF ! of if (planet_type.EQ."earth")
[5105]419
[5159]420
421
[5105]422  ! .....   FIN  de l'integration  de   q    .......
423
[5116]424      !WRITE(*,*) 'integrd 410'
[5105]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.