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