/[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 28 by guez, Fri Mar 26 18:33:04 2010 UTC revision 32 by guez, Tue Apr 6 17:52:58 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
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, nday, offline, &      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
19           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      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      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 integrd_m, only: integrd
28      USE logic, ONLY: iflag_phys, ok_guide      USE logic, ONLY: iflag_phys, ok_guide
29      USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1      USE paramet_m, ONLY: ip1jmp1
30      USE pression_m, ONLY: pression      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(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
36      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
37      REAL ps(ip1jmp1) ! pression au sol, en Pa      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature
38        REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa
39    
40      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
41      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
# Line 44  contains Line 47  contains
47      ! Variables dynamiques:      ! Variables dynamiques:
48    
49      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
50      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
51      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
52      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(ip1jmp1, llm) ! geopotential
53      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
54    
55      ! variables dynamiques intermediaire pour le transport      ! variables dynamiques intermediaire pour le transport
56      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse      REAL pbaru(ip1jmp1, llm), pbarv((iim + 1) * jjm, llm) !flux de masse
57    
58      ! variables dynamiques au pas - 1      ! variables dynamiques au pas - 1
59      REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)      REAL vcovm1((iim + 1) * jjm, llm), ucovm1(ip1jmp1, llm)
60      REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)      REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
61      REAL massem1(ip1jmp1, llm)      REAL massem1(ip1jmp1, llm)
62    
63      ! tendances dynamiques      ! tendances dynamiques
64      REAL dv(ip1jm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)
65      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
66    
67      ! tendances de la dissipation      ! tendances de la dissipation
68      REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
69      REAL dtetadis(ip1jmp1, llm)      REAL dtetadis(iim + 1, jjm + 1, llm)
70    
71      ! tendances physiques      ! tendances physiques
72      REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
73      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
74    
75      ! variables pour le fichier histoire      ! variables pour le fichier histoire
76    
     REAL tppn(iim), tpps(iim), tpn, tps  
   
77      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
78      INTEGER itaufin      INTEGER itaufin
79      INTEGER iday ! jour julien      INTEGER iday ! jour julien
80      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
81      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
82      LOGICAL:: lafin=.false.      LOGICAL:: lafin=.false.
83      INTEGER ij, l      INTEGER i, j, l
84    
85      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
86    
87      ! Variables test conservation energie      ! Variables test conservation energie
88      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
89      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
90      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
91      ! cree par la dissipation      ! cree par la dissipation
92      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
93      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
     logical forward, leapf  
     REAL dt  
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 107  contains Line 108  contains
108      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
109    
110      ! Début de l'integration temporelle :      ! Début de l'integration temporelle :
111      outer_loop:do      period_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           ! Calcul des tendances dynamiques:
127           CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
128           CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
129                MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
130                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                dp, vcov, ucov, teta, q, ps, masse, 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                dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
164                dtvr)
165    
166           CALL pression(ip1jmp1, ap, bp, ps, p3d)
167           CALL exner_hyb(ps, p3d, pks, pk, pkf)
168    
169           ! 3. Leapfrog:
170    
171         do         leapfrog_loop: do j = 1, iperiod - 1
172            ! Calcul des tendances dynamiques:            ! 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                 MOD(itau, iconser) == 0, 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'humidité)            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, dp, vcov, ucov, teta, q, ps, masse, &
187                 finvmaold, leapf, dt)                 finvmaold, .true., 2 * dtvr)
188    
189            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
190               ! calcul des tendances physiques:               ! calcul des tendances physiques:
# Line 158  contains Line 201  contains
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)  
206            ENDIF            ENDIF
207    
208            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
# Line 186  contains Line 228  contains
228               dtetadis=dtetadis + dtetaecdt               dtetadis=dtetadis + dtetaecdt
229               teta=teta + dtetadis               teta=teta + dtetadis
230    
231               ! Calcul de la valeur moyenne unique de h aux pôles               ! Calcul de la valeur moyenne 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
259           end do leapfrog_loop
260        end do period_loop
261    
262            IF (itau == itaufin) THEN      ! {itau == itaufin}
263               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
264                    itau=itau_dyn+itaufin)           itau=itau_dyn+itaufin)
265            ENDIF  
266        ! Calcul des tendances dynamiques:
267            ! gestion de l'integration temporelle:      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
268            IF (MOD(itau, iperiod) == 0) exit      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
269            IF (MOD(itau - 1, iperiod) == 0) THEN           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
270               IF (forward) THEN           time + iday - day_ini)
                 ! 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  
        end do  
     end do outer_loop  
271    
272    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
273    
Legend:
Removed from v.28  
changed lines
  Added in v.32
  ViewVC Help
Powered by ViewVC 1.1.21