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

Last change on this file since 1523 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
RevLine 
[764]1!
[1279]2! $Id: leapfrog_p.F 1446 2010-10-22 09:27:25Z idelkadi $
[630]3!
4c
5c
6
[1146]7      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
[630]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
[985]18       USE timer_filtre, ONLY : print_filtre_timer
[1146]19       USE infotrac
[1279]20       USE guide_p_mod, ONLY : guide_main
21       USE getparam
[1403]22       USE control_mod
[630]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"
[1403]68!#include "com_io_dyn.h"
[630]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
[774]80      REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
81      REAL :: teta(ip1jmp1,llm)                 ! temperature potentielle
[1146]82      REAL :: q(ip1jmp1,llm,nqtot)              ! champs advectes
[774]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
[630]92
93c variables dynamiques intermediaire pour le transport
[774]94      REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
[630]95
96c   variables dynamiques au pas -1
[774]97      REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
98      REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
99      REAL,SAVE :: massem1(ip1jmp1,llm)
[630]100
101c   tendances dynamiques
[774]102      REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
[1146]103      REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
104      REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
[630]105
106c   tendances de la dissipation
[774]107      REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
108      REAL,SAVE :: dtetadis(ip1jmp1,llm)
[630]109
110c   tendances physiques
[774]111      REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
112      REAL,SAVE :: dtetafi(ip1jmp1,llm)
[1146]113      REAL,SAVE :: dpfi(ip1jmp1)
114      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
[630]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
[1279]122!      INTEGER  iday ! jour julien
123      REAL       time
[630]124
125      REAL  SSUM
[774]126      REAL time_0
127      REAL,SAVE :: finvmaold(ip1jmp1,llm)
[630]128
129cym      LOGICAL  lafin
[774]130      LOGICAL :: lafin
[630]131      INTEGER ij,iq,l
132      INTEGER ik
133
134      real time_step, t_wrt, t_ops
135
[1279]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
[630]142      LOGICAL first,callinigrads
143
144      data callinigrads/.true./
145      character*10 string10
146
[774]147      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[960]148      REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
[630]149
150c+jld variables test conservation energie
[774]151      REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
[630]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
[774]155      REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
156      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
157      REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
[630]158      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
159      CHARACTER*15 ztit
[764]160!      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
161!      SAVE      ip_ebil_dyn
162!      DATA      ip_ebil_dyn/0/
[630]163c-jld
164
165      character*80 dynhist_file, dynhistave_file
166      character*20 modname
167      character*80 abort_message
168
169
[774]170      logical,PARAMETER :: dissip_conservative=.TRUE.
171 
[630]172      INTEGER testita
173      PARAMETER (testita = 9)
[1286]174
175      logical , parameter :: flag_verif = .false.
[630]176     
[764]177c declaration liees au parallelisme
[630]178      INTEGER :: ierr
[774]179      LOGICAL :: FirstCaldyn
180      LOGICAL :: FirstPhysic
[630]181      INTEGER :: ijb,ije,j,i
182      type(Request) :: TestRequest
183      type(Request) :: Request_Dissip
184      type(Request) :: Request_physic
[774]185      REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
[1146]186      REAL,SAVE :: dtetafi_tmp(iip1,llm)
187      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
[774]188      REAL,SAVE :: dpfi_tmp(iip1)
189
190      INTEGER :: true_itau
[630]191      LOGICAL :: verbose=.true.
[774]192      INTEGER :: iapptrac
193      INTEGER :: AdjustCount
[1279]194!      INTEGER :: var_time
[774]195      LOGICAL :: ok_start_timer=.FALSE.
[1146]196      LOGICAL, SAVE :: firstcall=.TRUE.
[774]197
198c$OMP MASTER
[630]199      ItCount=0
[774]200c$OMP END MASTER     
201      true_itau=0
202      FirstCaldyn=.TRUE.
203      FirstPhysic=.TRUE.
204      iapptrac=0
205      AdjustCount = 0
206      lafin=.false.
[630]207     
208      itaufin   = nday*day_step
209      itaufinp1 = itaufin +1
[1146]210      modname="leapfrog_p"
[630]211
212      itau = 0
[1279]213!      iday = day_ini+itau/day_step
[1403]214!      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]215!         IF(time.GT.1.) THEN
216!          time = time-1.
217!          iday = iday+1
218!         ENDIF
[630]219
[1146]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
[630]231c-----------------------------------------------------------------------
232c   On initialise la pression et la fonction d'Exner :
233c   --------------------------------------------------
234
[774]235c$OMP MASTER
[1446]236      dq(:,:,:)=0.
[630]237      CALL pression ( ip1jmp1, ap, bp, ps, p       )
238      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[774]239c$OMP END MASTER
[630]240c-----------------------------------------------------------------------
241c   Debut de l'integration temporelle:
242c   ----------------------------------
[764]243c et du parallelisme !!
[630]244
245   1  CONTINUE
246
[1279]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))
[985]250
[630]251
252#ifdef CPP_IOIPSL
[1279]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
[630]258      endif
259#endif
[1279]260
[630]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
[774]274c$OMP MASTER
[630]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 )
[774]283c$OMP END MASTER
284c$OMP BARRIER
[630]285       else
[1146]286! Save fields obtained at previous time step as '...m1'
[630]287         ijb=ij_begin
288         ije=ij_end
[774]289
290c$OMP MASTER           
[630]291         psm1     (ijb:ije) = ps    (ijb:ije)
[774]292c$OMP END MASTER
293
294c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
[949]295         DO l=1,llm     
[774]296           ije=ij_end
[949]297           ucovm1   (ijb:ije,l) = ucov  (ijb:ije,l)
[774]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)
[949]301                 
[774]302           if (pole_sud) ije=ij_end-iip1
303           vcovm1(ijb:ije,l) = vcov  (ijb:ije,l)
[630]304       
[774]305
306         ENDDO
[985]307c$OMP ENDDO 
[774]308
309
310          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
[630]311     .                    llm, -2,2, .TRUE., 1 )
312
[1146]313       endif ! of if (FirstCaldyn)
[630]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
[764]324cym  ne sert a rien
[630]325cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
326
327   2  CONTINUE
328
[774]329c$OMP MASTER
[630]330      ItCount=ItCount+1
[995]331      if (MOD(ItCount,1)==1) then
[630]332        debug=.true.
333      else
334        debug=.false.
335      endif
[774]336c$OMP END MASTER
[630]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.
[764]352c      idissip=1
[630]353      IF( purmats ) THEN
[1403]354      ! Purely Matsuno time stepping
[630]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
[1146]358     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
[630]359      ELSE
[1403]360      ! Leapfrog/Matsuno time stepping
[630]361         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
362         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
[1146]363         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
[630]364      END IF
365
[1403]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
[630]372cym    ---> Pour le moment     
373cym      apphys = .FALSE.
374      statcl = .FALSE.
[1403]375      conser = .FALSE. ! ie: no output of control variables to stdout in //
[630]376     
377      if (firstCaldyn) then
[774]378c$OMP MASTER
[630]379          call SetDistrib(jj_Nb_Caldyn)
[774]380c$OMP END MASTER
381c$OMP BARRIER
[630]382          firstCaldyn=.FALSE.
383cym          call InitTime
[774]384c$OMP MASTER
[630]385          call Init_timer
[774]386c$OMP END MASTER
[630]387      endif
[774]388
389c$OMP MASTER     
390      IF (ok_start_timer) THEN
391        CALL InitTime
[949]392        ok_start_timer=.FALSE.
[774]393      ENDIF     
394c$OMP END MASTER     
395     
[630]396      if (Adjust) then
[774]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
[630]401           AdjustCount=0
402           call allgather_timer_average
403
404        if (Verbose) then
405       
[949]406        print *,'*********************************'
407        print *,'******    TIMER CALDYN     ******'
408        do i=0,mpi_size-1
409          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
[630]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 *,'*********************************'
[949]416        print *,'******    TIMER VANLEER    ******'
417        do i=0,mpi_size-1
418          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
[630]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 *,'*********************************'
[949]425        print *,'******    TIMER DISSIP    ******'
426        do i=0,mpi_size-1
427          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
[630]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       
[949]433        if (mpi_rank==0) call WriteBands
434       
[630]435       endif
436       
[949]437         call AdjustBands_caldyn
[630]438         if (mpi_rank==0) call WriteBands
[949]439         
440         call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
[630]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)
[774]468         call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
[630]469     &                                jj_Nb_caldyn,0,0,TestRequest)
470         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
471     &                                jj_Nb_caldyn,0,0,TestRequest)
472 
[1146]473        do j=1,nqtot
[764]474         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
475     &                                jj_nb_caldyn,0,0,TestRequest)
476        enddo
477
[630]478         call SetDistrib(jj_nb_caldyn)
479         call SendRequest(TestRequest)
480         call WaitRequest(TestRequest)
481         
[949]482        call AdjustBands_dissip
483        call AdjustBands_physic
[774]484
[630]485      endif
[774]486c$OMP END MASTER 
[630]487      endif       
[774]488     
[630]489     
490     
491c-----------------------------------------------------------------------
492c   calcul des tendances dynamiques:
493c   --------------------------------
[774]494c$OMP BARRIER
495c$OMP MASTER
[630]496       call VTb(VThallo)
[985]497c$OMP END MASTER
498
[630]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       
[1146]508c       do j=1,nqtot
[630]509c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
510c     *                       TestRequest)
511c        enddo
512
513       call SendRequest(TestRequest)
[985]514c$OMP BARRIER
[630]515       call WaitRequest(TestRequest)
[985]516
517c$OMP MASTER
[630]518       call VTe(VThallo)
[985]519c$OMP END MASTER
520c$OMP BARRIER
[764]521     
[949]522      if (debug) then       
[985]523!$OMP BARRIER
524!$OMP MASTER
[630]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/)))
[949]530        call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
[630]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/)))
[1146]534        do j=1,nqtot
[764]535          call WriteField_p('q'//trim(int2str(j)),
536     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
[985]537        enddo
538!$OMP END MASTER       
539c$OMP BARRIER
[630]540      endif
[985]541
[630]542     
543      True_itau=True_itau+1
[774]544
545c$OMP MASTER
[1146]546      IF (prt_level>9) THEN
547        WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
548      ENDIF
[630]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)
[774]557c$OMP END MASTER
[1279]558!      var_time=time+iday-day_ini
[764]559
[985]560c$OMP BARRIER
561!      CALL FTRACE_REGION_BEGIN("caldyn")
[1279]562      time = jD_cur + jH_cur
[630]563      CALL caldyn_p
564     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]565     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[764]566
[985]567!      CALL FTRACE_REGION_END("caldyn")
[1279]568
[774]569c$OMP MASTER
[630]570      call VTe(VTcaldyn)
[774]571c$OMP END MASTER     
[985]572
573cc$OMP BARRIER
574cc$OMP MASTER
[1279]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/)))
[985]585cc$OMP END MASTER
[630]586
587c-----------------------------------------------------------------------
588c   calcul des tendances advection des traceurs (dont l'humidite)
589c   -------------------------------------------------------------
590
591      IF( forward. OR . leapf )  THEN
[774]592cc$OMP PARALLEL DEFAULT(SHARED)
[630]593c
[960]594         CALL caladvtrac_p(q,pbaru,pbarv,
595     *        p, masse, dq,  teta,
596     .        flxw,pk, iapptrac)
[985]597
[1438]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
[985]603
[1146]604      ENDIF ! of IF( forward. OR . leapf )
[774]605cc$OMP END PARALLEL
[630]606
607c-----------------------------------------------------------------------
608c   integrations dynamique et traceurs:
609c   ----------------------------------
[774]610
611c$OMP MASTER
[630]612       call VTb(VTintegre)
[774]613c$OMP END MASTER
[764]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/)))
[774]622cc$OMP PARALLEL DEFAULT(SHARED)
[985]623c$OMP BARRIER
624!       CALL FTRACE_REGION_BEGIN("integrd")
[1146]625
[630]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
[985]630!       CALL FTRACE_REGION_END("integrd")
631c$OMP BARRIER
632cc$OMP MASTER
[764]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/)))
[985]641c
642c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
[1146]643c      do j=1,nqtot
[985]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
[764]650
[985]651
[774]652c$OMP MASTER
[630]653       call VTe(VTintegre)
[774]654c$OMP END MASTER
[630]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
[774]667
668cc$OMP END PARALLEL
669
[630]670c
671c
672       IF( apphys )  THEN
673c
674c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
675c
[774]676cc$OMP PARALLEL DEFAULT(SHARED)
677cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
[764]678
679c$OMP MASTER
[630]680         call suspend_timer(timer_caldyn)
[1146]681
[1403]682        if (prt_level >= 10) then
[1146]683         write(lunout,*)
684     &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
[1403]685        endif
[764]686c$OMP END MASTER
687
[630]688         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )
[764]689
[985]690c$OMP BARRIER
[630]691         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
[985]692c$OMP BARRIER
[1279]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)
[630]698
699c rajout debug
700c       lafin = .true.
701
702
703c   Inbterface avec les routines de phylmd (phymars ... )
704c   -----------------------------------------------------
705
706c+jld
707
[764]708c  Diagnostique de conservation de l'energie : initialisation
[630]709      IF (ip_ebil_dyn.ge.1 ) THEN
710          ztit='bil dyn'
[1146]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
[630]716      ENDIF
717c-jld
[764]718c$OMP BARRIER
719c$OMP MASTER
[630]720        call VTb(VThallo)
[985]721c$OMP END MASTER
722
[630]723        call SetTag(Request_physic,800)
[949]724       
[630]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       
[1279]734        call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
735     *                               jj_Nb_physic,1,2,Request_physic)
736
[630]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       
[949]752c        call SetDistrib(jj_nb_vanleer)
[1146]753        do j=1,nqtot
[630]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
[960]758
[630]759        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
760     *                               jj_Nb_physic,2,2,Request_physic)
[949]761       
[630]762        call SendRequest(Request_Physic)
[985]763c$OMP BARRIER
[630]764        call WaitRequest(Request_Physic)       
[985]765
766c$OMP BARRIER
767c$OMP MASTER
768        call SetDistrib(jj_nb_Physic)
[949]769        call VTe(VThallo)
770       
771        call VTb(VTphysiq)
[764]772c$OMP END MASTER
773c$OMP BARRIER
774
[949]775cc$OMP MASTER       
[764]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
[985]783!        CALL FTRACE_REGION_BEGIN("calfis")
[1279]784        CALL calfis_p(lafin ,jD_cur, jH_cur,
[630]785     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]786     $               du,dv,dteta,dq,
[630]787     $               flxw,
788     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
[985]789!        CALL FTRACE_REGION_END("calfis")
[630]790        ijb=ij_begin
791        ije=ij_end 
792        if ( .not. pole_nord) then
[764]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
[630]803          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim) 
[764]804c$OMP END MASTER
[1146]805        endif ! of if ( .not. pole_nord)
[764]806
807c$OMP BARRIER
808c$OMP MASTER
[949]809        call SetDistrib(jj_nb_Physic_bis)
[630]810
811        call VTb(VThallo)
[985]812c$OMP END MASTER
813c$OMP BARRIER
[630]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
[1146]827        do j=1,nqtot
[630]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)
[985]833c$OMP BARRIER
[630]834        call WaitRequest(Request_Physic)
835             
[985]836c$OMP BARRIER
837c$OMP MASTER
[630]838        call VTe(VThallo)
839 
840        call SetDistrib(jj_nb_Physic)
[764]841c$OMP END MASTER
[949]842c$OMP BARRIER       
843                ijb=ij_begin
844        if (.not. pole_nord) then
845       
[764]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
[630]858          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
[764]859c$OMP END MASTER
[949]860         
[1146]861        endif ! of if (.not. pole_nord)
[764]862c$OMP BARRIER
[949]863cc$OMP MASTER       
[630]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/)))
[764]868cc$OMP END MASTER
[630]869c     
[1146]870c      do j=1,nqtot
[630]871c        call WriteField_p('dqfi'//trim(int2str(j)),
872c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
873c      enddo
874
875c      ajout des tendances physiques:
876c      ------------------------------
[1000]877         IF (ok_strato) THEN
[1279]878           CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1000]879         ENDIF
880       
[1146]881          CALL addfi_p( dtphys, leapf, forward   ,
[630]882     $                  ucov, vcov, teta , q   ,ps ,
883     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
884
[764]885c$OMP BARRIER
886c$OMP MASTER
[630]887        call VTe(VTphysiq)
888
889        call VTb(VThallo)
[985]890c$OMP END MASTER
[630]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       
[1279]902        call Register_SwapField(masse,masse,ip1jmp1,llm,
903     *                               jj_Nb_caldyn,Request_physic)
904
[630]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
[1146]920        do j=1,nqtot
[630]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)
[985]928c$OMP BARRIER
[630]929        call WaitRequest(Request_Physic)     
930
[985]931c$OMP BARRIER
932c$OMP MASTER
933       call VTe(VThallo)
934       call SetDistrib(jj_Nb_caldyn)
[764]935c$OMP END MASTER
936c$OMP BARRIER
[630]937c
[764]938c  Diagnostique de conservation de l'energie : difference
[630]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
[764]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
[630]953
[1146]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
[1437]965!   Academic case : Simple friction and Newtonan relaxation
966!   -------------------------------------------------------
[630]967       ijb=ij_begin
968       ije=ij_end
[1403]969!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
970       do l=1,llm
[1437]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
[1403]975!$OMP END DO
[630]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)
[1279]980c$OMP BARRIER
[630]981       call WaitRequest(Request_Physic)     
[1403]982c$OMP BARRIER
[1437]983       call friction_p(ucov,vcov,dtvr)
[1403]984!$OMP BARRIER
[1437]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)
[1146]1016      ENDIF ! of IF(iflag_phys.EQ.2)
[630]1017
1018
[985]1019        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
[774]1020c$OMP BARRIER
[630]1021        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[774]1022c$OMP BARRIER
[630]1023
[774]1024cc$OMP END PARALLEL
[630]1025
1026c-----------------------------------------------------------------------
1027c   dissipation horizontale et verticale  des petites echelles:
1028c   ----------------------------------------------------------
1029
1030      IF(apdiss) THEN
[774]1031cc$OMP  PARALLEL DEFAULT(SHARED)
1032cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
[764]1033c$OMP MASTER
[630]1034        call suspend_timer(timer_caldyn)
[949]1035       
1036c       print*,'Entree dans la dissipation : Iteration No ',true_itau
[630]1037c   calcul de l'energie cinetique avant dissipation
[949]1038c       print *,'Passage dans la dissipation'
[630]1039
1040        call VTb(VThallo)
[764]1041c$OMP END MASTER
[630]1042
[764]1043c$OMP BARRIER
[985]1044
[630]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)       
[985]1061c$OMP BARRIER
[630]1062        call WaitRequest(Request_dissip)       
[985]1063
1064c$OMP BARRIER
1065c$OMP MASTER
[630]1066        call SetDistrib(jj_Nb_dissip)
[949]1067        call VTe(VThallo)
[630]1068        call VTb(VTdissipation)
[949]1069        call start_timer(timer_dissip)
[764]1070c$OMP END MASTER
1071c$OMP BARRIER
1072
[630]1073        call covcont_p(llm,ucov,vcov,ucont,vcont)
1074        call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1075
1076c   dissipation
1077
[985]1078!        CALL FTRACE_REGION_BEGIN("dissip")
[630]1079        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
[985]1080!        CALL FTRACE_REGION_END("dissip")
[764]1081         
[949]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)
[764]1087        ENDDO
[949]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)
[764]1093        ENDDO
[949]1094c$OMP END DO NOWAIT       
[764]1095
[630]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
[764]1103c$OMP BARRIER
1104c$OMP MASTER
[630]1105            call suspend_timer(timer_dissip)
[949]1106            call VTb(VThallo)
[985]1107c$OMP END MASTER
[630]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)
[985]1111c$OMP BARRIER
[630]1112            call WaitRequest(Request_Dissip)
[985]1113c$OMP MASTER
[630]1114            call VTe(VThallo)
1115            call resume_timer(timer_dissip)
[764]1116c$OMP END MASTER
[949]1117c$OMP BARRIER           
[630]1118            call covcont_p(llm,ucov,vcov,ucont,vcont)
1119            call enercin_p(vcov,ucov,vcont,ucont,ecin)
1120           
[949]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)
[630]1128              enddo
[949]1129            enddo
1130c$OMP END DO NOWAIT           
[1403]1131       endif ! of if (dissip_conservative)
[630]1132
1133       ijb=ij_begin
1134       ije=ij_end
[949]1135c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)           
1136         do l=1,llm
1137           do ij=ijb,ije
[630]1138              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
[949]1139           enddo
1140         enddo
1141c$OMP END DO NOWAIT         
[630]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
[764]1153c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]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
[764]1164c$OMP END DO NOWAIT
1165
1166c$OMP MASTER               
[630]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
[764]1175c$OMP END MASTER
[630]1176        endif
1177       
1178        if (pole_sud) then
[764]1179c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]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
[764]1190c$OMP END DO NOWAIT
1191
1192c$OMP MASTER               
[630]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
[764]1201c$OMP END MASTER
[630]1202        endif
1203
[764]1204
1205c$OMP BARRIER
1206c$OMP MASTER
[630]1207        call VTe(VTdissipation)
1208
1209        call stop_timer(timer_dissip)
[949]1210       
[630]1211        call VTb(VThallo)
[985]1212c$OMP END MASTER
[630]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)       
[985]1229c$OMP BARRIER
[630]1230        call WaitRequest(Request_dissip)       
[985]1231
1232c$OMP BARRIER
1233c$OMP MASTER
[630]1234        call SetDistrib(jj_Nb_caldyn)
1235        call VTe(VThallo)
1236        call resume_timer(timer_caldyn)
[961]1237c        print *,'fin dissipation'
[764]1238c$OMP END MASTER
[985]1239c$OMP BARRIER
[1403]1240      END IF ! of IF(apdiss)
[630]1241
[774]1242cc$OMP END PARALLEL
1243
[630]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/)))
[774]1259c$OMP MASTER     
[630]1260      call stop_timer(timer_caldyn)
[774]1261c$OMP END MASTER
[630]1262      IF (itau==itaumax) then
[774]1263c$OMP MASTER
[630]1264            call allgather_timer_average
1265
1266      if (mpi_rank==0) then
1267       
1268        print *,'*********************************'
[949]1269        print *,'******    TIMER CALDYN     ******'
1270        do i=0,mpi_size-1
1271          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
[630]1272     &            '  : temps moyen :',
1273     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1274        enddo
1275     
1276        print *,'*********************************'
[949]1277        print *,'******    TIMER VANLEER    ******'
1278        do i=0,mpi_size-1
1279          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
[630]1280     &            '  : temps moyen :',
1281     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1282        enddo
1283     
1284        print *,'*********************************'
[949]1285        print *,'******    TIMER DISSIP    ******'
1286        do i=0,mpi_size-1
1287          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
[630]1288     &            '  : temps moyen :',
1289     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
1290        enddo
[949]1291       
1292        print *,'*********************************'
1293        print *,'******    TIMER PHYSIC    ******'
1294        do i=0,mpi_size-1
1295          print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
[630]1296     &            '  : temps moyen :',
1297     &             timer_average(jj_nb_physic(i),timer_physic,i)
1298        enddo
[949]1299       
[630]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()
[985]1306      CALL print_filtre_timer
[1279]1307      call fin_getparam
[630]1308        call finalize_parallel
[774]1309c$OMP END MASTER
[985]1310c$OMP BARRIER
[949]1311        RETURN
[630]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
[1279]1321!              iday= day_ini+itau/day_step
[1403]1322!              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]1323!                IF(time.GT.1.) THEN
1324!                  time = time-1.
1325!                  iday = iday+1
1326!                ENDIF
[630]1327            ENDIF
1328
1329
[1286]1330            IF( itau. EQ. itaufinp1 ) then
[630]1331
[1286]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
[774]1343c$OMP MASTER
[1279]1344              call fin_getparam
[764]1345              call finalize_parallel
[774]1346c$OMP END MASTER
[630]1347              abort_message = 'Simulation finished'
1348              call abort_gcm(modname,abort_message,0)
[985]1349              RETURN
[630]1350            ENDIF
1351c-----------------------------------------------------------------------
1352c   ecriture du fichier histoire moyenne:
1353c   -------------------------------------
1354
1355            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[774]1356c$OMP BARRIER
[630]1357               IF(itau.EQ.itaufin) THEN
1358                  iav=1
1359               ELSE
1360                  iav=0
1361               ENDIF
1362#ifdef CPP_IOIPSL
[1146]1363             IF (ok_dynzon) THEN
[774]1364             call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
[985]1365             call SendRequest(TestRequest)
1366c$OMP BARRIER
[630]1367              call WaitRequest(TestRequest)
[1000]1368c$OMP BARRIER
[985]1369c$OMP MASTER
[1279]1370!              CALL writedynav_p(histaveid, itau,vcov ,
1371!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
[1146]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)
[985]1376c$OMP END MASTER
[1146]1377              ENDIF !ok_dynzon
[630]1378#endif
[1403]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))
[630]1402
1403c-----------------------------------------------------------------------
1404c   ecriture de la bande histoire:
1405c   ------------------------------
1406
[1403]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
[774]1410c$OMP BARRIER
1411c$OMP MASTER
[1146]1412              nbetat = nbetatdem
1413              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[774]1414       
[949]1415cym        unat=0.
1416       
[1146]1417              ijb=ij_begin
1418              ije=ij_end
[949]1419       
[1146]1420              if (pole_nord) then
1421                ijb=ij_begin+iip1
1422                unat(1:iip1,:)=0.
1423              endif
[949]1424       
[1146]1425              if (pole_sud) then
1426                ije=ij_end-iip1
1427                unat(ij_end-iip1+1:ij_end,:)=0.
1428              endif
[949]1429           
[1146]1430              do l=1,llm
1431                unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1432              enddo
[630]1433
[1146]1434              ijb=ij_begin
1435              ije=ij_end
1436              if (pole_sud) ije=ij_end-iip1
[949]1437       
[1146]1438              do l=1,llm
1439                vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1440              enddo
[949]1441       
[630]1442#ifdef CPP_IOIPSL
[1403]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
[1279]1458!              CALL writehist_p(histid,histvid, itau,vcov,
1459!     &                         ucov,teta,phi,q,masse,ps,phis)
[1403]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)
[630]1466#endif
[1146]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)
[774]1481c$OMP END MASTER
[1403]1482             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]1483            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]1484
1485            IF(itau.EQ.itaufin) THEN
1486
[774]1487c$OMP BARRIER
1488c$OMP MASTER
[630]1489
[1446]1490!              if (planet_type.eq."earth") then
[1146]1491! Write an Earth-format restart file
1492                CALL dynredem1_p("restart.nc",0.0,
1493     &                           vcov,ucov,teta,q,masse,ps)
[1446]1494!              endif ! of if (planet_type.eq."earth")
[630]1495
[1279]1496!              CLOSE(99)
[774]1497c$OMP END MASTER
[1146]1498            ENDIF ! of IF (itau.EQ.itaufin)
[630]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
[1146]1529            END IF ! of IF (MOD(itau,iperiod).EQ.0)
1530                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
[630]1531
1532
[1146]1533      ELSE ! of IF (.not.purmats)
1534
[630]1535c       ........................................................
1536c       ..............       schema  matsuno        ...............
1537c       ........................................................
1538            IF( forward )  THEN
1539
1540             itau =  itau + 1
[1279]1541!             iday = day_ini+itau/day_step
[1403]1542!             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]1543!
1544!                  IF(time.GT.1.) THEN
1545!                   time = time-1.
1546!                   iday = iday+1
1547!                  ENDIF
[630]1548
1549               forward =  .FALSE.
1550               IF( itau. EQ. itaufinp1 ) then 
[774]1551c$OMP MASTER
[1279]1552                 call fin_getparam
[949]1553                 call finalize_parallel
[774]1554c$OMP END MASTER
[630]1555                 abort_message = 'Simulation finished'
1556                 call abort_gcm(modname,abort_message,0)
[985]1557                 RETURN
[630]1558               ENDIF
1559               GO TO 2
1560
[1403]1561            ELSE ! of IF(forward) i.e. backward step
[630]1562
[1146]1563              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[630]1564               IF(itau.EQ.itaufin) THEN
1565                  iav=1
1566               ELSE
1567                  iav=0
1568               ENDIF
1569#ifdef CPP_IOIPSL
[1146]1570               IF (ok_dynzon) THEN
[774]1571c$OMP BARRIER
[1146]1572               call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1573               call SendRequest(TestRequest)
[985]1574c$OMP BARRIER
[1146]1575               call WaitRequest(TestRequest)
[1000]1576c$OMP BARRIER
[985]1577c$OMP MASTER
[1279]1578!               CALL writedynav_p(histaveid, itau,vcov ,
1579!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
[1146]1580               CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
[985]1581     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[774]1582c$OMP END MASTER
[1146]1583               END IF !ok_dynzon
[630]1584#endif
[1403]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
[1146]1608              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[630]1609
[1146]1610
[1403]1611               IF(MOD(itau,iecri         ).EQ.0) THEN
1612c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
[774]1613c$OMP BARRIER
1614c$OMP MASTER
[1146]1615                nbetat = nbetatdem
1616                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[630]1617
1618cym        unat=0.
[1146]1619                ijb=ij_begin
1620                ije=ij_end
[949]1621       
[1146]1622                if (pole_nord) then
1623                  ijb=ij_begin+iip1
1624                  unat(1:iip1,:)=0.
1625                endif
[949]1626       
[1146]1627                if (pole_sud) then
1628                  ije=ij_end-iip1
1629                  unat(ij_end-iip1+1:ij_end,:)=0.
1630                endif
[949]1631           
[1146]1632                do l=1,llm
1633                  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1634                enddo
[630]1635
[1146]1636                ijb=ij_begin
1637                ije=ij_end
1638                if (pole_sud) ije=ij_end-iip1
[949]1639       
[1146]1640                do l=1,llm
1641                  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1642                enddo
[630]1643
1644#ifdef CPP_IOIPSL
[1403]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
[1279]1660!                CALL writehist_p(histid, histvid, itau,vcov ,
1661!     &                           ucov,teta,phi,q,masse,ps,phis)
[1403]1662              endif ! of if (ok_dyn_ins)
[1146]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
[630]1673c     
[1146]1674                  if (mpi_rank==0) then
1675#include "write_grads_dyn.h"
1676                  endif
1677                endif ! of if (output_grads_dyn)
[630]1678
[774]1679c$OMP END MASTER
[1403]1680              ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]1681
[1146]1682              IF(itau.EQ.itaufin) THEN
[1446]1683!                if (planet_type.eq."earth") then
[774]1684c$OMP MASTER
1685                   CALL dynredem1_p("restart.nc",0.0,
[1146]1686     .                               vcov,ucov,teta,q,masse,ps)
[774]1687c$OMP END MASTER
[1446]1688!                endif ! of if (planet_type.eq."earth")
[1146]1689              ENDIF ! of IF(itau.EQ.itaufin)
[630]1690
[1146]1691              forward = .TRUE.
1692              GO TO  1
[630]1693
[1146]1694            ENDIF ! of IF (forward)
1695
1696      END IF ! of IF(.not.purmats)
[774]1697c$OMP MASTER
[1279]1698      call fin_getparam
[1146]1699      call finalize_parallel
[774]1700c$OMP END MASTER
[1146]1701      RETURN
[630]1702      END
Note: See TracBrowser for help on using the repository browser.