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