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