source: LMDZ.3.3/trunk/libf/dyn3d/gcmtest.F @ 2429

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