/[lmdze]/trunk/dyn3d/leapfrog.f
ViewVC logotype

Diff of /trunk/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 30 by guez, Thu Apr 1 09:07:28 2010 UTC revision 47 by guez, Fri Jul 1 15:00:48 2011 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      ! 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 com_io_dyn, ONLY: histaveid
19      USE comconst, ONLY: daysec, dtphys, dtvr      USE comconst, ONLY: daysec, dtphys, dtvr
# Line 18  contains Line 22  contains
22      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
23           periodav           periodav
24      USE dimens_m, ONLY: iim, jjm, 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      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      use filtreg_m, only: filtreg
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
34      USE logic, ONLY: iflag_phys, ok_guide      USE logic, ONLY: iflag_phys, ok_guide
35      USE paramet_m, ONLY: ip1jmp1      USE paramet_m, ONLY: ip1jmp1
     USE pression_m, ONLY: pression  
36      USE pressure_var, ONLY: p3d      USE pressure_var, ONLY: p3d
37      USE temps, ONLY: itau_dyn      USE temps, ONLY: itau_dyn
38    
39      ! Variables dynamiques:      ! Variables dynamiques:
40      REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
41      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature      REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
     REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa  
42    
43        REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
44        ! potential temperature
45    
46        REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
47      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
48      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
49      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
50        REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
51        ! mass fractions of advected fields
52    
53      REAL, intent(in):: time_0      REAL, intent(in):: time_0
54    
55      ! Variables local to the procedure:      ! Variables local to the procedure:
# Line 47  contains Line 59  contains
59      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
60      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches      REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
61      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
62      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(iim + 1, jjm + 1, llm) ! geopotential
63      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
64    
65      ! variables dynamiques intermediaire pour le transport      ! variables dynamiques intermediaire pour le transport
# Line 59  contains Line 71  contains
71      REAL massem1(ip1jmp1, llm)      REAL massem1(ip1jmp1, llm)
72    
73      ! tendances dynamiques      ! tendances dynamiques
74      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)      REAL dv((iim + 1) * jjm, llm), dudyn(ip1jmp1, llm)
75      REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)      REAL dteta(iim + 1, jjm + 1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
76    
77      ! tendances de la dissipation      ! tendances de la dissipation
78      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)      REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
# Line 68  contains Line 80  contains
80    
81      ! tendances physiques      ! tendances physiques
82      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)      REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
83      REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)      REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
84    
85      ! variables pour le fichier histoire      ! variables pour le fichier histoire
86    
87      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
88      INTEGER itaufin      INTEGER itaufin
     INTEGER iday ! jour julien  
89      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
90      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
91      LOGICAL:: lafin=.false.      INTEGER l
     INTEGER i, j, l  
   
92      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
93    
94      ! Variables test conservation energie      ! Variables test conservation energie
95      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)
96      ! Tendance de la temp. potentiel d (theta) / d t due a la  
     ! tansformation d'energie cinetique en energie thermique  
     ! cree par la dissipation  
97      REAL dtetaecdt(iim + 1, jjm + 1, llm)      REAL dtetaecdt(iim + 1, jjm + 1, llm)
98        ! tendance de la température potentielle due à la tansformation
99        ! d'énergie cinétique en énergie thermique créée par la dissipation
100    
101      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)      REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
102      logical forward, leapf      logical leapf
103      REAL dt      real dt
104    
105      !---------------------------------------------------      !---------------------------------------------------
106    
# Line 99  contains Line 109  contains
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 "itaufin" is one too
111    
     itau = 0  
     iday = day_ini  
     time = time_0  
112      dq = 0.      dq = 0.
113    
114      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
115      CALL pression(ip1jmp1, ap, bp, ps, p3d)      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
116      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
117    
118      ! Début de l'integration temporelle :      time_integration: do itau = 0, itaufin - 1
119      outer_loop:do i = 1, itaufin / iperiod         leapf = mod(itau, iperiod) /= 0
120         ! {itau is a multiple of iperiod}         if (leapf) then
121              dt = 2 * dtvr
122         ! 1. Matsuno forward:         else
123              ! Matsuno
124         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &            dt = dtvr
125              call guide(itau, ucov, vcov, teta, q, masse, ps)            if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
126         vcovm1 = vcov                 call guide(itau, ucov, vcov, teta, q, masse, ps)
127         ucovm1 = ucov            vcovm1 = vcov
128         tetam1 = teta            ucovm1 = ucov
129         massem1 = masse            tetam1 = teta
130         psm1 = ps            massem1 = masse
131         finvmaold = masse            psm1 = ps
132         CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)            finvmaold = masse
133              CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
134           end if
135    
136         ! Calcul des tendances dynamiques:         ! Calcul des tendances dynamiques:
137         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
138         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
139              MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &              dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
140              time + iday - day_ini)              conser=MOD(itau, iconser)==0)
141    
142         ! Calcul des tendances advection des traceurs (dont l'humidité)         ! Calcul des tendances advection des traceurs (dont l'humidité)
143         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
144    
145         ! Stokage du flux de masse pour traceurs offline:         ! Stokage du flux de masse pour traceurs offline:
146         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
147              dtvr, itau)              dtvr, itau)
148    
149         ! integrations dynamique et traceurs:         ! Integrations dynamique et traceurs:
150         CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
151              dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., &              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
152              dtvr)              leapf)
153    
154         CALL pression(ip1jmp1, ap, bp, ps, p3d)         if (.not. leapf) then
155         CALL exner_hyb(ps, p3d, pks, pk, pkf)            ! Matsuno backward
156              forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
157         ! 2. Matsuno backward:            CALL exner_hyb(ps, p3d, pks, pk, pkf)
   
        itau = itau + 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  
   
        ! Calcul des tendances dynamiques:  
        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)  
        CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &  
             .false., du, dv, dteta, dp, w, pbaru, pbarv, time + iday - day_ini)  
   
        ! integrations dynamique et traceurs:  
        CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &  
             dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., &  
             dtvr)  
   
        CALL pression(ip1jmp1, ap, bp, ps, p3d)  
        CALL exner_hyb(ps, p3d, pks, pk, pkf)  
   
        ! 3. Leapfrog:  
158    
        do j = 1, iperiod - 1  
159            ! Calcul des tendances dynamiques:            ! Calcul des tendances dynamiques:
160            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
161            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
162                 .false., du, dv, dteta, dp, w, pbaru, pbarv, &                 phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
163                 time + iday - day_ini)                 conser=.false.)
   
           ! Calcul des tendances advection des traceurs (dont l'humidité)  
           CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)  
           ! Stokage du flux de masse pour traceurs off-line:  
           IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &  
                dtvr, itau)  
164    
165            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
166            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
167                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
168                 finvmaold, .true., 2 * dtvr)                 finvmaold, dtvr, leapf=.false.)
169           end if
           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) 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  
   
              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)  
           ENDIF  
170    
171            CALL pression(ip1jmp1, ap, bp, ps, p3d)         IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
172              ! calcul des tendances physiques:
173    
174              forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
175            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
176    
177            IF (MOD(itau + 1, idissip) == 0) THEN            rdaym_ini = itau * dtvr / daysec
178               ! dissipation horizontale et verticale des petites echelles:            rdayvrai = rdaym_ini + day_ini
179              time = REAL(mod(itau, day_step)) / day_step + time_0
180              IF (time > 1.) time = time - 1.
181    
182              CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
183                   phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
184                   lafin=itau+1==itaufin)
185    
186              ! ajout des tendances physiques:
187              CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
188                   dtetafi, dqfi, dpfi)
189           ENDIF
190    
191           forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
192           CALL exner_hyb(ps, p3d, pks, pk, pkf)
193    
194           IF (MOD(itau + 1, idissip) == 0) THEN
195              ! dissipation horizontale et verticale des petites echelles:
196    
197               ! calcul de l'energie cinetique avant dissipation            ! calcul de l'energie cinetique avant dissipation
198               call covcont(llm, ucov, vcov, ucont, vcont)            call covcont(llm, ucov, vcov, ucont, vcont)
199               call enercin(vcov, ucov, vcont, ucont, ecin0)            call enercin(vcov, ucov, vcont, ucont, ecin0)
200    
201               ! dissipation            ! dissipation
202               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
203               ucov=ucov + dudis            ucov=ucov + dudis
204               vcov=vcov + dvdis            vcov=vcov + dvdis
205    
206               ! On rajoute la tendance due à la transformation Ec -> E            ! On rajoute la tendance due à la transformation Ec -> E
207               ! thermique créée lors de la dissipation            ! thermique créée lors de la dissipation
208               call covcont(llm, ucov, vcov, ucont, vcont)            call covcont(llm, ucov, vcov, ucont, vcont)
209               call enercin(vcov, ucov, vcont, ucont, ecin)            call enercin(vcov, ucov, vcont, ucont, ecin)
210               dtetaecdt= (ecin0 - ecin) / pk            dtetaecdt= (ecin0 - ecin) / pk
211               dtetadis=dtetadis + dtetaecdt            dtetadis=dtetadis + dtetaecdt
212               teta=teta + dtetadis            teta=teta + dtetadis
213    
214               ! Calcul de la valeur moyenne unique de h aux pôles            ! Calcul de la valeur moyenne aux pôles :
215               forall (l = 1: llm)            forall (l = 1: llm)
216                  teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &               teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
217                       / apoln                    / apoln
218                  teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
219                       * teta(:iim, jjm + 1, l)) / apols                    * teta(:iim, jjm + 1, l)) / apols
220               END forall            END forall
221    
222               ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln            ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
223               ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &            ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
224                    / apols                 / apols
225            END IF         END IF
226    
227            itau = itau + 1         IF (MOD(itau + 1, iperiod) == 0) THEN
228            iday = day_ini + itau / day_step            ! Écriture du fichier histoire moyenne:
229            time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0            CALL writedynav(histaveid, nqmx, itau + 1, vcov, ucov, teta, pk, &
230            IF (time > 1.) THEN                 phi, q, masse, ps, phis)
231               time = time - 1.            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
232               iday = iday + 1                 q(:, :, :, 1), dt_app = dtvr * iperiod, &
233            ENDIF                 dt_cum = dtvr * day_step * periodav)
234           ENDIF
235            IF (MOD(itau, iperiod) == 0) THEN      end do time_integration
              ! 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  
        end do  
     end do outer_loop  
236    
     ! {itau == itaufin}  
237      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &      CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
238           itau=itau_dyn+itaufin)           itau=itau_dyn+itaufin)
239    
     vcovm1 = vcov  
     ucovm1 = ucov  
     tetam1 = teta  
     massem1 = masse  
     psm1 = ps  
     finvmaold = masse  
     CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)  
   
240      ! Calcul des tendances dynamiques:      ! Calcul des tendances dynamiques:
241      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
242      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
243           MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &           dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
244           time + iday - day_ini)           conser=MOD(itaufin, iconser)==0)
   
     ! Calcul des tendances advection des traceurs (dont l'humidité)  
     CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)  
     ! Stokage du flux de masse pour traceurs off-line:  
     IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
          itaufin)  
   
     ! integrations dynamique et traceurs:  
     CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, dq, &  
          dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold, .false., dtvr)  
       
     CALL pression(ip1jmp1, ap, bp, ps, p3d)  
     CALL exner_hyb(ps, p3d, pks, pk, pkf)  
   
245    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
246    
247  end module leapfrog_m  end module leapfrog_m
Legend:
Removed from v.30  
changed lines
  Added in v.47
  ViewVC Help
Powered by ViewVC 1.1.21