/[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 29 by guez, Tue Mar 30 10:44:42 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
# Line 23  contains Line 24  contains
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: ip1jmp1      USE paramet_m, ONLY: ip1jmp1
30      USE pression_m, ONLY: pression      USE pression_m, ONLY: pression
# Line 30  contains Line 32  contains
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):: vcov((iim + 1) * jjm, llm) ! vent covariant
36        REAL, intent(inout):: ucov(ip1jmp1, 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    
# Line 77  contains Line 80  contains
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 l      INTEGER i, j, l
84    
85      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
86    
# Line 88  contains Line 91  contains
91      ! cree par la dissipation      ! cree par la dissipation
92      REAL dtetaecdt(iim + 1, jjm + 1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
93      REAL vcont((iim + 1) * jjm, 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 105  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 113  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, jjm + 1, llm, - 2, 2, .TRUE., 1)         CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
125    
126         do         ! 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           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 183  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               forall (l = 1: llm)               forall (l = 1: llm)
233                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
234                       / apoln                       / apoln
# Line 196  contains Line 241  contains
241                    / apols                    / apols
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) CALL dynredem1("restart.nc", vcov, ucov, teta, &      ! {itau == itaufin}
263                 q, masse, ps, itau=itau_dyn+itaufin)      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
264             itau=itau_dyn+itaufin)
265            ! gestion de l'integration temporelle:  
266            IF (MOD(itau, iperiod) == 0) exit      ! Calcul des tendances dynamiques:
267            IF (MOD(itau - 1, iperiod) == 0) THEN      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
268               IF (forward) THEN      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
269                  ! fin du pas forward et debut du pas backward           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
270                  forward = .FALSE.           time + iday - day_ini)
                 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.29  
changed lines
  Added in v.32
  ViewVC Help
Powered by ViewVC 1.1.21