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

Diff of /trunk/Sources/dyn3d/leapfrog.f

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

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