/[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

revision 23 by guez, Mon Dec 14 15:25:16 2009 UTC revision 27 by guez, Thu Mar 25 14:29:07 2010 UTC
# Line 1  Line 1 
1  module leapfrog_m  module leapfrog_m
2    
   ! This module is clean: no C preprocessor directive, no include line.  
   
3    IMPLICIT NONE    IMPLICIT NONE
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, time_0)
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
     ! Version du 10/01/98, avec coordonnees verticales hybrides, avec  
     ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)  
11    
12      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      USE calfis_m, ONLY: calfis
13      ! Objet:      USE com_io_dyn, ONLY: histaveid
14      ! GCM LMD nouvelle grille      USE comconst, ONLY: daysec, dtphys, dtvr
15        USE comgeom, ONLY: aire, apoln, apols
16      ! ... Dans inigeom, nouveaux calculs pour les elongations cu, cv      USE comvert, ONLY: ap, bp
17      ! et possibilite d'appeler une fonction f(y) a derivee tangente      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
18      ! hyperbolique a la place de la fonction a derivee sinusoidale.           periodav
19        USE dimens_m, ONLY: iim, llm, nqmx
20      ! ... Possibilité de choisir le schéma pour l'advection de      USE dynetat0_m, ONLY: day_ini
21      ! q, en modifiant iadv dans "traceur.def" (10/02) .      use dynredem1_m, only: dynredem1
22        USE exner_hyb_m, ONLY: exner_hyb
23      ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99)      use filtreg_m, only: filtreg
24      ! Pour Van-Leer iadv=10      USE guide_m, ONLY: guide
25        use inidissip_m, only: idissip
26      use dimens_m, only: iim, llm, nqmx      USE logic, ONLY: iflag_phys, ok_guide
27      use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1      USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1
28      use comconst, only: dtvr, daysec, dtphys      USE pression_m, ONLY: pression
29      use comvert, only: ap, bp      USE pressure_var, ONLY: p3d
30      use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &      USE temps, ONLY: dt, itau_dyn
          offline, periodav  
     use logic, only: ok_guide, iflag_phys  
     use comgeom  
     use serre  
     use temps, only: itaufin, day_ini, dt  
     use iniprint, only: prt_level  
     use com_io_dyn  
     use ener  
     use calfis_m, only: calfis  
     use exner_hyb_m, only: exner_hyb  
     use guide_m, only: guide  
     use pression_m, only: pression  
     use pressure_var, only: p3d  
31    
32      ! Variables dynamiques:      ! Variables dynamiques:
33      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
34      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL teta(ip1jmp1, llm) ! temperature potentielle
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
35      REAL ps(ip1jmp1) ! pression au sol, en Pa      REAL ps(ip1jmp1) ! pression au sol, en Pa
36    
37      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
38      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
39        REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
40      REAL time_0      REAL, intent(in):: time_0
41    
42      ! Variables local to the procedure:      ! Variables local to the procedure:
43    
# Line 91  contains Line 74  contains
74      REAL tppn(iim), tpps(iim), tpn, tps      REAL tppn(iim), tpps(iim), tpn, tps
75    
76      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
77      integer itaufinp1      INTEGER itaufin
78      INTEGER iday ! jour julien      INTEGER iday ! jour julien
79      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
   
     REAL SSUM  
80      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
81        LOGICAL:: lafin=.false.
     LOGICAL :: lafin=.false.  
82      INTEGER ij, l      INTEGER ij, l
83    
84      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
85    
86      !+jld variables test conservation energie      ! Variables test conservation energie
87      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
88      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
89      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
90      ! cree par la dissipation      ! cree par la dissipation
91      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(ip1jmp1, llm)
92      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
93      CHARACTER*15 ztit      logical forward, leapf
     INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.  
   
     logical:: dissip_conservative = .true.  
     logical forward, leapf, apphys, conser, apdiss  
94    
95      !---------------------------------------------------      !---------------------------------------------------
96    
97      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
98    
99      itaufin = nday * day_step      itaufin = nday * day_step
     itaufinp1 = itaufin + 1  
   
100      itau = 0      itau = 0
101      iday = day_ini      iday = day_ini
102      time = time_0      time = time_0
103      IF (time > 1.) THEN      dq = 0.
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
   
104      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
     dq=0.  
105      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
106      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
107    
108      ! Debut de l'integration temporelle:      ! Début de l'integration temporelle :
109      outer_loop:do      outer_loop:do
110         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
111            call guide(itau, ucov, vcov, teta, q, masse, ps)              call guide(itau, ucov, vcov, teta, q, masse, ps)
112         else         vcovm1 = vcov
113            IF (prt_level > 9) print *, &         ucovm1 = ucov
114                 'Attention : on ne guide pas les 6 dernieres heures.'         tetam1 = teta
115         endif         massem1 = masse
116           psm1 = ps
        CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)  
        CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)  
        CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)  
        CALL SCOPY(ijp1llm, masse, 1, massem1, 1)  
        CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)  
   
117         forward = .TRUE.         forward = .TRUE.
118         leapf = .FALSE.         leapf = .FALSE.
119         dt = dtvr         dt = dtvr
120           finvmaold = masse
        CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)  
121         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)
122    
123         do         do
124            ! gestion des appels de la physique et des dissipations:            ! Calcul des tendances dynamiques:
   
           apphys = .FALSE.  
           conser = .FALSE.  
           apdiss = .FALSE.  
   
           IF (MOD(itau, iconser) == 0) conser = .TRUE.  
           IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.  
           IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.  
   
           ! calcul des tendances dynamiques:  
   
125            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
   
126            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
127                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &                 MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
128                 time + iday - day_ini)                 time + iday - day_ini)
129    
           ! calcul des tendances advection des traceurs (dont l'humidite)  
   
130            IF (forward .OR. leapf) THEN            IF (forward .OR. leapf) THEN
131                 ! Calcul des tendances advection des traceurs (dont l'humidité)
132               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
133               IF (offline) THEN               IF (offline) THEN
134                  !maf stokage du flux de masse pour traceurs OFF-LINE                  ! Stokage du flux de masse pour traceurs off-line
135                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
136                       itau)                       itau)
137               ENDIF               ENDIF
# Line 193  contains Line 142  contains
142                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
143                 finvmaold, leapf)                 finvmaold, leapf)
144    
145            ! calcul des tendances physiques:            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
146                 ! calcul des tendances physiques:
           IF (apphys) THEN  
147               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
148    
149               CALL pression(ip1jmp1, ap, bp, ps, p3d)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
# Line 204  contains Line 152  contains
152               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
153               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
154    
              ! Interface avec les routines de phylmd (phymars ...)  
   
              ! Diagnostique de conservation de l'énergie : initialisation  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit='bil dyn'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
   
155               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
156                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
157                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
# Line 221  contains Line 160  contains
160               CALL addfi(nqmx, dtphys, &               CALL addfi(nqmx, dtphys, &
161                    ucov, vcov, teta, q, ps, &                    ucov, vcov, teta, q, ps, &
162                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
   
              ! Diagnostique de conservation de l'énergie : difference  
              IF (ip_ebil_dyn >= 1) THEN  
                 ztit = 'bil phys'  
                 CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &  
                      teta, q(:, :, 1), q(:, :, 2))  
              ENDIF  
163            ENDIF            ENDIF
164    
165            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
166            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
167    
168            ! dissipation horizontale et verticale des petites echelles:            IF (MOD(itau + 1, idissip) == 0) THEN
169                 ! dissipation horizontale et verticale des petites echelles:
170    
           IF (apdiss) THEN  
171               ! calcul de l'energie cinetique avant dissipation               ! calcul de l'energie cinetique avant dissipation
172               call covcont(llm, ucov, vcov, ucont, vcont)               call covcont(llm, ucov, vcov, ucont, vcont)
173               call enercin(vcov, ucov, vcont, ucont, ecin0)               call enercin(vcov, ucov, vcont, ucont, ecin0)
# Line 245  contains Line 177  contains
177               ucov=ucov + dudis               ucov=ucov + dudis
178               vcov=vcov + dvdis               vcov=vcov + dvdis
179    
180               if (dissip_conservative) then               ! On rajoute la tendance due a la transform. Ec -> E
181                  ! On rajoute la tendance due a la transform. Ec -> E               ! therm. cree lors de la dissipation
182                  ! therm. cree lors de la dissipation               call covcont(llm, ucov, vcov, ucont, vcont)
183                  call covcont(llm, ucov, vcov, ucont, vcont)               call enercin(vcov, ucov, vcont, ucont, ecin)
184                  call enercin(vcov, ucov, vcont, ucont, ecin)               dtetaecdt= (ecin0 - ecin) / pk
185                  dtetaecdt= (ecin0 - ecin) / pk               dtetadis=dtetadis + dtetaecdt
                 dtetadis=dtetadis + dtetaecdt  
              endif  
186               teta=teta + dtetadis               teta=teta + dtetadis
187    
188               ! Calcul de la valeur moyenne, unique de h aux poles .....               ! Calcul de la valeur moyenne, unique de h aux poles .....
   
189               DO l = 1, llm               DO l = 1, llm
190                  DO ij = 1, iim                  DO ij = 1, iim
191                     tppn(ij) = aire(ij) * teta(ij, l)                     tppn(ij) = aire(ij) * teta(ij, l)
192                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
193                  ENDDO                  ENDDO
194                  tpn = SSUM(iim, tppn, 1) / apoln                  tpn = SUM(tppn) / apoln
195                  tps = SSUM(iim, tpps, 1) / apols                  tps = SUM(tpps) / apols
196    
197                  DO ij = 1, iip1                  DO ij = 1, iip1
198                     teta(ij, l) = tpn                     teta(ij, l) = tpn
# Line 275  contains Line 204  contains
204                  tppn(ij) = aire(ij) * ps(ij)                  tppn(ij) = aire(ij) * ps(ij)
205                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
206               ENDDO               ENDDO
207               tpn = SSUM(iim, tppn, 1) / apoln               tpn = SUM(tppn) / apoln
208               tps = SSUM(iim, tpps, 1) / apols               tps = SUM(tpps) / apols
209    
210               DO ij = 1, iip1               DO ij = 1, iip1
211                  ps(ij) = tpn                  ps(ij) = tpn
212                  ps(ij + ip1jm) = tps                  ps(ij + ip1jm) = tps
213               ENDDO               ENDDO
   
214            END IF            END IF
215    
216            ! fin de l'intégration dynamique et physique pour le pas "itau"            ! fin de l'intégration dynamique et physique pour le pas "itau"
# Line 300  contains Line 228  contains
228               ENDIF               ENDIF
229            ENDIF            ENDIF
230    
231            IF (itau == itaufinp1) exit outer_loop            IF (itau == itaufin + 1) exit outer_loop
   
           ! ecriture du fichier histoire moyenne:  
232    
           ! Comment out the following calls when you do not want the output  
           ! files "dyn_hist_ave.nc" and "dynzon.nc"  
233            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
234                 ! ecriture du fichier histoire moyenne:
235               CALL writedynav(histaveid, nqmx, itau, vcov, &               CALL writedynav(histaveid, nqmx, itau, vcov, &
236                    ucov, teta, pk, phi, q, masse, ps, phis)                    ucov, teta, pk, phi, q, masse, ps, phis)
237               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
# Line 314  contains Line 239  contains
239            ENDIF            ENDIF
240    
241            IF (itau == itaufin) THEN            IF (itau == itaufin) THEN
242               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
243               CLOSE(99)                    itau=itau_dyn+itaufin)
244            ENDIF            ENDIF
245    
246            ! gestion de l'integration temporelle:            ! gestion de l'integration temporelle:
   
247            IF (MOD(itau, iperiod) == 0) exit            IF (MOD(itau, iperiod) == 0) exit
248            IF (MOD(itau - 1, iperiod) == 0) THEN            IF (MOD(itau - 1, iperiod) == 0) THEN
249               IF (forward) THEN               IF (forward) THEN
Legend:
Removed from v.23  
changed lines
  Added in v.27
  ViewVC Help
Powered by ViewVC 1.1.21