/[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 36 by guez, Thu Dec 2 17:11:04 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 calfis_m, ONLY: calfis      USE calfis_m, ONLY: calfis
14      USE com_io_dyn, ONLY: histaveid      USE com_io_dyn, ONLY: histaveid
# Line 32  contains Line 32  contains
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)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
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              CALL pression(ip1jmp1, ap, bp, ps, p3d)
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.  
   
              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  
160    
161            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
162            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
163    
164            IF (MOD(itau + 1, idissip) == 0) THEN            rdaym_ini = itau * dtvr / daysec
165               ! dissipation horizontale et verticale des petites echelles:            rdayvrai = rdaym_ini + day_ini
166              time = REAL(mod(itau, day_step)) / day_step + time_0
167              IF (time > 1.) time = time - 1.
168    
169              CALL calfis(itau + 1 == itaufin, rdayvrai, time, ucov, vcov, &
170                   teta, q, masse, ps, pk, phis, phi, du, dv, dteta, dq, w, dufi, &
171                   dvfi, dtetafi, dqfi, dpfi)
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)
179           CALL exner_hyb(ps, p3d, pks, pk, pkf)
180    
181               ! calcul de l'energie cinetique avant dissipation         IF (MOD(itau + 1, idissip) == 0) THEN
182               call covcont(llm, ucov, vcov, ucont, vcont)            ! dissipation horizontale et verticale des petites echelles:
183               call enercin(vcov, ucov, vcont, ucont, ecin0)  
184              ! calcul de l'energie cinetique avant dissipation
185               ! dissipation            call covcont(llm, ucov, vcov, ucont, vcont)
186               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            call enercin(vcov, ucov, vcont, ucont, ecin0)
187               ucov=ucov + dudis  
188               vcov=vcov + dvdis            ! dissipation
189              CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
190               ! On rajoute la tendance due à la transformation Ec -> E            ucov=ucov + dudis
191               ! thermique créée lors de la dissipation            vcov=vcov + dvdis
192               call covcont(llm, ucov, vcov, ucont, vcont)  
193               call enercin(vcov, ucov, vcont, ucont, ecin)            ! On rajoute la tendance due à la transformation Ec -> E
194               dtetaecdt= (ecin0 - ecin) / pk            ! thermique créée lors de la dissipation
195               dtetadis=dtetadis + dtetaecdt            call covcont(llm, ucov, vcov, ucont, vcont)
196               teta=teta + dtetadis            call enercin(vcov, ucov, vcont, ucont, ecin)
197              dtetaecdt= (ecin0 - ecin) / pk
198               ! Calcul de la valeur moyenne aux pôles :            dtetadis=dtetadis + dtetaecdt
199               forall (l = 1: llm)            teta=teta + dtetadis
200                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &  
201                       / apoln            ! Calcul de la valeur moyenne aux pôles :
202                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &            forall (l = 1: llm)
203                       * teta(:iim, jjm + 1, l)) / apols               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
204               END forall                    / apoln
205                 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
206               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln                    * teta(:iim, jjm + 1, l)) / apols
207               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &            END forall
208                    / apols  
209            END IF            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)) &
211            itau = itau + 1                 / apols
212            iday = day_ini + itau / day_step         END IF
213            time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0  
214            IF (time > 1.) THEN         IF (MOD(itau + 1, iperiod) == 0) THEN
215               time = time - 1.            ! ecriture du fichier histoire moyenne:
216               iday = iday + 1            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
217            ENDIF                 phi, q, masse, ps, phis)
218              call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
219            IF (MOD(itau, iperiod) == 0) THEN                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
220               ! ecriture du fichier histoire moyenne:         ENDIF
221               CALL writedynav(histaveid, nqmx, itau, vcov, &      end do
                   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.36
  ViewVC Help
Powered by ViewVC 1.1.21