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

Diff of /trunk/phylmd/physiq.f

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

trunk/libf/phylmd/physiq.f90 revision 62 by guez, Thu Jul 26 14:37:37 2012 UTC trunk/phylmd/physiq.f revision 103 by guez, Fri Aug 29 13:00:05 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)
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
     USE calendar, ONLY: ymds2ju  
21      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
22      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &
23           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin
24      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
25           ok_orodr, ok_orolf, soil_model           ok_orodr, ok_orolf
26      USE clmain_m, ONLY: clmain      USE clmain_m, ONLY: clmain
27        use clouds_gno_m, only: clouds_gno
28      USE comgeomphy, ONLY: airephy, cuphy, cvphy      USE comgeomphy, ONLY: airephy, cuphy, cvphy
29      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
30      USE conf_gcm_m, ONLY: offline, raz_date      USE conf_gcm_m, ONLY: offline, raz_date
# Line 31  contains Line 34  contains
34      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
35      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
36      use diagphy_m, only: diagphy      use diagphy_m, only: diagphy
37      USE dimens_m, ONLY: iim, jjm, llm, nqmx      USE dimens_m, ONLY: llm, nqmx
38      USE dimphy, ONLY: klon, nbtr      USE dimphy, ONLY: klon
39      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
40      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
41      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
42        use fisrtilp_m, only: fisrtilp
43      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
     USE histsync_m, ONLY: histsync  
     USE histwrite_m, ONLY: histwrite  
44      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, &
45           nbsrf           nbsrf
     USE ini_histhf_m, ONLY: ini_histhf  
     USE ini_histday_m, ONLY: ini_histday  
46      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins
47      USE oasis_m, ONLY: ok_oasis      use newmicro_m, only: newmicro
48      USE orbite_m, ONLY: orbite, zenang      USE orbite_m, ONLY: orbite
49      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
50      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
51      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
# Line 53  contains Line 53  contains
53      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
54      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
55      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
56        use readsulfate_m, only: readsulfate
57      use sugwd_m, only: sugwd      use sugwd_m, only: sugwd
58      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
59      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE temps, ONLY: annee_ref, day_ref, itau_phy
60        use unit_nml_m, only: unit_nml
61        USE ymds2ju_m, ONLY: ymds2ju
62      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
63        use zenang_m, only: zenang
64    
65      ! Arguments:      logical, intent(in):: lafin ! dernier passage
66    
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)
     logical, intent(in):: lafin ! dernier passage  
72    
73      REAL, intent(in):: paprs(klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
74      ! (pression pour chaque inter-couche, en Pa)      ! pression pour chaque inter-couche, en Pa
75    
76      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
77      ! (input pression pour le mileu de chaque couche (en Pa))      ! pression pour le mileu de chaque couche (en Pa)
78    
79      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
80      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! géopotentiel de chaque couche (référence sol)
81    
82      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
83    
84      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
85      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m/s
86    
87      REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s
88      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
89    
90      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
91      ! (humidité spécifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
92    
93      REAL omega(klon, llm) ! input vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
94      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-2)
95      REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/s)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
96      REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)
     REAL d_qx(klon, llm, nqmx) ! output tendance physique de "qx" (kg/kg/s)  
     REAL d_ps(klon) ! output tendance physique de la pression au sol  
97    
98      LOGICAL:: firstcal = .true.      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
99        ! tendance physique de "qx" (s-1)
100    
101      INTEGER nbteta      ! Local:
     PARAMETER(nbteta = 3)  
102    
103      REAL PVteta(klon, nbteta)      LOGICAL:: firstcal = .true.
     ! (output vorticite potentielle a des thetas constantes)  
104    
105      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
106      PARAMETER (ok_gust = .FALSE.)      PARAMETER (ok_gust = .FALSE.)
107    
108      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL, PARAMETER:: check = .FALSE.
109      PARAMETER (check = .FALSE.)      ! Verifier la conservation du modele en eau
110    
111      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
112      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
113    
     ! Parametres lies au coupleur OASIS:  
     INTEGER, SAVE:: npas, nexca  
     logical rnpb  
     parameter(rnpb = .true.)  
   
     character(len = 6), save:: ocean  
     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")  
   
     logical ok_ocean  
     SAVE ok_ocean  
   
114      ! "slab" ocean      ! "slab" ocean
115      REAL, save:: tslab(klon) ! temperature of ocean slab      REAL, save:: tslab(klon) ! temperature of ocean slab
116      REAL, save:: seaice(klon) ! glace de mer (kg/m2)      REAL, save:: seaice(klon) ! glace de mer (kg/m2)
117      REAL fluxo(klon) ! flux turbulents ocean-glace de mer      REAL fluxo(klon) ! flux turbulents ocean-glace de mer
118      REAL fluxg(klon) ! flux turbulents ocean-atmosphere      REAL fluxg(klon) ! flux turbulents ocean-atmosphere
119    
120      ! Modele thermique du sol, a activer pour le cycle diurne:      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
121      logical, save:: ok_veget      ! sorties journalieres, mensuelles et instantanees dans les
122      LOGICAL, save:: ok_journe ! sortir le fichier journalier      ! fichiers histday, histmth et histins
   
     LOGICAL ok_mensuel ! sortir le fichier mensuel  
   
     LOGICAL ok_instan ! sortir le fichier instantane  
     save ok_instan  
123    
124      LOGICAL ok_region ! sortir le fichier regional      LOGICAL ok_region ! sortir le fichier regional
125      PARAMETER (ok_region = .FALSE.)      PARAMETER (ok_region = .FALSE.)
# Line 145  contains Line 129  contains
129      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
130      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
131    
132      INTEGER ivap ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
133      PARAMETER (ivap = 1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     INTEGER iliq ! indice de traceurs pour eau liquide  
     PARAMETER (iliq = 2)  
134    
135      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
136      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
# Line 158  contains Line 140  contains
140    
141      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
142    
143      !IM Amip2 PV a theta constante      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
144        REAL swup0(klon, llm + 1), swup(klon, llm + 1)
     CHARACTER(LEN = 3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
   
     !MI Amip2 PV a theta constante  
   
     INTEGER klevp1  
     PARAMETER(klevp1 = llm + 1)  
   
     REAL swdn0(klon, klevp1), swdn(klon, klevp1)  
     REAL swup0(klon, klevp1), swup(klon, klevp1)  
145      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
146    
147      REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
148      REAL lwup0(klon, klevp1), lwup(klon, klevp1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
149      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
150    
151      !IM Amip2      ! Amip2
152      ! variables a une pression donnee      ! variables a une pression donnee
153    
154      integer nlevSTD      integer nlevSTD
# Line 206  contains Line 176  contains
176      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)
177    
178      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
179      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./
180      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
181    
182      ! cldtopres pression au sommet des nuages      ! cldtopres pression au sommet des nuages
# Line 247  contains Line 217  contains
217           'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &           'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &
218           'pc= 680-800hPa, tau> 60.'/           'pc= 680-800hPa, tau> 60.'/
219    
220      !IM ISCCP simulator v3.4      ! ISCCP simulator v3.4
   
     integer nid_hf, nid_hf3d  
     save nid_hf, nid_hf3d  
221    
222      ! Variables propres a la physique      ! Variables propres a la physique
223    
# Line 268  contains Line 235  contains
235      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
236      ! soil temperature of surface fraction      ! soil temperature of surface fraction
237    
238      REAL fevap(klon, nbsrf)      REAL, save:: fevap(klon, nbsrf) ! evaporation
     SAVE fevap ! evaporation  
239      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
240      SAVE fluxlat      SAVE fluxlat
241    
242      REAL fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
243      SAVE fqsurf ! humidite de l'air au contact de la surface      ! humidite de l'air au contact de la surface
   
     REAL, save:: qsol(klon) ! hauteur d'eau dans le sol  
244    
245      REAL fsnow(klon, nbsrf)      REAL, save:: qsol(klon)
246      SAVE fsnow ! epaisseur neigeuse      ! column-density of water in soil, in kg m-2
247    
248      REAL falbe(klon, nbsrf)      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
249      SAVE falbe ! albedo par type de surface      REAL, save:: falbe(klon, nbsrf) ! albedo par type de surface
250      REAL falblw(klon, nbsrf)      REAL, save:: falblw(klon, nbsrf) ! albedo par type de surface
     SAVE falblw ! albedo par type de surface  
251    
252      ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
253      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
254      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
255      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 308  contains Line 271  contains
271      !KE43      !KE43
272      ! Variables liees a la convection de K. Emanuel (sb):      ! Variables liees a la convection de K. Emanuel (sb):
273    
     REAL bas, top ! cloud base and top levels  
     SAVE bas  
     SAVE top  
   
274      REAL Ma(klon, llm) ! undilute upward mass flux      REAL Ma(klon, llm) ! undilute upward mass flux
275      SAVE Ma      SAVE Ma
276      REAL qcondc(klon, llm) ! in-cld water content from convect      REAL qcondc(klon, llm) ! in-cld water content from convect
277      SAVE qcondc      SAVE qcondc
278      REAL ema_work1(klon, llm), ema_work2(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
279      SAVE ema_work1, ema_work2      REAL, save:: wd(klon)
   
     REAL wd(klon) ! sb  
     SAVE wd ! sb  
280    
281      ! Variables locales pour la couche limite (al1):      ! Variables locales pour la couche limite (al1):
282    
# Line 329  contains Line 285  contains
285      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
286      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
287    
288      !AA Pour phytrac      ! Pour phytrac :
289      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
290      REAL yu1(klon) ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
291      REAL yv1(klon) ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
# Line 348  contains Line 304  contains
304      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
305      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
306    
307      REAL, save:: rain_fall(klon) ! pluie      REAL, save:: rain_fall(klon)
308      REAL, save:: snow_fall(klon) ! neige      ! liquid water mass flux (kg/m2/s), positive down
309    
310        REAL, save:: snow_fall(klon)
311        ! solid water mass flux (kg/m2/s), positive down
312    
313      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
314    
315      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon), devap(klon) ! evaporation and its derivative
316      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
317      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
318      SAVE dlw      SAVE dlw
# Line 365  contains Line 324  contains
324      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
325      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
326    
327      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
328      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
329    
330      ! Conditions aux limites      ! Conditions aux limites
331    
332      INTEGER julien      INTEGER julien
   
333      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
334      REAL pctsrf(klon, nbsrf)      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
335      !IM      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
336      REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE      REAL, save:: albsol(klon) ! albedo du sol total
337        REAL, save:: albsollw(klon) ! albedo du sol total
     SAVE pctsrf ! sous-fraction du sol  
     REAL albsol(klon)  
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
338      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
339    
340      ! Declaration des procedures appelees      ! Declaration des procedures appelees
341    
     EXTERNAL alboc ! calculer l'albedo sur ocean  
     !KE43  
     EXTERNAL conema3 ! convect4.3  
     EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)  
342      EXTERNAL nuage ! calculer les proprietes radiatives      EXTERNAL nuage ! calculer les proprietes radiatives
343      EXTERNAL transp ! transport total de l'eau et de l'energie      EXTERNAL transp ! transport total de l'eau et de l'energie
344    
345      ! Variables locales      ! Variables locales
346    
347      real clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
348      real clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
   
     save rnebcon, clwcon  
349    
350      REAL rhcl(klon, llm) ! humiditi relative ciel clair      REAL rhcl(klon, llm) ! humiditi relative ciel clair
351      REAL dialiq(klon, llm) ! eau liquide nuageuse      REAL dialiq(klon, llm) ! eau liquide nuageuse
# Line 420  contains Line 365  contains
365      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
366      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
367    
368      ! Le rayonnement n'est pas calculé tous les pas, il faut donc que      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
369      ! les variables soient rémanentes.      ! les variables soient r\'emanentes.
370      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
371      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL heat0(klon, llm) ! chauffage solaire ciel clair
372      REAL, save:: cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
373      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
374      REAL, save:: topsw(klon), toplw(klon), solsw(klon), sollw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
375      real sollwdown(klon) ! downward LW flux at surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
376        real, save:: sollwdown(klon) ! downward LW flux at surface
377      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
378      REAL albpla(klon)      REAL albpla(klon)
379      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
380      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
381      SAVE albpla, sollwdown      SAVE albpla
382      SAVE heat0, cool0      SAVE heat0, cool0
383    
384      INTEGER itaprad      INTEGER itaprad
# Line 447  contains Line 393  contains
393      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
394    
395      REAL dist, rmu0(klon), fract(klon)      REAL dist, rmu0(klon), fract(klon)
     REAL zdtime ! pas de temps du rayonnement (s)  
396      real zlongi      real zlongi
   
397      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
398      REAL za, zb      REAL za, zb
399      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
400      real zqsat(klon, llm)      real zqsat(klon, llm)
401      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
402      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup = 234.0)  
   
403      REAL zphi(klon, llm)      REAL zphi(klon, llm)
404    
405      !IM cf. AM Variables locales pour la CLA (hbtm2)      ! cf. AM Variables locales pour la CLA (hbtm2)
406    
407      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
408      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
# Line 485  contains Line 425  contains
425      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
426      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
427      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
     REAL tvp(klon, llm) ! virtual temp of lifted parcel  
428      REAL cape(klon) ! CAPE      REAL cape(klon) ! CAPE
429      SAVE cape      SAVE cape
430    
     REAL pbase(klon) ! cloud base pressure  
     SAVE pbase  
     REAL bbase(klon) ! cloud base buoyancy  
     SAVE bbase  
     REAL rflag(klon) ! flag fonctionnement de convect  
431      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
     ! -- convect43:  
     INTEGER ntra ! nb traceurs pour convect4.3  
     REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)  
     REAL dplcldt(klon), dplcldr(klon)  
432    
433      ! Variables du changement      ! Variables du changement
434    
435      ! con: convection      ! con: convection
436      ! lsc: large scale condensation      ! lsc: large scale condensation
437      ! ajs: ajustement sec      ! ajs: ajustement sec
438      ! eva: évaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
439      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
440      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
441      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
# Line 514  contains Line 444  contains
444      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
445      REAL rneb(klon, llm)      REAL rneb(klon, llm)
446    
447      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
448      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
449      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
450      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
# Line 535  contains Line 465  contains
465      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
466      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
467    
468      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
469      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
470      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
471    
472      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
473      real, save:: fact_cldcon      real:: fact_cldcon = 0.375
474      real, save:: facttemps      real:: facttemps = 1.e-4
475      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
476      real facteur      real facteur
477    
478      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
479      logical ptconv(klon, llm)      logical ptconv(klon, llm)
480    
481      ! Variables locales pour effectuer les appels en série :      ! Variables locales pour effectuer les appels en s\'erie :
482    
483      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
484      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
485      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
486        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
487    
488      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
489    
# Line 567  contains Line 492  contains
492      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
493      REAL aam, torsfc      REAL aam, torsfc
494    
     REAL dudyn(iim + 1, jjm + 1, llm)  
   
495      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)  
496    
497      INTEGER, SAVE:: nid_day, nid_ins      INTEGER, SAVE:: nid_ins
498    
499      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
500      REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.      REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
# Line 580  contains Line 502  contains
502      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
503    
504      REAL zsto      REAL zsto
   
     character(len = 20) modname  
     character(len = 80) abort_message  
     logical ok_sync  
505      real date0      real date0
506    
507      ! Variables liées au bilan d'énergie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
508      REAL ztsol(klon)      REAL ztsol(klon)
509      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
510      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
511      REAL zero_v(klon)      REAL zero_v(klon)
512      CHARACTER(LEN = 15) tit      CHARACTER(LEN = 20) tit
513      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
514      INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
515    
516      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
517      REAL ZRCPD      REAL ZRCPD
518    
519      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
520      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
521      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
522      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
523      !jq Aerosol effects (Johannes Quaas, 27/11/2003)  
524      REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]      ! Aerosol effects:
525    
526        REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
527    
528      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
529      ! (SO4 aerosol concentration, in ug/m3, pre-industrial value)      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
530    
531      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
532      ! (Cloud optical thickness for pre-industrial (pi) aerosols)      ! cloud optical thickness for pre-industrial (pi) aerosols
533    
534      REAL re(klon, llm) ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
535      REAL fl(klon, llm) ! denominator of re      REAL fl(klon, llm) ! denominator of re
536    
537      ! Aerosol optical properties      ! Aerosol optical properties
538      REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)      REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
539      REAL cg_ae(klon, llm, 2)      REAL, save:: cg_ae(klon, llm, 2)
   
     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.  
     ! ok_ade = True -ADE = topswad-topsw  
540    
541        REAL topswad(klon), solswad(klon) ! aerosol direct effect
542      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  
543    
544      REAL aerindex(klon) ! POLDER aerosol index      REAL aerindex(klon) ! POLDER aerosol index
545    
546      ! Parameters      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
547      LOGICAL, save:: ok_ade ! apply aerosol direct effect      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
548      LOGICAL, save:: ok_aie ! Apply aerosol indirect effect  
549      REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)      REAL:: bl95_b0 = 2., bl95_b1 = 0.2
550        ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
551        ! B). They link cloud droplet number concentration to aerosol mass
552        ! concentration.
553    
     SAVE bl95_b0, bl95_b1  
554      SAVE u10m      SAVE u10m
555      SAVE v10m      SAVE v10m
556      SAVE t2m      SAVE t2m
557      SAVE q2m      SAVE q2m
558      SAVE ffonte      SAVE ffonte
559      SAVE fqcalving      SAVE fqcalving
     SAVE piz_ae  
     SAVE tau_ae  
     SAVE cg_ae  
560      SAVE rain_con      SAVE rain_con
561      SAVE snow_con      SAVE snow_con
562      SAVE topswai      SAVE topswai
# Line 652  contains Line 565  contains
565      SAVE solswad      SAVE solswad
566      SAVE d_u_con      SAVE d_u_con
567      SAVE d_v_con      SAVE d_v_con
     SAVE rnebcon0  
     SAVE clwcon0  
568    
569      real zmasse(klon, llm)      real zmasse(klon, llm)
570      ! (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)
571    
572      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
573    
574        namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
575             facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
576             ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals
577    
578      !----------------------------------------------------------------      !----------------------------------------------------------------
579    
580      modname = 'physiq'      IF (if_ebil >= 1) zero_v = 0.
581      IF (if_ebil >= 1) THEN      IF (nqmx < 2) CALL abort_gcm('physiq', &
582         DO i = 1, klon           'eaux vapeur et liquide sont indispensables', 1)
           zero_v(i) = 0.  
        END DO  
     END IF  
     ok_sync = .TRUE.  
     IF (nqmx < 2) THEN  
        abort_message = 'eaux vapeur et liquide sont indispensables'  
        CALL abort_gcm(modname, abort_message, 1)  
     ENDIF  
583    
584      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
585         ! initialiser         ! initialiser
# Line 685  contains Line 592  contains
592         piz_ae = 0.         piz_ae = 0.
593         tau_ae = 0.         tau_ae = 0.
594         cg_ae = 0.         cg_ae = 0.
595         rain_con(:) = 0.         rain_con = 0.
596         snow_con(:) = 0.         snow_con = 0.
597         bl95_b0 = 0.         topswai = 0.
598         bl95_b1 = 0.         topswad = 0.
599         topswai(:) = 0.         solswai = 0.
600         topswad(:) = 0.         solswad = 0.
601         solswai(:) = 0.  
602         solswad(:) = 0.         d_u_con = 0.
603           d_v_con = 0.
604         d_u_con = 0.0         rnebcon0 = 0.
605         d_v_con = 0.0         clwcon0 = 0.
606         rnebcon0 = 0.0         rnebcon = 0.
607         clwcon0 = 0.0         clwcon = 0.
        rnebcon = 0.0  
        clwcon = 0.0  
608    
609         pblh =0. ! Hauteur de couche limite         pblh =0. ! Hauteur de couche limite
610         plcl =0. ! Niveau de condensation de la CLA         plcl =0. ! Niveau de condensation de la CLA
# Line 714  contains Line 619  contains
619    
620         IF (if_ebil >= 1) d_h_vcol_phy = 0.         IF (if_ebil >= 1) d_h_vcol_phy = 0.
621    
622         ! Appel à la lecture du run.def physique         iflag_thermals = 0
623         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &         nsplit_thermals = 1
624              fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &         print *, "Enter namelist 'physiq_nml'."
625              ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, &         read(unit=*, nml=physiq_nml)
626              iflag_thermals, nsplit_thermals)         write(unit_nml, nml=physiq_nml)
627    
628           call conf_phys
629    
630         ! Initialiser les compteurs:         ! Initialiser les compteurs:
631    
632         frugs = 0.         frugs = 0.
633         itap = 0         itap = 0
634         itaprad = 0         itaprad = 0
635         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &         CALL phyetat0(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &
636              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              fsnow, falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollw, &
637              snow_fall, solsw, sollwdown, dlw, radsol, frugs, agesno, zmea, &              dlw, radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, &
638              zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &              zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
639              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)              run_off_lic_0, sig1, w01)
640    
641         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
642         q2 = 1.e-8         q2 = 1e-8
643    
644         radpas = NINT(86400. / dtphys / nbapp_rad)         radpas = NINT(86400. / dtphys / nbapp_rad)
645    
# Line 740  contains Line 647  contains
647         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
648    
649         PRINT *, 'cycle_diurne = ', cycle_diurne         PRINT *, 'cycle_diurne = ', cycle_diurne
650           CALL printflag(radpas, ok_journe, ok_instan, ok_region)
651    
652         IF(ocean.NE.'force ') THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
           ok_ocean = .TRUE.  
        ENDIF  
   
        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'  
653            print *, "Au minimum 4 appels par jour si cycle diurne"            print *, "Au minimum 4 appels par jour si cycle diurne"
654            abort_message = 'Nbre d appels au rayonnement insuffisant'            call abort_gcm('physiq', &
655            call abort_gcm(modname, abort_message, 1)                 "Nombre d'appels au rayonnement insuffisant", 1)
656         ENDIF         ENDIF
        print *, "Clef pour la convection, iflag_con = ", iflag_con  
657    
658         ! Initialisation pour la convection de K.E. (sb):         ! Initialisation pour le sch\'ema de convection d'Emanuel :
659         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
660            print *, "Convection de Kerry Emanuel 4.3"            ibas_con = 1
661              itop_con = 1
           DO i = 1, klon  
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
662         ENDIF         ENDIF
663    
664         IF (ok_orodr) THEN         IF (ok_orodr) THEN
# Line 782  contains Line 677  contains
677         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
678         ecrit_reg = NINT(ecrit_reg/dtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
679    
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
680         ! Initialisation des sorties         ! Initialisation des sorties
681    
        call ini_histhf(dtphys, nid_hf, nid_hf3d)  
        call ini_histday(dtphys, ok_journe, nid_day, nqmx)  
682         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
683         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
684         !XXXPB Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
685         WRITE(*, *) 'physiq date0: ', date0         print *, 'physiq date0: ', date0
686      ENDIF test_firstcal      ENDIF test_firstcal
687    
688      ! Mettre a zero des variables de sortie (pour securite)      ! We will modify variables *_seri and we will not touch variables
689        ! u, v, t, qx:
690      DO i = 1, klon      t_seri = t
691         d_ps(i) = 0.0      u_seri = u
692      ENDDO      v_seri = v
693      DO iq = 1, nqmx      q_seri = qx(:, :, ivap)
694         DO k = 1, llm      ql_seri = qx(:, :, iliq)
695            DO i = 1, klon      tr_seri = qx(:, :, 3: nqmx)
              d_qx(i, k, iq) = 0.0  
           ENDDO  
        ENDDO  
     ENDDO  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
     ! Ne pas affecter les valeurs entrées de u, v, h, et q :  
696    
697      DO k = 1, llm      ztsol = sum(ftsol * pctsrf, dim = 2)
        DO i = 1, klon  
           t_seri(i, k) = t(i, k)  
           u_seri(i, k) = u(i, k)  
           v_seri(i, k) = v(i, k)  
           q_seri(i, k) = qx(i, k, ivap)  
           ql_seri(i, k) = qx(i, k, iliq)  
           qs_seri(i, k) = 0.  
        ENDDO  
     ENDDO  
     IF (nqmx >= 3) THEN  
        tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)  
     ELSE  
        tr_seri(:, :, 1) = 0.  
     ENDIF  
   
     DO i = 1, klon  
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
698    
699      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
700         tit = 'after dynamics'         tit = 'after dynamics'
701         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
702              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
703              d_ql, d_qs, d_ec)         ! Comme les tendances de la physique sont ajout\'es dans la
        ! Comme les tendances de la physique sont ajoutés dans la  
704         !  dynamique, la variation d'enthalpie par la dynamique devrait         !  dynamique, la variation d'enthalpie par la dynamique devrait
705         !  être égale à la variation de la physique au pas de temps         !  \^etre \'egale \`a la variation de la physique au pas de temps
706         !  précédent.  Donc la somme de ces 2 variations devrait être         !  pr\'ec\'edent.  Donc la somme de ces 2 variations devrait \^etre
707         !  nulle.         !  nulle.
708         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
709              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, &
710              d_qt, 0., fs_bound, fq_bound)              d_qt, 0.)
711      END IF      END IF
712    
713      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
# Line 866  contains Line 721  contains
721      ELSE      ELSE
722         DO k = 1, llm         DO k = 1, llm
723            DO i = 1, klon            DO i = 1, klon
724               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
725               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
726            ENDDO            ENDDO
727         ENDDO         ENDDO
728         ancien_ok = .TRUE.         ancien_ok = .TRUE.
# Line 883  contains Line 738  contains
738      ! Check temperatures:      ! Check temperatures:
739      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
740    
741      ! Incrementer le compteur de la physique      ! Incrémenter le compteur de la physique
742      itap = itap + 1      itap = itap + 1
743      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(NINT(rdayvrai), 360)
744      if (julien == 0) julien = 360      if (julien == 0) julien = 360
745    
746      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
   
     ! Mettre en action les conditions aux limites (albedo, sst, etc.).  
747    
748      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone :
749      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
750    
751      ! Évaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
752      DO k = 1, llm      DO k = 1, llm
753         DO i = 1, klon         DO i = 1, klon
754            zb = MAX(0., ql_seri(i, k))            zb = MAX(0., ql_seri(i, k))
# Line 909  contains Line 762  contains
762      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
763         tit = 'after reevap'         tit = 'after reevap'
764         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
765              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
766         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
767              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)
             fs_bound, fq_bound)  
   
768      END IF      END IF
769    
770      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
771        zxrugs = sum(frugs * pctsrf, dim = 2)
772    
773      DO i = 1, klon      ! Calculs nécessaires au calcul de l'albedo dans l'interface
        zxrugs(i) = 0.0  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)  
        ENDDO  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
   
     ! calculs necessaires au calcul de l'albedo dans l'interface  
774    
775      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), zlongi, dist)
776      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
777         zdtime = dtphys * REAL(radpas)         CALL zenang(zlongi, time, dtphys * REAL(radpas), rmu0, fract)
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
778      ELSE      ELSE
779         rmu0 = -999.999         rmu0 = -999.999
780      ENDIF      ENDIF
781    
782      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
783      albsol(:) = 0.      albsol = sum(falbe * pctsrf, dim = 2)
784      albsollw(:) = 0.      albsollw = sum(falblw * pctsrf, dim = 2)
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)  
           albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
785    
786      ! Repartition sous maille des flux LW et SW      ! R\'epartition sous maille des flux longwave et shortwave
787      ! Repartition du longwave par sous-surface linearisee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
788    
789      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
790         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
791            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
792                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
793            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))      END forall
        ENDDO  
     ENDDO  
794    
795      fder = dlw      fder = dlw
796    
797      ! Couche limite:      ! Couche limite:
798    
799      CALL clmain(dtphys, itap, date0, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
800           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, &           v_seri, julien, rmu0, co2_ppm, ftsol, cdmmax, cdhmax, &
801           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
802           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &
803           rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder, rlon, rlat, &           fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &
804           cuphy, cvphy, frugs, firstcal, lafin, agesno, rugoro, d_t_vdf, &           d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &
805           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, &
806           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
807           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab)
          fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)  
808    
809      ! Incrémentation des flux      ! Incr\'ementation des flux
810    
811      zxfluxt = 0.      zxfluxt = 0.
812      zxfluxq = 0.      zxfluxq = 0.
# Line 988  contains Line 815  contains
815      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
816         DO k = 1, llm         DO k = 1, llm
817            DO i = 1, klon            DO i = 1, klon
818               zxfluxt(i, k) = zxfluxt(i, k) + &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
819                    fluxt(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
820               zxfluxq(i, k) = zxfluxq(i, k) + &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
821                    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)  
822            END DO            END DO
823         END DO         END DO
824      END DO      END DO
825      DO i = 1, klon      DO i = 1, klon
826         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
827         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
828         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
829      ENDDO      ENDDO
830    
# Line 1017  contains Line 840  contains
840      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
841         tit = 'after clmain'         tit = 'after clmain'
842         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
843              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
844         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
845              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)
             fs_bound, fq_bound)  
846      END IF      END IF
847    
848      ! Update surface temperature:      ! Update surface temperature:
849    
850      DO i = 1, klon      DO i = 1, klon
851         zxtsol(i) = 0.0         zxtsol(i) = 0.
852         zxfluxlat(i) = 0.0         zxfluxlat(i) = 0.
853    
854         zt2m(i) = 0.0         zt2m(i) = 0.
855         zq2m(i) = 0.0         zq2m(i) = 0.
856         zu10m(i) = 0.0         zu10m(i) = 0.
857         zv10m(i) = 0.0         zv10m(i) = 0.
858         zxffonte(i) = 0.0         zxffonte(i) = 0.
859         zxfqcalving(i) = 0.0         zxfqcalving(i) = 0.
860    
861         s_pblh(i) = 0.0         s_pblh(i) = 0.
862         s_lcl(i) = 0.0         s_lcl(i) = 0.
863         s_capCL(i) = 0.0         s_capCL(i) = 0.
864         s_oliqCL(i) = 0.0         s_oliqCL(i) = 0.
865         s_cteiCL(i) = 0.0         s_cteiCL(i) = 0.
866         s_pblT(i) = 0.0         s_pblT(i) = 0.
867         s_therm(i) = 0.0         s_therm(i) = 0.
868         s_trmb1(i) = 0.0         s_trmb1(i) = 0.
869         s_trmb2(i) = 0.0         s_trmb2(i) = 0.
870         s_trmb3(i) = 0.0         s_trmb3(i) = 0.
871    
872         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
873              pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.)  >  EPSFRA) &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
874              THEN              'physiq : probl\`eme sous surface au point ', i, &
875            WRITE(*, *) 'physiq : pb sous surface au point ', i, &              pctsrf(i, 1 : nbsrf)
                pctsrf(i, 1 : nbsrf)  
        ENDIF  
876      ENDDO      ENDDO
877      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
878         DO i = 1, klon         DO i = 1, klon
# Line 1081  contains Line 900  contains
900         ENDDO         ENDDO
901      ENDDO      ENDDO
902    
903      ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne      ! Si une sous-fraction n'existe pas, elle prend la température moyenne :
   
904      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
905         DO i = 1, klon         DO i = 1, klon
906            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
# Line 1107  contains Line 925  contains
925         ENDDO         ENDDO
926      ENDDO      ENDDO
927    
928      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
929    
930      DO i = 1, klon      DO i = 1, klon
931         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
932      ENDDO      ENDDO
933    
934      ! Appeler la convection (au choix)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
935    
936      DO k = 1, llm      ! Appeler la convection (au choix)
        DO i = 1, klon  
           conv_q(i, k) = d_q_dyn(i, k) &  
                + d_q_vdf(i, k)/dtphys  
           conv_t(i, k) = d_t_dyn(i, k) &  
                + d_t_vdf(i, k)/dtphys  
        ENDDO  
     ENDDO  
     IF (check) THEN  
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
     zx_ajustq = iflag_con == 2  
     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  
937    
938      select case (iflag_con)      if (iflag_con == 2) then
939      case (2)         conv_q = d_q_dyn + d_q_vdf / dtphys
940         CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &         conv_t = d_t_dyn + d_t_vdf / dtphys
941              zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
942              pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
943              pmflxs)              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
944                d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
945                mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
946                kdtop, pmflxr, pmflxs)
947         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
948         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
949         DO i = 1, klon         ibas_con = llm + 1 - kcbot
950            ibas_con(i) = llm + 1 - kcbot(i)         itop_con = llm + 1 - kctop
951            itop_con(i) = llm + 1 - kctop(i)      else
952         ENDDO         ! iflag_con >= 3
953      case (3:)  
954         ! number of tracers for the convection scheme of Kerry Emanuel:         da = 0.
955         ! la partie traceurs est faite dans phytrac         mp = 0.
956         ! on met ntra = 1 pour limiter les appels mais on peut         phi = 0.
957         ! supprimer les calculs / ftra.         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
958         ntra = 1              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &
959         ! Schéma de convection modularisé et vectorisé :              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &
960         ! (driver commun aux versions 3 et 4)              qcondc, wd, pmflxr, pmflxs, da, phi, mp)
   
        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)  
961         clwcon0 = qcondc         clwcon0 = qcondc
962         pmfu = upwd + dnwd         mfu = upwd + dnwd
963           IF (.NOT. ok_gust) wd = 0.
964    
965         IF (.NOT. ok_gust) THEN         IF (thermcep) THEN
966            do i = 1, klon            zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
967               wd(i) = 0.0            zqsat = zqsat / (1. - retv * zqsat)
968            enddo         ELSE
969              zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
970         ENDIF         ENDIF
971    
972         ! Calcul des propriétés des nuages convectifs         ! Properties of convective clouds
973           clwcon0 = fact_cldcon * clwcon0
        DO k = 1, llm  
           DO i = 1, klon  
              zx_t = t_seri(i, k)  
              IF (thermcep) THEN  
                 zdelta = MAX(0., SIGN(1., rtt-zx_t))  
                 zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)  
                 zx_qs = MIN(0.5, zx_qs)  
                 zcor = 1./(1.-retv*zx_qs)  
                 zx_qs = zx_qs*zcor  
              ELSE  
                 IF (zx_t < t_coup) THEN  
                    zx_qs = qsats(zx_t)/play(i, k)  
                 ELSE  
                    zx_qs = qsatl(zx_t)/play(i, k)  
                 ENDIF  
              ENDIF  
              zqsat(i, k) = zx_qs  
           ENDDO  
        ENDDO  
   
        ! calcul des proprietes des nuages convectifs  
        clwcon0 = fact_cldcon*clwcon0  
974         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
975              rnebcon0)              rnebcon0)
976      case default  
977         print *, "iflag_con non-prevu", iflag_con         mfd = 0.
978         stop 1         pen_u = 0.
979      END select         pen_d = 0.
980           pde_d = 0.
981           pde_u = 0.
982        END if
983    
984      DO k = 1, llm      DO k = 1, llm
985         DO i = 1, klon         DO i = 1, klon
# Line 1219  contains Line 993  contains
993      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
994         tit = 'after convect'         tit = 'after convect'
995         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
996              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
997         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
998              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)
             fs_bound, fq_bound)  
999      END IF      END IF
1000    
1001      IF (check) THEN      IF (check) THEN
1002         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1003         print *, "aprescon = ", za         print *, "aprescon = ", za
1004         zx_t = 0.0         zx_t = 0.
1005         za = 0.0         za = 0.
1006         DO i = 1, klon         DO i = 1, klon
1007            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1008            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
# Line 1239  contains Line 1011  contains
1011         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
1012         print *, "Precip = ", zx_t         print *, "Precip = ", zx_t
1013      ENDIF      ENDIF
1014      IF (zx_ajustq) THEN  
1015         DO i = 1, klon      IF (iflag_con == 2) THEN
1016            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1017         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  
1018         DO k = 1, llm         DO k = 1, llm
1019            DO i = 1, klon            DO i = 1, klon
1020               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 1261  contains Line 1023  contains
1023            ENDDO            ENDDO
1024         ENDDO         ENDDO
1025      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
1026    
1027      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
1028    
1029      d_t_ajs = 0.      d_t_ajs = 0.
1030      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1286  contains Line 1047  contains
1047      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1048         tit = 'after dry_adjust'         tit = 'after dry_adjust'
1049         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1050              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1051      END IF      END IF
1052    
1053      ! Caclul des ratqs      ! Caclul des ratqs
1054    
1055      ! 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
1056      ! on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
1057      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
1058         do k = 1, llm         do k = 1, llm
1059            do i = 1, klon            do i = 1, klon
1060               if(ptconv(i, k)) then               if(ptconv(i, k)) then
1061                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
1062                       +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)
1063               else               else
1064                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
1065               endif               endif
# Line 1310  contains Line 1070  contains
1070      ! ratqs stables      ! ratqs stables
1071      do k = 1, llm      do k = 1, llm
1072         do i = 1, klon         do i = 1, klon
1073            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1074                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1075         enddo         enddo
1076      enddo      enddo
1077    
1078      ! ratqs final      ! ratqs final
1079      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1080         ! les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
1081         ! ratqs final         ! ratqs final
1082         ! 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
1083         ! relaxation des ratqs         ! relaxation des ratqs
1084         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
1085         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
1086      else      else
1087         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
1088         ratqs = ratqss         ratqs = ratqss
1089      endif      endif
1090    
     ! Processus de condensation à grande echelle et processus de  
     ! précipitation :  
1091      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1092           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, &
1093           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
# Line 1348  contains Line 1105  contains
1105         ENDDO         ENDDO
1106      ENDDO      ENDDO
1107      IF (check) THEN      IF (check) THEN
1108         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1109         print *, "apresilp = ", za         print *, "apresilp = ", za
1110         zx_t = 0.0         zx_t = 0.
1111         za = 0.0         za = 0.
1112         DO i = 1, klon         DO i = 1, klon
1113            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1114            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
# Line 1364  contains Line 1121  contains
1121      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1122         tit = 'after fisrt'         tit = 'after fisrt'
1123         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1124              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1125         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1126              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)
             fs_bound, fq_bound)  
1127      END IF      END IF
1128    
1129      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
# Line 1393  contains Line 1148  contains
1148         endif         endif
1149    
1150         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1151         CALL diagcld1(paprs, play, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1152              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1153         DO k = 1, llm         DO k = 1, llm
1154            DO i = 1, klon            DO i = 1, klon
1155               IF (diafra(i, k) > cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
# Line 1405  contains Line 1159  contains
1159            ENDDO            ENDDO
1160         ENDDO         ENDDO
1161      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1162         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1163         ! 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
1164         ! facttemps         ! d'un facteur facttemps.
1165         facteur = dtphys *facttemps         facteur = dtphys * facttemps
1166         do k = 1, llm         do k = 1, llm
1167            do i = 1, klon            do i = 1, klon
1168               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1169               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1170                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1171                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1172                  clwcon(i, k) = clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1173               endif               endif
# Line 1446  contains Line 1200  contains
1200      ENDDO      ENDDO
1201    
1202      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1203           dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1204           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1205    
1206      ! Humidité relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
1207      DO k = 1, llm      DO k = 1, llm
1208         DO i = 1, klon         DO i = 1, klon
1209            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1210            IF (thermcep) THEN            IF (thermcep) THEN
1211               zdelta = MAX(0., SIGN(1., rtt-zx_t))               zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
              zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)  
1212               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1213               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1.-retv*zx_qs)
1214               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
# Line 1472  contains Line 1225  contains
1225      ENDDO      ENDDO
1226    
1227      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
     ! Johannes Quaas, 27/11/2003  
1228      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1229         ! Get sulfate aerosol distribution         ! Get sulfate aerosol distribution :
1230         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(rdayvrai, firstcal, sulfate)
1231         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1232    
        ! Calculate aerosol optical properties (Olivier Boucher)  
1233         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, &
1234              aerindex)              aerindex)
1235      ELSE      ELSE
# Line 1487  contains Line 1238  contains
1238         cg_ae = 0.         cg_ae = 0.
1239      ENDIF      ENDIF
1240    
1241      ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :      ! Param\`etres optiques des nuages et quelques param\`etres pour
1242        ! diagnostics :
1243      if (ok_newmicro) then      if (ok_newmicro) then
1244         CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1245              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1246              fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, &              sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
             re, fl)  
1247      else      else
1248         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1249              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
# Line 1533  contains Line 1284  contains
1284      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1285         tit = 'after rad'         tit = 'after rad'
1286         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1287              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1288         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1289              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)
             fs_bound, fq_bound)  
1290      END IF      END IF
1291    
1292      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
1293      DO i = 1, klon      DO i = 1, klon
1294         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1295         zxsnow(i) = 0.0         zxsnow(i) = 0.
1296      ENDDO      ENDDO
1297      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1298         DO i = 1, klon         DO i = 1, klon
# Line 1552  contains Line 1301  contains
1301         ENDDO         ENDDO
1302      ENDDO      ENDDO
1303    
1304      ! 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)
1305    
1306      DO i = 1, klon      DO i = 1, klon
1307         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1308      ENDDO      ENDDO
1309    
1310      ! Paramétrisation de l'orographie à l'échelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1311    
1312      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1313         ! selection des points pour lesquels le shema est actif:         ! selection des points pour lesquels le shema est actif:
1314         igwd = 0         igwd = 0
1315         DO i = 1, klon         DO i = 1, klon
1316            itest(i) = 0            itest(i) = 0
1317            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.0)) THEN            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN
1318               itest(i) = 1               itest(i) = 1
1319               igwd = igwd + 1               igwd = igwd + 1
1320               idx(igwd) = i               idx(igwd) = i
# Line 1587  contains Line 1336  contains
1336      ENDIF      ENDIF
1337    
1338      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1339         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
1340         igwd = 0         igwd = 0
1341         DO i = 1, klon         DO i = 1, klon
1342            itest(i) = 0            itest(i) = 0
# Line 1612  contains Line 1361  contains
1361         ENDDO         ENDDO
1362      ENDIF      ENDIF
1363    
1364      ! Stress nécessaires : toute la physique      ! Stress n\'ecessaires : toute la physique
1365    
1366      DO i = 1, klon      DO i = 1, klon
1367         zustrph(i) = 0.         zustrph(i) = 0.
# Line 1631  contains Line 1380  contains
1380           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1381    
1382      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1383           2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1384           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1385    
1386      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1387      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, u, t, &
1388           dtphys, u, t, paprs, play, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1389           ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, albsol, rhcl, &
1390           frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &           cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, &
1391           pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)           mp, upwd, dnwd, tr_seri, zmasse)
1392    
1393      IF (offline) THEN      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1394         call phystokenc(dtphys, rlon, rlat, t, pmfu, pmfd, pen_u, pde_u, &           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1395              pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
             pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)  
     ENDIF  
1396    
1397      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1398      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 1672  contains Line 1419  contains
1419      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1420         tit = 'after physic'         tit = 'after physic'
1421         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1422              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1423         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1424         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1425         ! 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.
1426         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1427         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1428              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)
             fs_bound, fq_bound)  
   
1429         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1430      END IF      END IF
1431    
1432      ! SORTIES      ! SORTIES
1433    
1434      !cc prw = eau precipitable      ! prw = eau precipitable
1435      DO i = 1, klon      DO i = 1, klon
1436         prw(i) = 0.         prw(i) = 0.
1437         DO k = 1, llm         DO k = 1, llm
# Line 1708  contains Line 1451  contains
1451         ENDDO         ENDDO
1452      ENDDO      ENDDO
1453    
1454      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1455         DO iq = 3, nqmx         DO k = 1, llm
1456            DO k = 1, llm            DO i = 1, klon
1457               DO i = 1, klon               d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys
                 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys  
              ENDDO  
1458            ENDDO            ENDDO
1459         ENDDO         ENDDO
1460      ENDIF      ENDDO
1461    
1462      ! Sauvegarder les valeurs de t et q a la fin de la physique:      ! Sauvegarder les valeurs de t et q a la fin de la physique:
1463      DO k = 1, llm      DO k = 1, llm
# Line 1727  contains Line 1468  contains
1468      ENDDO      ENDDO
1469    
1470      ! Ecriture des sorties      ! Ecriture des sorties
     call write_histhf  
     call write_histday  
1471      call write_histins      call write_histins
1472    
1473      ! 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 1736  contains Line 1475  contains
1475         itau_phy = itau_phy + itap         itau_phy = itau_phy + itap
1476         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &
1477              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &
1478              rain_fall, snow_fall, solsw, sollwdown, dlw, radsol, frugs, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
1479              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
1480              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
1481      ENDIF      ENDIF
1482    
1483      firstcal = .FALSE.      firstcal = .FALSE.
1484    
1485    contains    contains
1486    
     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  
   
     !***************************************************************  
   
1487      subroutine write_histins      subroutine write_histins
1488    
1489        ! 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
1490    
1491          use dimens_m, only: iim, jjm
1492          USE histsync_m, ONLY: histsync
1493          USE histwrite_m, ONLY: histwrite
1494    
1495        real zout        real zout
1496        integer itau_w ! pas de temps ecriture        integer itau_w ! pas de temps ecriture
1497          REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1498    
1499        !--------------------------------------------------        !--------------------------------------------------
1500    
# Line 2005  contains Line 1710  contains
1710           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)
1711           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1712    
1713           if (ok_sync) then           call histsync(nid_ins)
             call histsync(nid_ins)  
          endif  
1714        ENDIF        ENDIF
1715    
1716      end subroutine write_histins      end subroutine write_histins
1717    
     !****************************************************  
   
     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  
   
1718    END SUBROUTINE physiq    END SUBROUTINE physiq
1719    
1720  end module physiq_m  end module physiq_m
Legend:
Removed from v.62  
changed lines
  Added in v.103
  ViewVC Help
Powered by ViewVC 1.1.21