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

Diff of /trunk/libf/dyn3d/leapfrog.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 24 by guez, Wed Mar 3 13:23:49 2010 UTC revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC
# Line 7  contains Line 7  contains
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
   
     ! Version du 10/01/98, avec coordonnées verticales hybrides, avec  
     ! nouveaux opérateurs dissipation "*" (gradiv2, divgrad2, nxgraro2)  
   
10      ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin      ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin
     ! Objet: nouvelle grille  
11    
12      ! Dans "inigeom", nouveaux calculs pour les élongations cu, cv      USE dimens_m, ONLY : iim, llm, nqmx
13      ! et possibilité d'appeler une fonction f(y) à dérivée tangente      USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
14      ! hyperbolique à la place de la fonction à dérivée sinusoïdale.      USE comconst, ONLY : daysec, dtphys, dtvr
15        USE comvert, ONLY : ap, bp
16      ! Possibilité de choisir le schéma pour l'advection de      USE conf_gcm_m, ONLY : day_step, iconser, idissip, iperiod, iphysiq, &
17      ! q, en modifiant iadv dans "traceur.def".           nday, offline, periodav
18        USE logic, ONLY : iflag_phys, ok_guide
19      ! Pour Van-Leer + vapeur d'eau saturée, iadv(1)=4.      USE comgeom, ONLY : aire, apoln, apols
20      ! Pour Van-Leer iadv=10      USE temps, ONLY : dt, itaufin
21        USE dynetat0_m, ONLY : day_ini
22      use dimens_m, only: iim, llm, nqmx      USE iniprint, ONLY : prt_level
23      use paramet_m, only: ip1jmp1, ip1jm, ijp1llm, jjp1, iip1      USE com_io_dyn, ONLY : histaveid
24      use comconst, only: dtvr, daysec, dtphys      USE calfis_m, ONLY : calfis
25      use comvert, only: ap, bp      USE exner_hyb_m, ONLY : exner_hyb
26      use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &      USE guide_m, ONLY : guide
27           offline, periodav      USE pression_m, ONLY : pression
28      use logic, only: ok_guide, iflag_phys      USE pressure_var, ONLY : p3d
     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  
29    
30      ! Variables dynamiques:      ! Variables dynamiques:
31      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
32      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL teta(ip1jmp1, llm) ! temperature potentielle
     REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields  
33      REAL ps(ip1jmp1) ! pression au sol, en Pa      REAL ps(ip1jmp1) ! pression au sol, en Pa
34    
35      REAL masse(ip1jmp1, llm) ! masse d'air      REAL masse(ip1jmp1, llm) ! masse d'air
36      REAL phis(ip1jmp1) ! geopotentiel au sol      REAL phis(ip1jmp1) ! geopotentiel au sol
37        REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
38      REAL, intent(in):: time_0      REAL, intent(in):: time_0
39    
40      ! Variables local to the procedure:      ! Variables local to the procedure:
# Line 90  contains Line 74  contains
74      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
75      INTEGER iday ! jour julien      INTEGER iday ! jour julien
76      REAL time ! time of day, as a fraction of day length      REAL time ! time of day, as a fraction of day length
   
     REAL SSUM  
77      real finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
   
78      LOGICAL :: lafin=.false.      LOGICAL :: lafin=.false.
79      INTEGER ij, l      INTEGER ij, l
80    
# Line 120  contains Line 101  contains
101      itau = 0      itau = 0
102      iday = day_ini      iday = day_ini
103      time = time_0      time = time_0
     IF (time > 1.) THEN  
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
   
104      dq = 0.      dq = 0.
105      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
106      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
# Line 132  contains Line 108  contains
108    
109      ! Debut de l'integration temporelle:      ! Debut de l'integration temporelle:
110      outer_loop:do      outer_loop:do
111         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600.) then         if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
112            call guide(itau, ucov, vcov, teta, q, masse, ps)              call guide(itau, ucov, vcov, teta, q, masse, ps)
113         else         vcovm1 = vcov
114            IF (prt_level > 9) print *, &         ucovm1 = ucov
115                 'Attention : on ne guide pas les 6 dernières heures.'         tetam1 = teta
116         endif         massem1 = masse
117           psm1 = ps
        CALL SCOPY(ip1jm * llm, 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)  
   
118         forward = .TRUE.         forward = .TRUE.
119         leapf = .FALSE.         leapf = .FALSE.
120         dt = dtvr         dt = dtvr
121           finvmaold = masse
        CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)  
122         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)         CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)
123    
124         do         do
125            ! gestion des appels de la physique et des dissipations:            ! gestion des appels de la physique et des dissipations:
126              apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0
127            apphys = .FALSE.            conser = MOD(itau, iconser) == 0
128            conser = .FALSE.            apdiss = MOD(itau + 1, idissip) == 0
           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.  
129    
130            ! calcul des tendances dynamiques:            ! calcul des tendances dynamiques:
   
131            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, 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, &                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &
134                 time + iday - day_ini)                 time + iday - day_ini)
135    
           ! calcul des tendances advection des traceurs (dont l'humidite)  
   
136            IF (forward .OR. leapf) THEN            IF (forward .OR. leapf) THEN
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, dq, teta, pk)
139               IF (offline) THEN               IF (offline) THEN
140                  !maf stokage du flux de masse pour traceurs OFF-LINE                  ! Stokage du flux de masse pour traceurs off-line
141                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
142                       itau)                       itau)
143               ENDIF               ENDIF
# Line 187  contains Line 148  contains
148                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
149                 finvmaold, leapf)                 finvmaold, leapf)
150    
           ! calcul des tendances physiques:  
   
151            IF (apphys) THEN            IF (apphys) THEN
152                 ! calcul des tendances physiques:
153               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
154    
155               CALL pression(ip1jmp1, ap, bp, ps, p3d)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
# Line 198  contains Line 158  contains
158               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
159               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
160    
              ! Interface avec les routines de phylmd (phymars ...)  
   
161               ! Diagnostique de conservation de l'énergie : initialisation               ! Diagnostique de conservation de l'énergie : initialisation
162               IF (ip_ebil_dyn >= 1) THEN               IF (ip_ebil_dyn >= 1) THEN
163                  ztit='bil dyn'                  ztit='bil dyn'
# Line 227  contains Line 185  contains
185            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
186            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
187    
           ! dissipation horizontale et verticale des petites echelles:  
   
188            IF (apdiss) THEN            IF (apdiss) THEN
189                 ! dissipation horizontale et verticale des petites echelles:
190    
191               ! calcul de l'energie cinetique avant dissipation               ! calcul de l'energie cinetique avant dissipation
192               call covcont(llm, ucov, vcov, ucont, vcont)               call covcont(llm, ucov, vcov, ucont, vcont)
193               call enercin(vcov, ucov, vcont, ucont, ecin0)               call enercin(vcov, ucov, vcont, ucont, ecin0)
# Line 250  contains Line 208  contains
208               teta=teta + dtetadis               teta=teta + dtetadis
209    
210               ! Calcul de la valeur moyenne, unique de h aux poles .....               ! Calcul de la valeur moyenne, unique de h aux poles .....
   
211               DO l = 1, llm               DO l = 1, llm
212                  DO ij = 1, iim                  DO ij = 1, iim
213                     tppn(ij) = aire(ij) * teta(ij, l)                     tppn(ij) = aire(ij) * teta(ij, l)
214                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
215                  ENDDO                  ENDDO
216                  tpn = SSUM(iim, tppn, 1) / apoln                  tpn = SUM(tppn) / apoln
217                  tps = SSUM(iim, tpps, 1) / apols                  tps = SUM(tpps) / apols
218    
219                  DO ij = 1, iip1                  DO ij = 1, iip1
220                     teta(ij, l) = tpn                     teta(ij, l) = tpn
# Line 269  contains Line 226  contains
226                  tppn(ij) = aire(ij) * ps(ij)                  tppn(ij) = aire(ij) * ps(ij)
227                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
228               ENDDO               ENDDO
229               tpn = SSUM(iim, tppn, 1) / apoln               tpn = SUM(tppn) / apoln
230               tps = SSUM(iim, tpps, 1) / apols               tps = SUM(tpps) / apols
231    
232               DO ij = 1, iip1               DO ij = 1, iip1
233                  ps(ij) = tpn                  ps(ij) = tpn
234                  ps(ij + ip1jm) = tps                  ps(ij + ip1jm) = tps
235               ENDDO               ENDDO
   
236            END IF            END IF
237    
238            ! 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 296  contains Line 252  contains
252    
253            IF (itau == itaufin + 1) exit outer_loop            IF (itau == itaufin + 1) exit outer_loop
254    
           ! ecriture du fichier histoire moyenne:  
   
           ! Comment out the following calls when you do not want the output  
           ! files "dyn_hist_ave.nc" and "dynzon.nc"  
255            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
256                 ! ecriture du fichier histoire moyenne:
257               CALL writedynav(histaveid, nqmx, itau, vcov, &               CALL writedynav(histaveid, nqmx, itau, vcov, &
258                    ucov, teta, pk, phi, q, masse, ps, phis)                    ucov, teta, pk, phi, q, masse, ps, phis)
259               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
# Line 309  contains Line 262  contains
262    
263            IF (itau == itaufin) THEN            IF (itau == itaufin) THEN
264               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)
              CLOSE(99)  
265            ENDIF            ENDIF
266    
267            ! gestion de l'integration temporelle:            ! gestion de l'integration temporelle:
   
268            IF (MOD(itau, iperiod) == 0) exit            IF (MOD(itau, iperiod) == 0) exit
269            IF (MOD(itau - 1, iperiod) == 0) THEN            IF (MOD(itau - 1, iperiod) == 0) THEN
270               IF (forward) THEN               IF (forward) THEN
Legend:
Removed from v.24  
changed lines
  Added in v.25
  ViewVC Help
Powered by ViewVC 1.1.21