source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/dyn3d/leapfrog.F @ 5373

Last change on this file since 5373 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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