source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3dpar/leapfrog_p.F @ 3693

Last change on this file since 3693 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 54.8 KB
Line 
1!
2! $Id: leapfrog_p.F 1520 2011-05-23 11:37:09Z fhourdin $
3!
4c
5c
6
7      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
8     &                    time_0)
9
10       USE misc_mod
11       USE parallel
12       USE times
13       USE mod_hallo
14       USE Bands
15       USE Write_Field
16       USE Write_Field_p
17       USE vampir
18       USE timer_filtre, ONLY : print_filtre_timer
19       USE infotrac
20       USE guide_p_mod, ONLY : guide_main
21       USE getparam
22       USE control_mod
23
24      IMPLICIT NONE
25
26c      ......   Version  du 10/01/98    ..........
27
28c             avec  coordonnees  verticales hybrides
29c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30
31c=======================================================================
32c
33c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
34c   -------
35c
36c   Objet:
37c   ------
38c
39c   GCM LMD nouvelle grille
40c
41c=======================================================================
42c
43c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
44c      et possibilite d'appeler une fonction f(y)  a derivee tangente
45c      hyperbolique a la  place de la fonction a derivee sinusoidale.
46
47c  ... Possibilite de choisir le shema pour l'advection de
48c        q  , en modifiant iadv dans traceur.def  (10/02) .
49c
50c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51c      Pour Van-Leer iadv=10
52c
53c-----------------------------------------------------------------------
54c   Declarations:
55c   -------------
56
57#include "dimensions.h"
58#include "paramet.h"
59#include "comconst.h"
60#include "comdissnew.h"
61#include "comvert.h"
62#include "comgeom.h"
63#include "logic.h"
64#include "temps.h"
65#include "ener.h"
66#include "description.h"
67#include "serre.h"
68!#include "com_io_dyn.h"
69#include "iniprint.h"
70#include "academic.h"
71     
72      INTEGER         longcles
73      PARAMETER     ( longcles = 20 )
74      REAL  clesphy0( longcles )
75
76      real zqmin,zqmax
77      INTEGER nbetatmoy, nbetatdem,nbetat
78
79c   variables dynamiques
80      REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
81      REAL :: teta(ip1jmp1,llm)                 ! temperature potentielle
82      REAL :: q(ip1jmp1,llm,nqtot)              ! champs advectes
83      REAL :: ps(ip1jmp1)                       ! pression  au sol
84      REAL,SAVE :: p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
85      REAL,SAVE :: pks(ip1jmp1)                      ! exner au  sol
86      REAL,SAVE :: pk(ip1jmp1,llm)                   ! exner au milieu des couches
87      REAL,SAVE :: pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
88      REAL :: masse(ip1jmp1,llm)                ! masse d'air
89      REAL :: phis(ip1jmp1)                     ! geopotentiel au sol
90      REAL,SAVE :: phi(ip1jmp1,llm)                  ! geopotentiel
91      REAL,SAVE :: w(ip1jmp1,llm)                    ! vitesse verticale
92
93c variables dynamiques intermediaire pour le transport
94      REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
95
96c   variables dynamiques au pas -1
97      REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
98      REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
99      REAL,SAVE :: massem1(ip1jmp1,llm)
100
101c   tendances dynamiques
102      REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
103      REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
104      REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
105
106c   tendances de la dissipation
107      REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
108      REAL,SAVE :: dtetadis(ip1jmp1,llm)
109
110c   tendances physiques
111      REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
112      REAL,SAVE :: dtetafi(ip1jmp1,llm)
113      REAL,SAVE :: dpfi(ip1jmp1)
114      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
115
116c   variables pour le fichier histoire
117      REAL dtav      ! intervalle de temps elementaire
118
119      REAL tppn(iim),tpps(iim),tpn,tps
120c
121      INTEGER itau,itaufinp1,iav
122!      INTEGER  iday ! jour julien
123      REAL       time
124
125      REAL  SSUM
126      REAL time_0
127      REAL,SAVE :: finvmaold(ip1jmp1,llm)
128
129cym      LOGICAL  lafin
130      LOGICAL :: lafin
131      INTEGER ij,iq,l
132      INTEGER ik
133
134      real time_step, t_wrt, t_ops
135
136! jD_cur: jour julien courant
137! jH_cur: heure julienne courante
138      REAL :: jD_cur, jH_cur
139      INTEGER :: an, mois, jour
140      REAL :: secondes
141
142      LOGICAL first,callinigrads
143
144      data callinigrads/.true./
145      character*10 string10
146
147      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
148      REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
149
150c+jld variables test conservation energie
151      REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
152C     Tendance de la temp. potentiel d (theta)/ d t due a la
153C     tansformation d'energie cinetique en energie thermique
154C     cree par la dissipation
155      REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
156      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
157      REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
158      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
159      CHARACTER*15 ztit
160!      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
161!      SAVE      ip_ebil_dyn
162!      DATA      ip_ebil_dyn/0/
163c-jld
164
165      character*80 dynhist_file, dynhistave_file
166      character(len=*),parameter :: modname="leapfrog"
167      character*80 abort_message
168
169
170      logical,PARAMETER :: dissip_conservative=.TRUE.
171 
172      INTEGER testita
173      PARAMETER (testita = 9)
174
175      logical , parameter :: flag_verif = .false.
176     
177c declaration liees au parallelisme
178      INTEGER :: ierr
179      LOGICAL :: FirstCaldyn
180      LOGICAL :: FirstPhysic
181      INTEGER :: ijb,ije,j,i
182      type(Request) :: TestRequest
183      type(Request) :: Request_Dissip
184      type(Request) :: Request_physic
185      REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
186      REAL,SAVE :: dtetafi_tmp(iip1,llm)
187      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
188      REAL,SAVE :: dpfi_tmp(iip1)
189
190      INTEGER :: true_itau
191      LOGICAL :: verbose=.true.
192      INTEGER :: iapptrac
193      INTEGER :: AdjustCount
194!      INTEGER :: var_time
195      LOGICAL :: ok_start_timer=.FALSE.
196      LOGICAL, SAVE :: firstcall=.TRUE.
197
198c$OMP MASTER
199      ItCount=0
200c$OMP END MASTER     
201      true_itau=0
202      FirstCaldyn=.TRUE.
203      FirstPhysic=.TRUE.
204      iapptrac=0
205      AdjustCount = 0
206      lafin=.false.
207     
208      itaufin   = nday*day_step
209      itaufinp1 = itaufin +1
210
211      itau = 0
212!      iday = day_ini+itau/day_step
213!      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
214!         IF(time.GT.1.) THEN
215!          time = time-1.
216!          iday = iday+1
217!         ENDIF
218
219c Allocate variables depending on dynamic variable nqtot
220c$OMP MASTER
221         IF (firstcall) THEN
222            firstcall=.FALSE.
223            ALLOCATE(dq(ip1jmp1,llm,nqtot))
224            ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
225            ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
226         END IF
227c$OMP END MASTER     
228c$OMP BARRIER
229
230c-----------------------------------------------------------------------
231c   On initialise la pression et la fonction d'Exner :
232c   --------------------------------------------------
233
234c$OMP MASTER
235      dq(:,:,:)=0.
236      CALL pression ( ip1jmp1, ap, bp, ps, p       )
237      if (disvert_type==1) then
238        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
239      else ! we assume that we are in the disvert_type==2 case
240        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
241      endif
242c$OMP END MASTER
243c-----------------------------------------------------------------------
244c   Debut de l'integration temporelle:
245c   ----------------------------------
246c et du parallelisme !!
247
248   1  CONTINUE
249
250      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
251      jH_cur = jH_ref +                                                 &
252     &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
253
254
255#ifdef CPP_IOIPSL
256      if (ok_guide) then
257!$OMP MASTER
258        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
259!$OMP END MASTER
260!$OMP BARRIER
261      endif
262#endif
263
264c
265c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
266c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
267c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
268c     ENDIF
269c
270cym      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
271cym      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
272cym      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
273cym      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
274cym      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
275
276       if (FirstCaldyn) then
277c$OMP MASTER
278         ucovm1=ucov
279         vcovm1=vcov
280         tetam1= teta
281         massem1= masse
282         psm1= ps
283         
284         finvmaold = masse
285         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
286c$OMP END MASTER
287c$OMP BARRIER
288       else
289! Save fields obtained at previous time step as '...m1'
290         ijb=ij_begin
291         ije=ij_end
292
293c$OMP MASTER           
294         psm1     (ijb:ije) = ps    (ijb:ije)
295c$OMP END MASTER
296
297c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
298         DO l=1,llm     
299           ije=ij_end
300           ucovm1   (ijb:ije,l) = ucov  (ijb:ije,l)
301           tetam1   (ijb:ije,l) = teta  (ijb:ije,l)
302           massem1  (ijb:ije,l) = masse (ijb:ije,l)
303           finvmaold(ijb:ije,l)=masse(ijb:ije,l)
304                 
305           if (pole_sud) ije=ij_end-iip1
306           vcovm1(ijb:ije,l) = vcov  (ijb:ije,l)
307       
308
309         ENDDO
310c$OMP ENDDO 
311
312
313          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
314     .                    llm, -2,2, .TRUE., 1 )
315
316       endif ! of if (FirstCaldyn)
317       
318      forward = .TRUE.
319      leapf   = .FALSE.
320      dt      =  dtvr
321
322c   ...    P.Le Van .26/04/94  ....
323
324cym      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
325cym      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
326
327cym  ne sert a rien
328cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
329
330   2  CONTINUE
331
332c$OMP MASTER
333      ItCount=ItCount+1
334      if (MOD(ItCount,1)==1) then
335        debug=.true.
336      else
337        debug=.false.
338      endif
339c$OMP END MASTER
340c-----------------------------------------------------------------------
341
342c   date:
343c   -----
344
345
346c   gestion des appels de la physique et des dissipations:
347c   ------------------------------------------------------
348c
349c   ...    P.Le Van  ( 6/02/95 )  ....
350
351      apphys = .FALSE.
352      statcl = .FALSE.
353      conser = .FALSE.
354      apdiss = .FALSE.
355
356      IF( purmats ) THEN
357      ! Purely Matsuno time stepping
358         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
359         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
360     s        apdiss = .TRUE.
361         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
362     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
363      ELSE
364      ! Leapfrog/Matsuno time stepping
365         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
366         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
367     s        apdiss = .TRUE.
368         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
369      END IF
370
371! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
372!          supress dissipation step
373      if (llm.eq.1) then
374        apdiss=.false.
375      endif
376
377cym    ---> Pour le moment     
378cym      apphys = .FALSE.
379      statcl = .FALSE.
380      conser = .FALSE. ! ie: no output of control variables to stdout in //
381     
382      if (firstCaldyn) then
383c$OMP MASTER
384          call SetDistrib(jj_Nb_Caldyn)
385c$OMP END MASTER
386c$OMP BARRIER
387          firstCaldyn=.FALSE.
388cym          call InitTime
389c$OMP MASTER
390          call Init_timer
391c$OMP END MASTER
392      endif
393
394c$OMP MASTER     
395      IF (ok_start_timer) THEN
396        CALL InitTime
397        ok_start_timer=.FALSE.
398      ENDIF     
399c$OMP END MASTER     
400     
401      if (Adjust) then
402c$OMP MASTER
403        AdjustCount=AdjustCount+1
404        if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
405     &         .and. itau/iphysiq>2 .and. Adjustcount>30) then
406           AdjustCount=0
407           call allgather_timer_average
408
409        if (Verbose) then
410       
411        print *,'*********************************'
412        print *,'******    TIMER CALDYN     ******'
413        do i=0,mpi_size-1
414          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
415     &            '  : temps moyen :',
416     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i),
417     &            '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
418        enddo
419     
420        print *,'*********************************'
421        print *,'******    TIMER VANLEER    ******'
422        do i=0,mpi_size-1
423          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
424     &            '  : temps moyen :',
425     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i),
426     &            '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
427        enddo
428     
429        print *,'*********************************'
430        print *,'******    TIMER DISSIP    ******'
431        do i=0,mpi_size-1
432          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
433     &            '  : temps moyen :',
434     &             timer_average(jj_nb_dissip(i),timer_dissip,i),
435     &             '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
436        enddo
437       
438        if (mpi_rank==0) call WriteBands
439       
440       endif
441       
442         call AdjustBands_caldyn
443         if (mpi_rank==0) call WriteBands
444         
445         call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
446     &                                jj_Nb_caldyn,0,0,TestRequest)
447         call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
448     &                                jj_Nb_caldyn,0,0,TestRequest)
449         call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
450     &                                jj_Nb_caldyn,0,0,TestRequest)
451         call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
452     &                                jj_Nb_caldyn,0,0,TestRequest)
453         call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
454     &                                jj_Nb_caldyn,0,0,TestRequest)
455         call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
456     &                                jj_Nb_caldyn,0,0,TestRequest)
457         call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
458     &                                jj_Nb_caldyn,0,0,TestRequest)
459         call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
460     &                                jj_Nb_caldyn,0,0,TestRequest)
461         call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
462     &                                jj_Nb_caldyn,0,0,TestRequest)
463         call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
464     &                                jj_Nb_caldyn,0,0,TestRequest)
465         call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
466     &                                jj_Nb_caldyn,0,0,TestRequest)
467         call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
468     &                                jj_Nb_caldyn,0,0,TestRequest)
469         call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
470     &                                jj_Nb_caldyn,0,0,TestRequest)
471         call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
472     &                                jj_Nb_caldyn,0,0,TestRequest)
473         call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
474     &                                jj_Nb_caldyn,0,0,TestRequest)
475         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
476     &                                jj_Nb_caldyn,0,0,TestRequest)
477 
478        do j=1,nqtot
479         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
480     &                                jj_nb_caldyn,0,0,TestRequest)
481        enddo
482
483         call SetDistrib(jj_nb_caldyn)
484         call SendRequest(TestRequest)
485         call WaitRequest(TestRequest)
486         
487        call AdjustBands_dissip
488        call AdjustBands_physic
489
490      endif
491c$OMP END MASTER 
492      endif       
493     
494     
495     
496c-----------------------------------------------------------------------
497c   calcul des tendances dynamiques:
498c   --------------------------------
499c$OMP BARRIER
500c$OMP MASTER
501       call VTb(VThallo)
502c$OMP END MASTER
503
504       call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
505       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
506       call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
507       call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
508       call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
509       call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
510       call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
511       call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
512       
513c       do j=1,nqtot
514c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
515c     *                       TestRequest)
516c        enddo
517
518       call SendRequest(TestRequest)
519c$OMP BARRIER
520       call WaitRequest(TestRequest)
521
522c$OMP MASTER
523       call VTe(VThallo)
524c$OMP END MASTER
525c$OMP BARRIER
526     
527      if (debug) then       
528!$OMP BARRIER
529!$OMP MASTER
530        call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
531        call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
532        call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
533        call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
534        call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
535        call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
536        call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
537        call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
538        call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
539        do j=1,nqtot
540          call WriteField_p('q'//trim(int2str(j)),
541     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
542        enddo
543!$OMP END MASTER       
544c$OMP BARRIER
545      endif
546
547     
548      True_itau=True_itau+1
549
550c$OMP MASTER
551      IF (prt_level>9) THEN
552        WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
553      ENDIF
554
555
556      call start_timer(timer_caldyn)
557
558      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
559
560     
561      call VTb(VTcaldyn)
562c$OMP END MASTER
563!      var_time=time+iday-day_ini
564
565c$OMP BARRIER
566!      CALL FTRACE_REGION_BEGIN("caldyn")
567      time = jD_cur + jH_cur
568      CALL caldyn_p
569     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
570     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
571
572!      CALL FTRACE_REGION_END("caldyn")
573
574c$OMP MASTER
575      call VTe(VTcaldyn)
576c$OMP END MASTER     
577
578cc$OMP BARRIER
579cc$OMP MASTER
580!      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
581!      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
582!      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
583!      call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
584!      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
585!      call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
586!      call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
587!      call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
588!      call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
589!      call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
590cc$OMP END MASTER
591
592c-----------------------------------------------------------------------
593c   calcul des tendances advection des traceurs (dont l'humidite)
594c   -------------------------------------------------------------
595
596      IF( forward. OR . leapf )  THEN
597cc$OMP PARALLEL DEFAULT(SHARED)
598c
599         CALL caladvtrac_p(q,pbaru,pbarv,
600     *        p, masse, dq,  teta,
601     .        flxw,pk, iapptrac)
602
603C        Stokage du flux de masse pour traceurs OFF-LINE
604         IF (offline .AND. .NOT. adjust) THEN
605            CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
606     .           dtvr, itau)
607         ENDIF
608
609      ENDIF ! of IF( forward. OR . leapf )
610cc$OMP END PARALLEL
611
612c-----------------------------------------------------------------------
613c   integrations dynamique et traceurs:
614c   ----------------------------------
615
616c$OMP MASTER
617       call VTb(VTintegre)
618c$OMP END MASTER
619c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
620c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
621c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
622c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
623c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
624c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
625c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
626c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
627cc$OMP PARALLEL DEFAULT(SHARED)
628c$OMP BARRIER
629!       CALL FTRACE_REGION_BEGIN("integrd")
630
631       CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
632     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
633     $              finvmaold                                    )
634
635!       CALL FTRACE_REGION_END("integrd")
636c$OMP BARRIER
637cc$OMP MASTER
638c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
639c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
640c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
641c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
642c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
643c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
644c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
645c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
646c
647c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
648c      do j=1,nqtot
649c        call WriteField_p('q'//trim(int2str(j)),
650c     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
651c        call WriteField_p('dq'//trim(int2str(j)),
652c     .                reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
653c      enddo
654cc$OMP END MASTER
655
656
657c$OMP MASTER
658       call VTe(VTintegre)
659c$OMP END MASTER
660c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
661c
662c-----------------------------------------------------------------------
663c   calcul des tendances physiques:
664c   -------------------------------
665c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
666c
667       IF( purmats )  THEN
668          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
669       ELSE
670          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
671       ENDIF
672
673cc$OMP END PARALLEL
674
675c
676c
677       IF( apphys )  THEN
678c
679c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
680c
681cc$OMP PARALLEL DEFAULT(SHARED)
682cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
683
684c$OMP MASTER
685         call suspend_timer(timer_caldyn)
686
687        if (prt_level >= 10) then
688         write(lunout,*)
689     &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
690        endif
691c$OMP END MASTER
692
693         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )
694
695c$OMP BARRIER
696         if (disvert_type==1) then
697           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
698         else ! we assume that we are in the disvert_type==2 case
699           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
700         endif
701c$OMP BARRIER
702           jD_cur = jD_ref + day_ini - day_ref
703     $        + int (itau * dtvr / daysec)
704           jH_cur = jH_ref +                                            &
705     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
706!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
707
708c rajout debug
709c       lafin = .true.
710
711
712c   Inbterface avec les routines de phylmd (phymars ... )
713c   -----------------------------------------------------
714
715c+jld
716
717c  Diagnostique de conservation de l'energie : initialisation
718      IF (ip_ebil_dyn.ge.1 ) THEN
719          ztit='bil dyn'
720! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
721           IF (planet_type.eq."earth") THEN
722            CALL diagedyn(ztit,2,1,1,dtphys
723     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
724           ENDIF
725      ENDIF
726c-jld
727c$OMP BARRIER
728c$OMP MASTER
729        call VTb(VThallo)
730c$OMP END MASTER
731
732        call SetTag(Request_physic,800)
733       
734        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
735     *                               jj_Nb_physic,2,2,Request_physic)
736       
737        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
738     *                               jj_Nb_physic,2,2,Request_physic)
739       
740        call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
741     *                               jj_Nb_physic,2,2,Request_physic)
742       
743        call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
744     *                               jj_Nb_physic,1,2,Request_physic)
745
746        call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
747     *                               jj_Nb_physic,2,2,Request_physic)
748       
749        call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
750     *                               jj_Nb_physic,2,2,Request_physic)
751       
752        call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
753     *                               jj_Nb_physic,2,2,Request_physic)
754       
755        call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
756     *                               jj_Nb_physic,2,2,Request_physic)
757       
758        call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
759     *                               jj_Nb_physic,2,2,Request_physic)
760       
761c        call SetDistrib(jj_nb_vanleer)
762        do j=1,nqtot
763 
764          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
765     *                               jj_Nb_physic,2,2,Request_physic)
766        enddo
767
768        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
769     *                               jj_Nb_physic,2,2,Request_physic)
770       
771        call SendRequest(Request_Physic)
772c$OMP BARRIER
773        call WaitRequest(Request_Physic)       
774
775c$OMP BARRIER
776c$OMP MASTER
777        call SetDistrib(jj_nb_Physic)
778        call VTe(VThallo)
779       
780        call VTb(VTphysiq)
781c$OMP END MASTER
782c$OMP BARRIER
783
784cc$OMP MASTER       
785c      call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
786c      call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
787c      call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
788c      call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
789c      call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
790cc$OMP END MASTER
791cc$OMP BARRIER
792!        CALL FTRACE_REGION_BEGIN("calfis")
793        CALL calfis_p(lafin ,jD_cur, jH_cur,
794     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
795     $               du,dv,dteta,dq,
796     $               flxw,
797     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
798!        CALL FTRACE_REGION_END("calfis")
799        ijb=ij_begin
800        ije=ij_end 
801        if ( .not. pole_nord) then
802c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
803          DO l=1,llm
804          dufi_tmp(1:iip1,l)   = dufi(ijb:ijb+iim,l)
805          dvfi_tmp(1:iip1,l)   = dvfi(ijb:ijb+iim,l) 
806          dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l) 
807          dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:) 
808          ENDDO
809c$OMP END DO NOWAIT
810
811c$OMP MASTER
812          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim) 
813c$OMP END MASTER
814        endif ! of if ( .not. pole_nord)
815
816c$OMP BARRIER
817c$OMP MASTER
818        call SetDistrib(jj_nb_Physic_bis)
819
820        call VTb(VThallo)
821c$OMP END MASTER
822c$OMP BARRIER
823 
824        call Register_Hallo(dufi,ip1jmp1,llm,
825     *                      1,0,0,1,Request_physic)
826       
827        call Register_Hallo(dvfi,ip1jm,llm,
828     *                      1,0,0,1,Request_physic)
829       
830        call Register_Hallo(dtetafi,ip1jmp1,llm,
831     *                      1,0,0,1,Request_physic)
832
833        call Register_Hallo(dpfi,ip1jmp1,1,
834     *                      1,0,0,1,Request_physic)
835
836        do j=1,nqtot
837          call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
838     *                        1,0,0,1,Request_physic)
839        enddo
840       
841        call SendRequest(Request_Physic)
842c$OMP BARRIER
843        call WaitRequest(Request_Physic)
844             
845c$OMP BARRIER
846c$OMP MASTER
847        call VTe(VThallo)
848 
849        call SetDistrib(jj_nb_Physic)
850c$OMP END MASTER
851c$OMP BARRIER       
852                ijb=ij_begin
853        if (.not. pole_nord) then
854       
855c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
856          DO l=1,llm
857            dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
858            dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
859            dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
860     &                              +dtetafi_tmp(1:iip1,l)
861            dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
862     &                              + dqfi_tmp(1:iip1,l,:)
863          ENDDO
864c$OMP END DO NOWAIT
865
866c$OMP MASTER
867          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
868c$OMP END MASTER
869         
870        endif ! of if (.not. pole_nord)
871c$OMP BARRIER
872cc$OMP MASTER       
873c      call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
874c      call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
875c      call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
876c      call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
877cc$OMP END MASTER
878c     
879c      do j=1,nqtot
880c        call WriteField_p('dqfi'//trim(int2str(j)),
881c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
882c      enddo
883
884c      ajout des tendances physiques:
885c      ------------------------------
886         IF (ok_strato) THEN
887           CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
888         ENDIF
889       
890          CALL addfi_p( dtphys, leapf, forward   ,
891     $                  ucov, vcov, teta , q   ,ps ,
892     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
893
894c$OMP BARRIER
895c$OMP MASTER
896        call VTe(VTphysiq)
897
898        call VTb(VThallo)
899c$OMP END MASTER
900
901        call SetTag(Request_physic,800)
902        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
903     *                               jj_Nb_caldyn,Request_physic)
904       
905        call Register_SwapField(vcov,vcov,ip1jm,llm,
906     *                               jj_Nb_caldyn,Request_physic)
907       
908        call Register_SwapField(teta,teta,ip1jmp1,llm,
909     *                               jj_Nb_caldyn,Request_physic)
910       
911        call Register_SwapField(masse,masse,ip1jmp1,llm,
912     *                               jj_Nb_caldyn,Request_physic)
913
914        call Register_SwapField(p,p,ip1jmp1,llmp1,
915     *                               jj_Nb_caldyn,Request_physic)
916       
917        call Register_SwapField(pk,pk,ip1jmp1,llm,
918     *                               jj_Nb_caldyn,Request_physic)
919       
920        call Register_SwapField(phis,phis,ip1jmp1,1,
921     *                               jj_Nb_caldyn,Request_physic)
922       
923        call Register_SwapField(phi,phi,ip1jmp1,llm,
924     *                               jj_Nb_caldyn,Request_physic)
925       
926        call Register_SwapField(w,w,ip1jmp1,llm,
927     *                               jj_Nb_caldyn,Request_physic)
928
929        do j=1,nqtot
930       
931          call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
932     *                               jj_Nb_caldyn,Request_physic)
933       
934        enddo
935
936        call SendRequest(Request_Physic)
937c$OMP BARRIER
938        call WaitRequest(Request_Physic)     
939
940c$OMP BARRIER
941c$OMP MASTER
942       call VTe(VThallo)
943       call SetDistrib(jj_Nb_caldyn)
944c$OMP END MASTER
945c$OMP BARRIER
946c
947c  Diagnostique de conservation de l'energie : difference
948      IF (ip_ebil_dyn.ge.1 ) THEN
949          ztit='bil phys'
950          CALL diagedyn(ztit,2,1,1,dtphys
951     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
952      ENDIF
953
954cc$OMP MASTER     
955c      if (debug) then
956c       call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
957c       call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
958c       call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
959c      endif
960cc$OMP END MASTER
961
962
963c-jld
964c$OMP MASTER
965         call resume_timer(timer_caldyn)
966         if (FirstPhysic) then
967           ok_start_timer=.TRUE.
968           FirstPhysic=.false.
969         endif
970c$OMP END MASTER
971       ENDIF ! of IF( apphys )
972
973      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
974!   Academic case : Simple friction and Newtonan relaxation
975!   -------------------------------------------------------
976       ijb=ij_begin
977       ije=ij_end
978!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
979       do l=1,llm
980        teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
981     &         (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
982     &                  (knewt_g+knewt_t(l)*clat4(ijb:ije))
983       enddo ! of do l=1,llm
984!$OMP END DO
985
986!$OMP MASTER
987       if (planet_type.eq."giant") then
988         ! add an intrinsic heat flux at the base of the atmosphere
989         teta(ijb:ije,1) = teta(ijb:ije,1)
990     &        + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
991       endif
992!$OMP END MASTER
993!$OMP BARRIER
994
995       call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
996       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
997       call SendRequest(Request_Physic)
998c$OMP BARRIER
999       call WaitRequest(Request_Physic)     
1000c$OMP BARRIER
1001       call friction_p(ucov,vcov,dtvr)
1002!$OMP BARRIER
1003
1004        ! Sponge layer (if any)
1005        IF (ok_strato) THEN
1006          ! set dufi,dvfi,... to zero
1007          ijb=ij_begin
1008          ije=ij_end
1009!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1010          do l=1,llm
1011            dufi(ijb:ije,l)=0
1012            dtetafi(ijb:ije,l)=0
1013            dqfi(ijb:ije,l,1:nqtot)=0
1014          enddo
1015!$OMP END DO
1016!$OMP MASTER
1017          dpfi(ijb:ije)=0
1018!$OMP END MASTER
1019          ijb=ij_begin
1020          ije=ij_end
1021          if (pole_sud) ije=ije-iip1
1022!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1023          do l=1,llm
1024            dvfi(ijb:ije,l)=0
1025          enddo
1026!$OMP END DO
1027
1028          CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
1029          CALL addfi_p( dtvr, leapf, forward   ,
1030     $                  ucov, vcov, teta , q   ,ps ,
1031     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
1032!$OMP BARRIER
1033        ENDIF ! of IF (ok_strato)
1034      ENDIF ! of IF(iflag_phys.EQ.2)
1035
1036
1037        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
1038c$OMP BARRIER
1039        if (disvert_type==1) then
1040          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
1041        else ! we assume that we are in the disvert_type==2 case
1042          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
1043        endif
1044c$OMP BARRIER
1045
1046cc$OMP END PARALLEL
1047
1048c-----------------------------------------------------------------------
1049c   dissipation horizontale et verticale  des petites echelles:
1050c   ----------------------------------------------------------
1051
1052      IF(apdiss) THEN
1053cc$OMP  PARALLEL DEFAULT(SHARED)
1054cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1055c$OMP MASTER
1056        call suspend_timer(timer_caldyn)
1057       
1058c       print*,'Entree dans la dissipation : Iteration No ',true_itau
1059c   calcul de l'energie cinetique avant dissipation
1060c       print *,'Passage dans la dissipation'
1061
1062        call VTb(VThallo)
1063c$OMP END MASTER
1064
1065c$OMP BARRIER
1066
1067        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1068     *                          jj_Nb_dissip,1,1,Request_dissip)
1069
1070        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1071     *                          jj_Nb_dissip,1,1,Request_dissip)
1072
1073        call Register_SwapField(teta,teta,ip1jmp1,llm,
1074     *                          jj_Nb_dissip,Request_dissip)
1075
1076        call Register_SwapField(p,p,ip1jmp1,llmp1,
1077     *                          jj_Nb_dissip,Request_dissip)
1078
1079        call Register_SwapField(pk,pk,ip1jmp1,llm,
1080     *                          jj_Nb_dissip,Request_dissip)
1081
1082        call SendRequest(Request_dissip)       
1083c$OMP BARRIER
1084        call WaitRequest(Request_dissip)       
1085
1086c$OMP BARRIER
1087c$OMP MASTER
1088        call SetDistrib(jj_Nb_dissip)
1089        call VTe(VThallo)
1090        call VTb(VTdissipation)
1091        call start_timer(timer_dissip)
1092c$OMP END MASTER
1093c$OMP BARRIER
1094
1095        call covcont_p(llm,ucov,vcov,ucont,vcont)
1096        call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1097
1098c   dissipation
1099
1100!        CALL FTRACE_REGION_BEGIN("dissip")
1101        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1102!        CALL FTRACE_REGION_END("dissip")
1103         
1104        ijb=ij_begin
1105        ije=ij_end
1106c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
1107        DO l=1,llm
1108          ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1109        ENDDO
1110c$OMP END DO NOWAIT       
1111        if (pole_sud) ije=ije-iip1
1112c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
1113        DO l=1,llm
1114          vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1115        ENDDO
1116c$OMP END DO NOWAIT       
1117
1118c       teta=teta+dtetadis
1119
1120
1121c------------------------------------------------------------------------
1122        if (dissip_conservative) then
1123C       On rajoute la tendance due a la transform. Ec -> E therm. cree
1124C       lors de la dissipation
1125c$OMP BARRIER
1126c$OMP MASTER
1127            call suspend_timer(timer_dissip)
1128            call VTb(VThallo)
1129c$OMP END MASTER
1130            call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1131            call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1132            call SendRequest(Request_Dissip)
1133c$OMP BARRIER
1134            call WaitRequest(Request_Dissip)
1135c$OMP MASTER
1136            call VTe(VThallo)
1137            call resume_timer(timer_dissip)
1138c$OMP END MASTER
1139c$OMP BARRIER           
1140            call covcont_p(llm,ucov,vcov,ucont,vcont)
1141            call enercin_p(vcov,ucov,vcont,ucont,ecin)
1142           
1143            ijb=ij_begin
1144            ije=ij_end
1145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)           
1146            do l=1,llm
1147              do ij=ijb,ije
1148                dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1149                dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1150              enddo
1151            enddo
1152c$OMP END DO NOWAIT           
1153       endif ! of if (dissip_conservative)
1154
1155       ijb=ij_begin
1156       ije=ij_end
1157c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)           
1158         do l=1,llm
1159           do ij=ijb,ije
1160              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1161           enddo
1162         enddo
1163c$OMP END DO NOWAIT         
1164c------------------------------------------------------------------------
1165
1166
1167c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
1168c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
1169c
1170
1171        ijb=ij_begin
1172        ije=ij_end
1173         
1174        if (pole_nord) then
1175c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1176          DO l  =  1, llm
1177            DO ij =  1,iim
1178             tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
1179            ENDDO
1180             tpn  = SSUM(iim,tppn,1)/apoln
1181
1182            DO ij = 1, iip1
1183             teta(  ij    ,l) = tpn
1184            ENDDO
1185          ENDDO
1186c$OMP END DO NOWAIT
1187
1188c$OMP MASTER               
1189          DO ij =  1,iim
1190            tppn(ij)  = aire(  ij    ) * ps (  ij    )
1191          ENDDO
1192            tpn  = SSUM(iim,tppn,1)/apoln
1193 
1194          DO ij = 1, iip1
1195            ps(  ij    ) = tpn
1196          ENDDO
1197c$OMP END MASTER
1198        endif
1199       
1200        if (pole_sud) then
1201c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1202          DO l  =  1, llm
1203            DO ij =  1,iim
1204             tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1205            ENDDO
1206             tps  = SSUM(iim,tpps,1)/apols
1207
1208            DO ij = 1, iip1
1209             teta(ij+ip1jm,l) = tps
1210            ENDDO
1211          ENDDO
1212c$OMP END DO NOWAIT
1213
1214c$OMP MASTER               
1215          DO ij =  1,iim
1216            tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
1217          ENDDO
1218            tps  = SSUM(iim,tpps,1)/apols
1219 
1220          DO ij = 1, iip1
1221            ps(ij+ip1jm) = tps
1222          ENDDO
1223c$OMP END MASTER
1224        endif
1225
1226
1227c$OMP BARRIER
1228c$OMP MASTER
1229        call VTe(VTdissipation)
1230
1231        call stop_timer(timer_dissip)
1232       
1233        call VTb(VThallo)
1234c$OMP END MASTER
1235        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1236     *                          jj_Nb_caldyn,Request_dissip)
1237
1238        call Register_SwapField(vcov,vcov,ip1jm,llm,
1239     *                          jj_Nb_caldyn,Request_dissip)
1240
1241        call Register_SwapField(teta,teta,ip1jmp1,llm,
1242     *                          jj_Nb_caldyn,Request_dissip)
1243
1244        call Register_SwapField(p,p,ip1jmp1,llmp1,
1245     *                          jj_Nb_caldyn,Request_dissip)
1246
1247        call Register_SwapField(pk,pk,ip1jmp1,llm,
1248     *                          jj_Nb_caldyn,Request_dissip)
1249
1250        call SendRequest(Request_dissip)       
1251c$OMP BARRIER
1252        call WaitRequest(Request_dissip)       
1253
1254c$OMP BARRIER
1255c$OMP MASTER
1256        call SetDistrib(jj_Nb_caldyn)
1257        call VTe(VThallo)
1258        call resume_timer(timer_caldyn)
1259c        print *,'fin dissipation'
1260c$OMP END MASTER
1261c$OMP BARRIER
1262      END IF ! of IF(apdiss)
1263
1264cc$OMP END PARALLEL
1265
1266c ajout debug
1267c              IF( lafin ) then 
1268c                abort_message = 'Simulation finished'
1269c                call abort_gcm(modname,abort_message,0)
1270c              ENDIF
1271       
1272c   ********************************************************************
1273c   ********************************************************************
1274c   .... fin de l'integration dynamique  et physique pour le pas itau ..
1275c   ********************************************************************
1276c   ********************************************************************
1277
1278c   preparation du pas d'integration suivant  ......
1279cym      call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1280cym      call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1281c$OMP MASTER     
1282      call stop_timer(timer_caldyn)
1283c$OMP END MASTER
1284      IF (itau==itaumax) then
1285c$OMP MASTER
1286            call allgather_timer_average
1287
1288      if (mpi_rank==0) then
1289       
1290        print *,'*********************************'
1291        print *,'******    TIMER CALDYN     ******'
1292        do i=0,mpi_size-1
1293          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
1294     &            '  : temps moyen :',
1295     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1296        enddo
1297     
1298        print *,'*********************************'
1299        print *,'******    TIMER VANLEER    ******'
1300        do i=0,mpi_size-1
1301          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
1302     &            '  : temps moyen :',
1303     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1304        enddo
1305     
1306        print *,'*********************************'
1307        print *,'******    TIMER DISSIP    ******'
1308        do i=0,mpi_size-1
1309          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
1310     &            '  : temps moyen :',
1311     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
1312        enddo
1313       
1314        print *,'*********************************'
1315        print *,'******    TIMER PHYSIC    ******'
1316        do i=0,mpi_size-1
1317          print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
1318     &            '  : temps moyen :',
1319     &             timer_average(jj_nb_physic(i),timer_physic,i)
1320        enddo
1321       
1322      endif 
1323     
1324      print *,'Taille du Buffer MPI (REAL*8)',MaxBufferSize
1325      print *,'Taille du Buffer MPI utilise (REAL*8)',MaxBufferSize_Used
1326      print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1327      print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1328      CALL print_filtre_timer
1329      call fin_getparam
1330        call finalize_parallel
1331c$OMP END MASTER
1332c$OMP BARRIER
1333        RETURN
1334      ENDIF
1335     
1336      IF ( .NOT.purmats ) THEN
1337c       ........................................................
1338c       ..............  schema matsuno + leapfrog  ..............
1339c       ........................................................
1340
1341            IF(forward. OR. leapf) THEN
1342              itau= itau + 1
1343!              iday= day_ini+itau/day_step
1344!              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1345!                IF(time.GT.1.) THEN
1346!                  time = time-1.
1347!                  iday = iday+1
1348!                ENDIF
1349            ENDIF
1350
1351
1352            IF( itau. EQ. itaufinp1 ) then
1353
1354              if (flag_verif) then
1355                write(79,*) 'ucov',ucov
1356                write(80,*) 'vcov',vcov
1357                write(81,*) 'teta',teta
1358                write(82,*) 'ps',ps
1359                write(83,*) 'q',q
1360                WRITE(85,*) 'q1 = ',q(:,:,1)
1361                WRITE(86,*) 'q3 = ',q(:,:,3)
1362              endif
1363 
1364
1365c$OMP MASTER
1366              call fin_getparam
1367              call finalize_parallel
1368c$OMP END MASTER
1369              abort_message = 'Simulation finished'
1370              call abort_gcm(modname,abort_message,0)
1371              RETURN
1372            ENDIF
1373c-----------------------------------------------------------------------
1374c   ecriture du fichier histoire moyenne:
1375c   -------------------------------------
1376
1377            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1378c$OMP BARRIER
1379               IF(itau.EQ.itaufin) THEN
1380                  iav=1
1381               ELSE
1382                  iav=0
1383               ENDIF
1384#ifdef CPP_IOIPSL
1385             IF (ok_dynzon) THEN
1386             call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1387             call SendRequest(TestRequest)
1388c$OMP BARRIER
1389              call WaitRequest(TestRequest)
1390c$OMP BARRIER
1391c$OMP MASTER
1392!              CALL writedynav_p(histaveid, itau,vcov ,
1393!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1394
1395c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1396              CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1397     ,             ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1398c$OMP END MASTER
1399              ENDIF !ok_dynzon
1400#endif
1401               IF (ok_dyn_ave) THEN
1402!$OMP MASTER
1403#ifdef CPP_IOIPSL
1404! Ehouarn: Gather fields and make master send to output
1405                call Gather_Field(vcov,ip1jm,llm,0)
1406                call Gather_Field(ucov,ip1jmp1,llm,0)
1407                call Gather_Field(teta,ip1jmp1,llm,0)
1408                call Gather_Field(pk,ip1jmp1,llm,0)
1409                call Gather_Field(phi,ip1jmp1,llm,0)
1410                do iq=1,nqtot
1411                  call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1412                enddo
1413                call Gather_Field(masse,ip1jmp1,llm,0)
1414                call Gather_Field(ps,ip1jmp1,1,0)
1415                call Gather_Field(phis,ip1jmp1,1,0)
1416                if (mpi_rank==0) then
1417                 CALL writedynav(itau,vcov,
1418     &                 ucov,teta,pk,phi,q,masse,ps,phis)
1419                endif
1420#endif
1421!$OMP END MASTER
1422               ENDIF ! of IF (ok_dyn_ave)
1423            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1424
1425c-----------------------------------------------------------------------
1426c   ecriture de la bande histoire:
1427c   ------------------------------
1428
1429            IF( MOD(itau,iecri).EQ.0) THEN
1430             ! Ehouarn: output only during LF or Backward Matsuno
1431             if (leapf.or.(.not.leapf.and.(.not.forward))) then
1432c$OMP BARRIER
1433c$OMP MASTER
1434              nbetat = nbetatdem
1435              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1436       
1437cym        unat=0.
1438       
1439              ijb=ij_begin
1440              ije=ij_end
1441       
1442              if (pole_nord) then
1443                ijb=ij_begin+iip1
1444                unat(1:iip1,:)=0.
1445              endif
1446       
1447              if (pole_sud) then
1448                ije=ij_end-iip1
1449                unat(ij_end-iip1+1:ij_end,:)=0.
1450              endif
1451           
1452              do l=1,llm
1453                unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1454              enddo
1455
1456              ijb=ij_begin
1457              ije=ij_end
1458              if (pole_sud) ije=ij_end-iip1
1459       
1460              do l=1,llm
1461                vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1462              enddo
1463       
1464#ifdef CPP_IOIPSL
1465              if (ok_dyn_ins) then
1466! Ehouarn: Gather fields and make master write to output
1467                call Gather_Field(vcov,ip1jm,llm,0)
1468                call Gather_Field(ucov,ip1jmp1,llm,0)
1469                call Gather_Field(teta,ip1jmp1,llm,0)
1470                call Gather_Field(phi,ip1jmp1,llm,0)
1471                do iq=1,nqtot
1472                  call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1473                enddo
1474                call Gather_Field(masse,ip1jmp1,llm,0)
1475                call Gather_Field(ps,ip1jmp1,1,0)
1476                call Gather_Field(phis,ip1jmp1,1,0)
1477                if (mpi_rank==0) then
1478                 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1479                endif
1480!              CALL writehist_p(histid,histvid, itau,vcov,
1481!     &                         ucov,teta,phi,q,masse,ps,phis)
1482! or use writefield_p
1483!      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1484!      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1485!      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1486!      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1487              endif ! of if (ok_dyn_ins)
1488#endif
1489! For some Grads outputs of fields
1490              if (output_grads_dyn) then
1491! Ehouarn: hope this works the way I think it does:
1492                  call Gather_Field(unat,ip1jmp1,llm,0)
1493                  call Gather_Field(vnat,ip1jm,llm,0)
1494                  call Gather_Field(teta,ip1jmp1,llm,0)
1495                  call Gather_Field(ps,ip1jmp1,1,0)
1496                  do iq=1,nqtot
1497                    call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1498                  enddo
1499                  if (mpi_rank==0) then
1500#include "write_grads_dyn.h"
1501                  endif
1502              endif ! of if (output_grads_dyn)
1503c$OMP END MASTER
1504             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1505            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1506
1507            IF(itau.EQ.itaufin) THEN
1508
1509c$OMP BARRIER
1510c$OMP MASTER
1511
1512!              if (planet_type.eq."earth") then
1513! Write an Earth-format restart file
1514                CALL dynredem1_p("restart.nc",0.0,
1515     &                           vcov,ucov,teta,q,masse,ps)
1516!              endif ! of if (planet_type.eq."earth")
1517
1518!              CLOSE(99)
1519c$OMP END MASTER
1520            ENDIF ! of IF (itau.EQ.itaufin)
1521
1522c-----------------------------------------------------------------------
1523c   gestion de l'integration temporelle:
1524c   ------------------------------------
1525
1526            IF( MOD(itau,iperiod).EQ.0 )    THEN
1527                    GO TO 1
1528            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1529
1530                   IF( forward )  THEN
1531c      fin du pas forward et debut du pas backward
1532
1533                      forward = .FALSE.
1534                        leapf = .FALSE.
1535                           GO TO 2
1536
1537                   ELSE
1538c      fin du pas backward et debut du premier pas leapfrog
1539
1540                        leapf =  .TRUE.
1541                        dt  =  2.*dtvr
1542                        GO TO 2
1543                   END IF
1544            ELSE
1545
1546c      ......   pas leapfrog  .....
1547
1548                 leapf = .TRUE.
1549                 dt  = 2.*dtvr
1550                 GO TO 2
1551            END IF ! of IF (MOD(itau,iperiod).EQ.0)
1552                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
1553
1554
1555      ELSE ! of IF (.not.purmats)
1556
1557c       ........................................................
1558c       ..............       schema  matsuno        ...............
1559c       ........................................................
1560            IF( forward )  THEN
1561
1562             itau =  itau + 1
1563!             iday = day_ini+itau/day_step
1564!             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1565!
1566!                  IF(time.GT.1.) THEN
1567!                   time = time-1.
1568!                   iday = iday+1
1569!                  ENDIF
1570
1571               forward =  .FALSE.
1572               IF( itau. EQ. itaufinp1 ) then 
1573c$OMP MASTER
1574                 call fin_getparam
1575                 call finalize_parallel
1576c$OMP END MASTER
1577                 abort_message = 'Simulation finished'
1578                 call abort_gcm(modname,abort_message,0)
1579                 RETURN
1580               ENDIF
1581               GO TO 2
1582
1583            ELSE ! of IF(forward) i.e. backward step
1584
1585              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1586               IF(itau.EQ.itaufin) THEN
1587                  iav=1
1588               ELSE
1589                  iav=0
1590               ENDIF
1591#ifdef CPP_IOIPSL
1592               IF (ok_dynzon) THEN
1593c$OMP BARRIER
1594               call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1595               call SendRequest(TestRequest)
1596c$OMP BARRIER
1597               call WaitRequest(TestRequest)
1598c$OMP BARRIER
1599c$OMP MASTER
1600!               CALL writedynav_p(histaveid, itau,vcov ,
1601!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1602               CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1603     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1604c$OMP END MASTER
1605               END IF !ok_dynzon
1606#endif
1607               IF (ok_dyn_ave) THEN
1608!$OMP MASTER
1609#ifdef CPP_IOIPSL
1610! Ehouarn: Gather fields and make master send to output
1611                call Gather_Field(vcov,ip1jm,llm,0)
1612                call Gather_Field(ucov,ip1jmp1,llm,0)
1613                call Gather_Field(teta,ip1jmp1,llm,0)
1614                call Gather_Field(pk,ip1jmp1,llm,0)
1615                call Gather_Field(phi,ip1jmp1,llm,0)
1616                do iq=1,nqtot
1617                  call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1618                enddo
1619                call Gather_Field(masse,ip1jmp1,llm,0)
1620                call Gather_Field(ps,ip1jmp1,1,0)
1621                call Gather_Field(phis,ip1jmp1,1,0)
1622                if (mpi_rank==0) then
1623                 CALL writedynav(itau,vcov,
1624     &                 ucov,teta,pk,phi,q,masse,ps,phis)
1625                endif
1626#endif
1627!$OMP END MASTER
1628               ENDIF ! of IF (ok_dyn_ave)
1629
1630              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1631
1632
1633               IF(MOD(itau,iecri         ).EQ.0) THEN
1634c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
1635c$OMP BARRIER
1636c$OMP MASTER
1637                nbetat = nbetatdem
1638                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1639
1640cym        unat=0.
1641                ijb=ij_begin
1642                ije=ij_end
1643       
1644                if (pole_nord) then
1645                  ijb=ij_begin+iip1
1646                  unat(1:iip1,:)=0.
1647                endif
1648       
1649                if (pole_sud) then
1650                  ije=ij_end-iip1
1651                  unat(ij_end-iip1+1:ij_end,:)=0.
1652                endif
1653           
1654                do l=1,llm
1655                  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1656                enddo
1657
1658                ijb=ij_begin
1659                ije=ij_end
1660                if (pole_sud) ije=ij_end-iip1
1661       
1662                do l=1,llm
1663                  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1664                enddo
1665
1666#ifdef CPP_IOIPSL
1667              if (ok_dyn_ins) then
1668! Ehouarn: Gather fields and make master send to output
1669                call Gather_Field(vcov,ip1jm,llm,0)
1670                call Gather_Field(ucov,ip1jmp1,llm,0)
1671                call Gather_Field(teta,ip1jmp1,llm,0)
1672                call Gather_Field(phi,ip1jmp1,llm,0)
1673                do iq=1,nqtot
1674                  call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1675                enddo
1676                call Gather_Field(masse,ip1jmp1,llm,0)
1677                call Gather_Field(ps,ip1jmp1,1,0)
1678                call Gather_Field(phis,ip1jmp1,1,0)
1679                if (mpi_rank==0) then
1680                 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1681                endif
1682!                CALL writehist_p(histid, histvid, itau,vcov ,
1683!     &                           ucov,teta,phi,q,masse,ps,phis)
1684              endif ! of if (ok_dyn_ins)
1685#endif
1686! For some Grads output (but does it work?)
1687                if (output_grads_dyn) then
1688                  call Gather_Field(unat,ip1jmp1,llm,0)
1689                  call Gather_Field(vnat,ip1jm,llm,0)
1690                  call Gather_Field(teta,ip1jmp1,llm,0)
1691                  call Gather_Field(ps,ip1jmp1,1,0)
1692                  do iq=1,nqtot
1693                    call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1694                  enddo
1695c     
1696                  if (mpi_rank==0) then
1697#include "write_grads_dyn.h"
1698                  endif
1699                endif ! of if (output_grads_dyn)
1700
1701c$OMP END MASTER
1702              ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1703
1704              IF(itau.EQ.itaufin) THEN
1705!                if (planet_type.eq."earth") then
1706c$OMP MASTER
1707                   CALL dynredem1_p("restart.nc",0.0,
1708     .                               vcov,ucov,teta,q,masse,ps)
1709c$OMP END MASTER
1710!                endif ! of if (planet_type.eq."earth")
1711              ENDIF ! of IF(itau.EQ.itaufin)
1712
1713              forward = .TRUE.
1714              GO TO  1
1715
1716            ENDIF ! of IF (forward)
1717
1718      END IF ! of IF(.not.purmats)
1719c$OMP MASTER
1720      call fin_getparam
1721      call finalize_parallel
1722c$OMP END MASTER
1723      RETURN
1724      END
Note: See TracBrowser for help on using the repository browser.