/[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 5 by guez, Mon Mar 3 16:32:04 2008 UTC revision 24 by guez, Wed Mar 3 13:23:49 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    
11      ! Version du 10/01/98, avec coordonnees verticales hybrides, avec      ! Version du 10/01/98, avec coordonnées verticales hybrides, avec
12      ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)      ! nouveaux opérateurs dissipation "*" (gradiv2, divgrad2, nxgraro2)
13    
14      ! Auteur: P. Le Van /L. Fairhead/F.Hourdin      ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin
15      ! Objet:      ! Objet: nouvelle grille
     ! GCM LMD nouvelle grille  
16    
17      ! ... Dans inigeom, nouveaux calculs pour les elongations cu, cv      ! Dans "inigeom", nouveaux calculs pour les élongations cu, cv
18      ! et possibilite d'appeler une fonction f(y) a derivee tangente      ! et possibilité d'appeler une fonction f(y) à dérivée tangente
19      ! hyperbolique a la place de la fonction a derivee sinusoidale.      ! hyperbolique à la place de la fonction à dérivée sinusoïdale.
20    
21      ! ... Possibilite de choisir le shema pour l'advection de      ! Possibilité de choisir le schéma pour l'advection de
22      ! q, en modifiant iadv dans "traceur.def" (10/02) .      ! q, en modifiant iadv dans "traceur.def".
23    
24      ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99)      ! Pour Van-Leer + vapeur d'eau saturée, iadv(1)=4.
25      ! Pour Van-Leer iadv=10      ! Pour Van-Leer iadv=10
26    
27      use dimens_m, only: iim, llm, nqmx      use dimens_m, only: iim, llm, nqmx
28      use paramet_m, only: ip1jmp1, ip1jm, llmp1, ijmllm, ijp1llm, jjp1, iip1, &      use paramet_m, only: ip1jmp1, ip1jm, ijp1llm, jjp1, iip1
          iip2  
29      use comconst, only: dtvr, daysec, dtphys      use comconst, only: dtvr, daysec, dtphys
30      use comvert, only: ap, bp      use comvert, only: ap, bp
31      use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &      use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &
32           offline, periodav           offline, periodav
33      use logic, only: ok_guide, apdiss, apphys, conser, forward, iflag_phys, &      use logic, only: ok_guide, iflag_phys
          leapf, statcl  
34      use comgeom      use comgeom
35      use serre      use serre
36      use temps, only: itaufin, day_ini, dt      use temps, only: itaufin, day_ini, dt
37      use iniprint, only: prt_level      use iniprint, only: prt_level
38      use com_io_dyn      use com_io_dyn
     use abort_gcm_m, only: abort_gcm  
39      use ener      use ener
40      use calfis_m, only: calfis      use calfis_m, only: calfis
41      use exner_hyb_m, only: exner_hyb      use exner_hyb_m, only: exner_hyb
42      use guide_m, only: guide      use guide_m, only: guide
43      use pression_m, only: pression      use pression_m, only: pression
44        use pressure_var, only: p3d
45    
46      integer nq      ! Variables dynamiques:
   
     INTEGER longcles  
     PARAMETER (longcles = 20)  
     REAL clesphy0(longcles)  
   
     ! variables dynamiques  
47      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
48      REAL teta(ip1jmp1, llm) ! temperature potentielle      REAL teta(ip1jmp1, llm) ! temperature potentielle
49      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields      REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
50      REAL ps(ip1jmp1) ! pression au sol      REAL ps(ip1jmp1) ! pression au sol, en Pa
51      REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches      REAL masse(ip1jmp1, llm) ! masse d'air
52        REAL phis(ip1jmp1) ! geopotentiel au sol
53    
54        REAL, intent(in):: time_0
55    
56        ! Variables local to the procedure:
57    
58        ! Variables dynamiques:
59    
60      REAL pks(ip1jmp1) ! exner au sol      REAL pks(ip1jmp1) ! exner au sol
61      REAL pk(ip1jmp1, llm) ! exner au milieu des couches      REAL pk(ip1jmp1, llm) ! exner au milieu des couches
62      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  
63      REAL phi(ip1jmp1, llm) ! geopotential      REAL phi(ip1jmp1, llm) ! geopotential
64      REAL w(ip1jmp1, llm) ! vitesse verticale      REAL w(ip1jmp1, llm) ! vitesse verticale
65    
# Line 93  contains Line 87  contains
87    
88      REAL tppn(iim), tpps(iim), tpn, tps      REAL tppn(iim), tpps(iim), tpn, tps
89    
90      INTEGER itau, itaufinp1      INTEGER itau ! index of the time step of the dynamics, starts at 0
91      INTEGER iday ! jour julien      INTEGER iday ! jour julien
92      REAL time ! Heure de la journee en fraction d'1 jour      REAL time ! time of day, as a fraction of day length
93    
94      REAL SSUM      REAL SSUM
95      REAL time_0, finvmaold(ip1jmp1, llm)      real finvmaold(ip1jmp1, llm)
96    
97      LOGICAL :: lafin=.false.      LOGICAL :: lafin=.false.
98      INTEGER ij, l      INTEGER ij, l
99    
100      REAL rdayvrai, rdaym_ini      REAL rdayvrai, rdaym_ini
     LOGICAL callinigrads  
101    
102      data callinigrads/.true./      ! Variables test conservation energie
   
     !+jld variables test conservation energie  
103      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)      REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
104      ! Tendance de la temp. potentiel d (theta) / d t due a la      ! Tendance de la temp. potentiel d (theta) / d t due a la
105      ! tansformation d'energie cinetique en energie thermique      ! tansformation d'energie cinetique en energie thermique
# Line 116  contains Line 107  contains
107      REAL dtetaecdt(ip1jmp1, llm)      REAL dtetaecdt(ip1jmp1, llm)
108      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)      REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
109      CHARACTER*15 ztit      CHARACTER*15 ztit
110      INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.      INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
111      SAVE ip_ebil_dyn  
112      DATA ip_ebil_dyn /0/      logical:: dissip_conservative = .true.
113        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./  
114    
115      !---------------------------------------------------      !---------------------------------------------------
116    
117      print *, "Call sequence information: leapfrog"      print *, "Call sequence information: leapfrog"
118    
119      itaufin = nday * day_step      itaufin = nday * day_step
     itaufinp1 = itaufin + 1  
   
120      itau = 0      itau = 0
121      iday = day_ini      iday = day_ini
122      time = time_0      time = time_0
# Line 146  contains Line 125  contains
125         iday = iday + 1         iday = iday + 1
126      ENDIF      ENDIF
127    
128        dq = 0.
129      ! On initialise la pression et la fonction d'Exner :      ! On initialise la pression et la fonction d'Exner :
130      dq=0.      CALL pression(ip1jmp1, ap, bp, ps, p3d)
131      CALL pression(ip1jmp1, ap, bp, ps, p)      CALL exner_hyb(ps, p3d, pks, pk, pkf)
     CALL exner_hyb(ps, p, pks, pk, pkf)  
132    
133      ! Debut de l'integration temporelle:      ! Debut de l'integration temporelle:
134      do      outer_loop:do
135         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then         if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600.) then
136            call guide(itau, ucov, vcov, teta, q, masse, ps)            call guide(itau, ucov, vcov, teta, q, masse, ps)
137         else         else
138            IF (prt_level > 9) print *, &            IF (prt_level > 9) print *, &
139                 'Attention : on ne guide pas les 6 dernieres heures.'                 'Attention : on ne guide pas les 6 dernières heures.'
140         endif         endif
141    
142         CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)         CALL SCOPY(ip1jm * llm, vcov, 1, vcovm1, 1)
143         CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)         CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
144         CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)         CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
145         CALL SCOPY(ijp1llm, masse, 1, massem1, 1)         CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
# Line 177  contains Line 156  contains
156            ! gestion des appels de la physique et des dissipations:            ! gestion des appels de la physique et des dissipations:
157    
158            apphys = .FALSE.            apphys = .FALSE.
           statcl = .FALSE.  
159            conser = .FALSE.            conser = .FALSE.
160            apdiss = .FALSE.            apdiss = .FALSE.
161    
# Line 196  contains Line 174  contains
174            ! calcul des tendances advection des traceurs (dont l'humidite)            ! calcul des tendances advection des traceurs (dont l'humidite)
175    
176            IF (forward .OR. leapf) THEN            IF (forward .OR. leapf) THEN
177               CALL caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, pk)               CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
178               IF (offline) THEN               IF (offline) THEN
179                  !maf stokage du flux de masse pour traceurs OFF-LINE                  !maf stokage du flux de masse pour traceurs OFF-LINE
180                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &                  CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
# Line 206  contains Line 184  contains
184    
185            ! integrations dynamique et traceurs:            ! integrations dynamique et traceurs:
186            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &            CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
187                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold)                 dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
188                   finvmaold, leapf)
189    
190            ! calcul des tendances physiques:            ! calcul des tendances physiques:
191    
192            IF (apphys) THEN            IF (apphys) THEN
193               IF (itau + 1 == itaufin) lafin = .TRUE.               IF (itau + 1 == itaufin) lafin = .TRUE.
194    
195               CALL pression(ip1jmp1, ap, bp, ps, p)               CALL pression(ip1jmp1, ap, bp, ps, p3d)
196               CALL exner_hyb(ps, p, pks, pk, pkf)               CALL exner_hyb(ps, p3d, pks, pk, pkf)
197    
198               rdaym_ini = itau * dtvr / daysec               rdaym_ini = itau * dtvr / daysec
199               rdayvrai = rdaym_ini + day_ini               rdayvrai = rdaym_ini + day_ini
# Line 224  contains Line 203  contains
203               ! Diagnostique de conservation de l'énergie : initialisation               ! Diagnostique de conservation de l'énergie : initialisation
204               IF (ip_ebil_dyn >= 1) THEN               IF (ip_ebil_dyn >= 1) THEN
205                  ztit='bil dyn'                  ztit='bil dyn'
206                  CALL diagedyn(ztit, 2, 1, 1, dtphys &                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
207                       , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))                       teta, q(:, :, 1), q(:, :, 2))
208               ENDIF               ENDIF
209    
210               CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &               CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
211                    masse, ps, p, pk, phis, phi, du, dv, dteta, dq, w, &                    masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
212                    clesphy0, dufi, dvfi, dtetafi, dqfi, dpfi)                    dufi, dvfi, dtetafi, dqfi, dpfi)
213    
214               ! ajout des tendances physiques:               ! ajout des tendances physiques:
215               CALL addfi(nqmx, dtphys, &               CALL addfi(nqmx, dtphys, &
# Line 240  contains Line 219  contains
219               ! Diagnostique de conservation de l'énergie : difference               ! Diagnostique de conservation de l'énergie : difference
220               IF (ip_ebil_dyn >= 1) THEN               IF (ip_ebil_dyn >= 1) THEN
221                  ztit = 'bil phys'                  ztit = 'bil phys'
222                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p, pk, &                  CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
223                       teta, q(:, :, 1), q(:, :, 2))                       teta, q(:, :, 1), q(:, :, 2))
224               ENDIF               ENDIF
225            ENDIF            ENDIF
226    
227            CALL pression(ip1jmp1, ap, bp, ps, p)            CALL pression(ip1jmp1, ap, bp, ps, p3d)
228            CALL exner_hyb(ps, p, pks, pk, pkf)            CALL exner_hyb(ps, p3d, pks, pk, pkf)
229    
230            ! dissipation horizontale et verticale des petites echelles:            ! dissipation horizontale et verticale des petites echelles:
231    
# Line 256  contains Line 235  contains
235               call enercin(vcov, ucov, vcont, ucont, ecin0)               call enercin(vcov, ucov, vcont, ucont, ecin0)
236    
237               ! dissipation               ! dissipation
238               CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)               CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
239               ucov=ucov + dudis               ucov=ucov + dudis
240               vcov=vcov + dvdis               vcov=vcov + dvdis
241    
# Line 315  contains Line 294  contains
294               ENDIF               ENDIF
295            ENDIF            ENDIF
296    
297            IF (itau == itaufinp1) then            IF (itau == itaufin + 1) exit outer_loop
              abort_message = 'Simulation finished'  
              call abort_gcm(modname, abort_message, 0)  
           ENDIF  
298    
299            ! ecriture du fichier histoire moyenne:            ! ecriture du fichier histoire moyenne:
300    
# Line 332  contains Line 308  contains
308            ENDIF            ENDIF
309    
310            IF (itau == itaufin) THEN            IF (itau == itaufin) THEN
311               CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)               CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)
312               CLOSE(99)               CLOSE(99)
313            ENDIF            ENDIF
314    
# Line 350  contains Line 326  contains
326                  dt = 2. * dtvr                  dt = 2. * dtvr
327               END IF               END IF
328            ELSE            ELSE
329               ! ...... pas leapfrog .....               ! pas leapfrog
330               leapf = .TRUE.               leapf = .TRUE.
331               dt = 2. * dtvr               dt = 2. * dtvr
332            END IF            END IF
333         end do         end do
334      end do      end do outer_loop
335    
336    END SUBROUTINE leapfrog    END SUBROUTINE leapfrog
337    
Legend:
Removed from v.5  
changed lines
  Added in v.24
  ViewVC Help
Powered by ViewVC 1.1.21