/[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 36 by guez, Thu Dec 2 17:11:04 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      ! Matsuno-leapfrog scheme.
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):: ucov(ip1jmp1, llm) ! vent covariant
36      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL, intent(inout):: vcov((iim + 1) * jjm, 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      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
40      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
41        REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
42      REAL time_0      REAL, intent(in):: time_0
43    
44      ! Variables local to the procedure:      ! Variables local to the procedure:
45    
46      ! Variables dynamiques:      ! Variables dynamiques:
47    
48      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
49      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
50      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
51      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(ip1jmp1, llm) ! geopotential
52      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
53    
54      ! variables dynamiques intermediaire pour le transport      ! variables dynamiques intermediaire pour le transport
55      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse      REAL pbaru(ip1jmp1, llm), pbarv((iim + 1) * jjm, llm) !flux de masse
56    
57      ! variables dynamiques au pas - 1      ! variables dynamiques au pas - 1
58      REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)      REAL vcovm1((iim + 1) * jjm, llm), ucovm1(ip1jmp1, llm)
59      REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)      REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
60      REAL massem1(ip1jmp1, llm)      REAL massem1(ip1jmp1, llm)
61    
62      ! tendances dynamiques      ! tendances dynamiques
63      REAL dv(ip1jm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)
64      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
65    
66      ! tendances de la dissipation      ! tendances de la dissipation
67      REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
68      REAL dtetadis(ip1jmp1, llm)      REAL dtetadis(iim + 1, jjm + 1, llm)
69    
70      ! tendances physiques      ! tendances physiques
71      REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
72      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
73    
74      ! variables pour le fichier histoire      ! variables pour le fichier histoire
75    
76      REAL tppn(iim), tpps(iim), tpn, tps      INTEGER itau ! index of the time step of the dynamics, starts at 0
77        INTEGER itaufin
78      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  
79      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
80        INTEGER l
     LOGICAL :: lafin=.false.  
     INTEGER ij, l  
   
81      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL:: callinigrads = .true.  
82    
83      !+jld variables test conservation energie      ! Variables test conservation energie
84      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
85      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
86      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
87      ! cree par la dissipation      ! cree par la dissipation
88      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
89      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
90      CHARACTER*15 ztit      logical leapf
91      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  
92    
93      !---------------------------------------------------      !---------------------------------------------------
94    
95      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
96    
97      itaufin = nday * day_step      itaufin = nday * day_step
98      itaufinp1 = itaufin + 1      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
99    
100      itau = 0      dq = 0.
     iday = day_ini  
     time = time_0  
     IF (time > 1.) THEN  
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
101    
102      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
     dq=0.  
103      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
104      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
105    
106      ! Debut de l'integration temporelle:      ! Début de l'integration temporelle :
107      outer_loop:do      do itau = 0, itaufin - 1
108         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then         leapf = mod(itau, iperiod) /= 0
109            call guide(itau, ucov, vcov, teta, q, masse, ps)         if (leapf) then
110              dt = 2 * dtvr
111         else         else
112            IF (prt_level > 9) print *, &            ! Matsuno
113                 'Attention : on ne guide pas les 6 dernieres heures.'            dt = dtvr
114         endif            if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
115                   call guide(itau, ucov, vcov, teta, q, masse, ps)
116         CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)            vcovm1 = vcov
117         CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)            ucovm1 = ucov
118         CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)            tetam1 = teta
119         CALL SCOPY(ijp1llm, masse, 1, massem1, 1)            massem1 = masse
120         CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)            psm1 = ps
121              finvmaold = masse
122         forward = .TRUE.            CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
123         leapf = .FALSE.         end if
124         dt = dtvr  
125           ! Calcul des tendances dynamiques:
126         CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
127         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
128                MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
129         do              time_0)
130            ! gestion des appels de la physique et des dissipations:  
131           ! Calcul des tendances advection des traceurs (dont l'humidité)
132            apphys = .FALSE.         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
133            conser = .FALSE.  
134            apdiss = .FALSE.         ! Stokage du flux de masse pour traceurs offline:
135           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
136            IF (MOD(itau, iconser) == 0) conser = .TRUE.              dtvr, itau)
137            IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.  
138            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.         ! integrations dynamique et traceurs:
139           CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
140                dp, vcov, ucov, teta, q, ps, masse, finvmaold, leapf, dt)
141    
142            ! calcul des tendances dynamiques:         if (.not. leapf) then
143              ! Matsuno backward
144              CALL pression(ip1jmp1, ap, bp, ps, p3d)
145              CALL exner_hyb(ps, p3d, pks, pk, pkf)
146    
147              ! Calcul des tendances dynamiques:
148            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
149              CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
150            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  
151    
152            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
153            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
154                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
155                 finvmaold, leapf)                 dtvr)
156           end if
157    
158           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
159            ! calcul des tendances physiques:            ! calcul des tendances physiques:
160    
           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  
   
161            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
162            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
163    
164            ! dissipation horizontale et verticale des petites echelles:            rdaym_ini = itau * dtvr / daysec
165              rdayvrai = rdaym_ini + day_ini
166              time = REAL(mod(itau, day_step)) / day_step + time_0
167              IF (time > 1.) time = time - 1.
168    
169              CALL calfis(itau + 1 == itaufin, rdayvrai, time, ucov, vcov, &
170                   teta, q, masse, ps, pk, phis, phi, du, dv, dteta, dq, w, dufi, &
171                   dvfi, dtetafi, dqfi, dpfi)
172    
173              ! ajout des tendances physiques:
174              CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
175                   dtetafi, dqfi, dpfi)
176           ENDIF
177    
178            IF (apdiss) THEN         CALL pression(ip1jmp1, ap, bp, ps, p3d)
179               ! 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  
180    
181            IF (itau == itaufinp1) exit outer_loop         IF (MOD(itau + 1, idissip) == 0) THEN
182              ! dissipation horizontale et verticale des petites echelles:
183    
184            ! ecriture du fichier histoire moyenne:            ! calcul de l'energie cinetique avant dissipation
185              call covcont(llm, ucov, vcov, ucont, vcont)
186              call enercin(vcov, ucov, vcont, ucont, ecin0)
187    
188              ! dissipation
189              CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
190              ucov=ucov + dudis
191              vcov=vcov + dvdis
192    
193              ! On rajoute la tendance due à la transformation Ec -> E
194              ! thermique créée lors de la dissipation
195              call covcont(llm, ucov, vcov, ucont, vcont)
196              call enercin(vcov, ucov, vcont, ucont, ecin)
197              dtetaecdt= (ecin0 - ecin) / pk
198              dtetadis=dtetadis + dtetaecdt
199              teta=teta + dtetadis
200    
201              ! Calcul de la valeur moyenne aux pôles :
202              forall (l = 1: llm)
203                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
204                      / apoln
205                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
206                      * teta(:iim, jjm + 1, l)) / apols
207              END forall
208    
209              ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
210              ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
211                   / apols
212           END IF
213    
214            ! Comment out the following calls when you do not want the output         IF (MOD(itau + 1, iperiod) == 0) THEN
215            ! files "dyn_hist_ave.nc" and "dynzon.nc"            ! ecriture du fichier histoire moyenne:
216            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
217               CALL writedynav(histaveid, nqmx, itau, vcov, &                 phi, q, masse, ps, phis)
218                    ucov, teta, pk, phi, q, masse, ps, phis)            call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
219               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
220                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)         ENDIF
221            ENDIF      end do
222    
223            IF (itau == itaufin) THEN      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
224               CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)           itau=itau_dyn+itaufin)
225               CLOSE(99)  
226            ENDIF      ! Calcul des tendances dynamiques:
227        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
228            ! gestion de l'integration temporelle:      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
229             MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
230            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  
231    
232    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
233    
Legend:
Removed from v.18  
changed lines
  Added in v.36
  ViewVC Help
Powered by ViewVC 1.1.21