/[lmdze]/trunk/Sources/phylmd/physiq.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/physiq.f

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

trunk/libf/phylmd/physiq.f90 revision 61 by guez, Fri Apr 20 14:58:43 2012 UTC trunk/phylmd/physiq.f revision 90 by guez, Wed Mar 12 21:16:36 2014 UTC
# Line 5  module physiq_m Line 5  module physiq_m
5  contains  contains
6    
7    SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &    SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &
8         u, v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps, dudyn, PVteta)         u, v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps)
9    
10        ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11        ! (subversion revision 678)
12    
     ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28 (SVN revision 678)  
13      ! Author: Z.X. Li (LMD/CNRS) 1993      ! Author: Z.X. Li (LMD/CNRS) 1993
14    
15      ! This is the main procedure for the "physics" part of the program.      ! This is the main procedure for the "physics" part of the program.
16    
17      use aaam_bud_m, only: aaam_bud      use aaam_bud_m, only: aaam_bud
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
19        use aeropt_m, only: aeropt
20      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
21      USE calendar, ONLY: ymds2ju      USE calendar, ONLY: ymds2ju
22      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
# Line 22  contains Line 25  contains
25      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
26           ok_orodr, ok_orolf, soil_model           ok_orodr, ok_orolf, soil_model
27      USE clmain_m, ONLY: clmain      USE clmain_m, ONLY: clmain
28        use clouds_gno_m, only: clouds_gno
29      USE comgeomphy, ONLY: airephy, cuphy, cvphy      USE comgeomphy, ONLY: airephy, cuphy, cvphy
30      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
31      USE conf_gcm_m, ONLY: offline, raz_date      USE conf_gcm_m, ONLY: offline, raz_date
32      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
33        use conflx_m, only: conflx
34      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
35      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
36      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
37      USE dimens_m, ONLY: iim, jjm, llm, nqmx      use diagphy_m, only: diagphy
38        USE dimens_m, ONLY: llm, nqmx
39      USE dimphy, ONLY: klon, nbtr      USE dimphy, ONLY: klon, nbtr
40      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
41      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
42      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
43        use fisrtilp_m, only: fisrtilp
44      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
     USE histsync_m, ONLY: histsync  
     USE histwrite_m, ONLY: histwrite  
45      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
46           nbsrf           nbsrf
     USE ini_histhf_m, ONLY: ini_histhf  
     USE ini_histday_m, ONLY: ini_histday  
47      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins
48        use newmicro_m, only: newmicro
49      USE oasis_m, ONLY: ok_oasis      USE oasis_m, ONLY: ok_oasis
50      USE orbite_m, ONLY: orbite, zenang      USE orbite_m, ONLY: orbite, zenang
51      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
# Line 51  contains Line 55  contains
55      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
56      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
57      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
58        use readsulfate_m, only: readsulfate
59      use sugwd_m, only: sugwd      use sugwd_m, only: sugwd
60      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
61      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE temps, ONLY: annee_ref, day_ref, itau_phy
62        use unit_nml_m, only: unit_nml
63      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
64    
65      ! Arguments:      ! Arguments:
# Line 61  contains Line 67  contains
67      REAL, intent(in):: rdayvrai      REAL, intent(in):: rdayvrai
68      ! (elapsed time since January 1st 0h of the starting year, in days)      ! (elapsed time since January 1st 0h of the starting year, in days)
69    
70      REAL, intent(in):: time ! heure de la journée en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
71      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)
72      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
73    
# Line 83  contains Line 89  contains
89      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(klon, llm) ! input temperature (K)
90    
91      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(klon, llm, nqmx)
92      ! (humidité spécifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
93    
94      REAL omega(klon, llm) ! input vitesse verticale en Pa/s      REAL omega(klon, llm) ! input vitesse verticale en Pa/s
95      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)
# Line 97  contains Line 103  contains
103      INTEGER nbteta      INTEGER nbteta
104      PARAMETER(nbteta = 3)      PARAMETER(nbteta = 3)
105    
     REAL PVteta(klon, nbteta)  
     ! (output vorticite potentielle a des thetas constantes)  
   
     LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE  
     PARAMETER (ok_cvl = .TRUE.)  
106      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
107      PARAMETER (ok_gust = .FALSE.)      PARAMETER (ok_gust = .FALSE.)
108    
# Line 116  contains Line 117  contains
117      logical rnpb      logical rnpb
118      parameter(rnpb = .true.)      parameter(rnpb = .true.)
119    
120      character(len = 6), save:: ocean      character(len = 6):: ocean = 'force '
121      ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")      ! (type de mod\`ele oc\'ean \`a utiliser: "force" ou "slab" mais pas "couple")
   
     logical ok_ocean  
     SAVE ok_ocean  
122    
123      ! "slab" ocean      ! "slab" ocean
124      REAL, save:: tslab(klon) ! temperature of ocean slab      REAL, save:: tslab(klon) ! temperature of ocean slab
# Line 129  contains Line 127  contains
127      REAL fluxg(klon) ! flux turbulents ocean-atmosphere      REAL fluxg(klon) ! flux turbulents ocean-atmosphere
128    
129      ! Modele thermique du sol, a activer pour le cycle diurne:      ! Modele thermique du sol, a activer pour le cycle diurne:
130      logical, save:: ok_veget      logical:: ok_veget = .false. ! type de modele de vegetation utilise
     LOGICAL, save:: ok_journe ! sortir le fichier journalier  
   
     LOGICAL ok_mensuel ! sortir le fichier mensuel  
131    
132      LOGICAL ok_instan ! sortir le fichier instantane      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
133      save ok_instan      ! sorties journalieres, mensuelles et instantanees dans les
134        ! fichiers histday, histmth et histins
135    
136      LOGICAL ok_region ! sortir le fichier regional      LOGICAL ok_region ! sortir le fichier regional
137      PARAMETER (ok_region = .FALSE.)      PARAMETER (ok_region = .FALSE.)
# Line 167  contains Line 163  contains
163    
164      !MI Amip2 PV a theta constante      !MI Amip2 PV a theta constante
165    
166      INTEGER klevp1      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
167      PARAMETER(klevp1 = llm + 1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
   
     REAL swdn0(klon, klevp1), swdn(klon, klevp1)  
     REAL swup0(klon, klevp1), swup(klon, klevp1)  
168      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
169    
170      REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
171      REAL lwup0(klon, klevp1), lwup(klon, klevp1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
172      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
173    
174      !IM Amip2      !IM Amip2
# Line 206  contains Line 199  contains
199      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)
200    
201      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
202      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./      DATA zx_tau/0., 0.3, 1.3, 3.6, 9.4, 23., 60./
203      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
204    
205      ! cldtopres pression au sommet des nuages      ! cldtopres pression au sommet des nuages
# Line 268  contains Line 261  contains
261      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
262      ! soil temperature of surface fraction      ! soil temperature of surface fraction
263    
264      REAL fevap(klon, nbsrf)      REAL, save:: fevap(klon, nbsrf) ! evaporation
     SAVE fevap ! evaporation  
265      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
266      SAVE fluxlat      SAVE fluxlat
267    
# Line 286  contains Line 278  contains
278      REAL falblw(klon, nbsrf)      REAL falblw(klon, nbsrf)
279      SAVE falblw ! albedo par type de surface      SAVE falblw ! albedo par type de surface
280    
281      ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
282      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
283      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
284      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 316  contains Line 308  contains
308      SAVE Ma      SAVE Ma
309      REAL qcondc(klon, llm) ! in-cld water content from convect      REAL qcondc(klon, llm) ! in-cld water content from convect
310      SAVE qcondc      SAVE qcondc
311      REAL ema_work1(klon, llm), ema_work2(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
312      SAVE ema_work1, ema_work2      REAL, save:: wd(klon)
   
     REAL wd(klon) ! sb  
     SAVE wd ! sb  
313    
314      ! Variables locales pour la couche limite (al1):      ! Variables locales pour la couche limite (al1):
315    
# Line 329  contains Line 318  contains
318      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
319      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
320    
321      !AA Pour phytrac      ! Pour phytrac :
322      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
323      REAL yu1(klon) ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
324      REAL yv1(klon) ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
# Line 348  contains Line 337  contains
337      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
338      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
339    
340      !AA      REAL, save:: rain_fall(klon) ! pluie
341      REAL rain_fall(klon) ! pluie      REAL, save:: snow_fall(klon) ! neige
342      REAL snow_fall(klon) ! neige  
     save snow_fall, rain_fall  
     !IM cf FH pour Tiedtke 080604  
343      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
344    
345      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon), devap(klon) ! evaporation and its derivative
346      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
347      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
348      SAVE dlw      SAVE dlw
# Line 376  contains Line 363  contains
363      INTEGER julien      INTEGER julien
364    
365      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
366      REAL pctsrf(klon, nbsrf)      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
367      !IM      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE  
368    
     SAVE pctsrf ! sous-fraction du sol  
369      REAL albsol(klon)      REAL albsol(klon)
370      SAVE albsol ! albedo du sol total      SAVE albsol ! albedo du sol total
371      REAL albsollw(klon)      REAL albsollw(klon)
# Line 393  contains Line 378  contains
378      EXTERNAL alboc ! calculer l'albedo sur ocean      EXTERNAL alboc ! calculer l'albedo sur ocean
379      !KE43      !KE43
380      EXTERNAL conema3 ! convect4.3      EXTERNAL conema3 ! convect4.3
     EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)  
381      EXTERNAL nuage ! calculer les proprietes radiatives      EXTERNAL nuage ! calculer les proprietes radiatives
382      EXTERNAL transp ! transport total de l'eau et de l'energie      EXTERNAL transp ! transport total de l'eau et de l'energie
383    
384      ! Variables locales      ! Variables locales
385    
386      real clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
387      real clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
   
     save rnebcon, clwcon  
388    
389      REAL rhcl(klon, llm) ! humiditi relative ciel clair      REAL rhcl(klon, llm) ! humiditi relative ciel clair
390      REAL dialiq(klon, llm) ! eau liquide nuageuse      REAL dialiq(klon, llm) ! eau liquide nuageuse
# Line 422  contains Line 404  contains
404      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
405      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
406    
407      ! Le rayonnement n'est pas calcule tous les pas, il faut donc      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
408      ! que les variables soient rémanentes      ! les variables soient r\'emanentes.
409      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
410      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL heat0(klon, llm) ! chauffage solaire ciel clair
411      REAL cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
412      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
413      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
414      real sollwdown(klon) ! downward LW flux at surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
415      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      real, save:: sollwdown(klon) ! downward LW flux at surface
416        REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
417      REAL albpla(klon)      REAL albpla(klon)
418      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
419      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
420      SAVE cool, albpla, topsw, toplw, solsw, sollw, sollwdown      SAVE albpla
421      SAVE topsw0, toplw0, solsw0, sollw0, heat0, cool0      SAVE heat0, cool0
422    
423      INTEGER itaprad      INTEGER itaprad
424      SAVE itaprad      SAVE itaprad
# Line 451  contains Line 434  contains
434      REAL dist, rmu0(klon), fract(klon)      REAL dist, rmu0(klon), fract(klon)
435      REAL zdtime ! pas de temps du rayonnement (s)      REAL zdtime ! pas de temps du rayonnement (s)
436      real zlongi      real zlongi
   
437      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
438      REAL za, zb      REAL za, zb
439      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zdelta, zcor
440      real zqsat(klon, llm)      real zqsat(klon, llm)
441      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
442      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup = 234.0)  
   
443      REAL zphi(klon, llm)      REAL zphi(klon, llm)
444    
445      !IM cf. AM Variables locales pour la CLA (hbtm2)      !IM cf. AM Variables locales pour la CLA (hbtm2)
# Line 482  contains Line 460  contains
460      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
461      REAL s_trmb3(klon)      REAL s_trmb3(klon)
462    
463      ! Variables locales pour la convection de K. Emanuel (sb):      ! Variables locales pour la convection de K. Emanuel :
464    
465      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
466      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
# Line 498  contains Line 476  contains
476      REAL rflag(klon) ! flag fonctionnement de convect      REAL rflag(klon) ! flag fonctionnement de convect
477      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
478      ! -- convect43:      ! -- convect43:
     INTEGER ntra ! nb traceurs pour convect4.3  
479      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
480      REAL dplcldt(klon), dplcldr(klon)      REAL dplcldt(klon), dplcldr(klon)
481    
# Line 507  contains Line 484  contains
484      ! con: convection      ! con: convection
485      ! lsc: large scale condensation      ! lsc: large scale condensation
486      ! ajs: ajustement sec      ! ajs: ajustement sec
487      ! eva: évaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
488      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
489      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
490      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
# Line 516  contains Line 493  contains
493      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
494      REAL rneb(klon, llm)      REAL rneb(klon, llm)
495    
496      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
497      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
498      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
499      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
500      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
501      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
502    
503      INTEGER,save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
504    
505      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
506      REAL snow_con(klon), snow_lsc(klon)      REAL snow_con(klon), snow_lsc(klon)
# Line 537  contains Line 514  contains
514      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
515      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
516    
517      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
518      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
519      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
520    
521      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
522      real, save:: fact_cldcon      real:: fact_cldcon = 0.375
523      real, save:: facttemps      real:: facttemps = 1.e-4
524      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
525      real facteur      real facteur
526    
527      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
528      logical ptconv(klon, llm)      logical ptconv(klon, llm)
529    
530      ! Variables locales pour effectuer les appels en série :      ! Variables locales pour effectuer les appels en s\'erie :
531    
532      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
533      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm), qs_seri(klon, llm)
# Line 569  contains Line 543  contains
543      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
544      REAL aam, torsfc      REAL aam, torsfc
545    
     REAL dudyn(iim + 1, jjm + 1, llm)  
   
546      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)  
547    
548      INTEGER, SAVE:: nid_day, nid_ins      INTEGER, SAVE:: nid_day, nid_ins
549    
# Line 583  contains Line 554  contains
554    
555      REAL zsto      REAL zsto
556    
     character(len = 20) modname  
     character(len = 80) abort_message  
557      logical ok_sync      logical ok_sync
558      real date0      real date0
559    
560      ! Variables liées au bilan d'énergie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
561      REAL ztsol(klon)      REAL ztsol(klon)
562      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
563      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
564      REAL fs_bound, fq_bound      REAL fs_bound, fq_bound
565      REAL zero_v(klon)      REAL zero_v(klon)
566      CHARACTER(LEN = 15) ztit      CHARACTER(LEN = 15) tit
567      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
568      INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
569    
570      REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique      REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
571      REAL ZRCPD      REAL ZRCPD
572    
573      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
574      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
575      REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille      REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
576      REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille      REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
577      !jq Aerosol effects (Johannes Quaas, 27/11/2003)  
578      REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]      ! Aerosol effects:
579    
580        REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
581    
582      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
583      ! (SO4 aerosol concentration, in ug/m3, pre-industrial value)      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
584    
585      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
586      ! (Cloud optical thickness for pre-industrial (pi) aerosols)      ! cloud optical thickness for pre-industrial (pi) aerosols
587    
588      REAL re(klon, llm) ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
589      REAL fl(klon, llm) ! denominator of re      REAL fl(klon, llm) ! denominator of re
590    
591      ! Aerosol optical properties      ! Aerosol optical properties
592      REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)      REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
593      REAL cg_ae(klon, llm, 2)      REAL, save:: cg_ae(klon, llm, 2)
594    
595      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.      REAL topswad(klon), solswad(klon) ! aerosol direct effect
596      ! ok_ade = True -ADE = topswad-topsw      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
   
     REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.  
     ! ok_aie = True ->  
     ! ok_ade = True -AIE = topswai-topswad  
     ! ok_ade = F -AIE = topswai-topsw  
597    
598      REAL aerindex(klon) ! POLDER aerosol index      REAL aerindex(klon) ! POLDER aerosol index
599    
600      ! Parameters      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
601      LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
602      REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)  
603        REAL:: bl95_b0 = 2., bl95_b1 = 0.2
604        ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
605        ! B). They link cloud droplet number concentration to aerosol mass
606        ! concentration.
607    
     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1  
608      SAVE u10m      SAVE u10m
609      SAVE v10m      SAVE v10m
610      SAVE t2m      SAVE t2m
611      SAVE q2m      SAVE q2m
612      SAVE ffonte      SAVE ffonte
613      SAVE fqcalving      SAVE fqcalving
     SAVE piz_ae  
     SAVE tau_ae  
     SAVE cg_ae  
614      SAVE rain_con      SAVE rain_con
615      SAVE snow_con      SAVE snow_con
616      SAVE topswai      SAVE topswai
# Line 653  contains Line 619  contains
619      SAVE solswad      SAVE solswad
620      SAVE d_u_con      SAVE d_u_con
621      SAVE d_v_con      SAVE d_v_con
     SAVE rnebcon0  
     SAVE clwcon0  
622    
623      real zmasse(klon, llm)      real zmasse(klon, llm)
624      ! (column-density of mass of air in a cell, in kg m-2)      ! (column-density of mass of air in a cell, in kg m-2)
625    
626      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
627    
628        namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
629             fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &
630             ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &
631             nsplit_thermals
632    
633      !----------------------------------------------------------------      !----------------------------------------------------------------
634    
635      modname = 'physiq'      IF (if_ebil >= 1) zero_v = 0.
     IF (if_ebil >= 1) THEN  
        DO i = 1, klon  
           zero_v(i) = 0.  
        END DO  
     END IF  
636      ok_sync = .TRUE.      ok_sync = .TRUE.
637      IF (nqmx < 2) THEN      IF (nqmx < 2) CALL abort_gcm('physiq', &
638         abort_message = 'eaux vapeur et liquide sont indispensables'           'eaux vapeur et liquide sont indispensables', 1)
        CALL abort_gcm(modname, abort_message, 1)  
     ENDIF  
639    
640      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
641         ! initialiser         ! initialiser
# Line 688  contains Line 650  contains
650         cg_ae = 0.         cg_ae = 0.
651         rain_con(:) = 0.         rain_con(:) = 0.
652         snow_con(:) = 0.         snow_con(:) = 0.
        bl95_b0 = 0.  
        bl95_b1 = 0.  
653         topswai(:) = 0.         topswai(:) = 0.
654         topswad(:) = 0.         topswad(:) = 0.
655         solswai(:) = 0.         solswai(:) = 0.
656         solswad(:) = 0.         solswad(:) = 0.
657    
658         d_u_con = 0.0         d_u_con = 0.
659         d_v_con = 0.0         d_v_con = 0.
660         rnebcon0 = 0.0         rnebcon0 = 0.
661         clwcon0 = 0.0         clwcon0 = 0.
662         rnebcon = 0.0         rnebcon = 0.
663         clwcon = 0.0         clwcon = 0.
664    
665         pblh =0. ! Hauteur de couche limite         pblh =0. ! Hauteur de couche limite
666         plcl =0. ! Niveau de condensation de la CLA         plcl =0. ! Niveau de condensation de la CLA
# Line 715  contains Line 675  contains
675    
676         IF (if_ebil >= 1) d_h_vcol_phy = 0.         IF (if_ebil >= 1) d_h_vcol_phy = 0.
677    
678         ! appel a la lecture du run.def physique         iflag_thermals = 0
679           nsplit_thermals = 1
680           print *, "Enter namelist 'physiq_nml'."
681           read(unit=*, nml=physiq_nml)
682           write(unit_nml, nml=physiq_nml)
683    
684         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &         call conf_phys
             ok_instan, fact_cldcon, facttemps, ok_newmicro, &  
             iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &  
             ok_ade, ok_aie, &  
             bl95_b0, bl95_b1, &  
             iflag_thermals, nsplit_thermals)  
685    
686         ! Initialiser les compteurs:         ! Initialiser les compteurs:
687    
# Line 731  contains Line 690  contains
690         itaprad = 0         itaprad = 0
691         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &
692              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &
693              snow_fall, solsw, sollwdown, dlw, radsol, frugs, agesno, zmea, &              snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &
694              zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &              zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
695              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
696    
697         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
698         q2 = 1.e-8         q2 = 1e-8
699    
700         radpas = NINT(86400. / dtphys / nbapp_rad)         radpas = NINT(86400. / dtphys / nbapp_rad)
701    
# Line 744  contains Line 703  contains
703         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
704    
705         PRINT *, 'cycle_diurne = ', cycle_diurne         PRINT *, 'cycle_diurne = ', cycle_diurne
706           CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &
707                ok_instan, ok_region)
708    
709         IF(ocean.NE.'force ') THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
710            ok_ocean = .TRUE.            print *, "Au minimum 4 appels par jour si cycle diurne"
711         ENDIF            call abort_gcm('physiq', &
712                   "Nombre d'appels au rayonnement insuffisant", 1)
        CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &  
             ok_region)  
   
        IF (dtphys*REAL(radpas) > 21600..AND.cycle_diurne) THEN  
           print *,'Nbre d appels au rayonnement insuffisant'  
           print *,"Au minimum 4 appels par jour si cycle diurne"  
           abort_message = 'Nbre d appels au rayonnement insuffisant'  
           call abort_gcm(modname, abort_message, 1)  
713         ENDIF         ENDIF
        print *,"Clef pour la convection, iflag_con = ", iflag_con  
        print *,"Clef pour le driver de la convection, ok_cvl = ", &  
             ok_cvl  
714    
715         ! Initialisation pour la convection de K.E. (sb):         ! Initialisation pour le sch\'ema de convection d'Emanuel :
716         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
717              ibas_con = 1
718            print *,"*** Convection de Kerry Emanuel 4.3 "            itop_con = 1
   
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG  
           DO i = 1, klon  
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END  
   
719         ENDIF         ENDIF
720    
721         IF (ok_orodr) THEN         IF (ok_orodr) THEN
# Line 797  contains Line 739  contains
739         npas = 0         npas = 0
740         nexca = 0         nexca = 0
741    
        print *,'AVANT HIST IFLAG_CON = ', iflag_con  
   
742         ! Initialisation des sorties         ! Initialisation des sorties
743    
        call ini_histhf(dtphys, nid_hf, nid_hf3d)  
        call ini_histday(dtphys, ok_journe, nid_day, nqmx)  
744         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
745         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
746         !XXXPB Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
747         WRITE(*, *) 'physiq date0: ', date0         print *, 'physiq date0: ', date0
748      ENDIF test_firstcal      ENDIF test_firstcal
749    
750      ! Mettre a zero des variables de sortie (pour securite)      ! Mettre a zero des variables de sortie (pour securite)
751    
752      DO i = 1, klon      DO i = 1, klon
753         d_ps(i) = 0.0         d_ps(i) = 0.
754      ENDDO      ENDDO
755      DO iq = 1, nqmx      DO iq = 1, nqmx
756         DO k = 1, llm         DO k = 1, llm
757            DO i = 1, klon            DO i = 1, klon
758               d_qx(i, k, iq) = 0.0               d_qx(i, k, iq) = 0.
759            ENDDO            ENDDO
760         ENDDO         ENDDO
761      ENDDO      ENDDO
# Line 825  contains Line 763  contains
763      mp = 0.      mp = 0.
764      phi = 0.      phi = 0.
765    
766      ! Ne pas affecter les valeurs entrées de u, v, h, et q :      ! Ne pas affecter les valeurs entr\'ees de u, v, h, et q :
767    
768      DO k = 1, llm      DO k = 1, llm
769         DO i = 1, klon         DO i = 1, klon
# Line 853  contains Line 791  contains
791      ENDDO      ENDDO
792    
793      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
794         ztit = 'after dynamics'         tit = 'after dynamics'
795         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
796              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
797              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
798         ! Comme les tendances de la physique sont ajoutés dans la         ! Comme les tendances de la physique sont ajout\'es dans la
799         !  dynamique, la variation d'enthalpie par la dynamique devrait         !  dynamique, la variation d'enthalpie par la dynamique devrait
800         !  être égale à la variation de la physique au pas de temps         !  \^etre \'egale \`a la variation de la physique au pas de temps
801         !  précédent.  Donc la somme de ces 2 variations devrait être         !  pr\'ec\'edent.  Donc la somme de ces 2 variations devrait \^etre
802         !  nulle.         !  nulle.
803         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
804              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
805              d_qt, 0., fs_bound, fq_bound)              d_qt, 0., fs_bound, fq_bound)
806      END IF      END IF
# Line 878  contains Line 816  contains
816      ELSE      ELSE
817         DO k = 1, llm         DO k = 1, llm
818            DO i = 1, klon            DO i = 1, klon
819               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
820               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
821            ENDDO            ENDDO
822         ENDDO         ENDDO
823         ancien_ok = .TRUE.         ancien_ok = .TRUE.
# Line 902  contains Line 840  contains
840    
841      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
842    
843      ! Mettre en action les conditions aux limites (albedo, sst, etc.).      ! Mettre en action les conditions aux limites (albedo, sst etc.).
844    
845      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
846      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
847    
848      ! Évaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
849      DO k = 1, llm      DO k = 1, llm
850         DO i = 1, klon         DO i = 1, klon
851            zb = MAX(0., ql_seri(i, k))            zb = MAX(0., ql_seri(i, k))
# Line 919  contains Line 857  contains
857      ql_seri = 0.      ql_seri = 0.
858    
859      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
860         ztit = 'after reevap'         tit = 'after reevap'
861         CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
862              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
863              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
864         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
865              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
866              fs_bound, fq_bound)              fs_bound, fq_bound)
867    
# Line 932  contains Line 870  contains
870      ! Appeler la diffusion verticale (programme de couche limite)      ! Appeler la diffusion verticale (programme de couche limite)
871    
872      DO i = 1, klon      DO i = 1, klon
873         zxrugs(i) = 0.0         zxrugs(i) = 0.
874      ENDDO      ENDDO
875      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
876         DO i = 1, klon         DO i = 1, klon
# Line 965  contains Line 903  contains
903         ENDDO         ENDDO
904      ENDDO      ENDDO
905    
906      ! Repartition sous maille des flux LW et SW      ! R\'epartition sous maille des flux longwave et shortwave
907      ! Repartition du longwave par sous-surface linearisee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
908    
909      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
910         DO i = 1, klon         DO i = 1, klon
911            fsollw(i, nsrf) = sollw(i) &            fsollw(i, nsrf) = sollw(i) &
912                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))
913            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))
914         ENDDO         ENDDO
915      ENDDO      ENDDO
916    
# Line 980  contains Line 918  contains
918    
919      ! Couche limite:      ! Couche limite:
920    
921      CALL clmain(dtphys, itap, date0, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &
922           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, &           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &
923           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
924           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &
925           rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder, rlon, rlat, &           rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &
926           cuphy, cvphy, frugs, firstcal, lafin, agesno, rugoro, d_t_vdf, &           frugs, firstcal, agesno, rugoro, d_t_vdf, &
927           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &
928           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &
929           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
930           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)
931    
932      ! Incrémentation des flux      ! Incr\'ementation des flux
933    
934      zxfluxt = 0.      zxfluxt = 0.
935      zxfluxq = 0.      zxfluxq = 0.
# Line 1000  contains Line 938  contains
938      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
939         DO k = 1, llm         DO k = 1, llm
940            DO i = 1, klon            DO i = 1, klon
941               zxfluxt(i, k) = zxfluxt(i, k) + &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
942                    fluxt(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
943               zxfluxq(i, k) = zxfluxq(i, k) + &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
944                    fluxq(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
              zxfluxu(i, k) = zxfluxu(i, k) + &  
                   fluxu(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxv(i, k) = zxfluxv(i, k) + &  
                   fluxv(i, k, nsrf) * pctsrf(i, nsrf)  
945            END DO            END DO
946         END DO         END DO
947      END DO      END DO
948      DO i = 1, klon      DO i = 1, klon
949         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
950         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
951         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
952      ENDDO      ENDDO
953    
# Line 1027  contains Line 961  contains
961      ENDDO      ENDDO
962    
963      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
964         ztit = 'after clmain'         tit = 'after clmain'
965         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
966              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
967              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
968         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
969              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
970              fs_bound, fq_bound)              fs_bound, fq_bound)
971      END IF      END IF
# Line 1039  contains Line 973  contains
973      ! Update surface temperature:      ! Update surface temperature:
974    
975      DO i = 1, klon      DO i = 1, klon
976         zxtsol(i) = 0.0         zxtsol(i) = 0.
977         zxfluxlat(i) = 0.0         zxfluxlat(i) = 0.
978    
979         zt2m(i) = 0.0         zt2m(i) = 0.
980         zq2m(i) = 0.0         zq2m(i) = 0.
981         zu10m(i) = 0.0         zu10m(i) = 0.
982         zv10m(i) = 0.0         zv10m(i) = 0.
983         zxffonte(i) = 0.0         zxffonte(i) = 0.
984         zxfqcalving(i) = 0.0         zxfqcalving(i) = 0.
985    
986         s_pblh(i) = 0.0         s_pblh(i) = 0.
987         s_lcl(i) = 0.0         s_lcl(i) = 0.
988         s_capCL(i) = 0.0         s_capCL(i) = 0.
989         s_oliqCL(i) = 0.0         s_oliqCL(i) = 0.
990         s_cteiCL(i) = 0.0         s_cteiCL(i) = 0.
991         s_pblT(i) = 0.0         s_pblT(i) = 0.
992         s_therm(i) = 0.0         s_therm(i) = 0.
993         s_trmb1(i) = 0.0         s_trmb1(i) = 0.
994         s_trmb2(i) = 0.0         s_trmb2(i) = 0.
995         s_trmb3(i) = 0.0         s_trmb3(i) = 0.
996    
997         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
998              pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.)  >  EPSFRA) &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
999              THEN              'physiq : probl\`eme sous surface au point ', i, pctsrf(i, 1 : nbsrf)
           WRITE(*, *) 'physiq : pb sous surface au point ', i, &  
                pctsrf(i, 1 : nbsrf)  
        ENDIF  
1000      ENDDO      ENDDO
1001      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1002         DO i = 1, klon         DO i = 1, klon
# Line 1122  contains Line 1053  contains
1053      ! Calculer la derive du flux infrarouge      ! Calculer la derive du flux infrarouge
1054    
1055      DO i = 1, klon      DO i = 1, klon
1056         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
1057      ENDDO      ENDDO
1058    
1059      ! Appeler la convection (au choix)      ! Appeler la convection (au choix)
1060    
1061      DO k = 1, llm      DO k = 1, llm
1062         DO i = 1, klon         DO i = 1, klon
1063            conv_q(i, k) = d_q_dyn(i, k) &            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys
1064                 + d_q_vdf(i, k)/dtphys            conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys
           conv_t(i, k) = d_t_dyn(i, k) &  
                + d_t_vdf(i, k)/dtphys  
1065         ENDDO         ENDDO
1066      ENDDO      ENDDO
1067    
1068      IF (check) THEN      IF (check) THEN
1069         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1070         print *, "avantcon = ", za         print *, "avantcon = ", za
1071      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
     IF (iflag_con == 2) zx_ajustq = .TRUE.  
     IF (zx_ajustq) THEN  
        DO i = 1, klon  
           z_avant(i) = 0.0  
        ENDDO  
        DO k = 1, llm  
           DO i = 1, klon  
              z_avant(i) = z_avant(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
     ENDIF  
1072    
1073      select case (iflag_con)      if (iflag_con == 2) then
1074      case (1)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
1075         print *, 'Réactiver l''appel à "conlmd" dans "physiq.F".'         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
1076         stop 1              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
1077      case (2)              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
1078         CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
1079              zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &              kdtop, pmflxr, pmflxs)
             pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &  
             pmflxs)  
1080         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
1081         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
1082         DO i = 1, klon         ibas_con = llm + 1 - kcbot
1083            ibas_con(i) = llm + 1 - kcbot(i)         itop_con = llm + 1 - kctop
1084            itop_con(i) = llm + 1 - kctop(i)      else
1085         ENDDO         ! iflag_con >= 3
1086      case (3:)  
1087         ! number of tracers for the convection scheme of Kerry Emanuel:         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &
1088                v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &
1089                d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &
1090                itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &
1091                pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &
1092                wd, pmflxr, pmflxs, da, phi, mp, ntra=1)
1093           ! (number of tracers for the convection scheme of Kerry Emanuel:
1094         ! la partie traceurs est faite dans phytrac         ! la partie traceurs est faite dans phytrac
1095         ! on met ntra = 1 pour limiter les appels mais on peut         ! on met ntra = 1 pour limiter les appels mais on peut
1096         ! supprimer les calculs / ftra.         ! supprimer les calculs / ftra.)
        ntra = 1  
        ! Schéma de convection modularisé et vectorisé :  
        ! (driver commun aux versions 3 et 4)  
   
        IF (ok_cvl) THEN  
           ! new driver for convectL  
           CALL concvl(iflag_con, dtphys, paprs, play, t_seri, q_seri, &  
                u_seri, v_seri, tr_seri, ntra, ema_work1, ema_work2, d_t_con, &  
                d_q_con, d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &  
                itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, pbase, &  
                bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, pmflxr, &  
                pmflxs, da, phi, mp)  
           clwcon0 = qcondc  
           pmfu = upwd + dnwd  
        ELSE  
           ! conema3 ne contient pas les traceurs  
           CALL conema3(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, &  
                tr_seri, ntra, ema_work1, ema_work2, d_t_con, d_q_con, &  
                d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &  
                itop_con, upwd, dnwd, dnwd0, bas, top, Ma, cape, tvp, rflag, &  
                pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, clwcon0)  
        ENDIF  
1097    
1098         IF (.NOT. ok_gust) THEN         clwcon0 = qcondc
1099            do i = 1, klon         mfu = upwd + dnwd
1100               wd(i) = 0.0         IF (.NOT. ok_gust) wd = 0.
           enddo  
        ENDIF  
1101    
1102         ! Calcul des propriétés des nuages convectifs         ! Calcul des propri\'et\'es des nuages convectifs
1103    
1104         DO k = 1, llm         DO k = 1, llm
1105            DO i = 1, klon            DO i = 1, klon
1106               zx_t = t_seri(i, k)               zx_t = t_seri(i, k)
1107               IF (thermcep) THEN               IF (thermcep) THEN
1108                  zdelta = MAX(0., SIGN(1., rtt-zx_t))                  zdelta = MAX(0., SIGN(1., rtt-zx_t))
1109                  zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)                  zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)
1110                  zx_qs = MIN(0.5, zx_qs)                  zx_qs = MIN(0.5, zx_qs)
1111                  zcor = 1./(1.-retv*zx_qs)                  zcor = 1./(1.-retv*zx_qs)
1112                  zx_qs = zx_qs*zcor                  zx_qs = zx_qs*zcor
# Line 1225  contains Line 1122  contains
1122         ENDDO         ENDDO
1123    
1124         ! calcul des proprietes des nuages convectifs         ! calcul des proprietes des nuages convectifs
1125         clwcon0 = fact_cldcon*clwcon0         clwcon0 = fact_cldcon * clwcon0
1126         call clouds_gno &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
1127              (klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, rnebcon0)              rnebcon0)
1128      case default  
1129         print *, "iflag_con non-prevu", iflag_con         mfd = 0.
1130         stop 1         pen_u = 0.
1131      END select         pen_d = 0.
1132           pde_d = 0.
1133           pde_u = 0.
1134        END if
1135    
1136      DO k = 1, llm      DO k = 1, llm
1137         DO i = 1, klon         DO i = 1, klon
# Line 1243  contains Line 1143  contains
1143      ENDDO      ENDDO
1144    
1145      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1146         ztit = 'after convect'         tit = 'after convect'
1147         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1148              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1149              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1150         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1151              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &
1152              fs_bound, fq_bound)              fs_bound, fq_bound)
1153      END IF      END IF
1154    
1155      IF (check) THEN      IF (check) THEN
1156         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1157         print *,"aprescon = ", za         print *, "aprescon = ", za
1158         zx_t = 0.0         zx_t = 0.
1159         za = 0.0         za = 0.
1160         DO i = 1, klon         DO i = 1, klon
1161            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1162            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
1163                 snow_con(i))*airephy(i)/REAL(klon)                 snow_con(i))*airephy(i)/REAL(klon)
1164         ENDDO         ENDDO
1165         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
1166         print *,"Precip = ", zx_t         print *, "Precip = ", zx_t
1167      ENDIF      ENDIF
1168      IF (zx_ajustq) THEN  
1169         DO i = 1, klon      IF (iflag_con == 2) THEN
1170            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1171         ENDDO         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
        DO k = 1, llm  
           DO i = 1, klon  
              z_apres(i) = z_apres(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
        DO i = 1, klon  
           z_factor(i) = (z_avant(i)-(rain_con(i) + snow_con(i))*dtphys) &  
                /z_apres(i)  
        ENDDO  
1172         DO k = 1, llm         DO k = 1, llm
1173            DO i = 1, klon            DO i = 1, klon
1174               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
# Line 1287  contains Line 1177  contains
1177            ENDDO            ENDDO
1178         ENDDO         ENDDO
1179      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
1180    
1181      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
1182    
1183      d_t_ajs = 0.      d_t_ajs = 0.
1184      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1310  contains Line 1199  contains
1199      endif      endif
1200    
1201      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1202         ztit = 'after dry_adjust'         tit = 'after dry_adjust'
1203         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1204              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1205              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1206      END IF      END IF
1207    
1208      ! Caclul des ratqs      ! Caclul des ratqs
1209    
1210      ! ratqs convectifs a l'ancienne en fonction de q(z = 0)-q / q      ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
1211      ! on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
1212      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
1213         do k = 1, llm         do k = 1, llm
1214            do i = 1, klon            do i = 1, klon
1215               if(ptconv(i, k)) then               if(ptconv(i, k)) then
1216                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
1217                       +fact_cldcon*(q_seri(i, 1)-q_seri(i, k))/q_seri(i, k)                       * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
1218               else               else
1219                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
1220               endif               endif
# Line 1336  contains Line 1225  contains
1225      ! ratqs stables      ! ratqs stables
1226      do k = 1, llm      do k = 1, llm
1227         do i = 1, klon         do i = 1, klon
1228            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1229                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1230         enddo         enddo
1231      enddo      enddo
1232    
1233      ! ratqs final      ! ratqs final
1234      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1235         ! les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
1236         ! ratqs final         ! ratqs final
1237         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1238         ! relaxation des ratqs         ! relaxation des ratqs
1239         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
1240         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
1241      else      else
1242         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
1243         ratqs = ratqss         ratqs = ratqss
1244      endif      endif
1245    
     ! Processus de condensation à grande echelle et processus de  
     ! précipitation :  
1246      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1247           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
1248           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
# Line 1375  contains Line 1261  contains
1261      ENDDO      ENDDO
1262      IF (check) THEN      IF (check) THEN
1263         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1264         print *,"apresilp = ", za         print *, "apresilp = ", za
1265         zx_t = 0.0         zx_t = 0.
1266         za = 0.0         za = 0.
1267         DO i = 1, klon         DO i = 1, klon
1268            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1269            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
1270                 + snow_lsc(i))*airephy(i)/REAL(klon)                 + snow_lsc(i))*airephy(i)/REAL(klon)
1271         ENDDO         ENDDO
1272         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
1273         print *,"Precip = ", zx_t         print *, "Precip = ", zx_t
1274      ENDIF      ENDIF
1275    
1276      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1277         ztit = 'after fisrt'         tit = 'after fisrt'
1278         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1279              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1280              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1281         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1282              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &
1283              fs_bound, fq_bound)              fs_bound, fq_bound)
1284      END IF      END IF
# Line 1401  contains Line 1287  contains
1287    
1288      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1289    
1290      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke      IF (iflag_cldcon <= -1) THEN
1291           ! seulement pour Tiedtke
1292         snow_tiedtke = 0.         snow_tiedtke = 0.
1293         if (iflag_cldcon == -1) then         if (iflag_cldcon == -1) then
1294            rain_tiedtke = rain_con            rain_tiedtke = rain_con
# Line 1418  contains Line 1305  contains
1305         endif         endif
1306    
1307         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1308         CALL diagcld1(paprs, play, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1309              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1310         DO k = 1, llm         DO k = 1, llm
1311            DO i = 1, klon            DO i = 1, klon
1312               IF (diafra(i, k) > cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
# Line 1430  contains Line 1316  contains
1316            ENDDO            ENDDO
1317         ENDDO         ENDDO
1318      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1319         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1320         ! convection et du calcul du pas de temps précédent diminué d'un facteur         ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1321         ! facttemps         ! d'un facteur facttemps.
1322         facteur = dtphys *facttemps         facteur = dtphys * facttemps
1323         do k = 1, llm         do k = 1, llm
1324            do i = 1, klon            do i = 1, klon
1325               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1326               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1327                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1328                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1329                  clwcon(i, k) = clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1330               endif               endif
# Line 1465  contains Line 1351  contains
1351      ENDIF      ENDIF
1352    
1353      ! Precipitation totale      ! Precipitation totale
   
1354      DO i = 1, klon      DO i = 1, klon
1355         rain_fall(i) = rain_con(i) + rain_lsc(i)         rain_fall(i) = rain_con(i) + rain_lsc(i)
1356         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
1357      ENDDO      ENDDO
1358    
1359      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1360         ztit = "after diagcld"           dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1361         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
1362    
1363      ! Humidité relative pour diagnostic:      ! Humidit\'e relative pour diagnostic :
1364      DO k = 1, llm      DO k = 1, llm
1365         DO i = 1, klon         DO i = 1, klon
1366            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
# Line 1501  contains Line 1383  contains
1383      ENDDO      ENDDO
1384    
1385      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
     ! Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)  
1386      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1387         ! Get sulfate aerosol distribution         ! Get sulfate aerosol distribution :
1388         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(rdayvrai, firstcal, sulfate)
1389         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1390    
        ! Calculate aerosol optical properties (Olivier Boucher)  
1391         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1392              aerindex)              aerindex)
1393      ELSE      ELSE
# Line 1516  contains Line 1396  contains
1396         cg_ae = 0.         cg_ae = 0.
1397      ENDIF      ENDIF
1398    
1399      ! Paramètres optiques des nuages et quelques paramètres pour      ! Param\`etres optiques des nuages et quelques param\`etres pour diagnostics :
     ! diagnostics :  
1400      if (ok_newmicro) then      if (ok_newmicro) then
1401         CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1402              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1403              fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, &              sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
             re, fl)  
1404      else      else
1405         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1406              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
# Line 1541  contains Line 1419  contains
1419                 + falblw(i, is_ter) * pctsrf(i, is_ter) &                 + falblw(i, is_ter) * pctsrf(i, is_ter) &
1420                 + falblw(i, is_sic) * pctsrf(i, is_sic)                 + falblw(i, is_sic) * pctsrf(i, is_sic)
1421         ENDDO         ENDDO
1422         ! nouveau rayonnement (compatible Arpege-IFS):         ! Rayonnement (compatible Arpege-IFS) :
1423         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &
1424              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &
1425              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &
# Line 1561  contains Line 1439  contains
1439      ENDDO      ENDDO
1440    
1441      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1442         ztit = 'after rad'         tit = 'after rad'
1443         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1444              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1445              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1446         call diagphy(airephy, ztit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1447              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
1448              fs_bound, fq_bound)              fs_bound, fq_bound)
1449      END IF      END IF
1450    
1451      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
1452      DO i = 1, klon      DO i = 1, klon
1453         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1454         zxsnow(i) = 0.0         zxsnow(i) = 0.
1455      ENDDO      ENDDO
1456      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1457         DO i = 1, klon         DO i = 1, klon
# Line 1582  contains Line 1460  contains
1460         ENDDO         ENDDO
1461      ENDDO      ENDDO
1462    
1463      ! Calculer le bilan du sol et la dérive de température (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1464    
1465      DO i = 1, klon      DO i = 1, klon
1466         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1467      ENDDO      ENDDO
1468    
1469      ! Paramétrisation de l'orographie à l'échelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1470    
1471      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1472         ! selection des points pour lesquels le shema est actif:         ! selection des points pour lesquels le shema est actif:
1473         igwd = 0         igwd = 0
1474         DO i = 1, klon         DO i = 1, klon
1475            itest(i) = 0            itest(i) = 0
1476            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.0)) THEN            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN
1477               itest(i) = 1               itest(i) = 1
1478               igwd = igwd + 1               igwd = igwd + 1
1479               idx(igwd) = i               idx(igwd) = i
# Line 1617  contains Line 1495  contains
1495      ENDIF      ENDIF
1496    
1497      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1498         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
1499         igwd = 0         igwd = 0
1500         DO i = 1, klon         DO i = 1, klon
1501            itest(i) = 0            itest(i) = 0
# Line 1642  contains Line 1520  contains
1520         ENDDO         ENDDO
1521      ENDIF      ENDIF
1522    
1523      ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE      ! Stress n\'ecessaires : toute la physique
1524    
1525      DO i = 1, klon      DO i = 1, klon
1526         zustrph(i) = 0.         zustrph(i) = 0.
# Line 1650  contains Line 1528  contains
1528      ENDDO      ENDDO
1529      DO k = 1, llm      DO k = 1, llm
1530         DO i = 1, klon         DO i = 1, klon
1531            zustrph(i) = zustrph(i) + (u_seri(i, k)-u(i, k))/dtphys* zmasse(i, k)            zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1532            zvstrph(i) = zvstrph(i) + (v_seri(i, k)-v(i, k))/dtphys* zmasse(i, k)                 * zmasse(i, k)
1533              zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1534                   * zmasse(i, k)
1535         ENDDO         ENDDO
1536      ENDDO      ENDDO
1537    
1538      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &
1539           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1540    
1541      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1542         ztit = 'after orography'           2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1543         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
1544    
1545      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1546      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, &      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &
1547           nqmx-2, dtphys, u, t, paprs, play, pmfu, pmfd, pen_u, pde_u, &           dtphys, u, t, paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, &
1548           pen_d, pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &           entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, &
1549           frac_impa, frac_nucl, pphis, albsol, rhcl, cldfra, rneb, &           albsol, rhcl, cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, &
1550           diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, &           psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)
1551           tr_seri, zmasse)  
1552        IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1553      IF (offline) THEN           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1554         call phystokenc(dtphys, rlon, rlat, t, pmfu, pmfd, pen_u, pde_u, &           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
             pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &  
             pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)  
     ENDIF  
1555    
1556      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1557      CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &      CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &
# Line 1702  contains Line 1576  contains
1576      END DO      END DO
1577    
1578      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1579         ztit = 'after physic'         tit = 'after physic'
1580         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1581              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1582              d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1583         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1584         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1585         ! est egale a la variation de la physique au pas de temps precedent.         ! est egale a la variation de la physique au pas de temps precedent.
1586         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1587         call diagphy(airephy, ztit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1588              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &
1589              fs_bound, fq_bound)              fs_bound, fq_bound)
1590    
# Line 1720  contains Line 1594  contains
1594    
1595      ! SORTIES      ! SORTIES
1596    
1597      !cc prw = eau precipitable      ! prw = eau precipitable
1598      DO i = 1, klon      DO i = 1, klon
1599         prw(i) = 0.         prw(i) = 0.
1600         DO k = 1, llm         DO k = 1, llm
# Line 1759  contains Line 1633  contains
1633      ENDDO      ENDDO
1634    
1635      ! Ecriture des sorties      ! Ecriture des sorties
     call write_histhf  
     call write_histday  
1636      call write_histins      call write_histins
1637    
1638      ! Si c'est la fin, il faut conserver l'etat de redemarrage      ! Si c'est la fin, il faut conserver l'etat de redemarrage
# Line 1768  contains Line 1640  contains
1640         itau_phy = itau_phy + itap         itau_phy = itau_phy + itap
1641         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &
1642              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &
1643              rain_fall, snow_fall, solsw, sollwdown, dlw, radsol, frugs, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
1644              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
1645              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
1646      ENDIF      ENDIF
1647    
1648      firstcal = .FALSE.      firstcal = .FALSE.
1649    
1650    contains    contains
1651    
     subroutine write_histday  
   
       use gr_phy_write_3d_m, only: gr_phy_write_3d  
       integer itau_w ! pas de temps ecriture  
   
       !------------------------------------------------  
   
       if (ok_journe) THEN  
          itau_w = itau_phy + itap  
          if (nqmx <= 4) then  
             call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &  
                  gr_phy_write_3d(wo) * 1e3)  
             ! (convert "wo" from kDU to DU)  
          end if  
          if (ok_sync) then  
             call histsync(nid_day)  
          endif  
       ENDIF  
   
     End subroutine write_histday  
   
     !****************************  
   
     subroutine write_histhf  
   
       ! From phylmd/write_histhf.h, version 1.5 2005/05/25 13:10:09  
   
       !------------------------------------------------  
   
       call write_histhf3d  
   
       IF (ok_sync) THEN  
          call histsync(nid_hf)  
       ENDIF  
   
     end subroutine write_histhf  
   
     !***************************************************************  
   
1652      subroutine write_histins      subroutine write_histins
1653    
1654        ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09        ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09
1655    
1656          use dimens_m, only: iim, jjm
1657          USE histsync_m, ONLY: histsync
1658          USE histwrite_m, ONLY: histwrite
1659    
1660        real zout        real zout
1661        integer itau_w ! pas de temps ecriture        integer itau_w ! pas de temps ecriture
1662          REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1663    
1664        !--------------------------------------------------        !--------------------------------------------------
1665    
# Line 2044  contains Line 1882  contains
1882    
1883      end subroutine write_histins      end subroutine write_histins
1884    
     !****************************************************  
   
     subroutine write_histhf3d  
   
       ! From phylmd/write_histhf3d.h, version 1.2 2005/05/25 13:10:09  
   
       integer itau_w ! pas de temps ecriture  
   
       !-------------------------------------------------------  
   
       itau_w = itau_phy + itap  
   
       ! Champs 3D:  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, qx(1, 1, ivap), zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)  
   
       if (nbtr >= 3) then  
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, tr_seri(1, 1, 3), &  
               zx_tmp_3d)  
          CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)  
       end if  
   
       if (ok_sync) then  
          call histsync(nid_hf3d)  
       endif  
   
     end subroutine write_histhf3d  
   
1885    END SUBROUTINE physiq    END SUBROUTINE physiq
1886    
1887  end module physiq_m  end module physiq_m
Legend:
Removed from v.61  
changed lines
  Added in v.90
  ViewVC Help
Powered by ViewVC 1.1.21