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

Diff of /trunk/dyn3d/leapfrog.f

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

revision 32 by guez, Tue Apr 6 17:52:58 2010 UTC revision 33 by guez, Fri Apr 9 10:56:14 2010 UTC
# Line 76  contains Line 76  contains
76    
77      INTEGER itau ! index of the time step of the dynamics, starts at 0      INTEGER itau ! index of the time step of the dynamics, starts at 0
78      INTEGER itaufin      INTEGER itaufin
     INTEGER iday ! jour julien  
79      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
80      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
81      LOGICAL:: lafin=.false.      INTEGER l
     INTEGER i, j, l  
   
82      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
83    
84      ! Variables test conservation energie      ! Variables test conservation energie
# Line 91  contains Line 88  contains
88      ! cree par la dissipation      ! cree par la dissipation
89      REAL dtetaecdt(iim + 1, jjm + 1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
90      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
91        logical leapf
92        real dt
93    
94      !---------------------------------------------------      !---------------------------------------------------
95    
# Line 99  contains Line 98  contains
98      itaufin = nday * day_step      itaufin = nday * day_step
99      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
100    
     itau = 0  
     iday = day_ini  
     time = time_0  
101      dq = 0.      dq = 0.
102    
103      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
104      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
105      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
106    
107      ! Début de l'integration temporelle :      ! Début de l'integration temporelle :
108      period_loop:do i = 1, itaufin / iperiod      do itau = 0, itaufin - 1
109         ! {"itau" is a multiple of "iperiod"}         leapf = mod(itau, iperiod) /= 0
110           if (leapf) then
111         ! 1. Matsuno forward:            dt = 2 * dtvr
112           else
113         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &            ! Matsuno
114              call guide(itau, ucov, vcov, teta, q, masse, ps)            dt = dtvr
115         vcovm1 = vcov            if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
116         ucovm1 = ucov                 call guide(itau, ucov, vcov, teta, q, masse, ps)
117         tetam1 = teta            vcovm1 = vcov
118         massem1 = masse            ucovm1 = ucov
119         psm1 = ps            tetam1 = teta
120         finvmaold = masse            massem1 = masse
121         CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)            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:         ! Calcul des tendances dynamiques:
127         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
128         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
129              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
130              time + iday - day_ini)              time_0)
131    
132         ! Calcul des tendances advection des traceurs (dont l'humidité)         ! Calcul des tendances advection des traceurs (dont l'humidité)
133         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134    
135         ! Stokage du flux de masse pour traceurs offline:         ! Stokage du flux de masse pour traceurs offline:
136         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
137              dtvr, itau)              dtvr, itau)
138    
139         ! integrations dynamique et traceurs:         ! integrations dynamique et traceurs:
140         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
141              dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &              dp, vcov, ucov, teta, q, ps, masse, finvmaold, leapf, dt)
             dtvr)  
   
        CALL pression(ip1jmp1, ap, bp, ps, p3d)  
        CALL exner_hyb(ps, p3d, pks, pk, pkf)  
   
        ! 2. Matsuno backward:  
   
        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  
   
        ! Calcul des tendances dynamiques:  
        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)  
        CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &  
             .false., du, dv, dteta, dp, w, pbaru, pbarv, time + iday - day_ini)  
142    
143         ! integrations dynamique et traceurs:         if (.not. leapf) then
144         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &            ! Matsuno backward
145              dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &            CALL pression(ip1jmp1, ap, bp, ps, p3d)
146              dtvr)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
   
        CALL pression(ip1jmp1, ap, bp, ps, p3d)  
        CALL exner_hyb(ps, p3d, pks, pk, pkf)  
   
        ! 3. Leapfrog:  
147    
        leapfrog_loop: do j = 1, iperiod - 1  
148            ! Calcul des tendances dynamiques:            ! 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, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
151                 .false., du, dv, dteta, dp, w, pbaru, pbarv, &                 phi, .false., du, dv, dteta, dp, w, pbaru, pbarv, time_0)
                time + iday - day_ini)  
   
           ! Calcul des tendances advection des traceurs (dont l'humidité)  
           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)  
           ! Stokage du flux de masse pour traceurs off-line:  
           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &  
                dtvr, itau)  
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, dp, vcov, ucov, teta, q, ps, masse, &                 dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
156                 finvmaold, .true., 2 * dtvr)                 dtvr)
157           end if
158    
159            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
160               ! calcul des tendances physiques:            ! calcul des tendances physiques:
              IF (itau + 1 == itaufin) lafin = .TRUE.  
   
              CALL pression(ip1jmp1, ap, bp, ps, p3d)  
              CALL exner_hyb(ps, p3d, pks, pk, pkf)  
   
              rdaym_ini = itau * dtvr / daysec  
              rdayvrai = rdaym_ini + day_ini  
   
              CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
                   masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &  
                   dufi, dvfi, dtetafi, dqfi, dpfi)  
   
              ! ajout des tendances physiques:  
              CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &  
                   dtetafi, dqfi, dpfi)  
           ENDIF  
161    
162            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
163            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
164    
165            IF (MOD(itau + 1, idissip) == 0) THEN            rdaym_ini = itau * dtvr / daysec
166               ! dissipation horizontale et verticale des petites echelles:            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(nqmx, 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               ! calcul de l'energie cinetique avant dissipation         CALL pression(ip1jmp1, ap, bp, ps, p3d)
180               call covcont(llm, ucov, vcov, ucont, vcont)         CALL exner_hyb(ps, p3d, pks, pk, pkf)
181               call enercin(vcov, ucov, vcont, ucont, ecin0)  
182           IF (MOD(itau + 1, idissip) == 0) THEN
183               ! dissipation            ! dissipation horizontale et verticale des petites echelles:
184               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)  
185               ucov=ucov + dudis            ! calcul de l'energie cinetique avant dissipation
186               vcov=vcov + dvdis            call covcont(llm, ucov, vcov, ucont, vcont)
187              call enercin(vcov, ucov, vcont, ucont, ecin0)
188               ! On rajoute la tendance due à la transformation Ec -> E  
189               ! thermique créée lors de la dissipation            ! dissipation
190               call covcont(llm, ucov, vcov, ucont, vcont)            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
191               call enercin(vcov, ucov, vcont, ucont, ecin)            ucov=ucov + dudis
192               dtetaecdt= (ecin0 - ecin) / pk            vcov=vcov + dvdis
193               dtetadis=dtetadis + dtetaecdt  
194               teta=teta + dtetadis            ! On rajoute la tendance due à la transformation Ec -> E
195              ! thermique créée lors de la dissipation
196               ! Calcul de la valeur moyenne aux pôles :            call covcont(llm, ucov, vcov, ucont, vcont)
197               forall (l = 1: llm)            call enercin(vcov, ucov, vcont, ucont, ecin)
198                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &            dtetaecdt= (ecin0 - ecin) / pk
199                       / apoln            dtetadis=dtetadis + dtetaecdt
200                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &            teta=teta + dtetadis
201                       * teta(:iim, jjm + 1, l)) / apols  
202               END forall            ! Calcul de la valeur moyenne aux pôles :
203              forall (l = 1: llm)
204               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
205               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &                    / apoln
206                    / apols               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
207            END IF                    * teta(:iim, jjm + 1, l)) / apols
208              END forall
209            itau = itau + 1  
210            iday = day_ini + itau / day_step            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
211            time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
212            IF (time > 1.) THEN                 / apols
213               time = time - 1.         END IF
214               iday = iday + 1  
215            ENDIF         IF (MOD(itau + 1, iperiod) == 0) THEN
216              ! ecriture du fichier histoire moyenne:
217            IF (MOD(itau, iperiod) == 0) THEN            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
218               ! ecriture du fichier histoire moyenne:                 phi, q, masse, ps, phis)
219               CALL writedynav(histaveid, nqmx, itau, vcov, &            call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
220                    ucov, teta, pk, phi, q, masse, ps, phis)                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
221               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &         ENDIF
222                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)      end do
           ENDIF  
        end do leapfrog_loop  
     end do period_loop  
223    
     ! {itau == itaufin}  
224      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
225           itau=itau_dyn+itaufin)           itau=itau_dyn+itaufin)
226    
# Line 267  contains Line 228  contains
228      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
229      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
230           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
231           time + iday - day_ini)           time_0)
232    
233    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
234    
Legend:
Removed from v.32  
changed lines
  Added in v.33
  ViewVC Help
Powered by ViewVC 1.1.21