/[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 32 by guez, Tue Apr 6 17:52:58 2010 UTC revision 45 by guez, Wed Apr 27 13:00:12 2011 UTC
# Line 8  contains Line 8  contains
8    
9      ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34      ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34
10      ! Authors: P. Le Van, L. Fairhead, F. Hourdin      ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11      ! schema matsuno + leapfrog      ! Matsuno-leapfrog scheme.
12    
13        use addfi_m, only: addfi
14        use bilan_dyn_m, only: bilan_dyn
15        use caladvtrac_m, only: caladvtrac
16        use caldyn_m, only: caldyn
17      USE calfis_m, ONLY: calfis      USE calfis_m, ONLY: calfis
18      USE com_io_dyn, ONLY: histaveid      USE com_io_dyn, ONLY: histaveid
19      USE comconst, ONLY: daysec, dtphys, dtvr      USE comconst, ONLY: daysec, dtphys, dtvr
# Line 22  contains Line 26  contains
26      use dynredem1_m, only: dynredem1      use dynredem1_m, only: dynredem1
27      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
28      use filtreg_m, only: filtreg      use filtreg_m, only: filtreg
29        use geopot_m, only: geopot
30      USE guide_m, ONLY: guide      USE guide_m, ONLY: guide
31      use inidissip_m, only: idissip      use inidissip_m, only: idissip
32      use integrd_m, only: integrd      use integrd_m, only: integrd
33      USE logic, ONLY: iflag_phys, ok_guide      USE logic, ONLY: iflag_phys, ok_guide
34      USE paramet_m, ONLY: ip1jmp1      USE paramet_m, ONLY: ip1jmp1
     USE pression_m, ONLY: pression  
35      USE pressure_var, ONLY: p3d      USE pressure_var, ONLY: p3d
36      USE temps, ONLY: itau_dyn      USE temps, ONLY: itau_dyn
37    
38      ! Variables dynamiques:      ! Variables dynamiques:
     REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant  
39      REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant      REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
40      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature      REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
     REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa  
41    
42        REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
43        ! potential temperature
44    
45        REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
46      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
47      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
48      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
49        REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
50        ! mass fractions of advected fields
51    
52      REAL, intent(in):: time_0      REAL, intent(in):: time_0
53    
54      ! Variables local to the procedure:      ! Variables local to the procedure:
# Line 62  contains Line 71  contains
71    
72      ! tendances dynamiques      ! tendances dynamiques
73      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)
74      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(iim + 1, jjm + 1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
75    
76      ! tendances de la dissipation      ! tendances de la dissipation
77      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
# Line 70  contains Line 79  contains
79    
80      ! tendances physiques      ! tendances physiques
81      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
82      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
83    
84      ! variables pour le fichier histoire      ! variables pour le fichier histoire
85    
86      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
87      INTEGER itaufin      INTEGER itaufin
     INTEGER iday ! jour julien  
88      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
89      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
90      LOGICAL:: lafin=.false.      INTEGER l
     INTEGER i, j, l  
   
91      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
92    
93      ! Variables test conservation energie      ! Variables test conservation energie
94      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
95      ! Tendance de la temp. potentiel d (theta) / d t due a la  
     ! tansformation d'energie cinetique en energie thermique  
     ! cree par la dissipation  
96      REAL dtetaecdt(iim + 1, jjm + 1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
97        ! tendance de la température potentielle due à la tansformation
98        ! d'énergie cinétique en énergie thermique créée par la dissipation
99    
100      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
101        logical leapf
102        real dt
103    
104      !---------------------------------------------------      !---------------------------------------------------
105    
# Line 99  contains Line 108  contains
108      itaufin = nday * day_step      itaufin = nday * day_step
109      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too      ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
110    
     itau = 0  
     iday = day_ini  
     time = time_0  
111      dq = 0.      dq = 0.
112    
113      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
114      CALL pression(ip1jmp1, ap, bp, ps, p3d)      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
115      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
116    
117      ! Début de l'integration temporelle :      time_integration: do itau = 0, itaufin - 1
118      period_loop:do i = 1, itaufin / iperiod         leapf = mod(itau, iperiod) /= 0
119         ! {"itau" is a multiple of "iperiod"}         if (leapf) then
120              dt = 2 * dtvr
121         ! 1. Matsuno forward:         else
122              ! Matsuno
123         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &            dt = dtvr
124              call guide(itau, ucov, vcov, teta, q, masse, ps)            if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
125         vcovm1 = vcov                 call guide(itau, ucov, vcov, teta, q, masse, ps)
126         ucovm1 = ucov            vcovm1 = vcov
127         tetam1 = teta            ucovm1 = ucov
128         massem1 = masse            tetam1 = teta
129         psm1 = ps            massem1 = masse
130         finvmaold = masse            psm1 = ps
131         CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)            finvmaold = masse
132              CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
133           end if
134    
135         ! Calcul des tendances dynamiques:         ! Calcul des tendances dynamiques:
136         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
137         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
138              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
139              time + iday - day_ini)              time_0)
140    
141         ! Calcul des tendances advection des traceurs (dont l'humidité)         ! Calcul des tendances advection des traceurs (dont l'humidité)
142         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
143    
144         ! Stokage du flux de masse pour traceurs offline:         ! Stokage du flux de masse pour traceurs offline:
145         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
146              dtvr, itau)              dtvr, itau)
147    
148         ! integrations dynamique et traceurs:         ! integrations dynamique et traceurs:
149         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, dp, &
150              dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &              vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, leapf)
             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)  
   
        ! integrations dynamique et traceurs:  
        CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &  
             dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &  
             dtvr)  
   
        CALL pression(ip1jmp1, ap, bp, ps, p3d)  
        CALL exner_hyb(ps, p3d, pks, pk, pkf)  
151    
152         ! 3. Leapfrog:         if (.not. leapf) then
153              ! Matsuno backward
154              forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
155              CALL exner_hyb(ps, p3d, pks, pk, pkf)
156    
        leapfrog_loop: do j = 1, iperiod - 1  
157            ! Calcul des tendances dynamiques:            ! Calcul des tendances dynamiques:
158            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
159            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
160                 .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)  
161    
162            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
163            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
164                 dteta, dp, vcov, ucov, teta, q, ps, masse, &                 dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, &
165                 finvmaold, .true., 2 * dtvr)                 dtvr, leapf=.false.)
166           end if
167            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN  
168               ! calcul des tendances physiques:         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
169               IF (itau + 1 == itaufin) lafin = .TRUE.            ! calcul des tendances physiques:
   
              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  
170    
171            CALL pression(ip1jmp1, ap, bp, ps, p3d)            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
172            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
173    
174            IF (MOD(itau + 1, idissip) == 0) THEN            rdaym_ini = itau * dtvr / daysec
175               ! dissipation horizontale et verticale des petites echelles:            rdayvrai = rdaym_ini + day_ini
176              time = REAL(mod(itau, day_step)) / day_step + time_0
177              IF (time > 1.) time = time - 1.
178    
179              CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
180                   phis, phi, du, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
181                   lafin=itau+1==itaufin)
182    
183              ! ajout des tendances physiques:
184              CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
185                   dtetafi, dqfi, dpfi)
186           ENDIF
187    
188               ! calcul de l'energie cinetique avant dissipation         forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
189               call covcont(llm, ucov, vcov, ucont, vcont)         CALL exner_hyb(ps, p3d, pks, pk, pkf)
190               call enercin(vcov, ucov, vcont, ucont, ecin0)  
191           IF (MOD(itau + 1, idissip) == 0) THEN
192               ! dissipation            ! dissipation horizontale et verticale des petites echelles:
193               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)  
194               ucov=ucov + dudis            ! calcul de l'energie cinetique avant dissipation
195               vcov=vcov + dvdis            call covcont(llm, ucov, vcov, ucont, vcont)
196              call enercin(vcov, ucov, vcont, ucont, ecin0)
197               ! On rajoute la tendance due à la transformation Ec -> E  
198               ! thermique créée lors de la dissipation            ! dissipation
199               call covcont(llm, ucov, vcov, ucont, vcont)            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
200               call enercin(vcov, ucov, vcont, ucont, ecin)            ucov=ucov + dudis
201               dtetaecdt= (ecin0 - ecin) / pk            vcov=vcov + dvdis
202               dtetadis=dtetadis + dtetaecdt  
203               teta=teta + dtetadis            ! On rajoute la tendance due à la transformation Ec -> E
204              ! thermique créée lors de la dissipation
205               ! Calcul de la valeur moyenne aux pôles :            call covcont(llm, ucov, vcov, ucont, vcont)
206               forall (l = 1: llm)            call enercin(vcov, ucov, vcont, ucont, ecin)
207                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &            dtetaecdt= (ecin0 - ecin) / pk
208                       / apoln            dtetadis=dtetadis + dtetaecdt
209                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &            teta=teta + dtetadis
210                       * teta(:iim, jjm + 1, l)) / apols  
211               END forall            ! Calcul de la valeur moyenne aux pôles :
212              forall (l = 1: llm)
213               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
214               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &                    / apoln
215                    / apols               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
216            END IF                    * teta(:iim, jjm + 1, l)) / apols
217              END forall
218            itau = itau + 1  
219            iday = day_ini + itau / day_step            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
220            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)) &
221            IF (time > 1.) THEN                 / apols
222               time = time - 1.         END IF
223               iday = iday + 1  
224            ENDIF         IF (MOD(itau + 1, iperiod) == 0) THEN
225              ! Écriture du fichier histoire moyenne:
226            IF (MOD(itau, iperiod) == 0) THEN            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
227               ! ecriture du fichier histoire moyenne:                 phi, q, masse, ps, phis)
228               CALL writedynav(histaveid, nqmx, itau, vcov, &            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
229                    ucov, teta, pk, phi, q, masse, ps, phis)                 q(:, :, :, 1), dt_app = dtvr * iperiod, &
230               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &                 dt_cum = dtvr * day_step * periodav)
231                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)         ENDIF
232            ENDIF      end do time_integration
        end do leapfrog_loop  
     end do period_loop  
233    
     ! {itau == itaufin}  
234      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
235           itau=itau_dyn+itaufin)           itau=itau_dyn+itaufin)
236    
# Line 267  contains Line 238  contains
238      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
239      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
240           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
241           time + iday - day_ini)           time_0)
242    
243    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
244    
Legend:
Removed from v.32  
changed lines
  Added in v.45
  ViewVC Help
Powered by ViewVC 1.1.21