source: LMDZ4/trunk/libf/dyn3d/leapfrog.F @ 541

Last change on this file since 541 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

  • le offline
  • les thermiques
  • mellor & yamada dans la couche limite

LF

  • 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      character*80 dynhist_file, dynhistave_file
152      character*20 modname
153      character*80 abort_message
154
155C Calendrier
156      LOGICAL true_calendar
157      PARAMETER (true_calendar = .false.)
158
159      logical dissip_conservative
160      save dissip_conservative
161      data dissip_conservative/.true./
162
163      LOGICAL prem
164      save prem
165      DATA prem/.true./
166      INTEGER testita
167      PARAMETER (testita = 9)
168
169      itaufin   = nday*day_step
170      itaufinp1 = itaufin +1
171
172
173      itau = 0
174      iday = day_ini+itau/day_step
175      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
176         IF(time.GT.1.) THEN
177          time = time-1.
178          iday = iday+1
179         ENDIF
180
181
182c-----------------------------------------------------------------------
183c   On initialise la pression et la fonction d'Exner :
184c   --------------------------------------------------
185
186      dq=0.
187      CALL pression ( ip1jmp1, ap, bp, ps, p       )
188      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
189
190c-----------------------------------------------------------------------
191c   Debut de l'integration temporelle:
192c   ----------------------------------
193
194   1  CONTINUE
195
196
197#ifdef CPP_IOIPSL
198      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
199        call guide(itau,ucov,vcov,teta,q,masse,ps)
200      else
201        print*,'attention on ne guide pas les 6 dernieres heures'
202      endif
203#endif
204c
205c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
206c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
207c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
208c     ENDIF
209c
210      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
211      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
212      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
213      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
214      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
215
216      forward = .TRUE.
217      leapf   = .FALSE.
218      dt      =  dtvr
219
220c   ...    P.Le Van .26/04/94  ....
221
222      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
223      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
224
225      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
226
227   2  CONTINUE
228
229c-----------------------------------------------------------------------
230
231c   date:
232c   -----
233
234
235c   gestion des appels de la physique et des dissipations:
236c   ------------------------------------------------------
237c
238c   ...    P.Le Van  ( 6/02/95 )  ....
239
240      apphys = .FALSE.
241      statcl = .FALSE.
242      conser = .FALSE.
243      apdiss = .FALSE.
244
245      IF( purmats ) THEN
246         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
247         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
248         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
249     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
250      ELSE
251         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
252         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
253         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
254      END IF
255
256c-----------------------------------------------------------------------
257c   calcul des tendances dynamiques:
258c   --------------------------------
259
260      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
261
262      CALL caldyn
263     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
264     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
265
266c-----------------------------------------------------------------------
267c   calcul des tendances advection des traceurs (dont l'humidite)
268c   -------------------------------------------------------------
269
270      IF( forward. OR . leapf )  THEN
271
272c
273#ifdef INCA_CH4
274             CALL caladvtrac(q,pbaru,pbarv,
275     *                      p, masse, dq,  teta,
276     .             flxw,
277     .             pk,
278     .             mmt_adj,
279     .             hadv_flg)
280#else
281             CALL caladvtrac(q,pbaru,pbarv,
282     *                      p, masse, dq,  teta,
283     .             pk)
284#endif
285
286         IF (offline) THEN
287Cmaf stokage du flux de masse pour traceurs OFF-LINE
288
289#ifdef CPP_IOIPSL
290           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
291     .   dtvr, itau)
292#endif
293
294
295         ENDIF
296c
297      ENDIF
298
299
300c-----------------------------------------------------------------------
301c   integrations dynamique et traceurs:
302c   ----------------------------------
303
304
305       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
306     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
307     $              finvmaold                                    )
308
309
310c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
311c
312c-----------------------------------------------------------------------
313c   calcul des tendances physiques:
314c   -------------------------------
315c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
316c
317       IF( purmats )  THEN
318          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
319       ELSE
320          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
321       ENDIF
322c
323c
324       IF( apphys )  THEN
325c
326c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
327c
328
329         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
330         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
331
332           rdaym_ini  = itau * dtvr / daysec
333           rdayvrai   = rdaym_ini  + day_ini
334
335
336c rajout debug
337c       lafin = .true.
338
339
340c   Inbterface avec les routines de phylmd (phymars ... )
341c   -----------------------------------------------------
342
343#ifdef CPP_PHYS
344c+jld
345
346c  Diagnostique de conservation de l'énergie : initialisation
347      IF (ip_ebil_dyn.ge.1 ) THEN
348          ztit='bil dyn'
349          CALL diagedyn(ztit,2,1,1,dtphys
350     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
351      ENDIF
352c-jld
353
354        CALL calfis( nq, lafin ,rdayvrai,time  ,
355     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
356     $               du,dv,dteta,dq,w,
357#ifdef INCA_CH4
358     $               flxw,
359#endif
360     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
361
362c      ajout des tendances physiques:
363c      ------------------------------
364          CALL addfi( nqmx, dtphys, leapf, forward   ,
365     $                  ucov, vcov, teta , q   ,ps ,
366     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
367c
368c  Diagnostique de conservation de l'énergie : difference
369      IF (ip_ebil_dyn.ge.1 ) THEN
370          ztit='bil phys'
371          CALL diagedyn(ztit,2,1,1,dtphys
372     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
373      ENDIF
374#else
375
376c   Calcul academique de la physique = Rappel Newtonien + fritcion
377c   --------------------------------------------------------------
378       teta(:,:)=teta(:,:)
379     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
380       call friction(ucov,vcov,iphysiq*dtvr)
381
382#endif
383
384c-jld
385       ENDIF
386
387        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
388        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
389
390
391c-----------------------------------------------------------------------
392c   dissipation horizontale et verticale  des petites echelles:
393c   ----------------------------------------------------------
394
395      IF(apdiss) THEN
396
397
398c   calcul de l'energie cinetique avant dissipation
399        call covcont(llm,ucov,vcov,ucont,vcont)
400        call enercin(vcov,ucov,vcont,ucont,ecin0)
401
402c   dissipation
403        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
404        ucov=ucov+dudis
405        vcov=vcov+dvdis
406c       teta=teta+dtetadis
407
408
409c------------------------------------------------------------------------
410        if (dissip_conservative) then
411C       On rajoute la tendance due a la transform. Ec -> E therm. cree
412C       lors de la dissipation
413            call covcont(llm,ucov,vcov,ucont,vcont)
414            call enercin(vcov,ucov,vcont,ucont,ecin)
415            dtetaecdt= (ecin0-ecin)/ pk
416c           teta=teta+dtetaecdt
417            dtetadis=dtetadis+dtetaecdt
418        endif
419        teta=teta+dtetadis
420c------------------------------------------------------------------------
421
422
423c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
424c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
425c
426
427        DO l  =  1, llm
428          DO ij =  1,iim
429           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
430           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
431          ENDDO
432           tpn  = SSUM(iim,tppn,1)/apoln
433           tps  = SSUM(iim,tpps,1)/apols
434
435          DO ij = 1, iip1
436           teta(  ij    ,l) = tpn
437           teta(ij+ip1jm,l) = tps
438          ENDDO
439        ENDDO
440
441        DO ij =  1,iim
442          tppn(ij)  = aire(  ij    ) * ps (  ij    )
443          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
444        ENDDO
445          tpn  = SSUM(iim,tppn,1)/apoln
446          tps  = SSUM(iim,tpps,1)/apols
447
448        DO ij = 1, iip1
449          ps(  ij    ) = tpn
450          ps(ij+ip1jm) = tps
451        ENDDO
452
453
454      END IF
455
456c ajout debug
457c              IF( lafin ) then 
458c                abort_message = 'Simulation finished'
459c                call abort_gcm(modname,abort_message,0)
460c              ENDIF
461       
462c   ********************************************************************
463c   ********************************************************************
464c   .... fin de l'integration dynamique  et physique pour le pas itau ..
465c   ********************************************************************
466c   ********************************************************************
467
468c   preparation du pas d'integration suivant  ......
469
470      IF ( .NOT.purmats ) THEN
471c       ........................................................
472c       ..............  schema matsuno + leapfrog  ..............
473c       ........................................................
474
475            IF(forward. OR. leapf) THEN
476              itau= itau + 1
477              iday= day_ini+itau/day_step
478              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
479                IF(time.GT.1.) THEN
480                  time = time-1.
481                  iday = iday+1
482                ENDIF
483            ENDIF
484
485
486            IF( itau. EQ. itaufinp1 ) then 
487       write(79,*) 'ucov',ucov
488       write(80,*) 'vcov',vcov
489       write(81,*) 'teta',teta
490       write(82,*) 'ps',ps
491       write(83,*) 'q',q
492       WRITE(85,*) 'q1 = ',q(:,:,1)
493       WRITE(86,*) 'q3 = ',q(:,:,3)
494
495              abort_message = 'Simulation finished'
496
497              call abort_gcm(modname,abort_message,0)
498            ENDIF
499c-----------------------------------------------------------------------
500c   ecriture du fichier histoire moyenne:
501c   -------------------------------------
502
503            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
504               IF(itau.EQ.itaufin) THEN
505                  iav=1
506               ELSE
507                  iav=0
508               ENDIF
509#ifdef CPP_IOIPSL
510              CALL writedynav(histaveid, nqmx, itau,vcov ,
511     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
512               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
513     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
514#endif
515
516            ENDIF
517
518c-----------------------------------------------------------------------
519c   ecriture de la bande histoire:
520c   ------------------------------
521
522            IF( MOD(itau,iecri         ).EQ.0) THEN
523c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
524
525               nbetat = nbetatdem
526       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
527        unat=0.
528        do l=1,llm
529           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
530           vnat(:,l)=vcov(:,l)/cv(:)
531        enddo
532#ifdef CPP_IOIPSL
533c        CALL writehist(histid,histvid, nqmx,itau,vcov,
534c     s                       ucov,teta,phi,q,masse,ps,phis)
535#else
536#include "write_grads_dyn.h"
537#endif
538
539
540            ENDIF
541
542            IF(itau.EQ.itaufin) THEN
543
544
545#ifdef CPP_IOIPSL
546       CALL dynredem1("restart.nc",0.0,
547     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
548#endif
549
550              CLOSE(99)
551            ENDIF
552
553c-----------------------------------------------------------------------
554c   gestion de l'integration temporelle:
555c   ------------------------------------
556
557            IF( MOD(itau,iperiod).EQ.0 )    THEN
558                    GO TO 1
559            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
560
561                   IF( forward )  THEN
562c      fin du pas forward et debut du pas backward
563
564                      forward = .FALSE.
565                        leapf = .FALSE.
566                           GO TO 2
567
568                   ELSE
569c      fin du pas backward et debut du premier pas leapfrog
570
571                        leapf =  .TRUE.
572                        dt  =  2.*dtvr
573                        GO TO 2
574                   END IF
575            ELSE
576
577c      ......   pas leapfrog  .....
578
579                 leapf = .TRUE.
580                 dt  = 2.*dtvr
581                 GO TO 2
582            END IF
583
584      ELSE
585
586c       ........................................................
587c       ..............       schema  matsuno        ...............
588c       ........................................................
589            IF( forward )  THEN
590
591             itau =  itau + 1
592             iday = day_ini+itau/day_step
593             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
594
595                  IF(time.GT.1.) THEN
596                   time = time-1.
597                   iday = iday+1
598                  ENDIF
599
600               forward =  .FALSE.
601               IF( itau. EQ. itaufinp1 ) then 
602                 abort_message = 'Simulation finished'
603                 call abort_gcm(modname,abort_message,0)
604               ENDIF
605               GO TO 2
606
607            ELSE
608
609            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
610               IF(itau.EQ.itaufin) THEN
611                  iav=1
612               ELSE
613                  iav=0
614               ENDIF
615#ifdef CPP_IOIPSL
616              CALL writedynav(histaveid, nqmx, itau,vcov ,
617     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
618               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
619     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
620#endif
621
622            ENDIF
623
624               IF(MOD(itau,iecri         ).EQ.0) THEN
625c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
626                  nbetat = nbetatdem
627       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
628        unat=0.
629        do l=1,llm
630           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
631           vnat(:,l)=vcov(:,l)/cv(:)
632        enddo
633#ifdef CPP_IOIPSL
634c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
635c    ,                           ucov,teta,phi,q,masse,ps,phis)
636#else
637#include "write_grads_dyn.h"
638#endif
639
640
641               ENDIF
642
643#ifdef CPP_IOIPSL
644                 IF(itau.EQ.itaufin)
645     . CALL dynredem1("restart.nc",0.0,
646     .                     vcov,ucov,teta,q,nqmx,masse,ps)
647#endif
648
649                 forward = .TRUE.
650                 GO TO  1
651
652            ENDIF
653
654      END IF
655
656      STOP
657      END
Note: See TracBrowser for help on using the repository browser.