/[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

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