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

Diff of /trunk/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 161 by guez, Fri Jul 24 14:27:59 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 covcont_m, only: covcont
22      USE iniprint, ONLY : prt_level      USE disvert_m, ONLY: ap, bp
23      USE com_io_dyn, ONLY : histaveid      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
24      USE calfis_m, ONLY : calfis           iflag_phys, iecri
25      USE exner_hyb_m, ONLY : exner_hyb      USE conf_guide_m, ONLY: ok_guide
26      USE guide_m, ONLY : guide      USE dimens_m, ONLY: iim, jjm, llm, nqmx
27      USE pression_m, ONLY : pression      use dissip_m, only: dissip
28      USE pressure_var, ONLY : p3d      USE dynetat0_m, ONLY: day_ini
29        use dynredem1_m, only: dynredem1
30        USE exner_hyb_m, ONLY: exner_hyb
31        use filtreg_scal_m, only: filtreg_scal
32        use fluxstokenc_m, only: fluxstokenc
33        use geopot_m, only: geopot
34        USE guide_m, ONLY: guide
35        use inidissip_m, only: idissip
36        use integrd_m, only: integrd
37        use nr_util, only: assert
38        USE pressure_var, ONLY: p3d
39        USE temps, ONLY: itau_dyn
40        use writedynav_m, only: writedynav
41        use writehist_m, only: writehist
42    
43      ! Variables dynamiques:      ! Variables dynamiques:
44      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
45      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  
46    
47      ! Variables local to the procedure:      REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
48        ! potential temperature
49    
50      ! Variables dynamiques:      REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
51        REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
52        REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
53    
54      REAL pks(ip1jmp1) ! exner au sol      REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
55      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)  
56    
57      ! variables pour le fichier histoire      ! Local:
58    
59        ! Variables dynamiques:
60    
61      REAL tppn(iim), tpps(iim), tpn, tps      REAL pks(iim + 1, jjm + 1) ! exner au sol
62        REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
63        REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
64        REAL phi(iim + 1, jjm + 1, llm) ! geopotential
65        REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
66    
67        ! Variables dynamiques intermediaire pour le transport
68        ! Flux de masse :
69        REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
70    
71        ! Variables dynamiques au pas - 1
72        REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
73        REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
74        REAL massem1(iim + 1, jjm + 1, llm)
75    
76        ! Tendances dynamiques
77        REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
78        REAL dteta(iim + 1, jjm + 1, llm)
79        real dp((iim + 1) * (jjm + 1))
80    
81        ! Tendances de la dissipation :
82        REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
83        REAL dtetadis(iim + 1, jjm + 1, llm)
84    
85        ! Tendances physiques
86        REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
87        REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
88    
89        ! Variables pour le fichier histoire
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 iday ! jour julien      INTEGER itaufin
92      REAL time ! time of day, as a fraction of day length      INTEGER l
     real finvmaold(ip1jmp1, llm)  
     LOGICAL :: lafin=.false.  
     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".
     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         end if
129    
130            ! calcul des tendances dynamiques:         ! Calcul des tendances dynamiques:
131            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL geopot(teta, pk, pks, phis, phi)
132            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
133                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &              dudyn, dv, dteta, dp, w, pbaru, pbarv, &
134                 time + iday - day_ini)              conser = MOD(itau, iconser) == 0)
135    
136            IF (forward .OR. leapf) THEN         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
137               ! calcul des tendances advection des traceurs (dont l'humidite)  
138               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         ! Stokage du flux de masse pour traceurs offline:
139               IF (offline) THEN         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
140                  ! Stokage du flux de masse pour traceurs off-line              dtvr, itau)
141                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
142                       itau)         ! Int\'egrations dynamique et traceurs:
143               ENDIF         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
144            ENDIF              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dt, leapf)
145    
146           forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
147           CALL exner_hyb(ps, p3d, pks, pk)
148           pkf = pk
149           CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
150    
151           if (.not. leapf) then
152              ! Matsuno backward
153              ! Calcul des tendances dynamiques:
154              CALL geopot(teta, pk, pks, phis, phi)
155              CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
156                   phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
157    
158            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
159            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
160                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dtvr, &
161                 finvmaold, leapf)                 leapf=.false.)
162    
163            IF (apphys) THEN            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
164               ! calcul des tendances physiques:            CALL exner_hyb(ps, p3d, pks, pk)
165               IF (itau + 1 == itaufin) lafin = .TRUE.            pkf = pk
166              CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
167               CALL pression(ip1jmp1, ap, bp, ps, p3d)         end if
168               CALL exner_hyb(ps, p3d, pks, pk, pkf)  
169           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
170               rdaym_ini = itau * dtvr / daysec            CALL calfis(ucov, vcov, teta, q, pk, phis, phi, w, dufi, dvfi, &
171               rdayvrai = rdaym_ini + day_ini                 dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
172                   time = REAL(mod(itau, day_step)) / day_step, &
173               ! Diagnostique de conservation de l'énergie : initialisation                 lafin = itau + 1 == itaufin)
174               IF (ip_ebil_dyn >= 1) THEN  
175                  ztit='bil dyn'            CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
176                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &         ENDIF
177                       teta, q(:, :, 1), q(:, :, 2))  
178               ENDIF         IF (MOD(itau + 1, idissip) == 0) THEN
179              ! Dissipation horizontale et verticale des petites \'echelles
180               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
181                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &            ! calcul de l'\'energie cin\'etique avant dissipation
182                    dufi, dvfi, dtetafi, dqfi, dpfi)            call covcont(llm, ucov, vcov, ucont, vcont)
183              call enercin(vcov, ucov, vcont, ucont, ecin0)
184               ! ajout des tendances physiques:  
185               CALL addfi(nqmx, dtphys, &            ! dissipation
186                    ucov, vcov, teta, q, ps, &            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
187                    dufi, dvfi, dtetafi, dqfi, dpfi)            ucov = ucov + dudis
188              vcov = vcov + dvdis
189               ! Diagnostique de conservation de l'énergie : difference  
190               IF (ip_ebil_dyn >= 1) THEN            ! On ajoute la tendance due \`a la transformation \'energie
191                  ztit = 'bil phys'            ! cin\'etique en \'energie thermique par la dissipation
192                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &            call covcont(llm, ucov, vcov, ucont, vcont)
193                       teta, q(:, :, 1), q(:, :, 2))            call enercin(vcov, ucov, vcont, ucont, ecin)
194               ENDIF            dtetadis = dtetadis + (ecin0 - ecin) / pk
195            ENDIF            teta = teta + dtetadis
196    
197            CALL pression(ip1jmp1, ap, bp, ps, p3d)            ! Calcul de la valeur moyenne aux p\^oles :
198            CALL exner_hyb(ps, p3d, pks, pk, pkf)            forall (l = 1: llm)
199                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
200            IF (apdiss) THEN                    / apoln
201               ! dissipation horizontale et verticale des petites echelles:               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
202                      * teta(:iim, jjm + 1, l)) / apols
203               ! calcul de l'energie cinetique avant dissipation            END forall
204               call covcont(llm, ucov, vcov, ucont, vcont)         END IF
205               call enercin(vcov, ucov, vcont, ucont, ecin0)  
206           IF (MOD(itau + 1, iperiod) == 0) THEN
207               ! dissipation            ! \'Ecriture du fichier histoire moyenne:
208               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
209               ucov=ucov + dudis                 time = itau + 1)
210               vcov=vcov + dvdis            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
211                   q(:, :, :, 1))
212               if (dissip_conservative) then         ENDIF
213                  ! On rajoute la tendance due a la transform. Ec -> E  
214                  ! therm. cree lors de la dissipation         IF (MOD(itau + 1, iecri * day_step) == 0) THEN
215                  call covcont(llm, ucov, vcov, ucont, vcont)            CALL geopot(teta, pk, pks, phis, phi)
216                  call enercin(vcov, ucov, vcont, ucont, ecin)            CALL writehist(itau, vcov, ucov, teta, phi, masse, ps)
217                  dtetaecdt= (ecin0 - ecin) / pk         END IF
218                  dtetadis=dtetadis + dtetaecdt      end do time_integration
219               endif  
220               teta=teta + dtetadis      CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
221    
222               ! Calcul de la valeur moyenne, unique de h aux poles .....      ! Calcul des tendances dynamiques:
223               DO l = 1, llm      CALL geopot(teta, pk, pks, phis, phi)
224                  DO ij = 1, iim      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
225                     tppn(ij) = aire(ij) * teta(ij, l)           dudyn, dv, dteta, dp, w, pbaru, pbarv, &
226                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)           conser = MOD(itaufin, iconser) == 0)
                 ENDDO  
                 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  
227    
228    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
229    
Legend:
Removed from v.25  
changed lines
  Added in v.161
  ViewVC Help
Powered by ViewVC 1.1.21