source: LMDZ.3.3/branches/LF/libf/dyn3d/gcm.F

Last change on this file was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.1 KB
Line 
1
2      PROGRAM gcm
3
4      IMPLICIT NONE
5
6c      ......   Version  du 10/01/98    ..........
7
8c             avec  coordonnees  verticales hybrides
9c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
10
11c=======================================================================
12c
13c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
14c   -------
15c
16c   Objet:
17c   ------
18c
19c   GCM LMD nouvelle grille
20c
21c=======================================================================
22
23c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv
24c        et possibilite d'appeler une fonction f(y)  a derivee tangente
25c        hyperbolique a la  place de la fonction a derivee sinusoidale.         
26
27c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
28c         q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
29c
30c-----------------------------------------------------------------------
31c   Declarations:
32c   -------------
33
34#include "dimensions.h"
35#include "paramet.h"
36#include "comconst.h"
37#include "comdissnew.h"
38#include "comvert.h"
39#include "comgeom.h"
40#include "logic.h"
41#include "temps.h"
42#include "control.h"
43#include "ener.h"
44#include "netcdf.inc"
45#include "description.h"
46#include "serre.h"
47#include "tracstoke.h"
48
49
50      INTEGER         longcles
51      PARAMETER     ( longcles = 20 )
52      REAL  clesphy0( longcles )
53      SAVE  clesphy0
54
55      INTEGER*4  iday ! jour julien
56      REAL       time ! Heure de la journee en fraction d'1 jour
57      REAL zdtvr
58      INTEGER nbetatmoy, nbetatdem,nbetat
59
60c   variables dynamiques
61      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
62      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
63      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
64      REAL ps(ip1jmp1)                       ! pression  au sol
65      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
66      REAL pks(ip1jmp1)                      ! exner au  sol
67      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
68      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
69      REAL masse(ip1jmp1,llm)                ! masse d'air
70      REAL phis(ip1jmp1)                     ! geopotentiel au sol
71      REAL phi(ip1jmp1,llm)                  ! geopotentiel
72      REAL w(ip1jmp1,llm)                    ! vitesse verticale
73
74c variables dynamiques intermediaire pour le transport
75      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
76
77c   variables dynamiques au pas -1
78      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
79      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
80      REAL massem1(ip1jmp1,llm)
81
82c   tendances dynamiques
83      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
84      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
85
86c   tendances de la dissipation
87      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
88      REAL dhdis(ip1jmp1,llm)
89
90c   tendances physiques
91      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
92      REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
93
94c   variables pour le fichier histoire
95      REAL dtav      ! intervalle de temps elementaire
96
97      REAL tppn(iim),tpps(iim),tpn,tps
98c
99      INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
100
101      INTEGER itau,itaufinp1,iav
102
103      EXTERNAL caldyn, traceur
104      EXTERNAL dissip,geopot,iniconst,inifilr
105      EXTERNAL integrd,SCOPY
106      EXTERNAL iniav,writeav,writehist
107      EXTERNAL inigeom
108      EXTERNAL exner_hyb,addit
109      EXTERNAL defrun_new, test_period
110      REAL  SSUM
111      REAL time_0 , finvmaold(ip1jmp1,llm)
112
113      LOGICAL lafin
114      INTEGER ij,iq,l,numvanle,iapp_tracvl
115
116
117      integer histid, histvid, histaveid
118      real time_step, t_wrt, t_ops
119
120      REAL rdayvrai,rdaym_ini,rday_ecri
121      LOGICAL first
122      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
123
124      LOGICAL offline  ! Controle du stockage ds "fluxmass"
125      PARAMETER (offline=.false.)
126
127      character*80 dynhist_file, dynhistave_file
128      character*20 modname
129      character*80 abort_message
130
131c-----------------------------------------------------------------------
132c   Initialisations:
133c   ----------------
134
135      modname = 'gcm'
136      descript = 'Run GCM LMDZ'
137      lafin    = .FALSE.
138      dynhist_file = 'dyn_hist.nc'
139      dynhistave_file = 'dyn_hist_ave.nc'
140c-----------------------------------------------------------------------
141c
142c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
143c  ...................................................................
144c
145c     iadv = 1    shema  transport type "humidite specifique LMD" 
146c     iadv = 2    shema   amont
147c     iadv = 3    shema  Van-leer
148c
149c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
150c                                , iq = 2  pour l'eau liquide
151c         et eventuellement      , iq = 3, nqmx pour les autres traceurs
152c
153c     iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
154c
155      DO iq = 1, nqmx
156       iadv( iq ) = 3
157      ENDDO
158c
159      DO  iq = 1, nqmx
160       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
161     * ,' pour le traceur no ', iq
162       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
163     * ,' traceur no ', iq
164       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
165     * ,'le traceur no ', iq
166       IF( iadv(iq).LE.0.OR.iadv(iq).GT.3 )   THEN
167       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
168     * . '
169         STOP
170       ENDIF
171      ENDDO
172c
173      first = .TRUE.
174      numvanle = nqmx + 1
175      DO  iq = 1, nqmx
176        IF( iadv(iq).EQ.3.AND.first ) THEN
177          numvanle = iq
178          first    = .FALSE.
179        ENDIF
180      ENDDO
181c
182      DO  iq = 1, nqmx
183        IF( iadv(iq).NE.3.AND.iq.GT.numvanle )   THEN
184          PRINT *,' Il y a discontinuite dans le choix du shema de ',
185     *    'Van-leer pour les traceurs . Corriger et repasser . '
186           STOP
187        ENDIF
188        IF( iadv(iq).LT.1.OR.iadv(iq).GT.3 )    THEN
189          PRINT *,' Le choix de  iadv  est errone pour le traceur ',
190     *    iq
191         STOP
192        ENDIF
193      ENDDO
194c
195
196      CALL dynetat0("start.nc",nqmx,vcov,ucov,
197     .              teta,q,masse,ps,phis, time_0)
198
199
200c     OPEN( 99,file ='run.def',status='old',form='formatted')
201      CALL defrun_new( 99, .TRUE. , clesphy0 )
202c     CLOSE(99)
203
204c  on recalcule eventuellement le pas de temps
205
206      IF(MOD(day_step,iperiod).NE.0)
207     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
208
209      IF(MOD(day_step,iphysiq).NE.0)
210     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
211
212      zdtvr    = daysec/FLOAT(day_step)
213        IF(dtvr.NE.zdtvr) THEN
214         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
215        ENDIF
216
217c  nombre d'etats dans les fichiers demarrage et histoire
218      nbetatdem = nday / iecri
219      nbetatmoy = nday / periodav + 1
220
221      dtvr = zdtvr
222      CALL iniconst
223      CALL inigeom
224
225      CALL inifilr
226c
227c   ......    P.Le Van    ( modif  le 29/04/97 )   ......... 
228c
229      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
230     *                tetagdiv, tetagrot , tetatemp              )
231c
232
233      CALL pression ( ip1jmp1, ap, bp, ps, p       )
234      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
235c
236
237c  numero de stockage pour les fichiers de redemarrage:
238
239c-----------------------------------------------------------------------
240c   temps de depart et de fin:
241c   --------------------------
242
243      itau = 0
244      iday = day_ini+itau/day_step
245      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
246         IF(time.GT.1.) THEN
247          time = time-1.
248          iday = iday+1
249         ENDIF
250      itaufin   = nday*day_step
251      itaufinp1 = itaufin +1
252
253      day_end = day_ini + nday
254      PRINT 300, itau,itaufin,day_ini,day_end
255
256      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
257
258      ecripar = .TRUE.
259
260      time_step = zdtvr
261      t_ops = iecri * daysec
262      t_wrt = iecri * daysec
263      CALL inithist(dynhist_file,day_ini,anne_ini,time_step,
264     .              t_ops, t_wrt, nqmx, histid, histvid)
265
266      t_ops = iperiod * time_step
267      t_wrt = periodav * daysec
268      CALL initdynav(dynhistave_file,day_ini,anne_ini,time_step,
269     .              t_ops, t_wrt, nqmx, histaveid)
270
271      dtav = iperiod*dtvr/daysec
272
273
274c   Quelques initialisations pour les traceurs
275      call initial0(ijp1llm*nqmx,dq)
276      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
277      istphy=istdyn/iphysiq
278
279c-----------------------------------------------------------------------
280c   Debut de l'integration temporelle:
281c   ----------------------------------
282
283   1  CONTINUE
284c     call nudge(itau,ucov,vcov,teta,masse,ps)
285c
286c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
287        CALL  test_period ( ucov,vcov,teta,q,p,phis )
288        PRINT *,' ----   Test_period apres continue   OK ! -----', itau
289c     ENDIF
290c
291      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
292      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
293      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
294      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
295      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
296
297      forward = .TRUE.
298      leapf   = .FALSE.
299      dt      =  dtvr
300
301c   ...    P.Le Van .26/04/94  ....
302
303      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
304      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
305
306
307   2  CONTINUE
308
309c-----------------------------------------------------------------------
310
311c   date:
312c   -----
313
314
315c   gestion des appels de la physique et des dissipations:
316c   ------------------------------------------------------
317c
318c   ...    P.Le Van  ( 6/02/95 )  ....
319
320      apphys = .FALSE.
321      statcl = .FALSE.
322      conser = .FALSE.
323      apdiss = .FALSE.
324
325      IF( purmats ) THEN
326         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
327         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
328         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
329     $                              .AND.   physic    ) apphys = .TRUE.
330      ELSE
331         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
332         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
333         IF( MOD(itau+1,iphysiq).EQ.0. AND. physic    ) apphys = .TRUE.
334      END IF
335
336c-----------------------------------------------------------------------
337c   calcul des tendances dynamiques:
338c   --------------------------------
339
340      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
341c
342c
343      CALL caldyn
344     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
345     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
346
347c-----------------------------------------------------------------------
348c   calcul des tendances advection des traceurs (dont l'humidite)
349c   -------------------------------------------------------------
350
351      IF( forward. OR . leapf )  THEN
352c
353        DO 15  iq = 1, nqmx
354c
355         IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
356           CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
357
358         ELSE IF( iq.EQ. nqmx )   THEN
359c
360            iapp_tracvl = 5
361c
362cccc     iapp_tracvl est la frequence en pas du groupement des flux
363cccc      de masse pour  Van-Leer dans la routine  tracvl  .
364c
365            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
366     *                             p,  masse, dq                  )
367c
368         ENDIF
369c
370 15     CONTINUE
371C
372         IF (offline) THEN
373Cmaf stokage du flux de masse pour traceurs OFF-LINE
374           CALL fluxstoke(pbaru,pbarv,masse,teta,phi,phis)
375         ENDIF
376c
377      ENDIF
378
379
380c-----------------------------------------------------------------------
381c   integrations dynamique et traceurs:
382c   ----------------------------------
383
384
385       CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
386     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
387     $              finvmaold                                    )
388
389c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
390c
391c-----------------------------------------------------------------------
392c   calcul des tendances physiques:
393c   -------------------------------
394c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
395c
396       IF( purmats )  THEN
397          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
398       ELSE
399          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
400       ENDIF
401c
402c
403       IF( apphys )  THEN
404c
405c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
406c
407
408         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
409         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
410
411           rdaym_ini  = itau * dtvr / daysec
412           rdayvrai   = rdaym_ini  + day_ini
413
414           IF ( ecritphy.LT.1. )  THEN
415             rday_ecri = rdaym_ini
416           ELSE
417             rday_ecri = NINT( rdayvrai )
418           ENDIF
419c
420        CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
421     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
422     $     du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi  )
423
424c      ajout des tendances physiques:
425c      ------------------------------
426          CALL addfi( nqmx, dtphys, leapf, forward   ,
427     $                  ucov, vcov, teta , q   ,ps ,
428     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
429c
430       ENDIF
431
432        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
433        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
434
435c-----------------------------------------------------------------------
436c
437c
438c   dissipation horizontale et verticale  des petites echelles:
439c   ----------------------------------------------------------
440
441      IF(apdiss) THEN
442
443         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
444
445         CALL addit( ijp1llm,ucov ,dudis,ucov )
446         CALL addit( ijmllm ,vcov ,dvdis,vcov )
447         CALL addit( ijp1llm,teta ,dhdis,teta )
448
449
450c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
451c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
452c
453      IF( alphax.NE. 0. ) THEN
454
455        DO l  =  1, llm
456          DO ij =  1,iim
457           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
458           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
459          ENDDO
460           tpn  = SSUM(iim,tppn,1)/apoln
461           tps  = SSUM(iim,tpps,1)/apols
462
463          DO ij = 1, iip1
464           teta(  ij    ,l) = tpn
465           teta(ij+ip1jm,l) = tps
466          ENDDO
467        ENDDO
468
469      ENDIF
470
471      END IF
472       
473c   ********************************************************************
474c   ********************************************************************
475c   .... fin de l'integration dynamique  et physique pour le pas itau ..
476c   ********************************************************************
477c   ********************************************************************
478
479c   preparation du pas d'integration suivant  ......
480
481      IF ( .NOT.purmats ) THEN
482c       ........................................................
483c       ..............  schema matsuno + leapfrog  ..............
484c       ........................................................
485
486            IF(forward. OR. leapf) THEN
487              itau= itau + 1
488              iday= day_ini+itau/day_step
489              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
490                IF(time.GT.1.) THEN
491                  time = time-1.
492                  iday = iday+1
493                ENDIF
494            ENDIF
495
496
497            IF( itau. EQ. itaufinp1 ) then 
498              abort_message = 'Simulation finished'
499              call abort_gcm(modname,abort_message,0)
500            ENDIF
501c-----------------------------------------------------------------------
502c   ecriture du fichier histoire moyenne:
503c   -------------------------------------
504
505            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
506               IF(itau.EQ.itaufin) THEN
507                  iav=1
508               ELSE
509                  iav=0
510               ENDIF
511               CALL writedynav(histaveid, nqmx, itau,vcov ,
512     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
513            ENDIF
514
515c-----------------------------------------------------------------------
516c   ecriture de la bande histoire:
517c   ------------------------------
518
519            IF( MOD(itau,iecri*day_step).EQ.0) THEN
520
521               nbetat = nbetatdem
522       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
523       CALL writehist( histid, histvid, nqmx, itau,vcov ,
524     ,                          ucov,teta,phi,q,masse,ps,phis)
525
526
527            ENDIF
528
529            IF(itau.EQ.itaufin) THEN
530
531
532       PRINT *,' Appel test_period avant redem ', itau
533       CALL  test_period ( ucov,vcov,teta,q,p,phis )
534       CALL dynredem1("restart.nc",0.0,
535     .                     vcov,ucov,teta,q,nqmx,masse,ps)
536
537              CLOSE(99)
538            ENDIF
539
540c-----------------------------------------------------------------------
541c   gestion de l'integration temporelle:
542c   ------------------------------------
543
544            IF( MOD(itau,iperiod).EQ.0 )    THEN
545                    GO TO 1
546            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
547
548                   IF( forward )  THEN
549c      fin du pas forward et debut du pas backward
550
551                      forward = .FALSE.
552                        leapf = .FALSE.
553                           GO TO 2
554
555                   ELSE
556c      fin du pas backward et debut du premier pas leapfrog
557
558                        leapf =  .TRUE.
559                        dt  =  2.*dtvr
560                        GO TO 2
561                   END IF
562            ELSE
563
564c      ......   pas leapfrog  .....
565
566                 leapf = .TRUE.
567                 dt  = 2.*dtvr
568                 GO TO 2
569            END IF
570
571      ELSE
572
573c       ........................................................
574c       ..............       schema  matsuno        ...............
575c       ........................................................
576            IF( forward )  THEN
577
578             itau =  itau + 1
579             iday = day_ini+itau/day_step
580             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
581
582                  IF(time.GT.1.) THEN
583                   time = time-1.
584                   iday = iday+1
585                  ENDIF
586
587               forward =  .FALSE.
588               IF( itau. EQ. itaufinp1 ) then 
589                 abort_message = 'Simulation finished'
590                 call abort_gcm(modname,abort_message,0)
591               ENDIF
592               GO TO 2
593
594            ELSE
595
596            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
597               IF(itau.EQ.itaufin) THEN
598                  iav=1
599               ELSE
600                  iav=0
601               ENDIF
602               CALL writedynav(histaveid, nqmx, itau,vcov ,
603     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
604
605            ENDIF
606
607               IF(MOD(itau,iecri*day_step).EQ.0) THEN
608                  nbetat = nbetatdem
609       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
610       CALL writehist( histid, histvid, nqmx, itau,vcov ,
611     ,                          ucov,teta,phi,q,masse,ps,phis)
612               ENDIF
613
614                 IF(itau.EQ.itaufin)
615     . CALL dynredem1("restart.nc",0.0,
616     .                     vcov,ucov,teta,q,nqmx,masse,ps)
617
618                 forward = .TRUE.
619                 GO TO  1
620
621            ENDIF
622
623      END IF
624
625 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
626     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
627
628      STOP
629      END
Note: See TracBrowser for help on using the repository browser.