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

Diff of /trunk/dyn3d/leapfrog.f90

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

trunk/libf/dyn3d/leapfrog.f90 revision 10 by guez, Fri Apr 18 14:45:53 2008 UTC trunk/dyn3d/leapfrog.f revision 262 by guez, Wed Mar 7 13:46:18 2018 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)  
   
     ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34  
8    
9      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
10      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11    
12      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      ! Intégration temporelle du modèle : Matsuno-leapfrog scheme.
     ! 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.  
   
     ! ... Possibilité de choisir le schéma 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, jjm, llm, nqmx  
     use paramet_m, only: ip1jmp1, ip1jm, 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 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  
13    
14      integer nq      use addfi_m, only: addfi
15      REAL clesphy0(:)      use bilan_dyn_m, only: bilan_dyn
16        use caladvtrac_m, only: caladvtrac
17        use caldyn_m, only: caldyn
18        USE calfis_m, ONLY: calfis
19        USE comconst, ONLY: dtvr
20        USE comgeom, ONLY: aire_2d, apoln, apols
21        use covcont_m, only: covcont
22        USE disvert_m, ONLY: ap, bp
23        USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, &
24             iflag_phys, iecri
25        USE conf_guide_m, ONLY: ok_guide
26        USE dimens_m, ONLY: iim, jjm, llm, nqmx
27        use dissip_m, only: dissip
28        USE dynetat0_m, ONLY: day_ini
29        use dynredem1_m, only: dynredem1
30        use enercin_m, only: enercin
31        USE exner_hyb_m, ONLY: exner_hyb
32        use filtreg_scal_m, only: filtreg_scal
33        use geopot_m, only: geopot
34        USE guide_m, ONLY: guide
35        use inidissip_m, only: idissip
36        use integrd_m, only: integrd
37        use nr_util, only: assert
38        USE temps, ONLY: itau_dyn
39        use writehist_m, only: writehist
40    
41      ! Variables dynamiques:      ! Variables dynamiques:
42      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
43      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
     REAL ps(ip1jmp1) ! pression au sol, en Pa  
     REAL masse(ip1jmp1, llm) ! masse d'air  
     REAL phis(ip1jmp1) ! geopotentiel au sol  
44    
45      REAL time_0      REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
46        ! potential temperature
47    
48      ! Variables local to the procedure:      REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
49        REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
50        REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
51    
52      ! Variables dynamiques:      REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
53        ! mass fractions of advected fields
54    
55      REAL pks(ip1jmp1) ! exner au sol      ! Local:
56      REAL pk(ip1jmp1, llm) ! exner au milieu des couches  
57      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      ! Variables dynamiques:
     REAL phi(ip1jmp1, llm) ! geopotential  
     REAL w(ip1jmp1, llm) ! vitesse verticale  
   
     ! variables dynamiques intermediaire pour le transport  
     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse  
   
     ! variables dynamiques au pas - 1  
     REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)  
     REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)  
     REAL massem1(ip1jmp1, llm)  
   
     ! tendances dynamiques  
     REAL dv(ip1jm, llm), du(ip1jmp1, llm)  
     REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)  
   
     ! tendances de la dissipation  
     REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)  
     REAL dtetadis(ip1jmp1, llm)  
   
     ! tendances physiques  
     REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)  
     REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)  
   
     ! variables pour le fichier histoire  
   
     REAL tppn(iim), tpps(iim), tpn, tps  
   
     INTEGER itau, itaufinp1  
     INTEGER iday ! jour julien  
     REAL time ! Heure de la journee en fraction d'1 jour  
   
     REAL SSUM  
     real finvmaold(ip1jmp1, llm)  
   
     LOGICAL :: lafin=.false.  
     INTEGER ij, l  
   
     REAL rdayvrai, rdaym_ini  
     LOGICAL:: callinigrads = .true.  
   
     !+jld variables test conservation energie  
     REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)  
     ! Tendance de la temp. potentiel d (theta) / d t due a la  
     ! tansformation d'energie cinetique en energie thermique  
     ! cree par la dissipation  
     REAL dtetaecdt(ip1jmp1, llm)  
     REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)  
     CHARACTER*15 ztit  
     INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.  
58    
59      logical:: dissip_conservative = .true.      REAL pks(iim + 1, jjm + 1) ! exner au sol
60      LOGICAL:: prem = .true.      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
61        REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
62        REAL phi(iim + 1, jjm + 1, llm) ! geopotential
63        REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
64    
65        ! Variables dynamiques interm\'ediaires pour le transport
66        ! Flux de masse :
67        REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
68    
69        ! Variables dynamiques au pas - 1
70        REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
71        REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
72        REAL massem1(iim + 1, jjm + 1, llm)
73    
74        ! Tendances dynamiques
75        REAL dv((iim + 1) * jjm, llm), du(iim + 1, jjm + 1, llm)
76        REAL dteta(iim + 1, jjm + 1, llm)
77        real dp(iim + 1, jjm + 1)
78    
79        ! Tendances de la dissipation :
80        REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
81        REAL dtetadis(iim + 1, jjm + 1, llm)
82    
83        ! Tendances physiques
84        REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
85        REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
86    
87        ! Variables pour le fichier histoire
88        INTEGER itau ! index of the time step of the dynamics, starts at 0
89        INTEGER itaufin
90        INTEGER l
91    
92        ! Variables test conservation \'energie
93        REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
94    
95        REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
96        logical leapf
97        real dt ! time step, in s
98    
99        REAL p3d(iim + 1, jjm + 1, llm + 1) ! pressure at layer interfaces, in Pa
100        ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
101        ! for interface "l")
102    
103      !---------------------------------------------------      !---------------------------------------------------
104    
105      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
106        call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
107    
108      itaufin = nday * day_step      itaufin = nday * day_step
109      itaufinp1 = itaufin + 1      ! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
   
     itau = 0  
     iday = day_ini  
     time = time_0  
     IF (time > 1.) THEN  
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
110    
111      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
112      dq=0.      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
113      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL exner_hyb(ps, p3d, pks, pk)
114      CALL exner_hyb(ps, p3d, pks, pk, pkf)      pkf = pk
115        CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
116      ! Debut de l'integration temporelle:  
117      outer_loop:do      time_integration: do itau = 0, itaufin - 1
118         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then         leapf = mod(itau, iperiod) /= 0
119            call guide(itau, ucov, vcov, teta, q, masse, ps)         if (leapf) then
120              dt = 2 * dtvr
121         else         else
122            IF (prt_level > 9) print *, &            ! Matsuno
123                 'Attention : on ne guide pas les 6 dernieres heures.'            dt = dtvr
124         endif            if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
125              vcovm1 = vcov
126         CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)            ucovm1 = ucov
127         CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)            tetam1 = teta
128         CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)            massem1 = masse
129         CALL SCOPY(ijp1llm, masse, 1, massem1, 1)            psm1 = ps
130         CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)         end if
131    
132         forward = .TRUE.         ! Calcul des tendances dynamiques:
133         leapf = .FALSE.         CALL geopot(teta, pk, pks, phis, phi)
134         dt = dtvr         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
135                du, dv, dteta, dp, w, pbaru, pbarv, &
136         CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)              conser = MOD(itau, iconser) == 0)
137         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)  
138           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
139         do  
140            ! gestion des appels de la physique et des dissipations:         ! Int\'egrations dynamique et traceurs:
141           CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
142            apphys = .FALSE.              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dt, leapf)
143            statcl = .FALSE.  
144            conser = .FALSE.         forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
145            apdiss = .FALSE.         CALL exner_hyb(ps, p3d, pks, pk)
146           pkf = pk
147            IF (MOD(itau, iconser) == 0) conser = .TRUE.         CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
148            IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.  
149            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.         if (.not. leapf) then
150              ! Matsuno backward
151            ! calcul des tendances dynamiques:            ! Calcul des tendances dynamiques:
152              CALL geopot(teta, pk, pks, phis, phi)
153            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
154                   phi, du, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
           CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &  
                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  
155    
156            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
157            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
158                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold)                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dtvr, &
159                   leapf=.false.)
160            ! calcul des tendances physiques:  
161              forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
162            IF (apphys) THEN            CALL exner_hyb(ps, p3d, pks, pk)
163               IF (itau + 1 == itaufin) lafin = .TRUE.            pkf = pk
164              CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
165               CALL pression(ip1jmp1, ap, bp, ps, p3d)         end if
166               CALL exner_hyb(ps, p3d, pks, pk, pkf)  
167           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys) THEN
168               rdaym_ini = itau * dtvr / daysec            CALL calfis(ucov, vcov, teta, q, p3d, pk, phis, phi, w, dufi, dvfi, &
169               rdayvrai = rdaym_ini + day_ini                 dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
170                   time = REAL(mod(itau, day_step)) / day_step, &
171               ! Interface avec les routines de phylmd (phymars ...)                 lafin = itau + 1 == itaufin)
172    
173               ! Diagnostique de conservation de l'énergie : initialisation            CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
174               IF (ip_ebil_dyn >= 1) THEN         ENDIF
175                  ztit='bil dyn'  
176                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &         IF (MOD(itau + 1, idissip) == 0) THEN
177                       teta, q(:, :, 1), q(:, :, 2))            ! Dissipation horizontale et verticale des petites \'echelles
178               ENDIF  
179              ! calcul de l'\'energie cin\'etique avant dissipation
180               CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &            call covcont(llm, ucov, vcov, ucont, vcont)
181                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &            call enercin(vcov, ucov, vcont, ucont, ecin0)
182                    clesphy0, dufi, dvfi, dtetafi, dqfi, dpfi)  
183              ! dissipation
184               ! ajout des tendances physiques:            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
185               CALL addfi(nqmx, dtphys, &            ucov = ucov + dudis
186                    ucov, vcov, teta, q, ps, &            vcov = vcov + dvdis
187                    dufi, dvfi, dtetafi, dqfi, dpfi)  
188              ! On ajoute la tendance due \`a la transformation \'energie
189               ! Diagnostique de conservation de l'énergie : difference            ! cin\'etique en \'energie thermique par la dissipation
190               IF (ip_ebil_dyn >= 1) THEN            call covcont(llm, ucov, vcov, ucont, vcont)
191                  ztit = 'bil phys'            call enercin(vcov, ucov, vcont, ucont, ecin)
192                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &            dtetadis = dtetadis + (ecin0 - ecin) / pk
193                       teta, q(:, :, 1), q(:, :, 2))            teta = teta + dtetadis
194               ENDIF  
195            ENDIF            ! Calcul de la valeur moyenne aux p\^oles :
196              forall (l = 1: llm)
197            CALL pression(ip1jmp1, ap, bp, ps, p3d)               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
198            CALL exner_hyb(ps, p3d, pks, pk, pkf)                    / apoln
199                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm + 1) &
200            ! dissipation horizontale et verticale des petites echelles:                    * teta(:iim, jjm + 1, l)) / apols
201              END forall
202            IF (apdiss) THEN         END IF
203               ! calcul de l'energie cinetique avant dissipation  
204               call covcont(llm, ucov, vcov, ucont, vcont)         IF (MOD(itau + 1, iperiod) == 0) THEN
205               call enercin(vcov, ucov, vcont, ucont, ecin0)            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
206                   q(:, :, :, 1))
207               ! dissipation         ENDIF
208               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)  
209               ucov=ucov + dudis         IF (MOD(itau + 1, iecri * day_step) == 0) THEN
210               vcov=vcov + dvdis            CALL geopot(teta, pk, pks, phis, phi)
211              CALL writehist(vcov, ucov, teta, pk, phi, q, masse, ps, itau)
212               if (dissip_conservative) then         END IF
213                  ! On rajoute la tendance due a la transform. Ec -> E      end do time_integration
214                  ! therm. cree lors de la dissipation  
215                  call covcont(llm, ucov, vcov, ucont, vcont)      CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
216                  call enercin(vcov, ucov, vcont, ucont, ecin)  
217                  dtetaecdt= (ecin0 - ecin) / pk      ! Calcul des tendances dynamiques:
218                  dtetadis=dtetadis + dtetaecdt      CALL geopot(teta, pk, pks, phis, phi)
219               endif      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
220               teta=teta + dtetadis           du, dv, dteta, dp, w, pbaru, pbarv, &
221             conser = MOD(itaufin, iconser) == 0)
              ! 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  
   
           IF (itau == itaufinp1) exit outer_loop  
   
           ! 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., vcov, ucov, teta, q, 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 outer_loop  
222    
223    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
224    
Legend:
Removed from v.10  
changed lines
  Added in v.262
  ViewVC Help
Powered by ViewVC 1.1.21