/[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 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, clesphy0, &    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
8         time_0)  
9        ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34
10        ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11        ! Matsuno-leapfrog scheme.
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):: ucov(ip1jmp1, llm) ! vent covariant
36        REAL, intent(inout):: vcov((iim + 1) * jjm, 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        REAL masse(ip1jmp1, llm) ! masse d'air
40        REAL phis(ip1jmp1) ! geopotentiel au sol
41        REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
42        REAL, intent(in):: time_0
43    
44      ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34      ! Variables local to the procedure:
45    
46      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      ! Variables dynamiques:
     ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)  
47    
     ! Auteur: P. Le Van /L. Fairhead/F.Hourdin  
     ! 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.  
   
     ! ... 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  
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
     REAL ps(ip1jmp1) ! pression au sol  
     REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches  
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
     REAL masse(ip1jmp1, llm) ! masse d'air  
     REAL phis(ip1jmp1) ! geopotentiel au sol  
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
79      INTEGER iday ! jour julien      real finvmaold(ip1jmp1, llm)
80      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  
   
81      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
   
     data 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 ! 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./  
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 :
103      dq=0.      CALL pression(ip1jmp1, ap, bp, ps, p3d)
104      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.  
105    
106            ! calcul des tendances dynamiques:      ! Début de l'integration temporelle :
107        do itau = 0, itaufin - 1
108           leapf = mod(itau, iperiod) /= 0
109           if (leapf) then
110              dt = 2 * dtvr
111           else
112              ! Matsuno
113              dt = dtvr
114              if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
115                   call guide(itau, ucov, vcov, teta, q, masse, ps)
116              vcovm1 = vcov
117              ucovm1 = ucov
118              tetam1 = teta
119              massem1 = masse
120              psm1 = ps
121              finvmaold = masse
122              CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
123           end if
124    
125           ! Calcul des tendances dynamiques:
126           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
127           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                time_0)
130    
131           ! Calcul des tendances advection des traceurs (dont l'humidité)
132           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
133    
134           ! Stokage du flux de masse pour traceurs offline:
135           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
136                dtvr, itau)
137    
138           ! 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           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, 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  
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, finvmaold)                 dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
155                   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    
161            IF (apphys) THEN            CALL pression(ip1jmp1, ap, bp, ps, p3d)
162               IF (itau + 1 == itaufin) lafin = .TRUE.            CALL exner_hyb(ps, p3d, pks, pk, pkf)
163    
164               CALL pression(ip1jmp1, ap, bp, ps, p)            rdaym_ini = itau * dtvr / daysec
165               CALL exner_hyb(ps, p, pks, pk, pkf)            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               rdaym_ini = itau * dtvr / daysec         CALL pression(ip1jmp1, ap, bp, ps, p3d)
179               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)  
180    
181           IF (MOD(itau + 1, idissip) == 0) THEN
182            ! dissipation horizontale et verticale des petites echelles:            ! dissipation horizontale et verticale des petites echelles:
183    
184            IF (apdiss) THEN            ! calcul de l'energie cinetique avant dissipation
185               ! calcul de l'energie cinetique avant dissipation            call covcont(llm, ucov, vcov, ucont, vcont)
186               call covcont(llm, ucov, vcov, ucont, vcont)            call enercin(vcov, ucov, vcont, ucont, ecin0)
187               call enercin(vcov, ucov, vcont, ucont, ecin0)  
188              ! dissipation
189               ! dissipation            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
190               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)            ucov=ucov + dudis
191               ucov=ucov + dudis            vcov=vcov + dvdis
192               vcov=vcov + dvdis  
193              ! On rajoute la tendance due à la transformation Ec -> E
194               if (dissip_conservative) then            ! thermique créée lors de la dissipation
195                  ! On rajoute la tendance due a la transform. Ec -> E            call covcont(llm, ucov, vcov, ucont, vcont)
196                  ! therm. cree lors de la dissipation            call enercin(vcov, ucov, vcont, ucont, ecin)
197                  call covcont(llm, ucov, vcov, ucont, vcont)            dtetaecdt= (ecin0 - ecin) / pk
198                  call enercin(vcov, ucov, vcont, ucont, ecin)            dtetadis=dtetadis + dtetaecdt
199                  dtetaecdt= (ecin0 - ecin) / pk            teta=teta + dtetadis
200                  dtetadis=dtetadis + dtetaecdt  
201               endif            ! Calcul de la valeur moyenne aux pôles :
202               teta=teta + dtetadis            forall (l = 1: llm)
203                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
204               ! Calcul de la valeur moyenne, unique de h aux poles .....                    / apoln
205                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
206               DO l = 1, llm                    * teta(:iim, jjm + 1, l)) / apols
207                  DO ij = 1, iim            END forall
208                     tppn(ij) = aire(ij) * teta(ij, l)  
209                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
210                  ENDDO            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
211                  tpn = SSUM(iim, tppn, 1) / apoln                 / apols
212                  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  
213    
214           IF (MOD(itau + 1, iperiod) == 0) THEN
215            ! ecriture du fichier histoire moyenne:            ! ecriture du fichier histoire moyenne:
216              CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
217            ! Comment out the following calls when you do not want the output                 phi, q, masse, ps, phis)
218            ! files "dyn_hist_ave.nc" and "dynzon.nc"            call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
219            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
220               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  
221      end do      end do
222    
223        CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
224             itau=itau_dyn+itaufin)
225    
226        ! Calcul des tendances dynamiques:
227        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
228        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             time_0)
231    
232    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
233    
234  end module leapfrog_m  end module leapfrog_m
Legend:
Removed from v.5  
changed lines
  Added in v.36
  ViewVC Help
Powered by ViewVC 1.1.21