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

Last change on this file since 5248 was 5246, checked in by abarral, 20 hours ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

  • 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 mod_filtreg_p
10  USE write_field_loc
11  USE 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
18  IMPLICIT NONE
19
20
21  !=======================================================================
22  !
23  !   Auteur:  P. Le Van
24  !   -------
25  !
26  !   objet:
27  !   ------
28  !
29  !   Incrementation des tendances dynamiques
30  !
31  !=======================================================================
32  !-----------------------------------------------------------------------
33  !   Declarations:
34  !   -------------
35
36  include "dimensions.h"
37  include "paramet.h"
38  include "comgeom.h"
39  include "iniprint.h"
40
41  !   Arguments:
42  !   ----------
43
44  INTEGER,intent(in) :: nq ! number of tracers to handle in this routine
45
46  REAL,INTENT(INOUT) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
47  REAL,INTENT(INOUT) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
48  REAL,INTENT(INOUT) :: teta(ijb_u:ije_u,llm) ! potential temperature
49  REAL,INTENT(INOUT) :: q(ijb_u:ije_u,llm,nq) ! advected tracers
50  REAL,INTENT(INOUT) :: ps0(ijb_u:ije_u) ! surface pressure
51  REAL,INTENT(INOUT) :: masse(ijb_u:ije_u,llm) ! atmospheric mass
52  REAL,INTENT(INOUT) :: phis(ijb_u:ije_u) ! ground geopotential !!! unused
53  ! ! values at previous time step
54  REAL,INTENT(INOUT) :: vcovm1(ijb_v:ije_v,llm)
55  REAL,INTENT(INOUT) :: ucovm1(ijb_u:ije_u,llm)
56  REAL,INTENT(INOUT) :: tetam1(ijb_u:ije_u,llm)
57  REAL,INTENT(INOUT) :: psm1(ijb_u:ije_u)
58  REAL,INTENT(INOUT) :: massem1(ijb_u:ije_u,llm)
59  ! ! the tendencies to add
60  REAL,INTENT(INOUT) :: dv(ijb_v:ije_v,llm)
61  REAL,INTENT(INOUT) :: du(ijb_u:ije_u,llm)
62  REAL,INTENT(INOUT) :: dteta(ijb_u:ije_u,llm)
63  REAL,INTENT(INOUT) :: dp(ijb_u:ije_u)
64  REAL,INTENT(INOUT) :: dq(ijb_u:ije_u,llm,nq) !!! unused
65   ! REAL,INTENT(INOUT) ::finvmaold(ijb_u:ije_u,llm) !!! unused
66
67  !   Local:
68  !   ------
69
70  REAL :: vscr( ijb_v:ije_v ),uscr( ijb_u:ije_u )
71  REAL :: hscr( ijb_u:ije_u ),pscr(ijb_u:ije_u)
72  REAL :: massescr( ijb_u:ije_u,llm )
73   ! REAL finvmasse(ijb_u:ije_u,llm)
74  REAL :: tpn,tps,tppn(iim),tpps(iim)
75  REAL :: qpn,qps,qppn(iim),qpps(iim)
76
77  INTEGER :: l,ij,iq,i,j
78
79  REAL :: SSUM
80  EXTERNAL SSUM
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).LT.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.eq."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
42215   continue
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  RETURN
451END SUBROUTINE integrd_loc
Note: See TracBrowser for help on using the repository browser.