/[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 69 by guez, Mon Feb 18 16:33:12 2013 UTC trunk/phylmd/physiq.f revision 101 by guez, Mon Jul 7 17:45:21 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.
# Line 16  contains Line 18  contains
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
19      use aeropt_m, only: aeropt      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 32  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      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 newmicro_m, only: newmicro      use newmicro_m, only: newmicro
48      USE oasis_m, ONLY: ok_oasis      USE orbite_m, ONLY: orbite
     USE orbite_m, ONLY: orbite, zenang  
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 61  contains Line 58  contains
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      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):: ocean = 'force '  
     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")  
   
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    
     ! Modele thermique du sol, a activer pour le cycle diurne:  
     logical:: ok_veget = .false. ! type de modele de vegetation utilise  
   
120      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
121      ! sorties journalieres, mensuelles et instantanees dans les      ! sorties journalieres, mensuelles et instantanees dans les
122      ! fichiers histday, histmth et histins      ! fichiers histday, histmth et histins
# 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)
     SAVE ema_work1, ema_work2  
279      REAL, save:: wd(klon)      REAL, save:: wd(klon)
280    
281      ! Variables locales pour la couche limite (al1):      ! Variables locales pour la couche limite (al1):
# Line 346  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 363  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  
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 417  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 444  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)
398      REAL za, zb      REAL za, zb
# Line 454  contains Line 402  contains
402      REAL, PARAMETER:: t_coup = 234.      REAL, PARAMETER:: t_coup = 234.
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 477  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:  
     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 505  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 539  contains Line 478  contains
478      integer:: iflag_cldcon = 1      integer:: iflag_cldcon = 1
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 555  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 568  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
   
     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:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      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
# Line 634  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/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &      namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
575           fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &           facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
576           ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &           ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals
          nsplit_thermals  
577    
578      !----------------------------------------------------------------      !----------------------------------------------------------------
579    
580      IF (if_ebil >= 1) zero_v = 0.      IF (if_ebil >= 1) zero_v = 0.
     ok_sync = .TRUE.  
581      IF (nqmx < 2) CALL abort_gcm('physiq', &      IF (nqmx < 2) CALL abort_gcm('physiq', &
582           'eaux vapeur et liquide sont indispensables', 1)           'eaux vapeur et liquide sont indispensables', 1)
583    
# Line 665  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         topswai(:) = 0.         topswai = 0.
598         topswad(:) = 0.         topswad = 0.
599         solswai(:) = 0.         solswai = 0.
600         solswad(:) = 0.         solswad = 0.
601    
602         d_u_con = 0.0         d_u_con = 0.
603         d_v_con = 0.0         d_v_con = 0.
604         rnebcon0 = 0.0         rnebcon0 = 0.
605         clwcon0 = 0.0         clwcon0 = 0.
606         rnebcon = 0.0         rnebcon = 0.
607         clwcon = 0.0         clwcon = 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 705  contains Line 632  contains
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 = 1e-8         q2 = 1e-8
# Line 720  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, ocean /= 'force', ok_oasis, ok_journe, &         CALL printflag(radpas, ok_journe, ok_instan, ok_region)
             ok_instan, ok_region)  
651    
652         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
653            print *, "Au minimum 4 appels par jour si cycle diurne"            print *, "Au minimum 4 appels par jour si cycle diurne"
# Line 729  contains Line 655  contains
655                 "Nombre d'appels au rayonnement insuffisant", 1)                 "Nombre d'appels au rayonnement insuffisant", 1)
656         ENDIF         ENDIF
657    
658         ! Initialisation pour le schéma de convection d'Emanuel :         ! Initialisation pour le sch\'ema de convection d'Emanuel :
659         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
660            ibas_con = 1            ibas_con = 1
661            itop_con = 1            itop_con = 1
# Line 751  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         ! Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
685         print *, '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        t_seri = t
691        u_seri = u
692        v_seri = v
693        q_seri = qx(:, :, ivap)
694        ql_seri = qx(:, :, iliq)
695        tr_seri = qx(:, :, 3: nqmx)
696    
697      DO i = 1, klon      ztsol = sum(ftsol * pctsrf, dim = 2)
        d_ps(i) = 0.  
     ENDDO  
     DO iq = 1, nqmx  
        DO k = 1, llm  
           DO i = 1, klon  
              d_qx(i, k, iq) = 0.  
           ENDDO  
        ENDDO  
     ENDDO  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
     ! Ne pas affecter les valeurs entrées de u, v, h, et q :  
   
     DO k = 1, llm  
        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 835  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 852  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
747    
748      ! Mettre en action les conditions aux limites (albedo, sst etc.).      ! Prescrire l'ozone :
   
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.  
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 878  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)
     DO i = 1, klon  
        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  
772    
773      ! calculs necessaires au calcul de l'albedo dans l'interface      ! Calculs nécessaires 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. * 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 957  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'évaporation 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 986  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) + pctsrf(i, is_oce) &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
873              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
874              'physiq : problème sous surface au point ', i, pctsrf(i, 1 : nbsrf)              'physiq : probl\`eme sous surface au point ', i, &
875                pctsrf(i, 1 : nbsrf)
876      ENDDO      ENDDO
877      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
878         DO i = 1, klon         DO i = 1, klon
# Line 1047  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 1073  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. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
# Line 1083  contains Line 935  contains
935    
936      DO k = 1, llm      DO k = 1, llm
937         DO i = 1, klon         DO i = 1, klon
938            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k) / dtphys
939            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
940         ENDDO         ENDDO
941      ENDDO      ENDDO
942    
943      IF (check) THEN      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
944    
945      if (iflag_con == 2) then      if (iflag_con == 2) then
946         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
947         CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
948              zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
949              pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
950              pmflxs)              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
951                kdtop, pmflxr, pmflxs)
952         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
953         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
954         DO i = 1, klon         ibas_con = llm + 1 - kcbot
955            ibas_con(i) = llm + 1 - kcbot(i)         itop_con = llm + 1 - kctop
           itop_con(i) = llm + 1 - kctop(i)  
        ENDDO  
956      else      else
957         ! iflag_con >= 3         ! iflag_con >= 3
        CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &  
             v_seri, tr_seri, 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, ntra=1)  
        ! (number of tracers for the convection scheme of Kerry Emanuel:  
        ! la partie traceurs est faite dans phytrac  
        ! on met ntra = 1 pour limiter les appels mais on peut  
        ! supprimer les calculs / ftra.)  
958    
959           da = 0.
960           mp = 0.
961           phi = 0.
962           CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
963                w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &
964                ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &
965                qcondc, wd, pmflxr, pmflxs, da, phi, mp)
966         clwcon0 = qcondc         clwcon0 = qcondc
967         pmfu = upwd + dnwd         mfu = upwd + dnwd
968         IF (.NOT. ok_gust) wd = 0.         IF (.NOT. ok_gust) wd = 0.
969    
970         ! Calcul des propriétés des nuages convectifs         ! Calcul des propri\'et\'es des nuages convectifs
971    
972         DO k = 1, llm         DO k = 1, llm
973            DO i = 1, klon            DO i = 1, klon
              zx_t = t_seri(i, k)  
974               IF (thermcep) THEN               IF (thermcep) THEN
975                  zdelta = MAX(0., SIGN(1., rtt-zx_t))                  zdelta = MAX(0., SIGN(1., rtt - t_seri(i, k)))
976                  zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)                  zqsat(i, k) = r2es * FOEEW(t_seri(i, k), zdelta) / play(i, k)
977                  zx_qs = MIN(0.5, zx_qs)                  zqsat(i, k) = MIN(0.5, zqsat(i, k))
978                  zcor = 1./(1.-retv*zx_qs)                  zqsat(i, k) = zqsat(i, k) / (1.-retv*zqsat(i, k))
                 zx_qs = zx_qs*zcor  
979               ELSE               ELSE
980                  IF (zx_t < t_coup) THEN                  IF (t_seri(i, k) < t_coup) THEN
981                     zx_qs = qsats(zx_t)/play(i, k)                     zqsat(i, k) = qsats(t_seri(i, k))/play(i, k)
982                  ELSE                  ELSE
983                     zx_qs = qsatl(zx_t)/play(i, k)                     zqsat(i, k) = qsatl(t_seri(i, k))/play(i, k)
984                  ENDIF                  ENDIF
985               ENDIF               ENDIF
              zqsat(i, k) = zx_qs  
986            ENDDO            ENDDO
987         ENDDO         ENDDO
988    
989         ! calcul des proprietes des nuages convectifs         ! calcul des proprietes des nuages convectifs
990         clwcon0 = fact_cldcon*clwcon0         clwcon0 = fact_cldcon * clwcon0
991         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
992              rnebcon0)              rnebcon0)
993    
994           mfd = 0.
995           pen_u = 0.
996           pen_d = 0.
997           pde_d = 0.
998           pde_u = 0.
999      END if      END if
1000    
1001      DO k = 1, llm      DO k = 1, llm
# Line 1162  contains Line 1010  contains
1010      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1011         tit = 'after convect'         tit = 'after convect'
1012         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, &
1013              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)  
1014         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, &
1015              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)  
1016      END IF      END IF
1017    
1018      IF (check) THEN      IF (check) THEN
1019         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1020         print *, "aprescon = ", za         print *, "aprescon = ", za
1021         zx_t = 0.0         zx_t = 0.
1022         za = 0.0         za = 0.
1023         DO i = 1, klon         DO i = 1, klon
1024            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1025            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
# Line 1195  contains Line 1041  contains
1041         ENDDO         ENDDO
1042      ENDIF      ENDIF
1043    
1044      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
1045    
1046      d_t_ajs = 0.      d_t_ajs = 0.
1047      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1218  contains Line 1064  contains
1064      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1065         tit = 'after dry_adjust'         tit = 'after dry_adjust'
1066         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, &
1067              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)  
1068      END IF      END IF
1069    
1070      ! Caclul des ratqs      ! Caclul des ratqs
1071    
1072      ! 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
1073      ! on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
1074      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
1075         do k = 1, llm         do k = 1, llm
1076            do i = 1, klon            do i = 1, klon
1077               if(ptconv(i, k)) then               if(ptconv(i, k)) then
1078                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
1079                       +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)
1080               else               else
1081                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
1082               endif               endif
# Line 1242  contains Line 1087  contains
1087      ! ratqs stables      ! ratqs stables
1088      do k = 1, llm      do k = 1, llm
1089         do i = 1, klon         do i = 1, klon
1090            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1091                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1092         enddo         enddo
1093      enddo      enddo
1094    
# Line 1253  contains Line 1098  contains
1098         ! ratqs final         ! ratqs final
1099         ! 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
1100         ! relaxation des ratqs         ! relaxation des ratqs
1101         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
1102         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
1103      else      else
1104         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
1105         ratqs = ratqss         ratqs = ratqss
1106      endif      endif
1107    
     ! Processus de condensation à grande echelle et processus de  
     ! précipitation :  
1108      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1109           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, &
1110           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
# Line 1280  contains Line 1122  contains
1122         ENDDO         ENDDO
1123      ENDDO      ENDDO
1124      IF (check) THEN      IF (check) THEN
1125         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1126         print *, "apresilp = ", za         print *, "apresilp = ", za
1127         zx_t = 0.0         zx_t = 0.
1128         za = 0.0         za = 0.
1129         DO i = 1, klon         DO i = 1, klon
1130            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1131            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
# Line 1296  contains Line 1138  contains
1138      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1139         tit = 'after fisrt'         tit = 'after fisrt'
1140         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, &
1141              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)  
1142         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, &
1143              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)  
1144      END IF      END IF
1145    
1146      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
# Line 1336  contains Line 1176  contains
1176            ENDDO            ENDDO
1177         ENDDO         ENDDO
1178      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1179         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1180         ! 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
1181         ! facttemps         ! d'un facteur facttemps.
1182         facteur = dtphys *facttemps         facteur = dtphys * facttemps
1183         do k = 1, llm         do k = 1, llm
1184            do i = 1, klon            do i = 1, klon
1185               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1186               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1187                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1188                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1189                  clwcon(i, k) = clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1190               endif               endif
# Line 1377  contains Line 1217  contains
1217      ENDDO      ENDDO
1218    
1219      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, &
1220           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, &
1221           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1222    
1223      ! Humidité relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
1224      DO k = 1, llm      DO k = 1, llm
1225         DO i = 1, klon         DO i = 1, klon
1226            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
# Line 1416  contains Line 1256  contains
1256         cg_ae = 0.         cg_ae = 0.
1257      ENDIF      ENDIF
1258    
1259      ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :      ! Param\`etres optiques des nuages et quelques param\`etres pour
1260        ! diagnostics :
1261      if (ok_newmicro) then      if (ok_newmicro) then
1262         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1263              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
# Line 1461  contains Line 1302  contains
1302      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1303         tit = 'after rad'         tit = 'after rad'
1304         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, &
1305              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)  
1306         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1307              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)  
1308      END IF      END IF
1309    
1310      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
1311      DO i = 1, klon      DO i = 1, klon
1312         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1313         zxsnow(i) = 0.0         zxsnow(i) = 0.
1314      ENDDO      ENDDO
1315      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1316         DO i = 1, klon         DO i = 1, klon
# Line 1480  contains Line 1319  contains
1319         ENDDO         ENDDO
1320      ENDDO      ENDDO
1321    
1322      ! 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)
1323    
1324      DO i = 1, klon      DO i = 1, klon
1325         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1326      ENDDO      ENDDO
1327    
1328      ! Paramétrisation de l'orographie à l'échelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1329    
1330      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1331         ! selection des points pour lesquels le shema est actif:         ! selection des points pour lesquels le shema est actif:
1332         igwd = 0         igwd = 0
1333         DO i = 1, klon         DO i = 1, klon
1334            itest(i) = 0            itest(i) = 0
1335            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.0)) THEN            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN
1336               itest(i) = 1               itest(i) = 1
1337               igwd = igwd + 1               igwd = igwd + 1
1338               idx(igwd) = i               idx(igwd) = i
# Line 1515  contains Line 1354  contains
1354      ENDIF      ENDIF
1355    
1356      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1357         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
1358         igwd = 0         igwd = 0
1359         DO i = 1, klon         DO i = 1, klon
1360            itest(i) = 0            itest(i) = 0
# Line 1540  contains Line 1379  contains
1379         ENDDO         ENDDO
1380      ENDIF      ENDIF
1381    
1382      ! Stress nécessaires : toute la physique      ! Stress n\'ecessaires : toute la physique
1383    
1384      DO i = 1, klon      DO i = 1, klon
1385         zustrph(i) = 0.         zustrph(i) = 0.
# Line 1559  contains Line 1398  contains
1398           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1399    
1400      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1401           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, &
1402           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1403    
1404      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1405      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, u, t, &
1406           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, &
1407           ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, albsol, rhcl, &
1408           frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &           cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, &
1409           pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)           mp, upwd, dnwd, tr_seri, zmasse)
1410    
1411      IF (offline) THEN      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1412         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, &
1413              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  
1414    
1415      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1416      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 1600  contains Line 1437  contains
1437      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1438         tit = 'after physic'         tit = 'after physic'
1439         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, &
1440              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)  
1441         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1442         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1443         ! 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.
1444         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1445         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1446              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)  
   
1447         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1448      END IF      END IF
1449    
1450      ! SORTIES      ! SORTIES
# Line 1636  contains Line 1469  contains
1469         ENDDO         ENDDO
1470      ENDDO      ENDDO
1471    
1472      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1473         DO iq = 3, nqmx         DO k = 1, llm
1474            DO k = 1, llm            DO i = 1, klon
1475               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  
1476            ENDDO            ENDDO
1477         ENDDO         ENDDO
1478      ENDIF      ENDDO
1479    
1480      ! 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:
1481      DO k = 1, llm      DO k = 1, llm
# Line 1655  contains Line 1486  contains
1486      ENDDO      ENDDO
1487    
1488      ! Ecriture des sorties      ! Ecriture des sorties
     call write_histhf  
     call write_histday  
1489      call write_histins      call write_histins
1490    
1491      ! 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 1664  contains Line 1493  contains
1493         itau_phy = itau_phy + itap         itau_phy = itau_phy + itap
1494         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &
1495              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &
1496              rain_fall, snow_fall, solsw, sollwdown, dlw, radsol, frugs, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
1497              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
1498              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
1499      ENDIF      ENDIF
1500    
1501      firstcal = .FALSE.      firstcal = .FALSE.
1502    
1503    contains    contains
1504    
     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  
   
     !***************************************************************  
   
1505      subroutine write_histins      subroutine write_histins
1506    
1507        ! 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
1508    
1509          use dimens_m, only: iim, jjm
1510          USE histsync_m, ONLY: histsync
1511          USE histwrite_m, ONLY: histwrite
1512    
1513        real zout        real zout
1514        integer itau_w ! pas de temps ecriture        integer itau_w ! pas de temps ecriture
1515          REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1516    
1517        !--------------------------------------------------        !--------------------------------------------------
1518    
# Line 1933  contains Line 1728  contains
1728           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)
1729           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1730    
1731           if (ok_sync) then           call histsync(nid_ins)
             call histsync(nid_ins)  
          endif  
1732        ENDIF        ENDIF
1733    
1734      end subroutine write_histins      end subroutine write_histins
1735    
     !****************************************************  
   
     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  
   
1736    END SUBROUTINE physiq    END SUBROUTINE physiq
1737    
1738  end module physiq_m  end module physiq_m
Legend:
Removed from v.69  
changed lines
  Added in v.101
  ViewVC Help
Powered by ViewVC 1.1.21