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