source: trunk/LMDZ.MARS/libf/dyn3d/gcm.F @ 798

Last change on this file since 798 was 798, checked in by tnavarro, 12 years ago

modification that allows to have "1+1=2" with dynamical core only.

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