/[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 26 by guez, Tue Mar 9 15:27:15 2010 UTC revision 31 by guez, Thu Apr 1 14:59:19 2010 UTC
# Line 6  contains Line 6  contains
6    
7    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
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      ! Auteurs: P. Le Van, L. Fairhead, F. Hourdin      ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11        ! schema matsuno + leapfrog
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
15      USE comconst, ONLY: daysec, dtphys, dtvr      USE comconst, ONLY: daysec, dtphys, dtvr
16      USE comgeom, ONLY: aire, apoln, apols      USE comgeom, ONLY: aire_2d, apoln, apols
17      USE comvert, ONLY: ap, bp      USE comvert, ONLY: ap, bp
18      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, &      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
19           nday, offline, periodav           periodav
20      USE dimens_m, ONLY: iim, llm, nqmx      USE dimens_m, ONLY: iim, jjm, llm, nqmx
21      USE dynetat0_m, ONLY: day_ini      USE dynetat0_m, ONLY: day_ini
22        use dynredem1_m, only: dynredem1
23      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
24        use filtreg_m, only: filtreg
25      USE guide_m, ONLY: guide      USE guide_m, ONLY: guide
26      use inidissip_m, only: idissip      use inidissip_m, only: idissip
27      USE logic, ONLY: iflag_phys, ok_guide      USE logic, ONLY: iflag_phys, ok_guide
28      USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1      USE paramet_m, ONLY: ip1jmp1
29      USE pression_m, ONLY: pression      USE pression_m, ONLY: pression
30      USE pressure_var, ONLY: p3d      USE pressure_var, ONLY: p3d
31      USE temps, ONLY: dt, itaufin      USE temps, ONLY: itau_dyn
32    
33      ! Variables dynamiques:      ! Variables dynamiques:
34      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm) ! vents covariants
35      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature
36      REAL ps(ip1jmp1) ! pression au sol, en Pa      REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa
37    
38      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
39      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
# Line 42  contains Line 45  contains
45      ! Variables dynamiques:      ! Variables dynamiques:
46    
47      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
48      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
49      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
50      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(ip1jmp1, llm) ! geopotential
51      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
52    
53      ! variables dynamiques intermediaire pour le transport      ! variables dynamiques intermediaire pour le transport
54      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse      REAL pbaru(ip1jmp1, llm), pbarv((iim + 1) * jjm, llm) !flux de masse
55    
56      ! variables dynamiques au pas - 1      ! variables dynamiques au pas - 1
57      REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)      REAL vcovm1((iim + 1) * jjm, llm), ucovm1(ip1jmp1, llm)
58      REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)      REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
59      REAL massem1(ip1jmp1, llm)      REAL massem1(ip1jmp1, llm)
60    
61      ! tendances dynamiques      ! tendances dynamiques
62      REAL dv(ip1jm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)
63      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
64    
65      ! tendances de la dissipation      ! tendances de la dissipation
66      REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
67      REAL dtetadis(ip1jmp1, llm)      REAL dtetadis(iim + 1, jjm + 1, llm)
68    
69      ! tendances physiques      ! tendances physiques
70      REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
71      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
72    
73      ! variables pour le fichier histoire      ! variables pour le fichier histoire
74    
     REAL tppn(iim), tpps(iim), tpn, tps  
   
75      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
76        INTEGER itaufin
77      INTEGER iday ! jour julien      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.      LOGICAL:: lafin=.false.
81      INTEGER ij, l      INTEGER i, j, l
82    
83      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
84    
85      ! Variables test conservation energie      ! Variables test conservation energie
86      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
87      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
88      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
89      ! cree par la dissipation      ! cree par la dissipation
90      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
91      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
92      CHARACTER*15 ztit      logical forward, leapf
93      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.      REAL dt
   
     logical:: dissip_conservative = .true.  
     logical forward, leapf, apphys, conser, apdiss  
94    
95      !---------------------------------------------------      !---------------------------------------------------
96    
97      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
98    
99      itaufin = nday * day_step      itaufin = nday * day_step
100        ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
101    
102      itau = 0      itau = 0
103      iday = day_ini      iday = day_ini
104      time = time_0      time = time_0
# Line 106  contains Line 107  contains
107      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
108      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
109    
110      ! Debut de l'integration temporelle:      ! Début de l'integration temporelle :
111      outer_loop:do      outer_loop:do i = 1, itaufin / iperiod
112           ! {itau is a multiple of iperiod}
113    
114           ! 1. Matsuno forward:
115    
116         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
117              call guide(itau, ucov, vcov, teta, q, masse, ps)              call guide(itau, ucov, vcov, teta, q, masse, ps)
118         vcovm1 = vcov         vcovm1 = vcov
# Line 115  contains Line 120  contains
120         tetam1 = teta         tetam1 = teta
121         massem1 = masse         massem1 = masse
122         psm1 = ps         psm1 = ps
        forward = .TRUE.  
        leapf = .FALSE.  
        dt = dtvr  
123         finvmaold = masse         finvmaold = masse
124         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)         CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
125    
126         do         ! Calcul des tendances dynamiques:
127            ! gestion des appels de la physique et des dissipations:         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
128            apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
129            conser = MOD(itau, iconser) == 0              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
130            apdiss = MOD(itau + 1, idissip) == 0              time + iday - day_ini)
131    
132           ! Calcul des tendances advection des traceurs (dont l'humidité)
133           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134           ! Stokage du flux de masse pour traceurs offline:
135           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
136                dtvr, itau)
137    
138           ! integrations dynamique et traceurs:
139           CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
140                dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., &
141                dtvr)
142    
143           CALL pression(ip1jmp1, ap, bp, ps, p3d)
144           CALL exner_hyb(ps, p3d, pks, pk, pkf)
145    
146           ! 2. Matsuno backward:
147    
148           itau = itau + 1
149           iday = day_ini + itau / day_step
150           time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
151           IF (time > 1.) THEN
152              time = time - 1.
153              iday = iday + 1
154           ENDIF
155    
156           ! Calcul des tendances dynamiques:
157           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
158           CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
159                .false., du, dv, dteta, dp, w, pbaru, pbarv, time + iday - day_ini)
160    
161           ! integrations dynamique et traceurs:
162           CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
163                dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., &
164                dtvr)
165    
166            ! calcul des tendances dynamiques:         CALL pression(ip1jmp1, ap, bp, ps, p3d)
167           CALL exner_hyb(ps, p3d, pks, pk, pkf)
168    
169           ! 3. Leapfrog:
170    
171           do j = 1, iperiod - 1
172              ! Calcul des tendances dynamiques:
173            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
174            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
175                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &                 .false., du, dv, dteta, dp, w, pbaru, pbarv, &
176                 time + iday - day_ini)                 time + iday - day_ini)
177    
178            IF (forward .OR. leapf) THEN            ! Calcul des tendances advection des traceurs (dont l'humidité)
179               ! calcul des tendances advection des traceurs (dont l'humidite)            CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
180               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)            ! Stokage du flux de masse pour traceurs off-line:
181               IF (offline) THEN            IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
182                  ! Stokage du flux de masse pour traceurs off-line                 dtvr, itau)
                 CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
                      itau)  
              ENDIF  
           ENDIF  
183    
184            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
185            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
186                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
187                 finvmaold, leapf)                 finvmaold, .true., 2 * dtvr)
188    
189            IF (apphys) THEN            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
190               ! calcul des tendances physiques:               ! calcul des tendances physiques:
191               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
192    
# Line 158  contains Line 196  contains
196               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
197               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
198    
              ! Diagnostique de conservation de l'énergie : initialisation  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit='bil dyn'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
   
199               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
200                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
201                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
202    
203               ! ajout des tendances physiques:               ! ajout des tendances physiques:
204               CALL addfi(nqmx, dtphys, &               CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
205                    ucov, vcov, teta, q, ps, &                    dtetafi, dqfi, dpfi)
                   dufi, dvfi, dtetafi, dqfi, dpfi)  
   
              ! Diagnostique de conservation de l'énergie : difference  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit = 'bil phys'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
206            ENDIF            ENDIF
207    
208            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
209            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
210    
211            IF (apdiss) THEN            IF (MOD(itau + 1, idissip) == 0) THEN
212               ! dissipation horizontale et verticale des petites echelles:               ! dissipation horizontale et verticale des petites echelles:
213    
214               ! calcul de l'energie cinetique avant dissipation               ! calcul de l'energie cinetique avant dissipation
# Line 197  contains Line 220  contains
220               ucov=ucov + dudis               ucov=ucov + dudis
221               vcov=vcov + dvdis               vcov=vcov + dvdis
222    
223               if (dissip_conservative) then               ! On rajoute la tendance due à la transformation Ec -> E
224                  ! On rajoute la tendance due a la transform. Ec -> E               ! thermique créée lors de la dissipation
225                  ! therm. cree lors de la dissipation               call covcont(llm, ucov, vcov, ucont, vcont)
226                  call covcont(llm, ucov, vcov, ucont, vcont)               call enercin(vcov, ucov, vcont, ucont, ecin)
227                  call enercin(vcov, ucov, vcont, ucont, ecin)               dtetaecdt= (ecin0 - ecin) / pk
228                  dtetaecdt= (ecin0 - ecin) / pk               dtetadis=dtetadis + dtetaecdt
                 dtetadis=dtetadis + dtetaecdt  
              endif  
229               teta=teta + dtetadis               teta=teta + dtetadis
230    
231               ! Calcul de la valeur moyenne, unique de h aux poles .....               ! Calcul de la valeur moyenne unique de h aux pôles
232               DO l = 1, llm               forall (l = 1: llm)
233                  DO ij = 1, iim                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
234                     tppn(ij) = aire(ij) * teta(ij, l)                       / apoln
235                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
236                  ENDDO                       * teta(:iim, jjm + 1, l)) / apols
237                  tpn = SUM(tppn) / apoln               END forall
238                  tps = SUM(tpps) / apols  
239                 ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
240                  DO ij = 1, iip1               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
241                     teta(ij, l) = tpn                    / apols
                    teta(ij + ip1jm, l) = tps  
                 ENDDO  
              ENDDO  
   
              DO ij = 1, iim  
                 tppn(ij) = aire(ij) * ps(ij)  
                 tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)  
              ENDDO  
              tpn = SUM(tppn) / apoln  
              tps = SUM(tpps) / apols  
   
              DO ij = 1, iip1  
                 ps(ij) = tpn  
                 ps(ij + ip1jm) = tps  
              ENDDO  
242            END IF            END IF
243    
244            ! fin de l'intégration dynamique et physique pour le pas "itau"            itau = itau + 1
245            ! préparation du pas d'intégration suivant            iday = day_ini + itau / day_step
246              time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
247            ! schema matsuno + leapfrog            IF (time > 1.) THEN
248            IF (forward .OR. leapf) THEN               time = time - 1.
249               itau = itau + 1               iday = iday + 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  
250            ENDIF            ENDIF
251    
252            IF (itau == itaufin + 1) exit outer_loop            IF (MOD(itau, iperiod) == 0) THEN
   
           IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN  
253               ! ecriture du fichier histoire moyenne:               ! ecriture du fichier histoire moyenne:
254               CALL writedynav(histaveid, nqmx, itau, vcov, &               CALL writedynav(histaveid, nqmx, itau, vcov, &
255                    ucov, teta, pk, phi, q, masse, ps, phis)                    ucov, teta, pk, phi, q, masse, ps, phis)
256               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
257                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
258            ENDIF            ENDIF
   
           IF (itau == itaufin) THEN  
              CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)  
           ENDIF  
   
           ! gestion de l'integration temporelle:  
           IF (MOD(itau, iperiod) == 0) exit  
           IF (MOD(itau - 1, iperiod) == 0) THEN  
              IF (forward) THEN  
                 ! fin du pas forward et debut du pas backward  
                 forward = .FALSE.  
                 leapf = .FALSE.  
              ELSE  
                 ! fin du pas backward et debut du premier pas leapfrog  
                 leapf = .TRUE.  
                 dt = 2. * dtvr  
              END IF  
           ELSE  
              ! pas leapfrog  
              leapf = .TRUE.  
              dt = 2. * dtvr  
           END IF  
259         end do         end do
260      end do outer_loop      end do outer_loop
261    
262        ! {itau == itaufin}
263        CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
264             itau=itau_dyn+itaufin)
265    
266        vcovm1 = vcov
267        ucovm1 = ucov
268        tetam1 = teta
269        massem1 = masse
270        psm1 = ps
271        finvmaold = masse
272        CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
273    
274        ! Calcul des tendances dynamiques:
275        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
276        CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
277             MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
278             time + iday - day_ini)
279    
280        ! Calcul des tendances advection des traceurs (dont l'humidité)
281        CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
282        ! Stokage du flux de masse pour traceurs off-line:
283        IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
284             itaufin)
285    
286    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
287    
288  end module leapfrog_m  end module leapfrog_m
Legend:
Removed from v.26  
changed lines
  Added in v.31
  ViewVC Help
Powered by ViewVC 1.1.21