/[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 30 by guez, Thu Apr 1 09:07:28 2010 UTC revision 41 by guez, Tue Feb 22 15:09:57 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 24  contains Line 27  contains
27      use filtreg_m, only: filtreg      use filtreg_m, only: filtreg
28      USE guide_m, ONLY: guide      USE guide_m, ONLY: guide
29      use inidissip_m, only: idissip      use inidissip_m, only: idissip
30        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:
37      REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm) ! vents covariants      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, intent(inout):: 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 74  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 89  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 forward, leapf      logical leapf
96      REAL dt      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      outer_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              dq, dp, vcov, ucov, teta, q, ps, masse, phis, 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)  
145    
146         ! integrations dynamique et traceurs:         if (.not. leapf) then
147         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &            ! Matsuno backward
148              dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., &            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
149              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:  
150    
        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, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 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           forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
183           CALL exner_hyb(ps, p3d, pks, pk, pkf)
184    
185           IF (MOD(itau + 1, idissip) == 0) THEN
186              ! dissipation horizontale et verticale des petites echelles:
187    
188               ! calcul de l'energie cinetique avant dissipation            ! calcul de l'energie cinetique avant dissipation
189               call covcont(llm, ucov, vcov, ucont, vcont)            call covcont(llm, ucov, vcov, ucont, vcont)
190               call enercin(vcov, ucov, vcont, ucont, ecin0)            call enercin(vcov, ucov, vcont, ucont, ecin0)
191    
192               ! dissipation            ! dissipation
193               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
194               ucov=ucov + dudis            ucov=ucov + dudis
195               vcov=vcov + dvdis            vcov=vcov + dvdis
196    
197               ! On rajoute la tendance due à la transformation Ec -> E            ! On rajoute la tendance due à la transformation Ec -> E
198               ! thermique créée lors de la dissipation            ! thermique créée lors de la dissipation
199               call covcont(llm, ucov, vcov, ucont, vcont)            call covcont(llm, ucov, vcov, ucont, vcont)
200               call enercin(vcov, ucov, vcont, ucont, ecin)            call enercin(vcov, ucov, vcont, ucont, ecin)
201               dtetaecdt= (ecin0 - ecin) / pk            dtetaecdt= (ecin0 - ecin) / pk
202               dtetadis=dtetadis + dtetaecdt            dtetadis=dtetadis + dtetaecdt
203               teta=teta + dtetadis            teta=teta + dtetadis
204    
205               ! Calcul de la valeur moyenne unique de h aux pôles            ! Calcul de la valeur moyenne aux pôles :
206               forall (l = 1: llm)            forall (l = 1: llm)
207                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
208                       / apoln                    / apoln
209                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
210                       * teta(:iim, jjm + 1, l)) / apols                    * teta(:iim, jjm + 1, l)) / apols
211               END forall            END forall
212    
213               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
214               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
215                    / apols                 / apols
216            END IF         END IF
217    
218            itau = itau + 1         IF (MOD(itau + 1, iperiod) == 0) THEN
219            iday = day_ini + itau / day_step            ! Écriture du fichier histoire moyenne:
220            time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
221            IF (time > 1.) THEN                 phi, q, masse, ps, phis)
222               time = time - 1.            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
223               iday = iday + 1                 q(:, :, :, 1), dt_app = dtvr * iperiod, &
224            ENDIF                 dt_cum = dtvr * day_step * periodav)
225           ENDIF
226            IF (MOD(itau, iperiod) == 0) THEN      end do time_integration
              ! ecriture du fichier histoire moyenne:  
              CALL writedynav(histaveid, nqmx, itau, vcov, &  
                   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  
        end do  
     end do outer_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    
     vcovm1 = vcov  
     ucovm1 = ucov  
     tetam1 = teta  
     massem1 = masse  
     psm1 = ps  
     finvmaold = masse  
     CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)  
   
231      ! Calcul des tendances dynamiques:      ! Calcul des tendances dynamiques:
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)
   
     ! 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, &  
          itaufin)  
   
     ! integrations dynamique et traceurs:  
     CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, dq, &  
          dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., dtvr)  
       
     CALL pression(ip1jmp1, ap, bp, ps, p3d)  
     CALL exner_hyb(ps, p3d, pks, pk, pkf)  
236    
237    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
238    
Legend:
Removed from v.30  
changed lines
  Added in v.41
  ViewVC Help
Powered by ViewVC 1.1.21