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

Diff of /trunk/libf/dyn3d/leapfrog.f90

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

revision 5 by guez, Mon Mar 3 16:32:04 2008 UTC revision 32 by guez, Tue Apr 6 17:52:58 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
     INTEGER itau, itaufinp1  
79      INTEGER iday ! jour julien      INTEGER iday ! jour julien
80      REAL time ! Heure de la journee en fraction d'1 jour      REAL time ! time of day, as a fraction of day length
81        real finvmaold(ip1jmp1, llm)
82      REAL SSUM      LOGICAL:: lafin=.false.
83      REAL time_0, finvmaold(ip1jmp1, llm)      INTEGER i, j, l
   
     LOGICAL :: lafin=.false.  
     INTEGER ij, l  
84    
85      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
86    
87      data callinigrads/.true./      ! Variables test conservation energie
88        REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
     !+jld variables test conservation energie  
     REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)  
89      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
90      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
91      ! cree par la dissipation      ! cree par la dissipation
92      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
93      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
     CHARACTER*15 ztit  
     INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.  
     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./  
94    
95      !---------------------------------------------------      !---------------------------------------------------
96    
97      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
98    
99      itaufin = nday * day_step      itaufin = nday * day_step
100      itaufinp1 = itaufin + 1      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
101    
102      itau = 0      itau = 0
103      iday = day_ini      iday = day_ini
104      time = time_0      time = time_0
105      IF (time > 1.) THEN      dq = 0.
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
   
106      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
107      dq=0.      CALL pression(ip1jmp1, ap, bp, ps, p3d)
108      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.  
109    
110            ! calcul des tendances dynamiques:      ! Début de l'integration temporelle :
111        period_loop:do i = 1, itaufin / iperiod
112           ! {"itau" is a multiple of "iperiod"}
113    
114           ! 1. Matsuno forward:
115    
116           if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
117                call guide(itau, ucov, vcov, teta, q, masse, ps)
118           vcovm1 = vcov
119           ucovm1 = ucov
120           tetam1 = teta
121           massem1 = masse
122           psm1 = ps
123           finvmaold = masse
124           CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
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 + iday - day_ini)
131    
132           ! Calcul des tendances advection des traceurs (dont l'humidité)
133           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
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, .false., &
141                dtvr)
142    
143           CALL pression(ip1jmp1, ap, bp, ps, p3d)
144           CALL exner_hyb(ps, p3d, pks, pk, pkf)
145    
146           ! 2. Matsuno backward:
147    
148           itau = itau + 1
149           iday = day_ini + itau / day_step
150           time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
151           IF (time > 1.) THEN
152              time = time - 1.
153              iday = iday + 1
154           ENDIF
155    
156           ! Calcul des tendances dynamiques:
157           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
158           CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
159                .false., du, dv, dteta, dp, w, pbaru, pbarv, time + iday - day_ini)
160    
161           ! integrations dynamique et traceurs:
162           CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
163                dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
164                dtvr)
165    
166            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL pression(ip1jmp1, ap, bp, ps, p3d)
167           CALL exner_hyb(ps, p3d, pks, pk, pkf)
168    
169           ! 3. Leapfrog:
170    
171           leapfrog_loop: do j = 1, iperiod - 1
172              ! Calcul des tendances dynamiques:
173              CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
174            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
175                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &                 .false., du, dv, dteta, dp, w, pbaru, pbarv, &
176                 time + iday - day_ini)                 time + iday - day_ini)
177    
178            ! calcul des tendances advection des traceurs (dont l'humidite)            ! Calcul des tendances advection des traceurs (dont l'humidité)
179              CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
180            IF (forward .OR. leapf) THEN            ! Stokage du flux de masse pour traceurs off-line:
181               CALL caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, pk)            IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
182               IF (offline) THEN                 dtvr, itau)
                 !maf stokage du flux de masse pour traceurs OFF-LINE  
                 CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
                      itau)  
              ENDIF  
           ENDIF  
183    
184            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
185            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
186                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold)                 dteta, dp, vcov, ucov, teta, q, ps, masse, &
187                   finvmaold, .true., 2 * dtvr)
           ! calcul des tendances physiques:  
188    
189            IF (apphys) THEN            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
190                 ! calcul des tendances physiques:
191               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
192    
193               CALL pression(ip1jmp1, ap, bp, ps, p)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
194               CALL exner_hyb(ps, p, pks, pk, pkf)               CALL exner_hyb(ps, p3d, pks, pk, pkf)
195    
196               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
197               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
198    
199               ! Interface avec les routines de phylmd (phymars ...)               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
200                      masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
              ! 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, &  
201                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
202    
203               ! Diagnostique de conservation de l'énergie : difference               ! ajout des tendances physiques:
204               IF (ip_ebil_dyn >= 1) THEN               CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
205                  ztit = 'bil phys'                    dtetafi, dqfi, dpfi)
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
206            ENDIF            ENDIF
207    
208            CALL pression(ip1jmp1, ap, bp, ps, p)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
209            CALL exner_hyb(ps, p, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
210    
211            ! dissipation horizontale et verticale des petites echelles:            IF (MOD(itau + 1, idissip) == 0) THEN
212                 ! dissipation horizontale et verticale des petites echelles:
213    
           IF (apdiss) THEN  
214               ! calcul de l'energie cinetique avant dissipation               ! calcul de l'energie cinetique avant dissipation
215               call covcont(llm, ucov, vcov, ucont, vcont)               call covcont(llm, ucov, vcov, ucont, vcont)
216               call enercin(vcov, ucov, vcont, ucont, ecin0)               call enercin(vcov, ucov, vcont, ucont, ecin0)
217    
218               ! dissipation               ! dissipation
219               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
220               ucov=ucov + dudis               ucov=ucov + dudis
221               vcov=vcov + dvdis               vcov=vcov + dvdis
222    
223               if (dissip_conservative) then               ! On rajoute la tendance due à la transformation Ec -> E
224                  ! On rajoute la tendance due a la transform. Ec -> E               ! thermique créée lors de la dissipation
225                  ! therm. cree lors de la dissipation               call covcont(llm, ucov, vcov, ucont, vcont)
226                  call covcont(llm, ucov, vcov, ucont, vcont)               call enercin(vcov, ucov, vcont, ucont, ecin)
227                  call enercin(vcov, ucov, vcont, ucont, ecin)               dtetaecdt= (ecin0 - ecin) / pk
228                  dtetaecdt= (ecin0 - ecin) / pk               dtetadis=dtetadis + dtetaecdt
                 dtetadis=dtetadis + dtetaecdt  
              endif  
229               teta=teta + dtetadis               teta=teta + dtetadis
230    
231               ! Calcul de la valeur moyenne, unique de h aux poles .....               ! Calcul de la valeur moyenne aux pôles :
232                 forall (l = 1: llm)
233               DO l = 1, llm                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
234                  DO ij = 1, iim                       / apoln
235                     tppn(ij) = aire(ij) * teta(ij, l)                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
236                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)                       * teta(:iim, jjm + 1, l)) / apols
237                  ENDDO               END forall
238                  tpn = SSUM(iim, tppn, 1) / apoln  
239                  tps = SSUM(iim, tpps, 1) / apols               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
240                 ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
241                  DO ij = 1, iip1                    / apols
                    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  
   
242            END IF            END IF
243    
244            ! fin de l'intégration dynamique et physique pour le pas "itau"            itau = itau + 1
245            ! préparation du pas d'intégration suivant            iday = day_ini + itau / day_step
246              time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
247            ! schema matsuno + leapfrog            IF (time > 1.) THEN
248            IF (forward .OR. leapf) THEN               time = time - 1.
249               itau = itau + 1               iday = iday + 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)  
250            ENDIF            ENDIF
251    
252            ! ecriture du fichier histoire moyenne:            IF (MOD(itau, iperiod) == 0) THEN
253                 ! 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  
254               CALL writedynav(histaveid, nqmx, itau, vcov, &               CALL writedynav(histaveid, nqmx, itau, vcov, &
255                    ucov, teta, pk, phi, q, masse, ps, phis)                    ucov, teta, pk, phi, q, masse, ps, phis)
256               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
257                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
258            ENDIF            ENDIF
259           end do leapfrog_loop
260        end do period_loop
261    
262            IF (itau == itaufin) THEN      ! {itau == itaufin}
263               CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
264               CLOSE(99)           itau=itau_dyn+itaufin)
265            ENDIF  
266        ! Calcul des tendances dynamiques:
267            ! gestion de l'integration temporelle:      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
268        CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
269            IF (MOD(itau, iperiod) == 0) exit           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
270            IF (MOD(itau - 1, iperiod) == 0) THEN           time + iday - day_ini)
              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  
271    
272    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
273    
Legend:
Removed from v.5  
changed lines
  Added in v.32
  ViewVC Help
Powered by ViewVC 1.1.21