/[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 5 by guez, Mon Mar 3 16:32:04 2008 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, clesphy0, &    SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
        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      ! ... Possibilite de choisir le shema 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, llm, nqmx  
     use paramet_m, only: ip1jmp1, ip1jm, llmp1, 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, apdiss, apphys, conser, forward, iflag_phys, &  
          leapf, statcl  
     use comgeom  
     use serre  
     use temps, only: itaufin, day_ini, dt  
     use iniprint, only: prt_level  
     use com_io_dyn  
     use abort_gcm_m, only: abort_gcm  
     use ener  
     use calfis_m, only: calfis  
     use exner_hyb_m, only: exner_hyb  
     use guide_m, only: guide  
     use pression_m, only: pression  
   
     integer nq  
   
     INTEGER longcles  
     PARAMETER (longcles = 20)  
     REAL clesphy0(longcles)  
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
33        REAL ps(ip1jmp1) ! pression au sol, en Pa
34    
35        REAL masse(ip1jmp1, llm) ! masse d'air
36        REAL phis(ip1jmp1) ! geopotentiel au sol
37      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
38      REAL ps(ip1jmp1) ! pression au sol      REAL, intent(in):: time_0
39      REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches  
40        ! Variables local to the procedure:
41    
42        ! Variables dynamiques:
43    
44      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
45      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      REAL pk(ip1jmp1, llm) ! exner au milieu des couches
46      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches      REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
     REAL masse(ip1jmp1, llm) ! masse d'air  
     REAL phis(ip1jmp1) ! geopotentiel au sol  
47      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(ip1jmp1, llm) ! geopotential
48      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
49    
# Line 93  contains Line 71  contains
71    
72      REAL tppn(iim), tpps(iim), tpn, tps      REAL tppn(iim), tpps(iim), tpn, tps
73    
74      INTEGER itau, itaufinp1      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 ! Heure de la journee en fraction d'1 jour      REAL time ! time of day, as a fraction of day length
77        real finvmaold(ip1jmp1, llm)
78      REAL SSUM      LOGICAL:: lafin=.false.
     REAL time_0, finvmaold(ip1jmp1, llm)  
   
     LOGICAL :: lafin=.false.  
79      INTEGER ij, l      INTEGER ij, l
80    
81      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
   
     data 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 116  contains Line 88  contains
88      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(ip1jmp1, llm)
89      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
90      CHARACTER*15 ztit      CHARACTER*15 ztit
91      INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
92      SAVE ip_ebil_dyn  
93      DATA ip_ebil_dyn /0/      logical:: dissip_conservative = .true.
94        logical forward, leapf, apphys, conser, apdiss
     character(len=*), parameter:: modname = "leapfrog"  
     character*80 abort_message  
   
     logical dissip_conservative  
     save dissip_conservative  
     data dissip_conservative /.true./  
   
     LOGICAL prem  
     save prem  
     DATA prem /.true./  
95    
96      !---------------------------------------------------      !---------------------------------------------------
97    
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 :
106      dq=0.      CALL pression(ip1jmp1, ap, bp, ps, p3d)
107      CALL pression(ip1jmp1, ap, bp, ps, p)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
     CALL exner_hyb(ps, p, pks, pk, pkf)  
108    
109      ! Debut de l'integration temporelle:      ! Debut de l'integration temporelle:
110      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            statcl = .FALSE.            apdiss = MOD(itau + 1, idissip) == 0
           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.  
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               CALL caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, pk)               ! calcul des tendances advection des traceurs (dont l'humidite)
138                 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 206  contains Line 145  contains
145    
146            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
147            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
148                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold)                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
149                   finvmaold, leapf)
           ! calcul des tendances physiques:  
150    
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, p)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
156               CALL exner_hyb(ps, p, pks, pk, pkf)               CALL exner_hyb(ps, p3d, pks, pk, pkf)
157    
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'
164                  CALL diagedyn(ztit, 2, 1, 1, dtphys &                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
165                       , ucov, vcov, ps, p, pk, 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, p, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
170                    clesphy0, dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
171    
172               ! ajout des tendances physiques:               ! ajout des tendances physiques:
173               CALL addfi(nqmx, dtphys, &               CALL addfi(nqmx, dtphys, &
# Line 240  contains Line 177  contains
177               ! Diagnostique de conservation de l'énergie : difference               ! Diagnostique de conservation de l'énergie : difference
178               IF (ip_ebil_dyn >= 1) THEN               IF (ip_ebil_dyn >= 1) THEN
179                  ztit = 'bil phys'                  ztit = 'bil phys'
180                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p, pk, &                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
181                       teta, q(:, :, 1), q(:, :, 2))                       teta, q(:, :, 1), q(:, :, 2))
182               ENDIF               ENDIF
183            ENDIF            ENDIF
184    
185            CALL pression(ip1jmp1, ap, bp, ps, p)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
186            CALL exner_hyb(ps, p, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
   
           ! dissipation horizontale et verticale des petites echelles:  
187    
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)
194    
195               ! dissipation               ! dissipation
196               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
197               ucov=ucov + dudis               ucov=ucov + dudis
198               vcov=vcov + dvdis               vcov=vcov + dvdis
199    
# Line 271  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 290  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 315  contains Line 250  contains
250               ENDIF               ENDIF
251            ENDIF            ENDIF
252    
253            IF (itau == itaufinp1) then            IF (itau == itaufin + 1) exit outer_loop
              abort_message = 'Simulation finished'  
              call abort_gcm(modname, abort_message, 0)  
           ENDIF  
   
           ! 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 332  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
# Line 350  contains Line 277  contains
277                  dt = 2. * dtvr                  dt = 2. * dtvr
278               END IF               END IF
279            ELSE            ELSE
280               ! ...... pas leapfrog .....               ! pas leapfrog
281               leapf = .TRUE.               leapf = .TRUE.
282               dt = 2. * dtvr               dt = 2. * dtvr
283            END IF            END IF
284         end do         end do
285      end do      end do outer_loop
286    
287    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
288    
Legend:
Removed from v.5  
changed lines
  Added in v.26
  ViewVC Help
Powered by ViewVC 1.1.21