/[lmdze]/trunk/libf/dyn3d/leapfrog.f90
ViewVC logotype

Diff of /trunk/libf/dyn3d/leapfrog.f90

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

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