/[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 5 by guez, Mon Mar 3 16:32:04 2008 UTC revision 34 by guez, Wed Jun 2 11:01:12 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, 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  
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
10      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11        ! schema matsuno + leapfrog
12    
13        USE calfis_m, ONLY: calfis
14        USE com_io_dyn, ONLY: histaveid
15        USE comconst, ONLY: daysec, dtphys, dtvr
16        USE comgeom, ONLY: aire_2d, apoln, apols
17        USE comvert, ONLY: ap, bp
18        USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
19             periodav
20        USE dimens_m, ONLY: iim, jjm, llm, nqmx
21        USE dynetat0_m, ONLY: day_ini
22        use dynredem1_m, only: dynredem1
23        USE exner_hyb_m, ONLY: exner_hyb
24        use filtreg_m, only: filtreg
25        USE guide_m, ONLY: guide
26        use inidissip_m, only: idissip
27        use integrd_m, only: integrd
28        USE logic, ONLY: iflag_phys, ok_guide
29        USE paramet_m, ONLY: ip1jmp1
30        USE pression_m, ONLY: pression
31        USE pressure_var, ONLY: p3d
32        USE temps, ONLY: itau_dyn
33    
34        ! Variables dynamiques:
35        REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
36        REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
37        REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature
38        REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa
39    
40      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      REAL masse(ip1jmp1, llm) ! masse d'air
41      ! Objet:      REAL phis(ip1jmp1) ! geopotentiel au sol
     ! 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.  
   
     ! ... Possibilite de choisir le shema 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, llm, nqmx  
     use paramet_m, only: ip1jmp1, ip1jm, llmp1, 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 abort_gcm_m, only: abort_gcm  
     use ener  
     use calfis_m, only: calfis  
     use exner_hyb_m, only: exner_hyb  
     use guide_m, only: guide  
     use pression_m, only: pression  
   
     integer nq  
   
     INTEGER longcles  
     PARAMETER (longcles = 20)  
     REAL clesphy0(longcles)  
   
     ! variables dynamiques  
     REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants  
     REAL teta(ip1jmp1, llm) ! temperature potentielle  
42      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
43      REAL ps(ip1jmp1) ! pression au sol      REAL, intent(in):: time_0
44      REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches  
45        ! Variables local to the procedure:
46    
47        ! 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
     REAL masse(ip1jmp1, llm) ! masse d'air  
     REAL phis(ip1jmp1) ! geopotentiel au sol  
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
80      INTEGER iday ! jour julien      real finvmaold(ip1jmp1, llm)
81      REAL time ! Heure de la journee en fraction d'1 jour      INTEGER l
   
     REAL SSUM  
     REAL time_0, finvmaold(ip1jmp1, llm)  
   
     LOGICAL :: lafin=.false.  
     INTEGER ij, l  
   
82      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
   
     data 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 ! PRINT level for energy conserv. diag.      real dt
     SAVE ip_ebil_dyn  
     DATA ip_ebil_dyn /0/  
   
     character(len=*), parameter:: modname = "leapfrog"  
     character*80 abort_message  
   
     logical dissip_conservative  
     save dissip_conservative  
     data dissip_conservative /.true./  
   
     LOGICAL prem  
     save prem  
     DATA prem /.true./  
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 :
104      dq=0.      CALL pression(ip1jmp1, ap, bp, ps, p3d)
105      CALL pression(ip1jmp1, ap, bp, ps, p)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
     CALL exner_hyb(ps, p, pks, pk, pkf)  
   
     ! Debut de l'integration temporelle:  
     do  
        if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then  
           call guide(itau, ucov, vcov, teta, q, masse, ps)  
        else  
           IF (prt_level > 9) print *, &  
                'Attention : on ne guide pas les 6 dernieres heures.'  
        endif  
   
        CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)  
        CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)  
        CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)  
        CALL SCOPY(ijp1llm, masse, 1, massem1, 1)  
        CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)  
   
        forward = .TRUE.  
        leapf = .FALSE.  
        dt = dtvr  
   
        CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)  
        CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)  
   
        do  
           ! gestion des appels de la physique et des dissipations:  
   
           apphys = .FALSE.  
           statcl = .FALSE.  
           conser = .FALSE.  
           apdiss = .FALSE.  
   
           IF (MOD(itau, iconser) == 0) conser = .TRUE.  
           IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.  
           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.  
106    
107            ! calcul des tendances dynamiques:      ! Début de l'integration temporelle :
108        do itau = 0, itaufin - 1
109           leapf = mod(itau, iperiod) /= 0
110           if (leapf) then
111              dt = 2 * dtvr
112           else
113              ! Matsuno
114              dt = dtvr
115              if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
116                   call guide(itau, ucov, vcov, teta, q, masse, ps)
117              vcovm1 = vcov
118              ucovm1 = ucov
119              tetam1 = teta
120              massem1 = masse
121              psm1 = ps
122              finvmaold = masse
123              CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
124           end if
125    
126           ! Calcul des tendances dynamiques:
127           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
128           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                time_0)
131    
132           ! Calcul des tendances advection des traceurs (dont l'humidité)
133           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134    
135           ! Stokage du flux de masse pour traceurs offline:
136           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
137                dtvr, itau)
138    
139           ! 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           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, p, 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, finvmaold)                 dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
156                   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    
162            IF (apphys) THEN            CALL pression(ip1jmp1, ap, bp, ps, p3d)
163               IF (itau + 1 == itaufin) lafin = .TRUE.            CALL exner_hyb(ps, p3d, pks, pk, pkf)
164    
165               CALL pression(ip1jmp1, ap, bp, ps, p)            rdaym_ini = itau * dtvr / daysec
166               CALL exner_hyb(ps, p, pks, pk, pkf)            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(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               rdaym_ini = itau * dtvr / daysec         CALL pression(ip1jmp1, ap, bp, ps, p3d)
180               rdayvrai = rdaym_ini + day_ini         CALL exner_hyb(ps, p3d, pks, pk, pkf)
   
              ! 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, p, pk, teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
   
              CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
                   masse, ps, p, pk, phis, phi, du, dv, dteta, dq, w, &  
                   clesphy0, 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, p, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
           ENDIF  
   
           CALL pression(ip1jmp1, ap, bp, ps, p)  
           CALL exner_hyb(ps, p, pks, pk, pkf)  
181    
182           IF (MOD(itau + 1, idissip) == 0) THEN
183            ! dissipation horizontale et verticale des petites echelles:            ! dissipation horizontale et verticale des petites echelles:
184    
185            IF (apdiss) THEN            ! calcul de l'energie cinetique avant dissipation
186               ! calcul de l'energie cinetique avant dissipation            call covcont(llm, ucov, vcov, ucont, vcont)
187               call covcont(llm, ucov, vcov, ucont, vcont)            call enercin(vcov, ucov, vcont, ucont, ecin0)
188               call enercin(vcov, ucov, vcont, ucont, ecin0)  
189              ! dissipation
190               ! dissipation            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
191               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)            ucov=ucov + dudis
192               ucov=ucov + dudis            vcov=vcov + dvdis
193               vcov=vcov + dvdis  
194              ! On rajoute la tendance due à la transformation Ec -> E
195               if (dissip_conservative) then            ! thermique créée lors de la dissipation
196                  ! On rajoute la tendance due a la transform. Ec -> E            call covcont(llm, ucov, vcov, ucont, vcont)
197                  ! therm. cree lors de la dissipation            call enercin(vcov, ucov, vcont, ucont, ecin)
198                  call covcont(llm, ucov, vcov, ucont, vcont)            dtetaecdt= (ecin0 - ecin) / pk
199                  call enercin(vcov, ucov, vcont, ucont, ecin)            dtetadis=dtetadis + dtetaecdt
200                  dtetaecdt= (ecin0 - ecin) / pk            teta=teta + dtetadis
201                  dtetadis=dtetadis + dtetaecdt  
202               endif            ! Calcul de la valeur moyenne aux pôles :
203               teta=teta + dtetadis            forall (l = 1: llm)
204                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
205               ! Calcul de la valeur moyenne, unique de h aux poles .....                    / apoln
206                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
207               DO l = 1, llm                    * teta(:iim, jjm + 1, l)) / apols
208                  DO ij = 1, iim            END forall
209                     tppn(ij) = aire(ij) * teta(ij, l)  
210                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
211                  ENDDO            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
212                  tpn = SSUM(iim, tppn, 1) / apoln                 / apols
213                  tps = SSUM(iim, tpps, 1) / apols         END IF
   
                 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) then  
              abort_message = 'Simulation finished'  
              call abort_gcm(modname, abort_message, 0)  
           ENDIF  
214    
215           IF (MOD(itau + 1, iperiod) == 0) THEN
216            ! ecriture du fichier histoire moyenne:            ! ecriture du fichier histoire moyenne:
217              CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
218            ! Comment out the following calls when you do not want the output                 phi, q, masse, ps, phis)
219            ! files "dyn_hist_ave.nc" and "dynzon.nc"            call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
220            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
221               CALL writedynav(histaveid, nqmx, itau, vcov, &         ENDIF
                   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  
222      end do      end do
223    
224        CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
225             itau=itau_dyn+itaufin)
226    
227        ! Calcul des tendances dynamiques:
228        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
229        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             time_0)
232    
233    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
234    
235  end module leapfrog_m  end module leapfrog_m
Legend:
Removed from v.5  
changed lines
  Added in v.34
  ViewVC Help
Powered by ViewVC 1.1.21