/[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 22 by guez, Fri Jul 31 15:18:47 2009 UTC revision 26 by guez, Tue Mar 9 15:27:15 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, nq, 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        ! Auteurs: P. Le Van, L. Fairhead, F. Hourdin
11    
12      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      USE calfis_m, ONLY: calfis
13      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      USE com_io_dyn, ONLY: histaveid
14        USE comconst, ONLY: daysec, dtphys, dtvr
15      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      USE comgeom, ONLY: aire, apoln, apols
16      ! Objet:      USE comvert, ONLY: ap, bp
17      ! GCM LMD nouvelle grille      USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, &
18             nday, offline, periodav
19      ! ... Dans inigeom, nouveaux calculs pour les elongations cu, cv      USE dimens_m, ONLY: iim, llm, nqmx
20      ! et possibilite d'appeler une fonction f(y) a derivee tangente      USE dynetat0_m, ONLY: day_ini
21      ! hyperbolique a la place de la fonction a derivee sinusoidale.      USE exner_hyb_m, ONLY: exner_hyb
22        USE guide_m, ONLY: guide
23      ! ... Possibilité de choisir le schéma pour l'advection de      use inidissip_m, only: idissip
24      ! q, en modifiant iadv dans "traceur.def" (10/02) .      USE logic, ONLY: iflag_phys, ok_guide
25        USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1
26      ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99)      USE pression_m, ONLY: pression
27      ! Pour Van-Leer iadv=10      USE pressure_var, ONLY: p3d
28        USE temps, ONLY: dt, itaufin
     use dimens_m, only: iim, jjm, llm, nqmx  
     use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2  
     use comconst, only: dtvr, daysec, dtphys  
     use comvert, only: ap, bp  
     use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &  
          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  
   
     integer, intent(in):: nq  
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 time_0      REAL, intent(in):: time_0
39    
40      ! Variables local to the procedure:      ! Variables local to the procedure:
41    
# Line 93  contains Line 72  contains
72      REAL tppn(iim), tpps(iim), tpn, tps      REAL tppn(iim), tpps(iim), tpn, tps
73    
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
     integer itaufinp1  
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    
81      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL:: callinigrads = .true.  
82    
83      !+jld variables test conservation energie      ! Variables test conservation energie
84      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
85      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
86      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
# Line 117  contains Line 91  contains
91      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
92    
93      logical:: dissip_conservative = .true.      logical:: dissip_conservative = .true.
     LOGICAL:: prem = .true.  
94      logical forward, leapf, apphys, conser, apdiss      logical forward, leapf, apphys, conser, apdiss
95    
96      !---------------------------------------------------      !---------------------------------------------------
# Line 125  contains Line 98  contains
98      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
99    
100      itaufin = nday * day_step      itaufin = nday * day_step
     itaufinp1 = itaufin + 1  
   
101      itau = 0      itau = 0
102      iday = day_ini      iday = day_ini
103      time = time_0      time = time_0
104      IF (time > 1.) THEN      dq = 0.
        time = time - 1.  
        iday = iday + 1  
     ENDIF  
   
105      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
     dq=0.  
106      CALL pression(ip1jmp1, ap, bp, ps, p3d)      CALL pression(ip1jmp1, ap, bp, ps, p3d)
107      CALL exner_hyb(ps, p3d, pks, pk, pkf)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
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 dernieres heures.'         tetam1 = teta
116         endif         massem1 = masse
117           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)  
   
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 197  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 208  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 217  contains Line 165  contains
165                       teta, q(:, :, 1), q(:, :, 2))                       teta, q(:, :, 1), q(:, :, 2))
166               ENDIF               ENDIF
167    
168               CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
169                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
170                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
171    
# Line 237  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 260  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 279  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 304  contains Line 250  contains
250               ENDIF               ENDIF
251            ENDIF            ENDIF
252    
253            IF (itau == itaufinp1) exit outer_loop            IF (itau == itaufin + 1) exit outer_loop
   
           ! ecriture du fichier histoire moyenne:  
254    
           ! 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 318  contains Line 261  contains
261            ENDIF            ENDIF
262    
263            IF (itau == itaufin) THEN            IF (itau == itaufin) THEN
264               CALL dynredem1("restart.nc", 0., 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.22  
changed lines
  Added in v.26
  ViewVC Help
Powered by ViewVC 1.1.21