/[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 39 by guez, Tue Jan 25 15:11:05 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 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(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, dp, &
140              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)  
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(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
154                 dteta, dp, vcov, ucov, teta, q, ps, masse, &                 dp, vcov, ucov, teta, q(:, :, :2), ps, masse, finvmaold, dtvr, &
155                 finvmaold, .true., 2 * dtvr)                 leapf=.false.)
156           end if
157            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN  
158               ! calcul des tendances physiques:         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
159               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  
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)
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(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               ! calcul de l'energie cinetique avant dissipation         forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
179               call covcont(llm, ucov, vcov, ucont, vcont)         CALL exner_hyb(ps, p3d, pks, pk, pkf)
180               call enercin(vcov, ucov, vcont, ucont, ecin0)  
181           IF (MOD(itau + 1, idissip) == 0) THEN
182               ! dissipation            ! dissipation horizontale et verticale des petites echelles:
183               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)  
184               ucov=ucov + dudis            ! calcul de l'energie cinetique avant dissipation
185               vcov=vcov + dvdis            call covcont(llm, ucov, vcov, ucont, vcont)
186              call enercin(vcov, ucov, vcont, ucont, ecin0)
187               ! On rajoute la tendance due à la transformation Ec -> E  
188               ! thermique créée lors de la dissipation            ! dissipation
189               call covcont(llm, ucov, vcov, ucont, vcont)            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
190               call enercin(vcov, ucov, vcont, ucont, ecin)            ucov=ucov + dudis
191               dtetaecdt= (ecin0 - ecin) / pk            vcov=vcov + dvdis
192               dtetadis=dtetadis + dtetaecdt  
193               teta=teta + dtetadis            ! On rajoute la tendance due à la transformation Ec -> E
194              ! thermique créée lors de la dissipation
195               ! Calcul de la valeur moyenne aux pôles :            call covcont(llm, ucov, vcov, ucont, vcont)
196               forall (l = 1: llm)            call enercin(vcov, ucov, vcont, ucont, ecin)
197                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &            dtetaecdt= (ecin0 - ecin) / pk
198                       / apoln            dtetadis=dtetadis + dtetaecdt
199                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &            teta=teta + dtetadis
200                       * teta(:iim, jjm + 1, l)) / apols  
201               END forall            ! Calcul de la valeur moyenne aux pôles :
202              forall (l = 1: llm)
203               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
204               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &                    / apoln
205                    / apols               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
206            END IF                    * teta(:iim, jjm + 1, l)) / apols
207              END forall
208            itau = itau + 1  
209            iday = day_ini + itau / day_step            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
210            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)) &
211            IF (time > 1.) THEN                 / apols
212               time = time - 1.         END IF
213               iday = iday + 1  
214            ENDIF         IF (MOD(itau + 1, iperiod) == 0) THEN
215              ! ecriture du fichier histoire moyenne:
216            IF (MOD(itau, iperiod) == 0) THEN            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
217               ! ecriture du fichier histoire moyenne:                 phi, q, masse, ps, phis)
218               CALL writedynav(histaveid, nqmx, itau, vcov, &            call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, ps, &
219                    ucov, teta, pk, phi, q, masse, ps, phis)                 masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
220               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &         ENDIF
221                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)      end do
           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.39
  ViewVC Help
Powered by ViewVC 1.1.21