/[lmdze]/trunk/dyn3d/leapfrog.f
ViewVC logotype

Diff of /trunk/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 47 by guez, Fri Jul 1 15:00:48 2011 UTC
# Line 1  Line 1 
1  module leapfrog_m  module leapfrog_m
2    
   ! This module is clean: no C preprocessor directive, no include line.  
   
3    IMPLICIT NONE    IMPLICIT NONE
4    
5  contains  contains
6    
7    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, nq, q, clesphy0, &    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
        time_0)  
8    
9      ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34      ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34
10        ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11        ! Matsuno-leapfrog scheme.
12    
13        use addfi_m, only: addfi
14        use bilan_dyn_m, only: bilan_dyn
15        use caladvtrac_m, only: caladvtrac
16        use caldyn_m, only: caldyn
17        USE calfis_m, ONLY: calfis
18        USE com_io_dyn, ONLY: histaveid
19        USE comconst, ONLY: daysec, dtphys, dtvr
20        USE comgeom, ONLY: aire_2d, apoln, apols
21        USE comvert, ONLY: ap, bp
22        USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
23             periodav
24        USE dimens_m, ONLY: iim, jjm, llm, nqmx
25        use dissip_m, only: dissip
26        USE dynetat0_m, ONLY: day_ini
27        use dynredem1_m, only: dynredem1
28        USE exner_hyb_m, ONLY: exner_hyb
29        use filtreg_m, only: filtreg
30        use geopot_m, only: geopot
31        USE guide_m, ONLY: guide
32        use inidissip_m, only: idissip
33        use integrd_m, only: integrd
34        USE logic, ONLY: iflag_phys, ok_guide
35        USE paramet_m, ONLY: ip1jmp1
36        USE pressure_var, ONLY: p3d
37        USE temps, ONLY: itau_dyn
38    
39        ! Variables dynamiques:
40        REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
41        REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
42    
43      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
44      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      ! potential temperature
45    
46      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
     ! Objet:  
     ! GCM LMD nouvelle grille  
   
     ! ... Dans inigeom, nouveaux calculs pour les elongations cu, cv  
     ! et possibilite d'appeler une fonction f(y) a derivee tangente  
     ! hyperbolique a la place de la fonction a derivee sinusoidale.  
   
     ! ... Possibilite de choisir le shema pour l'advection de  
     ! q, en modifiant iadv dans "traceur.def" (10/02) .  
   
     ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99)  
     ! Pour Van-Leer iadv=10  
   
     use dimens_m, only: iim, llm, nqmx  
     use paramet_m, only: ip1jmp1, ip1jm, llmp1, ijmllm, ijp1llm, jjp1, iip1, &  
          iip2  
     use comconst, only: dtvr, daysec, dtphys  
     use comvert, only: ap, bp  
     use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &  
          offline, periodav  
     use logic, only: ok_guide, apdiss, apphys, conser, forward, iflag_phys, &  
          leapf, statcl  
     use comgeom  
     use serre  
     use temps, only: itaufin, day_ini, dt  
     use iniprint, only: prt_level  
     use com_io_dyn  
     use abort_gcm_m, only: abort_gcm  
     use ener  
     use calfis_m, only: calfis  
     use exner_hyb_m, only: exner_hyb  
     use guide_m, only: guide  
     use pression_m, only: pression  
   
     integer nq  
   
     INTEGER longcles  
     PARAMETER (longcles = 20)  
     REAL clesphy0(longcles)  
   
     ! variables dynamiques  
     REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants  
     REAL teta(ip1jmp1, llm) ! temperature potentielle  
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
     REAL ps(ip1jmp1) ! pression au sol  
     REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches  
     REAL pks(ip1jmp1) ! exner au sol  
     REAL pk(ip1jmp1, llm) ! exner au milieu des couches  
     REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches  
47      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
48      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
49      REAL phi(ip1jmp1, llm) ! geopotential  
50        REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
51        ! mass fractions of advected fields
52    
53        REAL, intent(in):: time_0
54    
55        ! Variables local to the procedure:
56    
57        ! Variables dynamiques:
58    
59        REAL pks(ip1jmp1) ! exner au sol
60        REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
61        REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
62        REAL phi(iim + 1, jjm + 1, llm) ! geopotential
63      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
64    
65      ! variables dynamiques intermediaire pour le transport      ! variables dynamiques intermediaire pour le transport
66      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse      REAL pbaru(ip1jmp1, llm), pbarv((iim + 1) * jjm, llm) !flux de masse
67    
68      ! variables dynamiques au pas - 1      ! variables dynamiques au pas - 1
69      REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)      REAL vcovm1((iim + 1) * jjm, llm), ucovm1(ip1jmp1, llm)
70      REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)      REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
71      REAL massem1(ip1jmp1, llm)      REAL massem1(ip1jmp1, llm)
72    
73      ! tendances dynamiques      ! tendances dynamiques
74      REAL dv(ip1jm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), dudyn(ip1jmp1, llm)
75      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(iim + 1, jjm + 1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
76    
77      ! tendances de la dissipation      ! tendances de la dissipation
78      REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
79      REAL dtetadis(ip1jmp1, llm)      REAL dtetadis(iim + 1, jjm + 1, llm)
80    
81      ! tendances physiques      ! tendances physiques
82      REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
83      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
84    
85      ! variables pour le fichier histoire      ! variables pour le fichier histoire
86    
87      REAL tppn(iim), tpps(iim), tpn, tps      INTEGER itau ! index of the time step of the dynamics, starts at 0
88        INTEGER itaufin
89      INTEGER itau, itaufinp1      REAL time ! time of day, as a fraction of day length
90      INTEGER iday ! jour julien      real finvmaold(ip1jmp1, llm)
91      REAL time ! Heure de la journee en fraction d'1 jour      INTEGER l
   
     REAL SSUM  
     REAL time_0, finvmaold(ip1jmp1, llm)  
   
     LOGICAL :: lafin=.false.  
     INTEGER ij, l  
   
92      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
93    
94      data callinigrads/.true./      ! Variables test conservation energie
95        REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
96    
97      !+jld variables test conservation energie      REAL dtetaecdt(iim + 1, jjm + 1, llm)
98      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      ! tendance de la température potentielle due à la tansformation
99      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! d'énergie cinétique en énergie thermique créée par la dissipation
100      ! tansformation d'energie cinetique en energie thermique  
101      ! cree par la dissipation      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
102      REAL dtetaecdt(ip1jmp1, llm)      logical leapf
103      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      real dt
     CHARACTER*15 ztit  
     INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.  
     SAVE ip_ebil_dyn  
     DATA ip_ebil_dyn /0/  
   
     character(len=*), parameter:: modname = "leapfrog"  
     character*80 abort_message  
   
     logical dissip_conservative  
     save dissip_conservative  
     data dissip_conservative /.true./  
   
     LOGICAL prem  
     save prem  
     DATA prem /.true./  
104    
105      !---------------------------------------------------      !---------------------------------------------------
106    
107      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
108    
109      itaufin = nday * day_step      itaufin = nday * day_step
110      itaufinp1 = itaufin + 1      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
111    
112      itau = 0      dq = 0.
     iday = day_ini  
     time = time_0  
     IF (time > 1.) THEN  
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
113    
114      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
115      dq=0.      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
116      CALL pression(ip1jmp1, ap, bp, ps, p)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
     CALL exner_hyb(ps, p, pks, pk, pkf)  
   
     ! Debut de l'integration temporelle:  
     do  
        if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then  
           call guide(itau, ucov, vcov, teta, q, masse, ps)  
        else  
           IF (prt_level > 9) print *, &  
                'Attention : on ne guide pas les 6 dernieres heures.'  
        endif  
   
        CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)  
        CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)  
        CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)  
        CALL SCOPY(ijp1llm, masse, 1, massem1, 1)  
        CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)  
   
        forward = .TRUE.  
        leapf = .FALSE.  
        dt = dtvr  
   
        CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)  
        CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)  
   
        do  
           ! gestion des appels de la physique et des dissipations:  
   
           apphys = .FALSE.  
           statcl = .FALSE.  
           conser = .FALSE.  
           apdiss = .FALSE.  
   
           IF (MOD(itau, iconser) == 0) conser = .TRUE.  
           IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.  
           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.  
117    
118            ! calcul des tendances dynamiques:      time_integration: do itau = 0, itaufin - 1
119           leapf = mod(itau, iperiod) /= 0
120           if (leapf) then
121              dt = 2 * dtvr
122           else
123              ! Matsuno
124              dt = dtvr
125              if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
126                   call guide(itau, ucov, vcov, teta, q, masse, ps)
127              vcovm1 = vcov
128              ucovm1 = ucov
129              tetam1 = teta
130              massem1 = masse
131              psm1 = ps
132              finvmaold = masse
133              CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
134           end if
135    
136           ! Calcul des tendances dynamiques:
137           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
138           CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
139                dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
140                conser=MOD(itau, iconser)==0)
141    
142           ! Calcul des tendances advection des traceurs (dont l'humidité)
143           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
144    
145           ! Stokage du flux de masse pour traceurs offline:
146           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
147                dtvr, itau)
148    
149           ! Integrations dynamique et traceurs:
150           CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
151                dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
152                leapf)
153    
154           if (.not. leapf) then
155              ! Matsuno backward
156              forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
157              CALL exner_hyb(ps, p3d, pks, pk, pkf)
158    
159              ! Calcul des tendances dynamiques:
160            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
161              CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
162            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &                 phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
163                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &                 conser=.false.)
                time + iday - day_ini)  
   
           ! calcul des tendances advection des traceurs (dont l'humidite)  
   
           IF (forward .OR. leapf) THEN  
              CALL caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, pk)  
              IF (offline) THEN  
                 !maf stokage du flux de masse pour traceurs OFF-LINE  
                 CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
                      itau)  
              ENDIF  
           ENDIF  
164    
165            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
166            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
167                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold)                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
168                   finvmaold, dtvr, leapf=.false.)
169           end if
170    
171           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
172            ! calcul des tendances physiques:            ! calcul des tendances physiques:
173    
174            IF (apphys) THEN            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
175               IF (itau + 1 == itaufin) lafin = .TRUE.            CALL exner_hyb(ps, p3d, pks, pk, pkf)
   
              CALL pression(ip1jmp1, ap, bp, ps, p)  
              CALL exner_hyb(ps, p, pks, pk, pkf)  
176    
177               rdaym_ini = itau * dtvr / daysec            rdaym_ini = itau * dtvr / daysec
178               rdayvrai = rdaym_ini + day_ini            rdayvrai = rdaym_ini + day_ini
179              time = REAL(mod(itau, day_step)) / day_step + time_0
180              IF (time > 1.) time = time - 1.
181    
182              CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
183                   phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
184                   lafin=itau+1==itaufin)
185    
186              ! ajout des tendances physiques:
187              CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
188                   dtetafi, dqfi, dpfi)
189           ENDIF
190    
191               ! Interface avec les routines de phylmd (phymars ...)         forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
192           CALL exner_hyb(ps, p3d, pks, pk, pkf)
              ! Diagnostique de conservation de l'énergie : initialisation  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit='bil dyn'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys &  
                      , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
   
              CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
                   masse, ps, p, pk, phis, phi, du, dv, dteta, dq, w, &  
                   clesphy0, dufi, dvfi, dtetafi, dqfi, dpfi)  
   
              ! ajout des tendances physiques:  
              CALL addfi(nqmx, dtphys, &  
                   ucov, vcov, teta, q, ps, &  
                   dufi, dvfi, dtetafi, dqfi, dpfi)  
   
              ! Diagnostique de conservation de l'énergie : difference  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit = 'bil phys'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
           ENDIF  
   
           CALL pression(ip1jmp1, ap, bp, ps, p)  
           CALL exner_hyb(ps, p, pks, pk, pkf)  
193    
194           IF (MOD(itau + 1, idissip) == 0) THEN
195            ! dissipation horizontale et verticale des petites echelles:            ! dissipation horizontale et verticale des petites echelles:
196    
197            IF (apdiss) THEN            ! calcul de l'energie cinetique avant dissipation
198               ! calcul de l'energie cinetique avant dissipation            call covcont(llm, ucov, vcov, ucont, vcont)
199               call covcont(llm, ucov, vcov, ucont, vcont)            call enercin(vcov, ucov, vcont, ucont, ecin0)
200               call enercin(vcov, ucov, vcont, ucont, ecin0)  
201              ! dissipation
202               ! dissipation            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
203               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)            ucov=ucov + dudis
204               ucov=ucov + dudis            vcov=vcov + dvdis
205               vcov=vcov + dvdis  
206              ! On rajoute la tendance due à la transformation Ec -> E
207               if (dissip_conservative) then            ! thermique créée lors de la dissipation
208                  ! On rajoute la tendance due a la transform. Ec -> E            call covcont(llm, ucov, vcov, ucont, vcont)
209                  ! therm. cree lors de la dissipation            call enercin(vcov, ucov, vcont, ucont, ecin)
210                  call covcont(llm, ucov, vcov, ucont, vcont)            dtetaecdt= (ecin0 - ecin) / pk
211                  call enercin(vcov, ucov, vcont, ucont, ecin)            dtetadis=dtetadis + dtetaecdt
212                  dtetaecdt= (ecin0 - ecin) / pk            teta=teta + dtetadis
213                  dtetadis=dtetadis + dtetaecdt  
214               endif            ! Calcul de la valeur moyenne aux pôles :
215               teta=teta + dtetadis            forall (l = 1: llm)
216                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
217               ! Calcul de la valeur moyenne, unique de h aux poles .....                    / apoln
218                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
219               DO l = 1, llm                    * teta(:iim, jjm + 1, l)) / apols
220                  DO ij = 1, iim            END forall
221                     tppn(ij) = aire(ij) * teta(ij, l)  
222                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
223                  ENDDO            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
224                  tpn = SSUM(iim, tppn, 1) / apoln                 / apols
225                  tps = SSUM(iim, tpps, 1) / apols         END IF
226    
227                  DO ij = 1, iip1         IF (MOD(itau + 1, iperiod) == 0) THEN
228                     teta(ij, l) = tpn            ! Écriture du fichier histoire moyenne:
229                     teta(ij + ip1jm, l) = tps            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
230                  ENDDO                 phi, q, masse, ps, phis)
231               ENDDO            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
232                   q(:, :, :, 1), dt_app = dtvr * iperiod, &
233               DO ij = 1, iim                 dt_cum = dtvr * day_step * periodav)
234                  tppn(ij) = aire(ij) * ps(ij)         ENDIF
235                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)      end do time_integration
236               ENDDO  
237               tpn = SSUM(iim, tppn, 1) / apoln      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
238               tps = SSUM(iim, tpps, 1) / apols           itau=itau_dyn+itaufin)
239    
240               DO ij = 1, iip1      ! Calcul des tendances dynamiques:
241                  ps(ij) = tpn      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
242                  ps(ij + ip1jm) = tps      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
243               ENDDO           dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
244             conser=MOD(itaufin, iconser)==0)
           END IF  
   
           ! fin de l'intégration dynamique et physique pour le pas "itau"  
           ! préparation du pas d'intégration suivant  
   
           ! schema matsuno + leapfrog  
           IF (forward .OR. leapf) THEN  
              itau = itau + 1  
              iday = day_ini + itau / day_step  
              time = REAL(itau - (iday - day_ini) * day_step) / day_step &  
                   + time_0  
              IF (time > 1.) THEN  
                 time = time - 1.  
                 iday = iday + 1  
              ENDIF  
           ENDIF  
   
           IF (itau == itaufinp1) then  
              abort_message = 'Simulation finished'  
              call abort_gcm(modname, abort_message, 0)  
           ENDIF  
   
           ! ecriture du fichier histoire moyenne:  
   
           ! Comment out the following calls when you do not want the output  
           ! files "dyn_hist_ave.nc" and "dynzon.nc"  
           IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN  
              CALL writedynav(histaveid, nqmx, itau, vcov, &  
                   ucov, teta, pk, phi, q, masse, ps, phis)  
              call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &  
                   ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)  
           ENDIF  
   
           IF (itau == itaufin) THEN  
              CALL dynredem1("restart.nc", 0.0, &  
                   vcov, ucov, teta, q, nqmx, masse, ps)  
              CLOSE(99)  
           ENDIF  
   
           ! gestion de l'integration temporelle:  
   
           IF (MOD(itau, iperiod) == 0) exit  
           IF (MOD(itau - 1, iperiod) == 0) THEN  
              IF (forward) THEN  
                 ! fin du pas forward et debut du pas backward  
                 forward = .FALSE.  
                 leapf = .FALSE.  
              ELSE  
                 ! fin du pas backward et debut du premier pas leapfrog  
                 leapf = .TRUE.  
                 dt = 2. * dtvr  
              END IF  
           ELSE  
              ! ...... pas leapfrog .....  
              leapf = .TRUE.  
              dt = 2. * dtvr  
           END IF  
        end do  
     end do  
   
245    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
246    
247  end module leapfrog_m  end module leapfrog_m
Legend:
Removed from v.3  
changed lines
  Added in v.47
  ViewVC Help
Powered by ViewVC 1.1.21