/[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

trunk/libf/dyn3d/leapfrog.f90 revision 26 by guez, Tue Mar 9 15:27:15 2010 UTC trunk/dyn3d/leapfrog.f revision 128 by guez, Thu Feb 12 16:23:33 2015 UTC
# Line 4  module leapfrog_m Line 4  module leapfrog_m
4    
5  contains  contains
6    
7    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q)
   
     ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34  
     ! Auteurs: P. Le Van, L. Fairhead, F. Hourdin  
8    
9        ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
10        ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11        ! Matsuno-leapfrog scheme.
12    
13        use addfi_m, only: addfi
14        use bilan_dyn_m, only: bilan_dyn
15        use caladvtrac_m, only: caladvtrac
16        use caldyn_m, only: caldyn
17      USE calfis_m, ONLY: calfis      USE calfis_m, ONLY: calfis
18      USE com_io_dyn, ONLY: histaveid      USE comconst, ONLY: daysec, dtvr
19      USE comconst, ONLY: daysec, dtphys, dtvr      USE comgeom, ONLY: aire_2d, apoln, apols
20      USE comgeom, ONLY: aire, apoln, apols      USE disvert_m, ONLY: ap, bp
21      USE comvert, ONLY: ap, bp      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
22      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, &           iflag_phys, iecri
23           nday, offline, periodav      USE conf_guide_m, ONLY: ok_guide
24      USE dimens_m, ONLY: iim, llm, nqmx      USE dimens_m, ONLY: iim, jjm, llm, nqmx
25        use dissip_m, only: dissip
26      USE dynetat0_m, ONLY: day_ini      USE dynetat0_m, ONLY: day_ini
27        use dynredem1_m, only: dynredem1
28      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
29        use filtreg_m, only: filtreg
30        use fluxstokenc_m, only: fluxstokenc
31        use geopot_m, only: geopot
32      USE guide_m, ONLY: guide      USE guide_m, ONLY: guide
33      use inidissip_m, only: idissip      use inidissip_m, only: idissip
34      USE logic, ONLY: iflag_phys, ok_guide      use integrd_m, only: integrd
35      USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1      use nr_util, only: assert
     USE pression_m, ONLY: pression  
36      USE pressure_var, ONLY: p3d      USE pressure_var, ONLY: p3d
37      USE temps, ONLY: dt, itaufin      USE temps, ONLY: itau_dyn
38        use writedynav_m, only: writedynav
39        use writehist_m, only: writehist
40    
41      ! Variables dynamiques:      ! Variables dynamiques:
42      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
43      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
     REAL ps(ip1jmp1) ! pression au sol, en Pa  
   
     REAL masse(ip1jmp1, llm) ! masse d'air  
     REAL phis(ip1jmp1) ! geopotentiel au sol  
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
     REAL, intent(in):: time_0  
44    
45      ! Variables local to the procedure:      REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
46        ! potential temperature
47    
48      ! Variables dynamiques:      REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
49        REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
50        REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
51    
52      REAL pks(ip1jmp1) ! exner au sol      REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
53      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      ! mass fractions of advected fields
     REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches  
     REAL phi(ip1jmp1, llm) ! geopotential  
     REAL w(ip1jmp1, llm) ! vitesse verticale  
   
     ! variables dynamiques intermediaire pour le transport  
     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse  
   
     ! variables dynamiques au pas - 1  
     REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)  
     REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)  
     REAL massem1(ip1jmp1, llm)  
   
     ! tendances dynamiques  
     REAL dv(ip1jm, llm), du(ip1jmp1, llm)  
     REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)  
   
     ! tendances de la dissipation  
     REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)  
     REAL dtetadis(ip1jmp1, llm)  
   
     ! tendances physiques  
     REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)  
     REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)  
54    
55      ! variables pour le fichier histoire      ! Local:
56    
57        ! Variables dynamiques:
58    
59      REAL tppn(iim), tpps(iim), tpn, tps      REAL pks(iim + 1, jjm + 1) ! exner au sol
60        REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
61        REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
62        REAL phi(iim + 1, jjm + 1, llm) ! geopotential
63        REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
64    
65        ! Variables dynamiques intermediaire pour le transport
66        ! Flux de masse :
67        REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
68    
69        ! Variables dynamiques au pas - 1
70        REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
71        REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
72        REAL massem1(iim + 1, jjm + 1, llm)
73    
74        ! Tendances dynamiques
75        REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
76        REAL dteta(iim + 1, jjm + 1, llm)
77        real dp((iim + 1) * (jjm + 1))
78    
79        ! Tendances de la dissipation :
80        REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
81        REAL dtetadis(iim + 1, jjm + 1, llm)
82    
83        ! Tendances physiques
84        REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
85        REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
86    
87        ! Variables pour le fichier histoire
88      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
89      INTEGER iday ! jour julien      INTEGER itaufin
90      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
91      real finvmaold(ip1jmp1, llm)      real finvmaold(iim + 1, jjm + 1, llm)
92      LOGICAL:: lafin=.false.      INTEGER l
     INTEGER ij, l  
   
     REAL rdayvrai, rdaym_ini  
   
     ! Variables test conservation energie  
     REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)  
     ! Tendance de la temp. potentiel d (theta) / d t due a la  
     ! tansformation d'energie cinetique en energie thermique  
     ! cree par la dissipation  
     REAL dtetaecdt(ip1jmp1, llm)  
     REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)  
     CHARACTER*15 ztit  
     INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.  
93    
94      logical:: dissip_conservative = .true.      ! Variables test conservation \'energie
95      logical forward, leapf, apphys, conser, apdiss      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
96    
97        REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
98        logical leapf
99        real dt ! time step, in s
100    
101      !---------------------------------------------------      !---------------------------------------------------
102    
103      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
104        call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
105    
106      itaufin = nday * day_step      itaufin = nday * day_step
107      itau = 0      ! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
108      iday = day_ini  
     time = time_0  
     dq = 0.  
109      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
110      CALL pression(ip1jmp1, ap, bp, ps, p3d)      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
111      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
112    
113      ! Debut de l'integration temporelle:      time_integration: do itau = 0, itaufin - 1
114      outer_loop:do         leapf = mod(itau, iperiod) /= 0
115         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &         if (leapf) then
116              call guide(itau, ucov, vcov, teta, q, masse, ps)            dt = 2 * dtvr
117         vcovm1 = vcov         else
118         ucovm1 = ucov            ! Matsuno
119         tetam1 = teta            dt = dtvr
120         massem1 = masse            if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
121         psm1 = ps            vcovm1 = vcov
122         forward = .TRUE.            ucovm1 = ucov
123         leapf = .FALSE.            tetam1 = teta
124         dt = dtvr            massem1 = masse
125         finvmaold = masse            psm1 = ps
126         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)            finvmaold = masse
127              CALL filtreg(finvmaold, direct = .false., intensive = .false.)
128         do         end if
129            ! gestion des appels de la physique et des dissipations:  
130            apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0         ! Calcul des tendances dynamiques:
131            conser = MOD(itau, iconser) == 0         CALL geopot(teta, pk, pks, phis, phi)
132            apdiss = MOD(itau + 1, idissip) == 0         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
133                dudyn, dv, dteta, dp, w, pbaru, pbarv, &
134            ! calcul des tendances dynamiques:              conser = MOD(itau, iconser) == 0)
135            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)  
136            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
137                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &  
138                 time + iday - day_ini)         ! Stokage du flux de masse pour traceurs offline:
139           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
140            IF (forward .OR. leapf) THEN              dtvr, itau)
141               ! calcul des tendances advection des traceurs (dont l'humidite)  
142               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         ! Int\'egrations dynamique et traceurs:
143               IF (offline) THEN         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
144                  ! Stokage du flux de masse pour traceurs off-line              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
145                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &              leapf)
146                       itau)  
147               ENDIF         forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
148            ENDIF         CALL exner_hyb(ps, p3d, pks, pk, pkf)
149    
150           if (.not. leapf) then
151              ! Matsuno backward
152              ! Calcul des tendances dynamiques:
153              CALL geopot(teta, pk, pks, phis, phi)
154              CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
155                   phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
156    
157            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
158            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
159                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
160                 finvmaold, leapf)                 finvmaold, dtvr, leapf=.false.)
   
           IF (apphys) 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)  
   
              rdaym_ini = itau * dtvr / daysec  
              rdayvrai = rdaym_ini + day_ini  
   
              ! 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  
   
              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)  
   
              ! 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  
           ENDIF  
161    
162            CALL pression(ip1jmp1, ap, bp, ps, p3d)            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
163            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
164           end if
165    
166            IF (apdiss) THEN         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
167               ! dissipation horizontale et verticale des petites echelles:            ! Calcul des tendances physiques :
168              time = REAL(mod(itau, day_step)) / day_step
169               ! calcul de l'energie cinetique avant dissipation            IF (time > 1.) time = time - 1.
170               call covcont(llm, ucov, vcov, ucont, vcont)            CALL calfis(itau * dtvr / daysec + day_ini, time, ucov, vcov, teta, &
171               call enercin(vcov, ucov, vcont, ucont, ecin0)                 q, pk, phis, phi, w, dufi, dvfi, dtetafi, dqfi, &
172                   lafin = itau + 1 == itaufin)
173               ! dissipation  
174               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
175               ucov=ucov + dudis         ENDIF
176               vcov=vcov + dvdis  
177           IF (MOD(itau + 1, idissip) == 0) THEN
178               if (dissip_conservative) then            ! Dissipation horizontale et verticale des petites \'echelles
179                  ! On rajoute la tendance due a la transform. Ec -> E  
180                  ! therm. cree lors de la dissipation            ! calcul de l'\'energie cin\'etique avant dissipation
181                  call covcont(llm, ucov, vcov, ucont, vcont)            call covcont(llm, ucov, vcov, ucont, vcont)
182                  call enercin(vcov, ucov, vcont, ucont, ecin)            call enercin(vcov, ucov, vcont, ucont, ecin0)
183                  dtetaecdt= (ecin0 - ecin) / pk  
184                  dtetadis=dtetadis + dtetaecdt            ! dissipation
185               endif            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
186               teta=teta + dtetadis            ucov = ucov + dudis
187              vcov = vcov + dvdis
188               ! Calcul de la valeur moyenne, unique de h aux poles .....  
189               DO l = 1, llm            ! On ajoute la tendance due \`a la transformation \'energie
190                  DO ij = 1, iim            ! cin\'etique en \'energie thermique par la dissipation
191                     tppn(ij) = aire(ij) * teta(ij, l)            call covcont(llm, ucov, vcov, ucont, vcont)
192                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)            call enercin(vcov, ucov, vcont, ucont, ecin)
193                  ENDDO            dtetadis = dtetadis + (ecin0 - ecin) / pk
194                  tpn = SUM(tppn) / apoln            teta = teta + dtetadis
195                  tps = SUM(tpps) / apols  
196              ! Calcul de la valeur moyenne aux p\^oles :
197                  DO ij = 1, iip1            forall (l = 1: llm)
198                     teta(ij, l) = tpn               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
199                     teta(ij + ip1jm, l) = tps                    / apoln
200                  ENDDO               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
201               ENDDO                    * teta(:iim, jjm + 1, l)) / apols
202              END forall
203               DO ij = 1, iim         END IF
204                  tppn(ij) = aire(ij) * ps(ij)  
205                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)         IF (MOD(itau + 1, iperiod) == 0) THEN
206               ENDDO            ! \'Ecriture du fichier histoire moyenne:
207               tpn = SUM(tppn) / apoln            CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
208               tps = SUM(tpps) / apols                 time = itau + 1)
209              call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
210               DO ij = 1, iip1                 q(:, :, :, 1))
211                  ps(ij) = tpn         ENDIF
212                  ps(ij + ip1jm) = tps  
213               ENDDO         IF (MOD(itau + 1, iecri * day_step) == 0) THEN
214            END IF            CALL geopot(teta, pk, pks, phis, phi)
215              CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps)
216            ! fin de l'intégration dynamique et physique pour le pas "itau"         END IF
217            ! préparation du pas d'intégration suivant      end do time_integration
218    
219            ! schema matsuno + leapfrog      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
220            IF (forward .OR. leapf) THEN           itau = itau_dyn + itaufin)
221               itau = itau + 1  
222               iday = day_ini + itau / day_step      ! Calcul des tendances dynamiques:
223               time = REAL(itau - (iday - day_ini) * day_step) / day_step &      CALL geopot(teta, pk, pks, phis, phi)
224                    + time_0      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
225               IF (time > 1.) THEN           dudyn, dv, dteta, dp, w, pbaru, pbarv, &
226                  time = time - 1.           conser = MOD(itaufin, iconser) == 0)
                 iday = iday + 1  
              ENDIF  
           ENDIF  
   
           IF (itau == itaufin + 1) exit outer_loop  
   
           IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) 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  
   
           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  
        end do  
     end do outer_loop  
227    
228    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
229    
Legend:
Removed from v.26  
changed lines
  Added in v.128
  ViewVC Help
Powered by ViewVC 1.1.21