source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/leapfrog_p.F @ 3524

Last change on this file since 3524 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

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