/[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 13 by guez, Fri Jul 25 19:59:34 2008 UTC revision 28 by guez, Fri Mar 26 18:33:04 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)
   
     ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34  
8    
9      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34
10      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      ! Authors: P. Le Van, L. Fairhead, F. Hourdin
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, jjm, llm, nqmx      USE logic, ONLY: iflag_phys, ok_guide
27      use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2      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: 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  
   
     integer nq  
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 92  contains Line 73  contains
73    
74      REAL tppn(iim), tpps(iim), tpn, tps      REAL tppn(iim), tpps(iim), tpn, tps
75    
76      INTEGER itau, itaufinp1      INTEGER itau ! index of the time step of the dynamics, starts at 0
77        INTEGER itaufin
78      INTEGER iday ! jour julien      INTEGER iday ! jour julien
79      REAL time ! Heure de la journee en fraction d'1 jour      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
     LOGICAL:: callinigrads = .true.  
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
94      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.      REAL dt
   
     logical:: dissip_conservative = .true.  
     LOGICAL:: prem = .true.  
     logical forward, leapf, apphys, conser, apdiss  
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 :
     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:      ! Début 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:            ! 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:  
   
126            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)            CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
   
127            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &            CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
128                 conser, du, dv, dteta, dp, w, pbaru, pbarv, &                 MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
129                 time + iday - day_ini)                 time + iday - day_ini)
130    
           ! calcul des tendances advection des traceurs (dont l'humidite)  
   
131            IF (forward .OR. leapf) THEN            IF (forward .OR. leapf) THEN
132                 ! Calcul des tendances advection des traceurs (dont l'humidité)
133               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134               IF (offline) THEN               IF (offline) THEN
135                  !maf stokage du flux de masse pour traceurs OFF-LINE                  ! Stokage du flux de masse pour traceurs off-line
136                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
137                       itau)                       itau)
138               ENDIF               ENDIF
# Line 194  contains Line 141  contains
141            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
142            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
143                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
144                 finvmaold, leapf)                 finvmaold, leapf, dt)
145    
146            ! calcul des tendances physiques:            IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
147                 ! calcul des tendances physiques:
           IF (apphys) THEN  
148               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
149    
150               CALL pression(ip1jmp1, ap, bp, ps, p3d)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
# Line 207  contains Line 153  contains
153               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
154               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
155    
156               ! Interface avec les routines de phylmd (phymars ...)               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
   
              ! 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  
   
              CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &  
157                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
158                    dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
159    
# Line 224  contains Line 161  contains
161               CALL addfi(nqmx, dtphys, &               CALL addfi(nqmx, dtphys, &
162                    ucov, vcov, teta, q, ps, &                    ucov, vcov, teta, q, ps, &
163                    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  
164            ENDIF            ENDIF
165    
166            CALL pression(ip1jmp1, ap, bp, ps, p3d)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
167            CALL exner_hyb(ps, p3d, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
168    
169            ! dissipation horizontale et verticale des petites echelles:            IF (MOD(itau + 1, idissip) == 0) THEN
170                 ! dissipation horizontale et verticale des petites echelles:
171    
           IF (apdiss) THEN  
172               ! calcul de l'energie cinetique avant dissipation               ! calcul de l'energie cinetique avant dissipation
173               call covcont(llm, ucov, vcov, ucont, vcont)               call covcont(llm, ucov, vcov, ucont, vcont)
174               call enercin(vcov, ucov, vcont, ucont, ecin0)               call enercin(vcov, ucov, vcont, ucont, ecin0)
# Line 248  contains Line 178  contains
178               ucov=ucov + dudis               ucov=ucov + dudis
179               vcov=vcov + dvdis               vcov=vcov + dvdis
180    
181               if (dissip_conservative) then               ! On rajoute la tendance due à la transformation Ec -> E
182                  ! On rajoute la tendance due a la transform. Ec -> E               ! thermique créée lors de la dissipation
183                  ! therm. cree lors de la dissipation               call covcont(llm, ucov, vcov, ucont, vcont)
184                  call covcont(llm, ucov, vcov, ucont, vcont)               call enercin(vcov, ucov, vcont, ucont, ecin)
185                  call enercin(vcov, ucov, vcont, ucont, ecin)               dtetaecdt= (ecin0 - ecin) / pk
186                  dtetaecdt= (ecin0 - ecin) / pk               dtetadis=dtetadis + dtetaecdt
                 dtetadis=dtetadis + dtetaecdt  
              endif  
187               teta=teta + dtetadis               teta=teta + dtetadis
188    
189               ! Calcul de la valeur moyenne, unique de h aux poles .....               ! Calcul de la valeur moyenne unique de h aux pôles
   
190               DO l = 1, llm               DO l = 1, llm
191                  DO ij = 1, iim                  DO ij = 1, iim
192                     tppn(ij) = aire(ij) * teta(ij, l)                     tppn(ij) = aire(ij) * teta(ij, l)
193                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)                     tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
194                  ENDDO                  ENDDO
195                  tpn = SSUM(iim, tppn, 1) / apoln                  tpn = SUM(tppn) / apoln
196                  tps = SSUM(iim, tpps, 1) / apols                  tps = SUM(tpps) / apols
197    
198                  DO ij = 1, iip1                  DO ij = 1, iip1
199                     teta(ij, l) = tpn                     teta(ij, l) = tpn
# Line 278  contains Line 205  contains
205                  tppn(ij) = aire(ij) * ps(ij)                  tppn(ij) = aire(ij) * ps(ij)
206                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)                  tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
207               ENDDO               ENDDO
208               tpn = SSUM(iim, tppn, 1) / apoln               tpn = SUM(tppn) / apoln
209               tps = SSUM(iim, tpps, 1) / apols               tps = SUM(tpps) / apols
210    
211               DO ij = 1, iip1               DO ij = 1, iip1
212                  ps(ij) = tpn                  ps(ij) = tpn
213                  ps(ij + ip1jm) = tps                  ps(ij + ip1jm) = tps
214               ENDDO               ENDDO
   
215            END IF            END IF
216    
217            ! 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 303  contains Line 229  contains
229               ENDIF               ENDIF
230            ENDIF            ENDIF
231    
232            IF (itau == itaufinp1) exit outer_loop            IF (itau == itaufin + 1) exit outer_loop
   
           ! ecriture du fichier histoire moyenne:  
233    
           ! Comment out the following calls when you do not want the output  
           ! files "dyn_hist_ave.nc" and "dynzon.nc"  
234            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN            IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
235                 ! ecriture du fichier histoire moyenne:
236               CALL writedynav(histaveid, nqmx, itau, vcov, &               CALL writedynav(histaveid, nqmx, itau, vcov, &
237                    ucov, teta, pk, phi, q, masse, ps, phis)                    ucov, teta, pk, phi, q, masse, ps, phis)
238               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &               call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
# Line 317  contains Line 240  contains
240            ENDIF            ENDIF
241    
242            IF (itau == itaufin) THEN            IF (itau == itaufin) THEN
243               CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
244               CLOSE(99)                    itau=itau_dyn+itaufin)
245            ENDIF            ENDIF
246    
247            ! gestion de l'integration temporelle:            ! gestion de l'integration temporelle:
   
248            IF (MOD(itau, iperiod) == 0) exit            IF (MOD(itau, iperiod) == 0) exit
249            IF (MOD(itau - 1, iperiod) == 0) THEN            IF (MOD(itau - 1, iperiod) == 0) THEN
250               IF (forward) THEN               IF (forward) THEN
# Line 335  contains Line 257  contains
257                  dt = 2. * dtvr                  dt = 2. * dtvr
258               END IF               END IF
259            ELSE            ELSE
260               ! ...... pas leapfrog .....               ! pas leapfrog
261               leapf = .TRUE.               leapf = .TRUE.
262               dt = 2. * dtvr               dt = 2. * dtvr
263            END IF            END IF
Legend:
Removed from v.13  
changed lines
  Added in v.28
  ViewVC Help
Powered by ViewVC 1.1.21