/[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 212 by guez, Thu Jan 12 12:31:31 2017 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: 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 enercin_m, only: enercin
31        USE exner_hyb_m, ONLY: exner_hyb
32        use filtreg_scal_m, only: filtreg_scal
33        use fluxstokenc_m, only: fluxstokenc
34        use geopot_m, only: geopot
35        USE guide_m, ONLY: guide
36        use inidissip_m, only: idissip
37        use integrd_m, only: integrd
38        use nr_util, only: assert
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        REAL p3d(iim + 1, jjm + 1, llm + 1) ! pressure at layer interfaces, in Pa
102        ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
103        ! for interface "l")
104    
105      !---------------------------------------------------      !---------------------------------------------------
106    
107      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
108        call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
109    
110      itaufin = nday * day_step      itaufin = nday * day_step
111      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)  
112    
113      ! Debut de l'integration temporelle:      ! On initialise la pression et la fonction d'Exner :
114      outer_loop:do      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
115         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &      CALL exner_hyb(ps, p3d, pks, pk)
116              call guide(itau, ucov, vcov, teta, q, masse, ps)      pkf = pk
117         vcovm1 = vcov      CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
118         ucovm1 = ucov  
119         tetam1 = teta      time_integration: do itau = 0, itaufin - 1
120         massem1 = masse         leapf = mod(itau, iperiod) /= 0
121         psm1 = ps         if (leapf) then
122         forward = .TRUE.            dt = 2 * dtvr
123         leapf = .FALSE.         else
124         dt = dtvr            ! Matsuno
125         finvmaold = masse            dt = dtvr
126         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)            if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
127              vcovm1 = vcov
128         do            ucovm1 = ucov
129            ! gestion des appels de la physique et des dissipations:            tetam1 = teta
130            apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0            massem1 = masse
131            conser = MOD(itau, iconser) == 0            psm1 = ps
132            apdiss = MOD(itau + 1, idissip) == 0         end if
133    
134            ! calcul des tendances dynamiques:         ! Calcul des tendances dynamiques:
135            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)         CALL geopot(teta, pk, pks, phis, phi)
136            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &         CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
137                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &              dudyn, dv, dteta, dp, w, pbaru, pbarv, &
138                 time + iday - day_ini)              conser = MOD(itau, iconser) == 0)
139    
140            IF (forward .OR. leapf) THEN         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
141               ! calcul des tendances advection des traceurs (dont l'humidite)  
142               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         ! Stokage du flux de masse pour traceurs offline:
143               IF (offline) THEN         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
144                  ! Stokage du flux de masse pour traceurs off-line              dtvr, itau)
145                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
146                       itau)         ! Int\'egrations dynamique et traceurs:
147               ENDIF         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
148            ENDIF              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dt, leapf)
149    
150           forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
151           CALL exner_hyb(ps, p3d, pks, pk)
152           pkf = pk
153           CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
154    
155           if (.not. leapf) then
156              ! Matsuno backward
157              ! Calcul des tendances dynamiques:
158              CALL geopot(teta, pk, pks, phis, phi)
159              CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
160                   phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
161    
162            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
163            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
164                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dtvr, &
165                 finvmaold, leapf)                 leapf=.false.)
166    
167            IF (apphys) THEN            forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
168               ! calcul des tendances physiques:            CALL exner_hyb(ps, p3d, pks, pk)
169               IF (itau + 1 == itaufin) lafin = .TRUE.            pkf = pk
170              CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
171               CALL pression(ip1jmp1, ap, bp, ps, p3d)         end if
172               CALL exner_hyb(ps, p3d, pks, pk, pkf)  
173           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys) THEN
174               rdaym_ini = itau * dtvr / daysec            CALL calfis(ucov, vcov, teta, q, p3d, pk, phis, phi, w, dufi, dvfi, &
175               rdayvrai = rdaym_ini + day_ini                 dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
176                   time = REAL(mod(itau, day_step)) / day_step, &
177               ! Diagnostique de conservation de l'énergie : initialisation                 lafin = itau + 1 == itaufin)
178               IF (ip_ebil_dyn >= 1) THEN  
179                  ztit='bil dyn'            CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
180                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &         ENDIF
181                       teta, q(:, :, 1), q(:, :, 2))  
182               ENDIF         IF (MOD(itau + 1, idissip) == 0) THEN
183              ! Dissipation horizontale et verticale des petites \'echelles
184               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
185                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &            ! calcul de l'\'energie cin\'etique avant dissipation
186                    dufi, dvfi, dtetafi, dqfi, dpfi)            call covcont(llm, ucov, vcov, ucont, vcont)
187              call enercin(vcov, ucov, vcont, ucont, ecin0)
188               ! ajout des tendances physiques:  
189               CALL addfi(nqmx, dtphys, &            ! dissipation
190                    ucov, vcov, teta, q, ps, &            CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
191                    dufi, dvfi, dtetafi, dqfi, dpfi)            ucov = ucov + dudis
192              vcov = vcov + dvdis
193               ! Diagnostique de conservation de l'énergie : difference  
194               IF (ip_ebil_dyn >= 1) THEN            ! On ajoute la tendance due \`a la transformation \'energie
195                  ztit = 'bil phys'            ! cin\'etique en \'energie thermique par la dissipation
196                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &            call covcont(llm, ucov, vcov, ucont, vcont)
197                       teta, q(:, :, 1), q(:, :, 2))            call enercin(vcov, ucov, vcont, ucont, ecin)
198               ENDIF            dtetadis = dtetadis + (ecin0 - ecin) / pk
199            ENDIF            teta = teta + dtetadis
200    
201            CALL pression(ip1jmp1, ap, bp, ps, p3d)            ! Calcul de la valeur moyenne aux p\^oles :
202            CALL exner_hyb(ps, p3d, pks, pk, pkf)            forall (l = 1: llm)
203                 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
204            IF (apdiss) THEN                    / apoln
205               ! dissipation horizontale et verticale des petites echelles:               teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm + 1) &
206                      * teta(:iim, jjm + 1, l)) / apols
207               ! calcul de l'energie cinetique avant dissipation            END forall
208               call covcont(llm, ucov, vcov, ucont, vcont)         END IF
209               call enercin(vcov, ucov, vcont, ucont, ecin0)  
210           IF (MOD(itau + 1, iperiod) == 0) THEN
211               ! dissipation            ! \'Ecriture du fichier histoire moyenne:
212               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)            CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
213               ucov=ucov + dudis                 time = itau + 1)
214               vcov=vcov + dvdis            call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
215                   q(:, :, :, 1))
216               if (dissip_conservative) then         ENDIF
217                  ! On rajoute la tendance due a la transform. Ec -> E  
218                  ! therm. cree lors de la dissipation         IF (MOD(itau + 1, iecri * day_step) == 0) THEN
219                  call covcont(llm, ucov, vcov, ucont, vcont)            CALL geopot(teta, pk, pks, phis, phi)
220                  call enercin(vcov, ucov, vcont, ucont, ecin)            CALL writehist(itau, vcov, ucov, teta, phi, masse, ps)
221                  dtetaecdt= (ecin0 - ecin) / pk         END IF
222                  dtetadis=dtetadis + dtetaecdt      end do time_integration
223               endif  
224               teta=teta + dtetadis      CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
225    
226               ! Calcul de la valeur moyenne, unique de h aux poles .....      ! Calcul des tendances dynamiques:
227               DO l = 1, llm      CALL geopot(teta, pk, pks, phis, phi)
228                  DO ij = 1, iim      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
229                     tppn(ij) = aire(ij) * teta(ij, l)           dudyn, dv, dteta, dp, w, pbaru, pbarv, &
230                     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  
231    
232    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
233    
Legend:
Removed from v.25  
changed lines
  Added in v.212
  ViewVC Help
Powered by ViewVC 1.1.21