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

Diff of /trunk/dyn3d/leapfrog.f90

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 178 by guez, Fri Mar 11 18:47:26 2016 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 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 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        ! Variables dynamiques:
59    
60      REAL tppn(iim), tpps(iim), tpn, tps      REAL pks(iim + 1, jjm + 1) ! exner au sol
61        REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
62        REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
63        REAL phi(iim + 1, jjm + 1, llm) ! geopotential
64        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 au pas - 1
71        REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
72        REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
73        REAL massem1(iim + 1, jjm + 1, llm)
74    
75        ! Tendances dynamiques
76        REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
77        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        ! Variables pour le fichier histoire
89      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
90      INTEGER iday ! jour julien      INTEGER itaufin
91      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.  
92    
93      logical:: dissip_conservative = .true.      ! Variables test conservation \'energie
94      logical forward, leapf, apphys, conser, apdiss      REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
95    
96        REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
97        logical leapf
98        real dt ! time step, in s
99    
100        REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
101        ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
102        ! for interface "l")
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      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)  
111    
112      ! Debut de l'integration temporelle:      ! On initialise la pression et la fonction d'Exner :
113      outer_loop:do      forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
114         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &      CALL exner_hyb(ps, p3d, pks, pk)
115              call guide(itau, ucov, vcov, teta, q, masse, ps)      pkf = pk
116         vcovm1 = vcov      CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
117         ucovm1 = ucov  
118         tetam1 = teta      time_integration: do itau = 0, itaufin - 1
119         massem1 = masse         leapf = mod(itau, iperiod) /= 0
120         psm1 = ps         if (leapf) then
121         forward = .TRUE.            dt = 2 * dtvr
122         leapf = .FALSE.         else
123         dt = dtvr            ! Matsuno
124         finvmaold = masse            dt = dtvr
125         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)            if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
126              vcovm1 = vcov
127         do            ucovm1 = ucov
128            ! gestion des appels de la physique et des dissipations:            tetam1 = teta
129            apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0            massem1 = masse
130            conser = MOD(itau, iconser) == 0            psm1 = ps
131            apdiss = MOD(itau + 1, idissip) == 0         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                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &              dudyn, dv, dteta, dp, w, pbaru, pbarv, &
137                 time + iday - day_ini)              conser = MOD(itau, iconser) == 0)
138    
139            IF (forward .OR. leapf) THEN         CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
140               ! calcul des tendances advection des traceurs (dont l'humidite)  
141               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)         ! Stokage du flux de masse pour traceurs offline:
142               IF (offline) THEN         IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
143                  ! Stokage du flux de masse pour traceurs off-line              dtvr, itau)
144                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &  
145                       itau)         ! Int\'egrations dynamique et traceurs:
146               ENDIF         CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
147            ENDIF              dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dt, 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, dtvr, &
164                 finvmaold, leapf)                 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, p3d, 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(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
224    
225               ! Calcul de la valeur moyenne, unique de h aux poles .....      ! Calcul des tendances dynamiques:
226               DO l = 1, llm      CALL geopot(teta, pk, pks, phis, phi)
227                  DO ij = 1, iim      CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
228                     tppn(ij) = aire(ij) * teta(ij, l)           dudyn, dv, dteta, dp, w, pbaru, pbarv, &
229                     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  
230    
231    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
232    
Legend:
Removed from v.25  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21