source: LMDZ4/branches/LMDZ4_par_0/libf/dyn3dpar/leapfrog_p.F

Last change on this file was 630, checked in by Laurent Fairhead, 19 years ago

Import d'une version parallele de la dynamique YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 40.1 KB
Line 
1!
2! $Header$
3!
4c
5c
6#define IO_DEBUG
7
8#undef CPP_IOIPSL
9
10      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
11     &                    time_0)
12
13       USE misc_mod
14       USE parallel
15       USE times
16       USE mod_hallo
17       USE Bands
18       USE Write_Field
19       USE Write_Field_p
20       USE vampir
21       
22#ifdef INCA
23      USE transport_controls, ONLY : hadv_flg, mmt_adj
24#endif
25
26      IMPLICIT NONE
27
28c      ......   Version  du 10/01/98    ..........
29
30c             avec  coordonnees  verticales hybrides
31c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
32
33c=======================================================================
34c
35c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
36c   -------
37c
38c   Objet:
39c   ------
40c
41c   GCM LMD nouvelle grille
42c
43c=======================================================================
44c
45c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
46c      et possibilite d'appeler une fonction f(y)  a derivee tangente
47c      hyperbolique a la  place de la fonction a derivee sinusoidale.
48
49c  ... Possibilite de choisir le shema pour l'advection de
50c        q  , en modifiant iadv dans traceur.def  (10/02) .
51c
52c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
53c      Pour Van-Leer iadv=10
54c
55c-----------------------------------------------------------------------
56c   Declarations:
57c   -------------
58
59#include "dimensions.h"
60#include "paramet.h"
61#include "comconst.h"
62#include "comdissnew.h"
63#include "comvert.h"
64#include "comgeom.h"
65#include "logic.h"
66#include "temps.h"
67#include "control.h"
68#include "ener.h"
69#include "description.h"
70#include "serre.h"
71#include "com_io_dyn.h"
72#include "iniprint.h"
73
74c#include "tracstoke.h"
75
76#include "academic.h"
77     
78      include 'mpif.h'
79      integer nq
80
81      INTEGER         longcles
82      PARAMETER     ( longcles = 20 )
83      REAL  clesphy0( longcles )
84
85      real zqmin,zqmax
86      INTEGER nbetatmoy, nbetatdem,nbetat
87
88c   variables dynamiques
89      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
90      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
91      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
92      REAL ps(ip1jmp1)                       ! pression  au sol
93      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
94      REAL pks(ip1jmp1)                      ! exner au  sol
95      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
96      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
97      REAL masse(ip1jmp1,llm)                ! masse d'air
98      REAL phis(ip1jmp1)                     ! geopotentiel au sol
99      REAL phi(ip1jmp1,llm)                  ! geopotentiel
100      REAL w(ip1jmp1,llm)                    ! vitesse verticale
101
102c variables dynamiques intermediaire pour le transport
103      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
104
105c   variables dynamiques au pas -1
106      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
107      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
108      REAL massem1(ip1jmp1,llm)
109
110c   tendances dynamiques
111      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
112      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
113
114c   tendances de la dissipation
115      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
116      REAL dtetadis(ip1jmp1,llm)
117
118c   tendances physiques
119      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
120      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
121
122c   variables pour le fichier histoire
123      REAL dtav      ! intervalle de temps elementaire
124
125      REAL tppn(iim),tpps(iim),tpn,tps
126c
127      INTEGER itau,itaufinp1,iav
128      INTEGER*4  iday ! jour julien
129      REAL       time ! Heure de la journee en fraction d'1 jour
130
131      REAL  SSUM
132      REAL time_0 , finvmaold(ip1jmp1,llm)
133
134cym      LOGICAL  lafin
135      LOGICAL :: lafin=.false.
136      INTEGER ij,iq,l
137      INTEGER ik
138
139      real time_step, t_wrt, t_ops
140
141      REAL rdayvrai,rdaym_ini
142      LOGICAL first,callinigrads
143
144      data callinigrads/.true./
145      character*10 string10
146
147      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
148#ifdef INCA_CH4
149      REAL :: flxw(ip1jmp1,llm)
150#endif
151
152c+jld variables test conservation energie
153      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
154C     Tendance de la temp. potentiel d (theta)/ d t due a la
155C     tansformation d'energie cinetique en energie thermique
156C     cree par la dissipation
157      REAL dtetaecdt(ip1jmp1,llm)
158      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
159      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
160      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
161      CHARACTER*15 ztit
162      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
163      SAVE      ip_ebil_dyn
164      DATA      ip_ebil_dyn/0/
165c-jld
166
167      character*80 dynhist_file, dynhistave_file
168      character*20 modname
169      character*80 abort_message
170
171C Calendrier
172      LOGICAL true_calendar
173      PARAMETER (true_calendar = .false.)
174
175      logical dissip_conservative
176      save dissip_conservative
177      data dissip_conservative/.true./
178
179      LOGICAL prem
180      save prem
181      DATA prem/.true./
182      INTEGER testita
183      PARAMETER (testita = 9)
184     
185c declaration liées au parallelisme
186      INTEGER :: ierr
187      LOGICAL :: FirstCaldyn=.TRUE.
188      LOGICAL :: FirstPhysic=.TRUE.
189      INTEGER :: ijb,ije,j,i
190      type(Request) :: TestRequest
191      type(Request) :: Request_Dissip
192      type(Request) :: Request_physic
193      REAL dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
194      REAL dtetafi_tmp(iip1,llm),dqfi_tmp(iip1,llm,nqmx)
195      REAL dpfi_tmp(iip1)
196      INTEGER :: true_itau=0
197      LOGICAL :: verbose=.true.
198      INTEGER :: iapptrac = 0
199      INTEGER :: AdjustCount = 0
200     
201      ItCount=0
202     
203      itaufin   = nday*day_step
204      itaufinp1 = itaufin +1
205
206
207      itau = 0
208      iday = day_ini+itau/day_step
209      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
210         IF(time.GT.1.) THEN
211          time = time-1.
212          iday = iday+1
213         ENDIF
214
215
216c-----------------------------------------------------------------------
217c   On initialise la pression et la fonction d'Exner :
218c   --------------------------------------------------
219
220      dq=0.
221      CALL pression ( ip1jmp1, ap, bp, ps, p       )
222      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
223
224c-----------------------------------------------------------------------
225c   Debut de l'integration temporelle:
226c   ----------------------------------
227c et du parallélisme !!
228
229   1  CONTINUE
230
231      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
232
233#ifdef CPP_IOIPSL
234      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
235        call guide(itau,ucov,vcov,teta,q,masse,ps)
236      else
237        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
238     .    '6 dernieres heures'
239      endif
240#endif
241c
242c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
243c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
244c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
245c     ENDIF
246c
247cym      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
248cym      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
249cym      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
250cym      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
251cym      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
252
253       if (FirstCaldyn) then
254         ucovm1=ucov
255         vcovm1=vcov
256         tetam1= teta
257         massem1= masse
258         psm1= ps
259         
260         finvmaold = masse
261         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
262
263       else
264         
265         ijb=ij_begin
266         ije=ij_end
267               
268         ucovm1   (ijb:ije,:) = ucov  (ijb:ije,:)
269         tetam1   (ijb:ije,:) = teta  (ijb:ije,:)
270         massem1  (ijb:ije,:) = masse (ijb:ije,:)
271         psm1     (ijb:ije) = ps    (ijb:ije)
272     
273         if (pole_sud) ije=ij_end-iip1
274         vcovm1(ijb:ije,:) = vcov  (ijb:ije,:)
275       
276         finvmaold(ij_begin:ij_end,1:llm)=masse(ij_begin:ij_end,1:llm)
277         CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
278     .                    llm, -2,2, .TRUE., 1 )
279
280       endif
281       
282      forward = .TRUE.
283      leapf   = .FALSE.
284      dt      =  dtvr
285
286c   ...    P.Le Van .26/04/94  ....
287
288cym      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
289cym      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
290
291cym  ne sert à rien
292cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
293
294   2  CONTINUE
295
296      ItCount=ItCount+1
297      if (MOD(ItCount,10000)==0) then
298        debug=.true.
299      else
300        debug=.false.
301      endif
302c-----------------------------------------------------------------------
303
304c   date:
305c   -----
306
307
308c   gestion des appels de la physique et des dissipations:
309c   ------------------------------------------------------
310c
311c   ...    P.Le Van  ( 6/02/95 )  ....
312
313      apphys = .FALSE.
314      statcl = .FALSE.
315      conser = .FALSE.
316      apdiss = .FALSE.
317
318      IF( purmats ) THEN
319         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
320         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
321         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
322     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
323      ELSE
324         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
325         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
326         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
327      END IF
328
329cym    ---> Pour le moment     
330cym      apphys = .FALSE.
331      statcl = .FALSE.
332      conser = .FALSE.
333     
334      if (firstCaldyn) then
335
336          call SetDistrib(jj_Nb_Caldyn)
337          firstCaldyn=.FALSE.
338cym          call InitTime
339          call Init_timer
340      endif
341     
342     
343     
344     
345      if (Adjust) then
346          AdjustCount=AdjustCount+1
347c         if (MOD(itau+1,20)==0) then
348          if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
349     &         .and. .not.FirstPhysic .and. Adjustcount>30) then
350           AdjustCount=0
351           call allgather_timer_average
352
353        if (Verbose) then
354       
355        print *,'*********************************'
356        print *,'******    TIMER CALDYN     ******'
357        do i=0,mpi_size-1
358          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
359     &            '  : temps moyen :',
360     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i),
361     &            '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
362        enddo
363     
364        print *,'*********************************'
365        print *,'******    TIMER VANLEER    ******'
366        do i=0,mpi_size-1
367          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
368     &            '  : temps moyen :',
369     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i),
370     &            '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
371        enddo
372     
373        print *,'*********************************'
374        print *,'******    TIMER DISSIP    ******'
375        do i=0,mpi_size-1
376          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
377     &            '  : temps moyen :',
378     &             timer_average(jj_nb_dissip(i),timer_dissip,i),
379     &             '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
380        enddo
381       
382        if (mpi_rank==0) call WriteBands
383       
384       endif
385       
386         call AdjustBands_caldyn
387         if (mpi_rank==0) call WriteBands
388         
389         call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
390     &                                jj_Nb_caldyn,0,0,TestRequest)
391         call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
392     &                                jj_Nb_caldyn,0,0,TestRequest)
393         call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
394     &                                jj_Nb_caldyn,0,0,TestRequest)
395         call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
396     &                                jj_Nb_caldyn,0,0,TestRequest)
397         call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
398     &                                jj_Nb_caldyn,0,0,TestRequest)
399         call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
400     &                                jj_Nb_caldyn,0,0,TestRequest)
401         call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
402     &                                jj_Nb_caldyn,0,0,TestRequest)
403         call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
404     &                                jj_Nb_caldyn,0,0,TestRequest)
405         call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
406     &                                jj_Nb_caldyn,0,0,TestRequest)
407         call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
408     &                                jj_Nb_caldyn,0,0,TestRequest)
409         call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
410     &                                jj_Nb_caldyn,0,0,TestRequest)
411         call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
412     &                                jj_Nb_caldyn,0,0,TestRequest)
413         call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
414     &                                jj_Nb_caldyn,0,0,TestRequest)
415         call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
416     &                                jj_Nb_caldyn,0,0,TestRequest)
417        call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
418     &                                jj_Nb_caldyn,0,0,TestRequest)
419         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
420     &                                jj_Nb_caldyn,0,0,TestRequest)
421 
422         call SetDistrib(jj_nb_caldyn)
423         call SendRequest(TestRequest)
424         call WaitRequest(TestRequest)
425         
426c       call SetDistrib(jj_Nb_vanleer)
427c        call AdjustBands_vanleer
428c       do j=1,nqmx
429c         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
430c     *                                jj_nb_vanleer,0,0,TestRequest)
431c        enddo
432c
433c         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
434c     &                                jj_Nb_vanleer,0,0,TestRequest)
435       
436c        call SendRequest(TestRequest)
437c        call WaitRequest(TestRequest)
438
439        call AdjustBands_dissip
440        call AdjustBands_physic
441c        call SetDistrib(jj_nb_caldyn)
442      endif
443     
444      endif       
445       
446     
447     
448c-----------------------------------------------------------------------
449c   calcul des tendances dynamiques:
450c   --------------------------------
451
452       call VTb(VThallo)
453       call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
454       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
455       call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
456       call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
457       call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
458       call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
459       call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
460       call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
461       
462c       do j=1,nqmx
463c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
464c     *                       TestRequest)
465c        enddo
466
467       call SendRequest(TestRequest)
468       call WaitRequest(TestRequest)
469       call VTe(VThallo)
470
471      if (debug) then   
472             
473        call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
474        call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
475        call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
476        call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
477        call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
478        call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
479        call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
480        call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
481c        do j=1,nqmx
482c          call WriteField_p('q'//trim(int2str(j)),
483c     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
484c        enddo       
485      endif
486     
487
488     
489      True_itau=True_itau+1
490      print*,"Iteration No",True_itau
491
492
493      call start_timer(timer_caldyn)
494
495      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
496
497     
498      call VTb(VTcaldyn)
499     
500      CALL caldyn_p
501     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
502     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
503     
504      call VTe(VTcaldyn)
505c      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
506c      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
507c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
508c      call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
509c      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
510c      call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
511c      call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
512
513c-----------------------------------------------------------------------
514c   calcul des tendances advection des traceurs (dont l'humidite)
515c   -------------------------------------------------------------
516
517      IF( forward. OR . leapf )  THEN
518
519c
520#ifdef INCA_CH4
521             CALL caladvtrac_p(q,pbaru,pbarv,
522     *                      p, masse, dq,  teta,
523     .             flxw,
524     .             pk,
525     .             mmt_adj,
526     .             hadv_flg,iapptrac)
527#else
528             CALL caladvtrac_p(q,pbaru,pbarv,
529     *                      p, masse, dq,  teta,
530     .             pk,iapptrac)
531#endif
532c      do j=1,nqmx
533c        call WriteField_p('q'//trim(int2str(j)),
534c    .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
535c        call WriteField_p('dq'//trim(int2str(j)),
536c     .                reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
537c      enddo
538
539         IF (offline) THEN
540Cmaf stokage du flux de masse pour traceurs OFF-LINE
541
542#ifdef CPP_IOIPSL
543           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
544     .   dtvr, itau)
545#endif
546
547
548         ENDIF
549c
550      ENDIF
551
552
553c-----------------------------------------------------------------------
554c   integrations dynamique et traceurs:
555c   ----------------------------------
556 
557       call VTb(VTintegre)
558       CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
559     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
560     $              finvmaold                                    )
561
562       call VTe(VTintegre)
563c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
564c
565c-----------------------------------------------------------------------
566c   calcul des tendances physiques:
567c   -------------------------------
568c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
569c
570       IF( purmats )  THEN
571          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
572       ELSE
573          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
574       ENDIF
575c
576c
577       IF( apphys )  THEN
578c
579c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
580c
581         call suspend_timer(timer_caldyn)
582         print*,'Entree dans la physique : Iteration No ',true_itau
583         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )
584         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
585
586           rdaym_ini  = itau * dtvr / daysec
587           rdayvrai   = rdaym_ini  + day_ini
588
589
590c rajout debug
591c       lafin = .true.
592
593
594c   Inbterface avec les routines de phylmd (phymars ... )
595c   -----------------------------------------------------
596
597#ifdef CPP_PHYS
598c+jld
599
600c  Diagnostique de conservation de l'énergie : initialisation
601      IF (ip_ebil_dyn.ge.1 ) THEN
602          ztit='bil dyn'
603          CALL diagedyn(ztit,2,1,1,dtphys
604     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
605      ENDIF
606c-jld
607        call VTb(VThallo)
608        call SetTag(Request_physic,800)
609       
610        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
611     *                               jj_Nb_physic,2,2,Request_physic)
612       
613        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
614     *                               jj_Nb_physic,2,2,Request_physic)
615       
616        call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
617     *                               jj_Nb_physic,2,2,Request_physic)
618       
619        call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
620     *                               jj_Nb_physic,2,2,Request_physic)
621       
622        call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
623     *                               jj_Nb_physic,2,2,Request_physic)
624       
625        call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
626     *                               jj_Nb_physic,2,2,Request_physic)
627       
628        call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
629     *                               jj_Nb_physic,2,2,Request_physic)
630       
631        call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
632     *                               jj_Nb_physic,2,2,Request_physic)
633       
634c       call SetDistrib(jj_nb_vanleer)
635        do j=1,nqmx
636 
637          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
638     *                               jj_Nb_physic,2,2,Request_physic)
639        enddo
640#ifdef INCA_CH4
641        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
642     *                               jj_Nb_physic,2,2,Request_physic)
643#endif
644        call SetDistrib(jj_nb_Physic)
645       
646        call SendRequest(Request_Physic)
647        call WaitRequest(Request_Physic)       
648       
649        call VTe(VThallo)
650       
651        call VTb(VTphysiq)
652        CALL calfis_p( nq, lafin ,rdayvrai,time  ,
653     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
654     $               du,dv,dteta,dq,w,
655#ifdef INCA_CH4
656     $               flxw,
657#endif
658     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
659     
660        ijb=ij_begin
661        ije=ij_end 
662        if ( .not. pole_nord) then
663          dufi_tmp(1:iip1,:)   = dufi(ijb:ijb+iim,:)
664          dvfi_tmp(1:iip1,:)   = dvfi(ijb:ijb+iim,:) 
665          dtetafi_tmp(1:iip1,:)= dtetafi(ijb:ijb+iim,:) 
666          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim) 
667          dqfi_tmp(1:iip1,:,:) = dqfi(ijb:ijb+iim,:,:) 
668        endif
669       
670        call SetDistrib(jj_nb_Physic_bis)
671
672        call VTb(VThallo)
673 
674        call Register_Hallo(dufi,ip1jmp1,llm,
675     *                      1,0,0,1,Request_physic)
676       
677        call Register_Hallo(dvfi,ip1jm,llm,
678     *                      1,0,0,1,Request_physic)
679       
680        call Register_Hallo(dtetafi,ip1jmp1,llm,
681     *                      1,0,0,1,Request_physic)
682
683        call Register_Hallo(dpfi,ip1jmp1,1,
684     *                      1,0,0,1,Request_physic)
685
686        do j=1,nqmx
687          call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
688     *                        1,0,0,1,Request_physic)
689        enddo
690       
691        call SendRequest(Request_Physic)
692        call WaitRequest(Request_Physic)
693             
694        call VTe(VThallo)
695 
696        call SetDistrib(jj_nb_Physic)
697       
698                ijb=ij_begin
699        if (.not. pole_nord) then
700          dufi(ijb:ijb+iim,:) = dufi(ijb:ijb+iim,:)+dufi_tmp(1:iip1,:)
701          dvfi(ijb:ijb+iim,:) = dvfi(ijb:ijb+iim,:)+dvfi_tmp(1:iip1,:)
702          dtetafi(ijb:ijb+iim,:) = dtetafi(ijb:ijb+iim,:)
703     &                           +dtetafi_tmp(1:iip1,:)
704          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
705          dqfi(ijb:ijb+iim,:,:) = dqfi(ijb:ijb+iim,:,:)
706     &                           + dqfi_tmp(1:iip1,:,:)
707        endif
708       
709c      call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
710c      call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
711c      call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
712c      call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
713c     
714c      do j=1,nqmx
715c        call WriteField_p('dqfi'//trim(int2str(j)),
716c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
717c      enddo
718
719c      ajout des tendances physiques:
720c      ------------------------------
721          CALL addfi_p( nqmx, dtphys, leapf, forward   ,
722     $                  ucov, vcov, teta , q   ,ps ,
723     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
724
725        call VTe(VTphysiq)
726
727        call VTb(VThallo)
728
729        call SetTag(Request_physic,800)
730        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
731     *                               jj_Nb_caldyn,Request_physic)
732       
733        call Register_SwapField(vcov,vcov,ip1jm,llm,
734     *                               jj_Nb_caldyn,Request_physic)
735       
736        call Register_SwapField(teta,teta,ip1jmp1,llm,
737     *                               jj_Nb_caldyn,Request_physic)
738       
739        call Register_SwapField(p,p,ip1jmp1,llmp1,
740     *                               jj_Nb_caldyn,Request_physic)
741       
742        call Register_SwapField(pk,pk,ip1jmp1,llm,
743     *                               jj_Nb_caldyn,Request_physic)
744       
745        call Register_SwapField(phis,phis,ip1jmp1,1,
746     *                               jj_Nb_caldyn,Request_physic)
747       
748        call Register_SwapField(phi,phi,ip1jmp1,llm,
749     *                               jj_Nb_caldyn,Request_physic)
750       
751        call Register_SwapField(w,w,ip1jmp1,llm,
752     *                               jj_Nb_caldyn,Request_physic)
753
754        do j=1,nqmx
755       
756          call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
757     *                               jj_Nb_caldyn,Request_physic)
758       
759        enddo
760
761        call SendRequest(Request_Physic)
762        call WaitRequest(Request_Physic)     
763
764        call VTe(VThallo)
765
766        call SetDistrib(jj_Nb_caldyn)
767c
768c  Diagnostique de conservation de l'énergie : difference
769      IF (ip_ebil_dyn.ge.1 ) THEN
770          ztit='bil phys'
771          CALL diagedyn(ztit,2,1,1,dtphys
772     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
773      ENDIF
774     
775      if (debug) then
776       call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
777       call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
778       call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
779      endif
780#else
781
782c   Calcul academique de la physique = Rappel Newtonien + fritcion
783c   --------------------------------------------------------------
784cym       teta(:,:)=teta(:,:)
785cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
786       ijb=ij_begin
787       ije=ij_end
788       teta(ijb:ije,:)=teta(ijb:ije,:)
789     s  -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
790
791       call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
792       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
793       call SendRequest(Request_Physic)
794       call WaitRequest(Request_Physic)     
795
796       call friction_p(ucov,vcov,iphysiq*dtvr)
797
798#endif
799
800c-jld
801         call resume_timer(timer_caldyn)
802         if (FirstPhysic) then
803           call InitTime
804           FirstPhysic=.false.
805         endif
806       ENDIF
807
808        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
809        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
810
811
812c-----------------------------------------------------------------------
813c   dissipation horizontale et verticale  des petites echelles:
814c   ----------------------------------------------------------
815
816      IF(apdiss) THEN
817        call suspend_timer(timer_caldyn)
818       
819        print*,'Entree dans la dissipation : Iteration No ',true_itau
820c   calcul de l'energie cinetique avant dissipation
821        print *,'Passage dans la dissipation'
822
823        call VTb(VThallo)
824
825        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
826     *                          jj_Nb_dissip,1,1,Request_dissip)
827
828        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
829     *                          jj_Nb_dissip,1,1,Request_dissip)
830
831        call Register_SwapField(teta,teta,ip1jmp1,llm,
832     *                          jj_Nb_dissip,Request_dissip)
833
834        call Register_SwapField(p,p,ip1jmp1,llmp1,
835     *                          jj_Nb_dissip,Request_dissip)
836
837        call Register_SwapField(pk,pk,ip1jmp1,llm,
838     *                          jj_Nb_dissip,Request_dissip)
839
840        call SendRequest(Request_dissip)       
841        call WaitRequest(Request_dissip)       
842        call SetDistrib(jj_Nb_dissip)
843       
844        call VTe(VThallo)
845
846        call VTb(VTdissipation)
847       
848        call start_timer(timer_dissip)
849        call covcont_p(llm,ucov,vcov,ucont,vcont)
850        call enercin_p(vcov,ucov,vcont,ucont,ecin0)
851
852c   dissipation
853
854        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
855       
856        ijb=ij_begin
857        ije=ij_end
858       
859        ucov(ijb:ije,1:llm)=ucov(ijb:ije,1:llm)+dudis(ijb:ije,1:llm)
860       
861        if (pole_sud) ije=ije-iip1
862        vcov(ijb:ije,1:llm)=vcov(ijb:ije,1:llm)+dvdis(ijb:ije,1:llm)
863c       teta=teta+dtetadis
864
865
866c------------------------------------------------------------------------
867        if (dissip_conservative) then
868C       On rajoute la tendance due a la transform. Ec -> E therm. cree
869C       lors de la dissipation
870            call suspend_timer(timer_dissip)
871            call VTb(VThallo)
872
873            call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
874            call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
875            call SendRequest(Request_Dissip)
876            call WaitRequest(Request_Dissip)
877            call VTe(VThallo)
878            call resume_timer(timer_dissip)
879           
880            call covcont_p(llm,ucov,vcov,ucont,vcont)
881            call enercin_p(vcov,ucov,vcont,ucont,ecin)
882           
883            ijb=ij_begin
884            ije=ij_end
885           
886            do l=1,llm
887              do ij=ijb,ije
888                dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
889                dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
890              enddo
891            enddo
892           
893       endif
894
895       ijb=ij_begin
896       ije=ij_end
897           
898         do l=1,llm
899           do ij=ijb,ije
900              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
901           enddo
902         enddo
903         
904c------------------------------------------------------------------------
905
906
907c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
908c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
909c
910
911        ijb=ij_begin
912        ije=ij_end
913         
914        if (pole_nord) then
915          DO l  =  1, llm
916            DO ij =  1,iim
917             tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
918            ENDDO
919             tpn  = SSUM(iim,tppn,1)/apoln
920
921            DO ij = 1, iip1
922             teta(  ij    ,l) = tpn
923            ENDDO
924          ENDDO
925               
926          DO ij =  1,iim
927            tppn(ij)  = aire(  ij    ) * ps (  ij    )
928          ENDDO
929            tpn  = SSUM(iim,tppn,1)/apoln
930 
931          DO ij = 1, iip1
932            ps(  ij    ) = tpn
933          ENDDO
934        endif
935       
936        if (pole_sud) then
937          DO l  =  1, llm
938            DO ij =  1,iim
939             tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
940            ENDDO
941             tps  = SSUM(iim,tpps,1)/apols
942
943            DO ij = 1, iip1
944             teta(ij+ip1jm,l) = tps
945            ENDDO
946          ENDDO
947               
948          DO ij =  1,iim
949            tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
950          ENDDO
951            tps  = SSUM(iim,tpps,1)/apols
952 
953          DO ij = 1, iip1
954            ps(ij+ip1jm) = tps
955          ENDDO
956        endif
957
958        call VTe(VTdissipation)
959
960        call stop_timer(timer_dissip)
961       
962        call VTb(VThallo)
963
964        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
965     *                          jj_Nb_caldyn,Request_dissip)
966
967        call Register_SwapField(vcov,vcov,ip1jm,llm,
968     *                          jj_Nb_caldyn,Request_dissip)
969
970        call Register_SwapField(teta,teta,ip1jmp1,llm,
971     *                          jj_Nb_caldyn,Request_dissip)
972
973        call Register_SwapField(p,p,ip1jmp1,llmp1,
974     *                          jj_Nb_caldyn,Request_dissip)
975
976        call Register_SwapField(pk,pk,ip1jmp1,llm,
977     *                          jj_Nb_caldyn,Request_dissip)
978
979        call SendRequest(Request_dissip)       
980        call WaitRequest(Request_dissip)       
981        call SetDistrib(jj_Nb_caldyn)
982        call VTe(VThallo)
983        call resume_timer(timer_caldyn)
984        print *,'fin dissipation'
985      END IF
986
987c ajout debug
988c              IF( lafin ) then 
989c                abort_message = 'Simulation finished'
990c                call abort_gcm(modname,abort_message,0)
991c              ENDIF
992       
993c   ********************************************************************
994c   ********************************************************************
995c   .... fin de l'integration dynamique  et physique pour le pas itau ..
996c   ********************************************************************
997c   ********************************************************************
998
999c   preparation du pas d'integration suivant  ......
1000cym      call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1001cym      call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1002     
1003      call stop_timer(timer_caldyn)
1004      IF (itau==itaumax) then
1005            call allgather_timer_average
1006
1007      if (mpi_rank==0) then
1008       
1009        print *,'*********************************'
1010        print *,'******    TIMER CALDYN     ******'
1011        do i=0,mpi_size-1
1012          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
1013     &            '  : temps moyen :',
1014     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1015        enddo
1016     
1017        print *,'*********************************'
1018        print *,'******    TIMER VANLEER    ******'
1019        do i=0,mpi_size-1
1020          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
1021     &            '  : temps moyen :',
1022     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1023        enddo
1024     
1025        print *,'*********************************'
1026        print *,'******    TIMER DISSIP    ******'
1027        do i=0,mpi_size-1
1028          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
1029     &            '  : temps moyen :',
1030     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
1031        enddo
1032       
1033        print *,'*********************************'
1034        print *,'******    TIMER PHYSIC    ******'
1035        do i=0,mpi_size-1
1036          print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
1037     &            '  : temps moyen :',
1038     &             timer_average(jj_nb_physic(i),timer_physic,i)
1039        enddo
1040       
1041      endif 
1042     
1043      print *,'Taille du Buffer MPI (REAL*8)',MaxBufferSize
1044      print *,'Taille du Buffer MPI utilise (REAL*8)',MaxBufferSize_Used
1045      print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1046      print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1047        call finalize_parallel
1048        STOP
1049      ENDIF
1050     
1051      IF ( .NOT.purmats ) THEN
1052c       ........................................................
1053c       ..............  schema matsuno + leapfrog  ..............
1054c       ........................................................
1055
1056            IF(forward. OR. leapf) THEN
1057              itau= itau + 1
1058              iday= day_ini+itau/day_step
1059              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1060                IF(time.GT.1.) THEN
1061                  time = time-1.
1062                  iday = iday+1
1063                ENDIF
1064            ENDIF
1065
1066
1067            IF( itau. EQ. itaufinp1 ) then 
1068c$$$       write(79,*) 'ucov',ucov
1069c$$$       write(80,*) 'vcov',vcov
1070c$$$       write(81,*) 'teta',teta
1071c$$$       write(82,*) 'ps',ps
1072c$$$       write(83,*) 'q',q
1073c$$$       WRITE(85,*) 'q1 = ',q(:,:,1)
1074c$$$       WRITE(86,*) 'q3 = ',q(:,:,3)
1075
1076              abort_message = 'Simulation finished'
1077
1078              call abort_gcm(modname,abort_message,0)
1079            ENDIF
1080c-----------------------------------------------------------------------
1081c   ecriture du fichier histoire moyenne:
1082c   -------------------------------------
1083
1084            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1085               IF(itau.EQ.itaufin) THEN
1086                  iav=1
1087               ELSE
1088                  iav=0
1089               ENDIF
1090#ifdef CPP_IOIPSL
1091              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1092              call SendRequest(TestRequest)
1093              call WaitRequest(TestRequest)
1094
1095              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
1096     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1097c               call bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1098c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1099#endif
1100
1101            ENDIF
1102
1103c-----------------------------------------------------------------------
1104c   ecriture de la bande histoire:
1105c   ------------------------------
1106
1107c      IF( MOD(itau,iecri         ).EQ.0) THEN
1108           IF( MOD(itau,iecri*day_step).EQ.0) THEN
1109
1110               nbetat = nbetatdem
1111        CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi )
1112       
1113cym     unat=0.
1114       
1115        ijb=ij_begin
1116        ije=ij_end
1117       
1118        if (pole_nord) then
1119          ijb=ij_begin+iip1
1120          unat(1:iip1,:)=0.
1121        endif
1122       
1123        if (pole_sud) then
1124          ije=ij_end-iip1
1125          unat(ij_end-iip1+1:ij_end,:)=0.
1126        endif
1127           
1128        do l=1,llm
1129           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1130        enddo
1131
1132        ijb=ij_begin
1133        ije=ij_end
1134        if (pole_sud) ije=ij_end-iip1
1135       
1136        do l=1,llm
1137           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1138        enddo
1139       
1140#ifdef CPP_IOIPSL
1141 
1142        CALL writehist_p(histid,histvid, nqmx,itau,vcov,
1143     s                       ucov,teta,phi,q,masse,ps,phis)
1144c#else
1145c       call Gather_Field(unat,ip1jmp1,llm,0)
1146c       call Gather_Field(vnat,ip1jm,llm,0)
1147c       call Gather_Field(teta,ip1jmp1,llm,0)
1148c       call Gather_Field(ps,ip1jmp1,1,0)
1149c       do iq=1,nqmx
1150c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1151c       enddo
1152c     
1153c       if (mpi_rank==0) then
1154c#include "write_grads_dyn.h"
1155c       endif
1156#endif
1157
1158            ENDIF
1159
1160            IF(itau.EQ.itaufin) THEN
1161
1162
1163#ifdef CPP_IOIPSL
1164       CALL dynredem1_p("restart.nc",0.0,
1165     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
1166#endif
1167
1168              CLOSE(99)
1169            ENDIF
1170
1171c-----------------------------------------------------------------------
1172c   gestion de l'integration temporelle:
1173c   ------------------------------------
1174
1175            IF( MOD(itau,iperiod).EQ.0 )    THEN
1176                    GO TO 1
1177            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1178
1179                   IF( forward )  THEN
1180c      fin du pas forward et debut du pas backward
1181
1182                      forward = .FALSE.
1183                        leapf = .FALSE.
1184                           GO TO 2
1185
1186                   ELSE
1187c      fin du pas backward et debut du premier pas leapfrog
1188
1189                        leapf =  .TRUE.
1190                        dt  =  2.*dtvr
1191                        GO TO 2
1192                   END IF
1193            ELSE
1194
1195c      ......   pas leapfrog  .....
1196
1197                 leapf = .TRUE.
1198                 dt  = 2.*dtvr
1199                 GO TO 2
1200            END IF
1201
1202      ELSE
1203
1204c       ........................................................
1205c       ..............       schema  matsuno        ...............
1206c       ........................................................
1207            IF( forward )  THEN
1208
1209             itau =  itau + 1
1210             iday = day_ini+itau/day_step
1211             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1212
1213                  IF(time.GT.1.) THEN
1214                   time = time-1.
1215                   iday = iday+1
1216                  ENDIF
1217
1218               forward =  .FALSE.
1219               IF( itau. EQ. itaufinp1 ) then 
1220                 abort_message = 'Simulation finished'
1221                 call abort_gcm(modname,abort_message,0)
1222               ENDIF
1223               GO TO 2
1224
1225            ELSE
1226
1227            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1228               IF(itau.EQ.itaufin) THEN
1229                  iav=1
1230               ELSE
1231                  iav=0
1232               ENDIF
1233#ifdef CPP_IOIPSL
1234              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1235              call SendRequest(TestRequest)
1236              call WaitRequest(TestRequest)
1237
1238              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
1239     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1240c               call bilan_dyn_p (2,dtvr*iperiod,dtvr*day_step*periodav,
1241c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1242#endif
1243
1244            ENDIF
1245
1246c               IF(MOD(itau,iecri         ).EQ.0) THEN
1247              IF(MOD(itau,iecri*day_step).EQ.0) THEN
1248                  nbetat = nbetatdem
1249       CALL geopot_p( ip1jmp1, teta  , pk , pks,  phis  , phi   )
1250
1251cym        unat=0.
1252        ijb=ij_begin
1253        ije=ij_end
1254       
1255        if (pole_nord) then
1256          ijb=ij_begin+iip1
1257          unat(1:iip1,:)=0.
1258        endif
1259       
1260        if (pole_sud) then
1261          ije=ij_end-iip1
1262          unat(ij_end-iip1+1:ij_end,:)=0.
1263        endif
1264           
1265        do l=1,llm
1266           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1267        enddo
1268
1269        ijb=ij_begin
1270        ije=ij_end
1271        if (pole_sud) ije=ij_end-iip1
1272       
1273        do l=1,llm
1274           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1275        enddo
1276
1277#ifdef CPP_IOIPSL
1278
1279       CALL writehist_p( histid, histvid, nqmx, itau,vcov ,
1280     ,                           ucov,teta,phi,q,masse,ps,phis)
1281c#else
1282c      call Gather_Field(unat,ip1jmp1,llm,0)
1283c      call Gather_Field(vnat,ip1jm,llm,0)
1284c      call Gather_Field(teta,ip1jmp1,llm,0)
1285c      call Gather_Field(ps,ip1jmp1,1,0)
1286c      do iq=1,nqmx
1287c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1288c      enddo
1289c     
1290c      if (mpi_rank==0) then
1291c#include "write_grads_dyn.h"
1292c      endif
1293#endif
1294
1295
1296               ENDIF
1297
1298#ifdef CPP_IOIPSL
1299                 IF(itau.EQ.itaufin)
1300     . CALL dynredem1_p("restart.nc",0.0,
1301     .                     vcov,ucov,teta,q,nqmx,masse,ps)
1302#endif
1303
1304                 forward = .TRUE.
1305                 GO TO  1
1306
1307            ENDIF
1308
1309      END IF
1310
1311      STOP
1312      END
Note: See TracBrowser for help on using the repository browser.