source: LMDZ4/branches/LMDZ4_V2_patch/libf/dyn3d/leapfrog.F @ 1350

Last change on this file since 1350 was 665, checked in by (none), 19 years ago

This commit was manufactured by cvs2svn to create branch
'LMDZ4_V2_patch'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.1 KB
Line 
1!
2! $Header$
3!
4c
5c
6      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
7     &                    time_0)
8
9#ifdef INCA
10      USE transport_controls, ONLY : hadv_flg, mmt_adj
11#endif
12
13      IMPLICIT NONE
14
15c      ......   Version  du 10/01/98    ..........
16
17c             avec  coordonnees  verticales hybrides
18c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
19
20c=======================================================================
21c
22c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
23c   -------
24c
25c   Objet:
26c   ------
27c
28c   GCM LMD nouvelle grille
29c
30c=======================================================================
31c
32c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
33c      et possibilite d'appeler une fonction f(y)  a derivee tangente
34c      hyperbolique a la  place de la fonction a derivee sinusoidale.
35
36c  ... Possibilite de choisir le shema pour l'advection de
37c        q  , en modifiant iadv dans traceur.def  (10/02) .
38c
39c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
40c      Pour Van-Leer iadv=10
41c
42c-----------------------------------------------------------------------
43c   Declarations:
44c   -------------
45
46#include "dimensions.h"
47#include "paramet.h"
48#include "comconst.h"
49#include "comdissnew.h"
50#include "comvert.h"
51#include "comgeom.h"
52#include "logic.h"
53#include "temps.h"
54#include "control.h"
55#include "ener.h"
56#include "description.h"
57#include "serre.h"
58#include "com_io_dyn.h"
59#include "iniprint.h"
60
61c#include "tracstoke.h"
62
63#include "academic.h"
64
65      integer nq
66
67      INTEGER         longcles
68      PARAMETER     ( longcles = 20 )
69      REAL  clesphy0( longcles )
70
71      real zqmin,zqmax
72      INTEGER nbetatmoy, nbetatdem,nbetat
73
74c   variables dynamiques
75      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
76      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
77      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
78      REAL ps(ip1jmp1)                       ! pression  au sol
79      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
80      REAL pks(ip1jmp1)                      ! exner au  sol
81      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
82      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
83      REAL masse(ip1jmp1,llm)                ! masse d'air
84      REAL phis(ip1jmp1)                     ! geopotentiel au sol
85      REAL phi(ip1jmp1,llm)                  ! geopotentiel
86      REAL w(ip1jmp1,llm)                    ! vitesse verticale
87
88c variables dynamiques intermediaire pour le transport
89      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
90
91c   variables dynamiques au pas -1
92      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
93      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
94      REAL massem1(ip1jmp1,llm)
95
96c   tendances dynamiques
97      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
98      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
99
100c   tendances de la dissipation
101      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
102      REAL dtetadis(ip1jmp1,llm)
103
104c   tendances physiques
105      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
106      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
107
108c   variables pour le fichier histoire
109      REAL dtav      ! intervalle de temps elementaire
110
111      REAL tppn(iim),tpps(iim),tpn,tps
112c
113      INTEGER itau,itaufinp1,iav
114      INTEGER*4  iday ! jour julien
115      REAL       time ! Heure de la journee en fraction d'1 jour
116
117      REAL  SSUM
118      REAL time_0 , finvmaold(ip1jmp1,llm)
119
120cym      LOGICAL  lafin
121      LOGICAL :: lafin=.false.
122      INTEGER ij,iq,l
123      INTEGER ik
124
125      real time_step, t_wrt, t_ops
126
127      REAL rdayvrai,rdaym_ini
128      LOGICAL first,callinigrads
129
130      data callinigrads/.true./
131      character*10 string10
132
133      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
134#ifdef INCA
135      REAL :: flxw(ip1jmp1,llm)
136#endif
137
138c+jld variables test conservation energie
139      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
140C     Tendance de la temp. potentiel d (theta)/ d t due a la
141C     tansformation d'energie cinetique en energie thermique
142C     cree par la dissipation
143      REAL dtetaecdt(ip1jmp1,llm)
144      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
145      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
146      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
147      CHARACTER*15 ztit
148      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
149      SAVE      ip_ebil_dyn
150      DATA      ip_ebil_dyn/0/
151c-jld
152
153      character*80 dynhist_file, dynhistave_file
154      character*20 modname
155      character*80 abort_message
156
157C Calendrier
158      LOGICAL true_calendar
159      PARAMETER (true_calendar = .false.)
160
161      logical dissip_conservative
162      save dissip_conservative
163      data dissip_conservative/.true./
164
165      LOGICAL prem
166      save prem
167      DATA prem/.true./
168      INTEGER testita
169      PARAMETER (testita = 9)
170
171      itaufin   = nday*day_step
172      itaufinp1 = itaufin +1
173
174
175      itau = 0
176      iday = day_ini+itau/day_step
177      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
178         IF(time.GT.1.) THEN
179          time = time-1.
180          iday = iday+1
181         ENDIF
182
183
184c-----------------------------------------------------------------------
185c   On initialise la pression et la fonction d'Exner :
186c   --------------------------------------------------
187
188      dq=0.
189      CALL pression ( ip1jmp1, ap, bp, ps, p       )
190      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
191
192c-----------------------------------------------------------------------
193c   Debut de l'integration temporelle:
194c   ----------------------------------
195
196   1  CONTINUE
197
198
199#ifdef CPP_IOIPSL
200      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
201        call guide(itau,ucov,vcov,teta,q,masse,ps)
202      else
203        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
204     .    '6 dernieres heures'
205      endif
206#endif
207c
208c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
209c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
210c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
211c     ENDIF
212c
213      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
214      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
215      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
216      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
217      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
218
219      forward = .TRUE.
220      leapf   = .FALSE.
221      dt      =  dtvr
222
223c   ...    P.Le Van .26/04/94  ....
224
225      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
226      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
227
228      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
229
230   2  CONTINUE
231
232c-----------------------------------------------------------------------
233
234c   date:
235c   -----
236
237
238c   gestion des appels de la physique et des dissipations:
239c   ------------------------------------------------------
240c
241c   ...    P.Le Van  ( 6/02/95 )  ....
242
243      apphys = .FALSE.
244      statcl = .FALSE.
245      conser = .FALSE.
246      apdiss = .FALSE.
247
248      IF( purmats ) THEN
249         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
250         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
251         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
252     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
253      ELSE
254         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
255         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
256         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
257      END IF
258
259c-----------------------------------------------------------------------
260c   calcul des tendances dynamiques:
261c   --------------------------------
262
263      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
264
265      CALL caldyn
266     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
267     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
268
269c-----------------------------------------------------------------------
270c   calcul des tendances advection des traceurs (dont l'humidite)
271c   -------------------------------------------------------------
272
273      IF( forward. OR . leapf )  THEN
274
275c
276#ifdef INCA
277             CALL caladvtrac(q,pbaru,pbarv,
278     *                      p, masse, dq,  teta,
279     .             flxw,
280     .             pk,
281     .             mmt_adj,
282     .             hadv_flg)
283#else
284             CALL caladvtrac(q,pbaru,pbarv,
285     *                      p, masse, dq,  teta,
286     .             pk)
287#endif
288
289         IF (offline) THEN
290Cmaf stokage du flux de masse pour traceurs OFF-LINE
291
292#ifdef CPP_IOIPSL
293           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
294     .   dtvr, itau)
295#endif
296
297
298         ENDIF
299c
300      ENDIF
301
302
303c-----------------------------------------------------------------------
304c   integrations dynamique et traceurs:
305c   ----------------------------------
306
307
308       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
309     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
310     $              finvmaold                                    )
311
312
313c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
314c
315c-----------------------------------------------------------------------
316c   calcul des tendances physiques:
317c   -------------------------------
318c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
319c
320       IF( purmats )  THEN
321          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
322       ELSE
323          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
324       ENDIF
325c
326c
327       IF( apphys )  THEN
328c
329c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
330c
331
332         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
333         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
334
335           rdaym_ini  = itau * dtvr / daysec
336           rdayvrai   = rdaym_ini  + day_ini
337
338
339c rajout debug
340c       lafin = .true.
341
342
343c   Inbterface avec les routines de phylmd (phymars ... )
344c   -----------------------------------------------------
345
346#ifdef CPP_PHYS
347c+jld
348
349c  Diagnostique de conservation de l'énergie : initialisation
350      IF (ip_ebil_dyn.ge.1 ) THEN
351          ztit='bil dyn'
352          CALL diagedyn(ztit,2,1,1,dtphys
353     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
354      ENDIF
355c-jld
356
357        CALL calfis( nq, lafin ,rdayvrai,time  ,
358     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
359     $               du,dv,dteta,dq,w,
360#ifdef INCA
361     $               flxw,
362#endif
363     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
364
365c      ajout des tendances physiques:
366c      ------------------------------
367          CALL addfi( nqmx, dtphys, leapf, forward   ,
368     $                  ucov, vcov, teta , q   ,ps ,
369     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
370c
371c  Diagnostique de conservation de l'énergie : difference
372      IF (ip_ebil_dyn.ge.1 ) THEN
373          ztit='bil phys'
374          CALL diagedyn(ztit,2,1,1,dtphys
375     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
376      ENDIF
377#else
378
379c   Calcul academique de la physique = Rappel Newtonien + fritcion
380c   --------------------------------------------------------------
381       teta(:,:)=teta(:,:)
382     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
383       call friction(ucov,vcov,iphysiq*dtvr)
384
385#endif
386
387c-jld
388       ENDIF
389
390        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
391        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
392
393
394c-----------------------------------------------------------------------
395c   dissipation horizontale et verticale  des petites echelles:
396c   ----------------------------------------------------------
397
398      IF(apdiss) THEN
399
400
401c   calcul de l'energie cinetique avant dissipation
402        call covcont(llm,ucov,vcov,ucont,vcont)
403        call enercin(vcov,ucov,vcont,ucont,ecin0)
404
405c   dissipation
406        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
407        ucov=ucov+dudis
408        vcov=vcov+dvdis
409c       teta=teta+dtetadis
410
411
412c------------------------------------------------------------------------
413        if (dissip_conservative) then
414C       On rajoute la tendance due a la transform. Ec -> E therm. cree
415C       lors de la dissipation
416            call covcont(llm,ucov,vcov,ucont,vcont)
417            call enercin(vcov,ucov,vcont,ucont,ecin)
418            dtetaecdt= (ecin0-ecin)/ pk
419c           teta=teta+dtetaecdt
420            dtetadis=dtetadis+dtetaecdt
421        endif
422        teta=teta+dtetadis
423c------------------------------------------------------------------------
424
425
426c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
427c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
428c
429
430        DO l  =  1, llm
431          DO ij =  1,iim
432           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
433           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
434          ENDDO
435           tpn  = SSUM(iim,tppn,1)/apoln
436           tps  = SSUM(iim,tpps,1)/apols
437
438          DO ij = 1, iip1
439           teta(  ij    ,l) = tpn
440           teta(ij+ip1jm,l) = tps
441          ENDDO
442        ENDDO
443
444        DO ij =  1,iim
445          tppn(ij)  = aire(  ij    ) * ps (  ij    )
446          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
447        ENDDO
448          tpn  = SSUM(iim,tppn,1)/apoln
449          tps  = SSUM(iim,tpps,1)/apols
450
451        DO ij = 1, iip1
452          ps(  ij    ) = tpn
453          ps(ij+ip1jm) = tps
454        ENDDO
455
456
457      END IF
458
459c ajout debug
460c              IF( lafin ) then 
461c                abort_message = 'Simulation finished'
462c                call abort_gcm(modname,abort_message,0)
463c              ENDIF
464       
465c   ********************************************************************
466c   ********************************************************************
467c   .... fin de l'integration dynamique  et physique pour le pas itau ..
468c   ********************************************************************
469c   ********************************************************************
470
471c   preparation du pas d'integration suivant  ......
472
473      IF ( .NOT.purmats ) THEN
474c       ........................................................
475c       ..............  schema matsuno + leapfrog  ..............
476c       ........................................................
477
478            IF(forward. OR. leapf) THEN
479              itau= itau + 1
480              iday= day_ini+itau/day_step
481              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
482                IF(time.GT.1.) THEN
483                  time = time-1.
484                  iday = iday+1
485                ENDIF
486            ENDIF
487
488
489            IF( itau. EQ. itaufinp1 ) then 
490c$$$       write(79,*) 'ucov',ucov
491c$$$       write(80,*) 'vcov',vcov
492c$$$       write(81,*) 'teta',teta
493c$$$       write(82,*) 'ps',ps
494c$$$       write(83,*) 'q',q
495c$$$       WRITE(85,*) 'q1 = ',q(:,:,1)
496c$$$       WRITE(86,*) 'q3 = ',q(:,:,3)
497
498              abort_message = 'Simulation finished'
499
500              call abort_gcm(modname,abort_message,0)
501            ENDIF
502c-----------------------------------------------------------------------
503c   ecriture du fichier histoire moyenne:
504c   -------------------------------------
505
506            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
507               IF(itau.EQ.itaufin) THEN
508                  iav=1
509               ELSE
510                  iav=0
511               ENDIF
512#ifdef CPP_IOIPSL
513              CALL writedynav(histaveid, nqmx, itau,vcov ,
514     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
515               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
516     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
517#endif
518
519            ENDIF
520
521c-----------------------------------------------------------------------
522c   ecriture de la bande histoire:
523c   ------------------------------
524
525            IF( MOD(itau,iecri         ).EQ.0) THEN
526c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
527
528               nbetat = nbetatdem
529       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
530        unat=0.
531        do l=1,llm
532           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
533           vnat(:,l)=vcov(:,l)/cv(:)
534        enddo
535#ifdef CPP_IOIPSL
536c        CALL writehist(histid,histvid, nqmx,itau,vcov,
537c     s                       ucov,teta,phi,q,masse,ps,phis)
538#else
539#include "write_grads_dyn.h"
540#endif
541
542
543            ENDIF
544
545            IF(itau.EQ.itaufin) THEN
546
547
548#ifdef CPP_IOIPSL
549       CALL dynredem1("restart.nc",0.0,
550     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
551#endif
552
553              CLOSE(99)
554            ENDIF
555
556c-----------------------------------------------------------------------
557c   gestion de l'integration temporelle:
558c   ------------------------------------
559
560            IF( MOD(itau,iperiod).EQ.0 )    THEN
561                    GO TO 1
562            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
563
564                   IF( forward )  THEN
565c      fin du pas forward et debut du pas backward
566
567                      forward = .FALSE.
568                        leapf = .FALSE.
569                           GO TO 2
570
571                   ELSE
572c      fin du pas backward et debut du premier pas leapfrog
573
574                        leapf =  .TRUE.
575                        dt  =  2.*dtvr
576                        GO TO 2
577                   END IF
578            ELSE
579
580c      ......   pas leapfrog  .....
581
582                 leapf = .TRUE.
583                 dt  = 2.*dtvr
584                 GO TO 2
585            END IF
586
587      ELSE
588
589c       ........................................................
590c       ..............       schema  matsuno        ...............
591c       ........................................................
592            IF( forward )  THEN
593
594             itau =  itau + 1
595             iday = day_ini+itau/day_step
596             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
597
598                  IF(time.GT.1.) THEN
599                   time = time-1.
600                   iday = iday+1
601                  ENDIF
602
603               forward =  .FALSE.
604               IF( itau. EQ. itaufinp1 ) then 
605                 abort_message = 'Simulation finished'
606                 call abort_gcm(modname,abort_message,0)
607               ENDIF
608               GO TO 2
609
610            ELSE
611
612            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
613               IF(itau.EQ.itaufin) THEN
614                  iav=1
615               ELSE
616                  iav=0
617               ENDIF
618#ifdef CPP_IOIPSL
619              CALL writedynav(histaveid, nqmx, itau,vcov ,
620     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
621               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
622     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
623#endif
624
625            ENDIF
626
627               IF(MOD(itau,iecri         ).EQ.0) THEN
628c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
629                  nbetat = nbetatdem
630       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
631        unat=0.
632        do l=1,llm
633           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
634           vnat(:,l)=vcov(:,l)/cv(:)
635        enddo
636#ifdef CPP_IOIPSL
637c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
638c    ,                           ucov,teta,phi,q,masse,ps,phis)
639#else
640#include "write_grads_dyn.h"
641#endif
642
643
644               ENDIF
645
646#ifdef CPP_IOIPSL
647                 IF(itau.EQ.itaufin)
648     . CALL dynredem1("restart.nc",0.0,
649     .                     vcov,ucov,teta,q,nqmx,masse,ps)
650#endif
651
652                 forward = .TRUE.
653                 GO TO  1
654
655            ENDIF
656
657      END IF
658
659      STOP
660      END
Note: See TracBrowser for help on using the repository browser.