/[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 31 by guez, Thu Apr 1 14:59:19 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 24  contains Line 25  contains
25      use filtreg_m, only: filtreg      use filtreg_m, only: filtreg
26      USE guide_m, ONLY: guide      USE guide_m, ONLY: guide
27      use inidissip_m, only: idissip      use inidissip_m, only: idissip
28        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:
35      REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm) ! vents covariants      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 74  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 89  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 forward, leapf      logical leapf
91      REAL dt      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      outer_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              dq, dp, vcov, ucov, teta, q, ps, masse, phis, 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, &  
             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)  
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    
        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, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
155                 finvmaold, .true., 2 * dtvr)                 dtvr)
156           end if
           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN  
              ! 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)  
157    
158               rdaym_ini = itau * dtvr / daysec         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
159               rdayvrai = rdaym_ini + day_ini            ! calcul des tendances physiques:
160    
161               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
162                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &            CALL exner_hyb(ps, p3d, pks, pk, pkf)
                   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 unique de h 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  
     end do outer_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    
     vcovm1 = vcov  
     ucovm1 = ucov  
     tetam1 = teta  
     massem1 = masse  
     psm1 = ps  
     finvmaold = masse  
     CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)  
   
226      ! Calcul des tendances dynamiques:      ! Calcul des tendances dynamiques:
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)
   
     ! 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)  
231    
232    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
233    
Legend:
Removed from v.31  
changed lines
  Added in v.37
  ViewVC Help
Powered by ViewVC 1.1.21