/[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/phylmd/physiq.f revision 101 by guez, Mon Jul 7 17:45:21 2014 UTC trunk/Sources/phylmd/physiq.f revision 200 by guez, Thu Jun 2 15:40:30 2016 UTC
# Line 4  module physiq_m Line 4  module physiq_m
4    
5  contains  contains
6    
7    SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &    SUBROUTINE physiq(lafin, dayvrai, time, paprs, play, pphi, pphis, u, v, t, &
8         u, v, t, qx, omega, d_u, d_v, d_t, d_qx)         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      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11      ! (subversion revision 678)      ! (subversion revision 678)
12    
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
     use aeropt_m, only: aeropt  
19      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
20      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
21      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &      USE clesphys, ONLY: cdhmax, cdmmax, ecrit_hf, ecrit_ins, ecrit_mth, &
22           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin           ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin, ok_instan
23      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &      USE clesphys2, ONLY: cycle_diurne, conv_emanuel, nbapp_rad, new_oliq, &
24           ok_orodr, ok_orolf           ok_orodr, ok_orolf
25      USE clmain_m, ONLY: clmain      USE clmain_m, ONLY: clmain
26      use clouds_gno_m, only: clouds_gno      use clouds_gno_m, only: clouds_gno
27      USE comgeomphy, ONLY: airephy, cuphy, cvphy      use comconst, only: dtphys
28        USE comgeomphy, ONLY: airephy
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, day_step, iphysiq
31      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
32      use conflx_m, only: conflx      use conflx_m, only: conflx
33      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
# Line 38  contains Line 38  contains
38      USE dimphy, ONLY: klon      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 dynetat0_m, only: day_ref, annee_ref
42      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
43      use fisrtilp_m, only: fisrtilp      use fisrtilp_m, only: fisrtilp
44      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
45        USE histsync_m, ONLY: histsync
46        USE histwrite_phy_m, ONLY: histwrite_phy
47      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, &
48           nbsrf           nbsrf
49      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins, nid_ins
50        use netcdf95, only: NF95_CLOSE
51      use newmicro_m, only: newmicro      use newmicro_m, only: newmicro
52        use nr_util, only: assert
53        use nuage_m, only: nuage
54      USE orbite_m, ONLY: orbite      USE orbite_m, ONLY: orbite
55      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
56      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
57      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
58        USE phyredem0_m, ONLY: phyredem0
59      USE phystokenc_m, ONLY: phystokenc      USE phystokenc_m, ONLY: phystokenc
60      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
61      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
62      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
63      use readsulfate_m, only: readsulfate      use yoegwd, only: sugwd
64      use sugwd_m, only: sugwd      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt
65      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      use time_phylmdz, only: itap, increment_itap
66      USE temps, ONLY: annee_ref, day_ref, itau_phy      use transp_m, only: transp
67        use transp_lay_m, only: transp_lay
68      use unit_nml_m, only: unit_nml      use unit_nml_m, only: unit_nml
69      USE ymds2ju_m, ONLY: ymds2ju      USE ymds2ju_m, ONLY: ymds2ju
70      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
# Line 64  contains Line 72  contains
72    
73      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
74    
75      REAL, intent(in):: rdayvrai      integer, intent(in):: dayvrai
76      ! (elapsed time since January 1st 0h of the starting year, in days)      ! current day number, based at value 1 on January 1st of annee_ref
77    
78      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
     REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)  
79    
80      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
81      ! pression pour chaque inter-couche, en Pa      ! pression pour chaque inter-couche, en Pa
# Line 76  contains Line 83  contains
83      REAL, intent(in):: play(:, :) ! (klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
84      ! pression pour le mileu de chaque couche (en Pa)      ! pression pour le mileu de chaque couche (en Pa)
85    
86      REAL, intent(in):: pphi(:, :) ! (klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
87      ! géopotentiel de chaque couche (référence sol)      ! géopotentiel de chaque couche (référence sol)
88    
89      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
# Line 102  contains Line 109  contains
109    
110      LOGICAL:: firstcal = .true.      LOGICAL:: firstcal = .true.
111    
112      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL, PARAMETER:: check = .FALSE.
     PARAMETER (ok_gust = .FALSE.)  
   
     LOGICAL, PARAMETER:: check = .FALSE.  
113      ! Verifier la conservation du modele en eau      ! Verifier la conservation du modele en eau
114    
115      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
116      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
117    
118      ! "slab" ocean      ! pour phystoke avec thermiques
     REAL, save:: tslab(klon) ! temperature of ocean slab  
     REAL, save:: seaice(klon) ! glace de mer (kg/m2)  
     REAL fluxo(klon) ! flux turbulents ocean-glace de mer  
     REAL fluxg(klon) ! flux turbulents ocean-atmosphere  
   
     logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.  
     ! sorties journalieres, mensuelles et instantanees dans les  
     ! fichiers histday, histmth et histins  
   
     LOGICAL ok_region ! sortir le fichier regional  
     PARAMETER (ok_region = .FALSE.)  
   
     ! pour phsystoke avec thermiques  
119      REAL fm_therm(klon, llm + 1)      REAL fm_therm(klon, llm + 1)
120      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
121      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
# Line 146  contains Line 137  contains
137    
138      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
139      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
140      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
   
     ! Amip2  
     ! variables a une pression donnee  
   
     integer nlevSTD  
     PARAMETER(nlevSTD = 17)  
     real rlevSTD(nlevSTD)  
     DATA rlevSTD/100000., 92500., 85000., 70000., &  
          60000., 50000., 40000., 30000., 25000., 20000., &  
          15000., 10000., 7000., 5000., 3000., 2000., 1000./  
     CHARACTER(LEN = 4) clevSTD(nlevSTD)  
     DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &  
          '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &  
          '70 ', '50 ', '30 ', '20 ', '10 '/  
141    
142      ! prw: precipitable water      ! prw: precipitable water
143      real prw(klon)      real prw(klon)
# Line 170  contains Line 147  contains
147      REAL flwp(klon), fiwp(klon)      REAL flwp(klon), fiwp(klon)
148      REAL flwc(klon, llm), fiwc(klon, llm)      REAL flwc(klon, llm), fiwc(klon, llm)
149    
     INTEGER kmax, lmax  
     PARAMETER(kmax = 8, lmax = 8)  
     INTEGER kmaxm1, lmaxm1  
     PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)  
   
     REAL zx_tau(kmaxm1), zx_pc(lmaxm1)  
     DATA zx_tau/0., 0.3, 1.3, 3.6, 9.4, 23., 60./  
     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./  
   
     ! cldtopres pression au sommet des nuages  
     REAL cldtopres(lmaxm1)  
     DATA cldtopres/50., 180., 310., 440., 560., 680., 800./  
   
     ! taulev: numero du niveau de tau dans les sorties ISCCP  
     CHARACTER(LEN = 4) taulev(kmaxm1)  
   
     DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/  
     CHARACTER(LEN = 3) pclev(lmaxm1)  
     DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/  
   
     CHARACTER(LEN = 28) cnameisccp(lmaxm1, kmaxm1)  
     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &  
          'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &  
          'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &  
          'pc= 680-800hPa, tau< 0.3', 'pc< 50hPa, tau= 0.3-1.3', &  
          'pc= 50-180hPa, tau= 0.3-1.3', 'pc= 180-310hPa, tau= 0.3-1.3', &  
          'pc= 310-440hPa, tau= 0.3-1.3', 'pc= 440-560hPa, tau= 0.3-1.3', &  
          'pc= 560-680hPa, tau= 0.3-1.3', 'pc= 680-800hPa, tau= 0.3-1.3', &  
          'pc< 50hPa, tau= 1.3-3.6', 'pc= 50-180hPa, tau= 1.3-3.6', &  
          'pc= 180-310hPa, tau= 1.3-3.6', 'pc= 310-440hPa, tau= 1.3-3.6', &  
          'pc= 440-560hPa, tau= 1.3-3.6', 'pc= 560-680hPa, tau= 1.3-3.6', &  
          'pc= 680-800hPa, tau= 1.3-3.6', 'pc< 50hPa, tau= 3.6-9.4', &  
          'pc= 50-180hPa, tau= 3.6-9.4', 'pc= 180-310hPa, tau= 3.6-9.4', &  
          'pc= 310-440hPa, tau= 3.6-9.4', 'pc= 440-560hPa, tau= 3.6-9.4', &  
          'pc= 560-680hPa, tau= 3.6-9.4', 'pc= 680-800hPa, tau= 3.6-9.4', &  
          'pc< 50hPa, tau= 9.4-23', 'pc= 50-180hPa, tau= 9.4-23', &  
          'pc= 180-310hPa, tau= 9.4-23', 'pc= 310-440hPa, tau= 9.4-23', &  
          'pc= 440-560hPa, tau= 9.4-23', 'pc= 560-680hPa, tau= 9.4-23', &  
          'pc= 680-800hPa, tau= 9.4-23', 'pc< 50hPa, tau= 23-60', &  
          'pc= 50-180hPa, tau= 23-60', 'pc= 180-310hPa, tau= 23-60', &  
          'pc= 310-440hPa, tau= 23-60', 'pc= 440-560hPa, tau= 23-60', &  
          'pc= 560-680hPa, tau= 23-60', 'pc= 680-800hPa, tau= 23-60', &  
          'pc< 50hPa, tau> 60.', 'pc= 50-180hPa, tau> 60.', &  
          'pc= 180-310hPa, tau> 60.', 'pc= 310-440hPa, tau> 60.', &  
          'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &  
          'pc= 680-800hPa, tau> 60.'/  
   
     ! ISCCP simulator v3.4  
   
150      ! Variables propres a la physique      ! Variables propres a la physique
151    
152      INTEGER, save:: radpas      INTEGER, save:: radpas
153      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
154      ! "physiq".)      ! "physiq".
155    
156      REAL radsol(klon)      REAL radsol(klon)
157      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
158    
     INTEGER, SAVE:: itap ! number of calls to "physiq"  
   
159      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
160    
161      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
# Line 246  contains Line 172  contains
172      ! column-density of water in soil, in kg m-2      ! column-density of water in soil, in kg m-2
173    
174      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
175      REAL, save:: falbe(klon, nbsrf) ! albedo par type de surface      REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
     REAL, save:: falblw(klon, nbsrf) ! albedo par type de surface  
176    
177      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
178      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
# Line 258  contains Line 183  contains
183      REAL, save:: zpic(klon) ! Maximum de l'OESM      REAL, save:: zpic(klon) ! Maximum de l'OESM
184      REAL, save:: zval(klon) ! Minimum de l'OESM      REAL, save:: zval(klon) ! Minimum de l'OESM
185      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
   
186      REAL zulow(klon), zvlow(klon)      REAL zulow(klon), zvlow(klon)
187        INTEGER igwd, itest(klon)
188    
189      INTEGER igwd, idx(klon), itest(klon)      REAL, save:: agesno(klon, nbsrf) ! age de la neige
190        REAL, save:: run_off_lic_0(klon)
191    
192      REAL agesno(klon, nbsrf)      ! Variables li\'ees \`a la convection d'Emanuel :
193      SAVE agesno ! age de la neige      REAL, save:: Ma(klon, llm) ! undilute upward mass flux
194        REAL, save:: qcondc(klon, llm) ! in-cld water content from convect
     REAL run_off_lic_0(klon)  
     SAVE run_off_lic_0  
     !KE43  
     ! Variables liees a la convection de K. Emanuel (sb):  
   
     REAL Ma(klon, llm) ! undilute upward mass flux  
     SAVE Ma  
     REAL qcondc(klon, llm) ! in-cld water content from convect  
     SAVE qcondc  
195      REAL, save:: sig1(klon, llm), w01(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
     REAL, save:: wd(klon)  
   
     ! Variables locales pour la couche limite (al1):  
   
     ! Variables locales:  
196    
197        ! Variables pour la couche limite (Alain Lahellec) :
198      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
199      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
200    
# Line 289  contains Line 202  contains
202      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
203      REAL yu1(klon) ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
204      REAL yv1(klon) ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
205      REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige      REAL ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
206      REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface  
207      ! !et necessaire pour limiter la      REAL fqcalving(klon, nbsrf)
208      ! !hauteur de neige, en kg/m2/s      ! flux d'eau "perdue" par la surface et necessaire pour limiter la
209        ! hauteur de neige, en kg/m2/s
210    
211      REAL zxffonte(klon), zxfqcalving(klon)      REAL zxffonte(klon), zxfqcalving(klon)
212    
213      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
# Line 317  contains Line 232  contains
232      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
233      SAVE dlw      SAVE dlw
234      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
235      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
     save fder  
236      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
237      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
238      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
# Line 333  contains Line 247  contains
247      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
248      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
249      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
250      REAL, save:: albsol(klon) ! albedo du sol total      REAL, save:: albsol(klon) ! albedo du sol total visible
     REAL, save:: albsollw(klon) ! albedo du sol total  
251      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
252    
     ! Declaration des procedures appelees  
   
     EXTERNAL nuage ! calculer les proprietes radiatives  
     EXTERNAL transp ! transport total de l'eau et de l'energie  
   
     ! Variables locales  
   
253      real, save:: clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
254      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
255    
# Line 368  contains Line 274  contains
274      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
275      ! les variables soient r\'emanentes.      ! les variables soient r\'emanentes.
276      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
277      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
278      REAL, save:: cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
279      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
280      REAL, save:: topsw(klon), toplw(klon), solsw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
281      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
282      real, save:: sollwdown(klon) ! downward LW flux at surface      real, save:: sollwdown(klon) ! downward LW flux at surface
283      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
284      REAL albpla(klon)      REAL, save:: albpla(klon)
285      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
286      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
     SAVE albpla  
     SAVE heat0, cool0  
   
     INTEGER itaprad  
     SAVE itaprad  
287    
288      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
289      REAL conv_t(klon, llm) ! convergence of temperature (K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K/s)
290    
291      REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut      REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
292      REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree      REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
293    
294      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
295    
296      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
297      real zlongi      real longi
298      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
299      REAL za, zb      REAL za, zb
300      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
301      real zqsat(klon, llm)      real zqsat(klon, llm)
302      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
303      REAL, PARAMETER:: t_coup = 234.      REAL, PARAMETER:: t_coup = 234.
304      REAL zphi(klon, llm)      REAL zphi(klon, llm)
305    
306      ! cf. AM Variables locales pour la CLA (hbtm2)      ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
307    
308      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
309      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
# Line 412  contains Line 313  contains
313      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
314      REAL, SAVE:: therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
315      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
316      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
317      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
318      ! Grdeurs de sorties      ! Grandeurs de sorties
319      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
320      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
321      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
322      REAL s_trmb3(klon)      REAL s_trmb3(klon)
323    
324      ! Variables locales pour la convection de K. Emanuel :      ! Variables pour la convection de K. Emanuel :
325    
326      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
327      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
# Line 452  contains Line 353  contains
353      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
354    
355      INTEGER, save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
356        real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
357    
358      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
359      REAL snow_con(klon), snow_lsc(klon)      REAL, save:: snow_con(klon) ! neige (mm / s)
360        real snow_lsc(klon)
361      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
362    
363      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
# Line 478  contains Line 381  contains
381      integer:: iflag_cldcon = 1      integer:: iflag_cldcon = 1
382      logical ptconv(klon, llm)      logical ptconv(klon, llm)
383    
384      ! Variables locales pour effectuer les appels en s\'erie :      ! Variables pour effectuer les appels en s\'erie :
385    
386      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
387      REAL ql_seri(klon, llm)      REAL ql_seri(klon, llm)
# Line 492  contains Line 395  contains
395      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
396      REAL aam, torsfc      REAL aam, torsfc
397    
     REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique  
   
     INTEGER, SAVE:: nid_ins  
   
398      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.
399      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.
400      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
401      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.
402    
     REAL zsto  
403      real date0      real date0
404    
405      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
# Line 511  contains Line 409  contains
409      REAL zero_v(klon)      REAL zero_v(klon)
410      CHARACTER(LEN = 20) tit      CHARACTER(LEN = 20) tit
411      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
412      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
413    
414        REAL d_t_ec(klon, llm)
415        ! tendance due \`a la conversion Ec en énergie thermique
416    
     REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique  
417      REAL ZRCPD      REAL ZRCPD
418    
419      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 526  contains Line 426  contains
426      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
427    
428      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
429      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value      ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value
430    
431      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
432      ! cloud optical thickness for pre-industrial (pi) aerosols      ! cloud optical thickness for pre-industrial aerosols
433    
434      REAL re(klon, llm) ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
435      REAL fl(klon, llm) ! denominator of re      REAL fl(klon, llm) ! denominator of re
# Line 541  contains Line 441  contains
441      REAL topswad(klon), solswad(klon) ! aerosol direct effect      REAL topswad(klon), solswad(klon) ! aerosol direct effect
442      REAL topswai(klon), solswai(klon) ! aerosol indirect effect      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
443    
     REAL aerindex(klon) ! POLDER aerosol index  
   
444      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
445      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
446    
# Line 558  contains Line 456  contains
456      SAVE ffonte      SAVE ffonte
457      SAVE fqcalving      SAVE fqcalving
458      SAVE rain_con      SAVE rain_con
     SAVE snow_con  
459      SAVE topswai      SAVE topswai
460      SAVE topswad      SAVE topswad
461      SAVE solswai      SAVE solswai
# Line 566  contains Line 463  contains
463      SAVE d_u_con      SAVE d_u_con
464      SAVE d_v_con      SAVE d_v_con
465    
466      real zmasse(klon, llm)      real zmasse(klon, llm)
467      ! (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)
468    
469      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      integer, save:: ncid_startphy
470    
471      namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &      namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, &
472           facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &           iflag_cldcon, ratqsbas, ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, &
473           ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals           bl95_b1, iflag_thermals, nsplit_thermals
474    
475      !----------------------------------------------------------------      !----------------------------------------------------------------
476    
477      IF (if_ebil >= 1) zero_v = 0.      IF (if_ebil >= 1) zero_v = 0.
478      IF (nqmx < 2) CALL abort_gcm('physiq', &      IF (nqmx < 2) CALL abort_gcm('physiq', &
479           'eaux vapeur et liquide sont indispensables', 1)           'eaux vapeur et liquide sont indispensables')
480    
481      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
482         ! initialiser         ! initialiser
# Line 614  contains Line 511  contains
511         pblt =0. ! T a la Hauteur de couche limite         pblt =0. ! T a la Hauteur de couche limite
512         therm =0.         therm =0.
513         trmb1 =0. ! deep_cape         trmb1 =0. ! deep_cape
514         trmb2 =0. ! inhibition         trmb2 =0. ! inhibition
515         trmb3 =0. ! Point Omega         trmb3 =0. ! Point Omega
516    
517         IF (if_ebil >= 1) d_h_vcol_phy = 0.         IF (if_ebil >= 1) d_h_vcol_phy = 0.
# Line 630  contains Line 527  contains
527         ! Initialiser les compteurs:         ! Initialiser les compteurs:
528    
529         frugs = 0.         frugs = 0.
530         itap = 0         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
531         itaprad = 0              fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
532         CALL phyetat0(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
533              fsnow, falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollw, &              q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
534              dlw, radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, &              w01, ncid_startphy)
             zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &  
             run_off_lic_0, sig1, w01)  
535    
536         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
537         q2 = 1e-8         q2 = 1e-8
538    
539         radpas = NINT(86400. / dtphys / nbapp_rad)         lmt_pas = day_step / iphysiq
540           print *, 'Number of time steps of "physics" per day: ', lmt_pas
        ! on remet le calendrier a zero  
        IF (raz_date) itau_phy = 0  
   
        PRINT *, 'cycle_diurne = ', cycle_diurne  
        CALL printflag(radpas, ok_journe, ok_instan, ok_region)  
541    
542         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN         radpas = lmt_pas / nbapp_rad
543            print *, "Au minimum 4 appels par jour si cycle diurne"         print *, "radpas = ", radpas
           call abort_gcm('physiq', &  
                "Nombre d'appels au rayonnement insuffisant", 1)  
        ENDIF  
544    
545         ! Initialisation pour le sch\'ema de convection d'Emanuel :         ! Initialisation pour le sch\'ema de convection d'Emanuel :
546         IF (iflag_con >= 3) THEN         IF (conv_emanuel) THEN
547            ibas_con = 1            ibas_con = 1
548            itop_con = 1            itop_con = 1
549         ENDIF         ENDIF
# Line 668  contains Line 555  contains
555            rugoro = 0.            rugoro = 0.
556         ENDIF         ENDIF
557    
        lmt_pas = NINT(86400. / dtphys) ! tous les jours  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
558         ecrit_ins = NINT(ecrit_ins/dtphys)         ecrit_ins = NINT(ecrit_ins/dtphys)
559         ecrit_hf = NINT(ecrit_hf/dtphys)         ecrit_hf = NINT(ecrit_hf/dtphys)
560         ecrit_mth = NINT(ecrit_mth/dtphys)         ecrit_mth = NINT(ecrit_mth/dtphys)
# Line 679  contains Line 563  contains
563    
564         ! Initialisation des sorties         ! Initialisation des sorties
565    
566         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys)
567         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
568         ! Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
569         print *, 'physiq date0: ', date0         print *, 'physiq date0: ', date0
570           CALL phyredem0(lmt_pas)
571      ENDIF test_firstcal      ENDIF test_firstcal
572    
573      ! We will modify variables *_seri and we will not touch variables      ! We will modify variables *_seri and we will not touch variables
# Line 692  contains Line 577  contains
577      v_seri = v      v_seri = v
578      q_seri = qx(:, :, ivap)      q_seri = qx(:, :, ivap)
579      ql_seri = qx(:, :, iliq)      ql_seri = qx(:, :, iliq)
580      tr_seri = qx(:, :, 3: nqmx)      tr_seri = qx(:, :, 3:nqmx)
581    
582      ztsol = sum(ftsol * pctsrf, dim = 2)      ztsol = sum(ftsol * pctsrf, dim = 2)
583    
584      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
585         tit = 'after dynamics'         tit = 'after dynamics'
586         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, &
587              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
588         ! Comme les tendances de la physique sont ajout\'es dans la         ! Comme les tendances de la physique sont ajout\'es dans la
589         !  dynamique, la variation d'enthalpie par la dynamique devrait         ! dynamique, la variation d'enthalpie par la dynamique devrait
590         !  \^etre \'egale \`a la variation de la physique au pas de temps         ! \^etre \'egale \`a la variation de la physique au pas de temps
591         !  pr\'ec\'edent.  Donc la somme de ces 2 variations devrait \^etre         ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
592         !  nulle.         ! nulle.
593         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, &
594              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, &
595              d_qt, 0.)              d_qt, 0.)
# Line 738  contains Line 623  contains
623      ! Check temperatures:      ! Check temperatures:
624      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
625    
626      ! Incrémenter le compteur de la physique      call increment_itap
627      itap = itap + 1      julien = MOD(dayvrai, 360)
     julien = MOD(NINT(rdayvrai), 360)  
628      if (julien == 0) julien = 360      if (julien == 0) julien = 360
629    
630      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
631    
632      ! Prescrire l'ozone :      ! Prescrire l'ozone :
633      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
# Line 759  contains Line 643  contains
643      ENDDO      ENDDO
644      ql_seri = 0.      ql_seri = 0.
645    
646      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
647         tit = 'after reevap'         tit = 'after reevap'
648         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, &
649              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 770  contains Line 654  contains
654      frugs = MAX(frugs, 0.000015)      frugs = MAX(frugs, 0.000015)
655      zxrugs = sum(frugs * pctsrf, dim = 2)      zxrugs = sum(frugs * pctsrf, dim = 2)
656    
657      ! Calculs nécessaires au calcul de l'albedo dans l'interface      ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
658        ! la surface.
659    
660      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), longi, dist)
661      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
662         CALL zenang(zlongi, time, dtphys * REAL(radpas), rmu0, fract)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
663      ELSE      ELSE
664         rmu0 = -999.999         mu0 = - 999.999
665      ENDIF      ENDIF
666    
667      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
668      albsol = sum(falbe * pctsrf, dim = 2)      albsol = sum(falbe * pctsrf, dim = 2)
     albsollw = sum(falblw * pctsrf, dim = 2)  
669    
670      ! R\'epartition sous maille des flux longwave et shortwave      ! R\'epartition sous maille des flux longwave et shortwave
671      ! R\'epartition du longwave par sous-surface lin\'earis\'ee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
# Line 796  contains Line 680  contains
680    
681      ! Couche limite:      ! Couche limite:
682    
683      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &      CALL clmain(dtphys, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, v_seri, &
684           v_seri, julien, rmu0, co2_ppm, ftsol, cdmmax, cdhmax, &           julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, &
685           ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &           ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, &
686           fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &           rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, firstcal, &
687           fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &           agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, &
688           d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &           fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, &
689           q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, &           yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, &
690           capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
          fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab)  
691    
692      ! Incr\'ementation des flux      ! Incr\'ementation des flux
693    
# Line 837  contains Line 720  contains
720         ENDDO         ENDDO
721      ENDDO      ENDDO
722    
723      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
724         tit = 'after clmain'         tit = 'after clmain'
725         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, &
726              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 858  contains Line 741  contains
741         zxffonte(i) = 0.         zxffonte(i) = 0.
742         zxfqcalving(i) = 0.         zxfqcalving(i) = 0.
743    
744         s_pblh(i) = 0.         s_pblh(i) = 0.
745         s_lcl(i) = 0.         s_lcl(i) = 0.
746         s_capCL(i) = 0.         s_capCL(i) = 0.
747         s_oliqCL(i) = 0.         s_oliqCL(i) = 0.
748         s_cteiCL(i) = 0.         s_cteiCL(i) = 0.
# Line 868  contains Line 751  contains
751         s_trmb1(i) = 0.         s_trmb1(i) = 0.
752         s_trmb2(i) = 0.         s_trmb2(i) = 0.
753         s_trmb3(i) = 0.         s_trmb3(i) = 0.
   
        IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &  
             + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &  
             'physiq : probl\`eme sous surface au point ', i, &  
             pctsrf(i, 1 : nbsrf)  
754      ENDDO      ENDDO
755    
756        call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
757    
758      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
759         DO i = 1, klon         DO i = 1, klon
760            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
# Line 928  contains Line 809  contains
809      ! Calculer la dérive du flux infrarouge      ! Calculer la dérive du flux infrarouge
810    
811      DO i = 1, klon      DO i = 1, klon
812         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
     ENDDO  
   
     ! Appeler la convection (au choix)  
   
     DO k = 1, llm  
        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  
813      ENDDO      ENDDO
814    
815      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
816    
817      if (iflag_con == 2) then      ! Appeler la convection
        z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)  
        CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &  
             q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &  
             d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &  
             mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &  
             kdtop, pmflxr, pmflxs)  
        WHERE (rain_con < 0.) rain_con = 0.  
        WHERE (snow_con < 0.) snow_con = 0.  
        ibas_con = llm + 1 - kcbot  
        itop_con = llm + 1 - kctop  
     else  
        ! iflag_con >= 3  
818    
819        if (conv_emanuel) then
820         da = 0.         da = 0.
821         mp = 0.         mp = 0.
822         phi = 0.         phi = 0.
823         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &         CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
824              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &              d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
825              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &              upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, da, phi, mp)
826              qcondc, wd, pmflxr, pmflxs, da, phi, mp)         snow_con = 0.
827         clwcon0 = qcondc         clwcon0 = qcondc
828         mfu = upwd + dnwd         mfu = upwd + dnwd
        IF (.NOT. ok_gust) wd = 0.  
   
        ! Calcul des propri\'et\'es des nuages convectifs  
829    
830         DO k = 1, llm         IF (thermcep) THEN
831            DO i = 1, klon            zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
832               IF (thermcep) THEN            zqsat = zqsat / (1. - retv * zqsat)
833                  zdelta = MAX(0., SIGN(1., rtt - t_seri(i, k)))         ELSE
834                  zqsat(i, k) = r2es * FOEEW(t_seri(i, k), zdelta) / play(i, k)            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
835                  zqsat(i, k) = MIN(0.5, zqsat(i, k))         ENDIF
                 zqsat(i, k) = zqsat(i, k) / (1.-retv*zqsat(i, k))  
              ELSE  
                 IF (t_seri(i, k) < t_coup) THEN  
                    zqsat(i, k) = qsats(t_seri(i, k))/play(i, k)  
                 ELSE  
                    zqsat(i, k) = qsatl(t_seri(i, k))/play(i, k)  
                 ENDIF  
              ENDIF  
           ENDDO  
        ENDDO  
836    
837         ! calcul des proprietes des nuages convectifs         ! Properties of convective clouds
838         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
839         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
840              rnebcon0)              rnebcon0)
841    
842           forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
843         mfd = 0.         mfd = 0.
844         pen_u = 0.         pen_u = 0.
845         pen_d = 0.         pen_d = 0.
846         pde_d = 0.         pde_d = 0.
847         pde_u = 0.         pde_u = 0.
848        else
849           conv_q = d_q_dyn + d_q_vdf / dtphys
850           conv_t = d_t_dyn + d_t_vdf / dtphys
851           z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
852           CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
853                q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
854                d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
855                mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
856                kdtop, pmflxr, pmflxs)
857           WHERE (rain_con < 0.) rain_con = 0.
858           WHERE (snow_con < 0.) snow_con = 0.
859           ibas_con = llm + 1 - kcbot
860           itop_con = llm + 1 - kctop
861      END if      END if
862    
863      DO k = 1, llm      DO k = 1, llm
# Line 1007  contains Line 869  contains
869         ENDDO         ENDDO
870      ENDDO      ENDDO
871    
872      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
873         tit = 'after convect'         tit = 'after convect'
874         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, &
875              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 1029  contains Line 891  contains
891         print *, "Precip = ", zx_t         print *, "Precip = ", zx_t
892      ENDIF      ENDIF
893    
894      IF (iflag_con == 2) THEN      IF (.not. conv_emanuel) THEN
895         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
896         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
897         DO k = 1, llm         DO k = 1, llm
# Line 1056  contains Line 918  contains
918         t_seri = t_seri + d_t_ajs         t_seri = t_seri + d_t_ajs
919         q_seri = q_seri + d_q_ajs         q_seri = q_seri + d_q_ajs
920      else      else
        ! Thermiques  
921         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
922              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
923      endif      endif
924    
925      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
926         tit = 'after dry_adjust'         tit = 'after dry_adjust'
927         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, &
928              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 1088  contains Line 949  contains
949      do k = 1, llm      do k = 1, llm
950         do i = 1, klon         do i = 1, klon
951            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
952                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
953         enddo         enddo
954      enddo      enddo
955    
# Line 1135  contains Line 996  contains
996         print *, "Precip = ", zx_t         print *, "Precip = ", zx_t
997      ENDIF      ENDIF
998    
999      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1000         tit = 'after fisrt'         tit = 'after fisrt'
1001         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, &
1002              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 1147  contains Line 1008  contains
1008    
1009      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1010    
1011      IF (iflag_cldcon <= -1) THEN      IF (iflag_cldcon <= - 1) THEN
1012         ! seulement pour Tiedtke         ! seulement pour Tiedtke
1013         snow_tiedtke = 0.         snow_tiedtke = 0.
1014         if (iflag_cldcon == -1) then         if (iflag_cldcon == - 1) then
1015            rain_tiedtke = rain_con            rain_tiedtke = rain_con
1016         else         else
1017            rain_tiedtke = 0.            rain_tiedtke = 0.
1018            do k = 1, llm            do k = 1, llm
1019               do i = 1, klon               do i = 1, klon
1020                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1021                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1022                          *zmasse(i, k)                          *zmasse(i, k)
1023                  endif                  endif
1024               enddo               enddo
# Line 1225  contains Line 1086  contains
1086         DO i = 1, klon         DO i = 1, klon
1087            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1088            IF (thermcep) THEN            IF (thermcep) THEN
1089               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)  
1090               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1091               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1. - retv*zx_qs)
1092               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
1093            ELSE            ELSE
1094               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
# Line 1243  contains Line 1103  contains
1103      ENDDO      ENDDO
1104    
1105      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
1106      IF (ok_ade .OR. ok_aie) THEN      tau_ae = 0.
1107         ! Get sulfate aerosol distribution :      piz_ae = 0.
1108         CALL readsulfate(rdayvrai, firstcal, sulfate)      cg_ae = 0.
        CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)  
   
        CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &  
             aerindex)  
     ELSE  
        tau_ae = 0.  
        piz_ae = 0.  
        cg_ae = 0.  
     ENDIF  
1109    
1110      ! Param\`etres optiques des nuages et quelques param\`etres pour      ! Param\`etres optiques des nuages et quelques param\`etres pour
1111      ! diagnostics :      ! diagnostics :
# Line 1268  contains Line 1119  contains
1119              bl95_b1, cldtaupi, re, fl)              bl95_b1, cldtaupi, re, fl)
1120      endif      endif
1121    
1122      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
1123      IF (MOD(itaprad, radpas) == 0) THEN         ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1124         DO i = 1, klon         ! Calcul de l'abedo moyen par maille
1125            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &         albsol = sum(falbe * pctsrf, dim = 2)
1126                 + falbe(i, is_lic) * pctsrf(i, is_lic) &  
                + falbe(i, is_ter) * pctsrf(i, is_ter) &  
                + falbe(i, is_sic) * pctsrf(i, is_sic)  
           albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &  
                + falblw(i, is_lic) * pctsrf(i, is_lic) &  
                + falblw(i, is_ter) * pctsrf(i, is_ter) &  
                + falblw(i, is_sic) * pctsrf(i, is_sic)  
        ENDDO  
1127         ! Rayonnement (compatible Arpege-IFS) :         ! Rayonnement (compatible Arpege-IFS) :
1128         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1129              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1130              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1131              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1132              lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &              swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1133              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)              solswad, cldtaupi, topswai, solswai)
        itaprad = 0  
1134      ENDIF      ENDIF
     itaprad = itaprad + 1  
1135    
1136      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
1137    
1138      DO k = 1, llm      DO k = 1, llm
1139         DO i = 1, klon         DO i = 1, klon
1140            t_seri(i, k) = t_seri(i, k) + (heat(i, k)-cool(i, k)) * dtphys/86400.            t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.
1141         ENDDO         ENDDO
1142      ENDDO      ENDDO
1143    
1144      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1145         tit = 'after rad'         tit = 'after rad'
1146         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, &
1147              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
# Line 1328  contains Line 1170  contains
1170      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1171    
1172      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1173         ! selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
1174         igwd = 0         igwd = 0
1175         DO i = 1, klon         DO i = 1, klon
1176            itest(i) = 0            itest(i) = 0
1177            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1178               itest(i) = 1               itest(i) = 1
1179               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1180            ENDIF            ENDIF
1181         ENDDO         ENDDO
1182    
1183         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1184              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1185              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1186    
1187         ! ajout des tendances         ! ajout des tendances
1188         DO k = 1, llm         DO k = 1, llm
# Line 1358  contains Line 1199  contains
1199         igwd = 0         igwd = 0
1200         DO i = 1, klon         DO i = 1, klon
1201            itest(i) = 0            itest(i) = 0
1202            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
1203               itest(i) = 1               itest(i) = 1
1204               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1205            ENDIF            ENDIF
1206         ENDDO         ENDDO
1207    
# Line 1394  contains Line 1234  contains
1234         ENDDO         ENDDO
1235      ENDDO      ENDDO
1236    
1237      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1238           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1239    
1240      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1241           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1242           d_qt, d_ec)           d_qt, d_ec)
1243    
1244      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1245      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, u, t, &      call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &
1246           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &           play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &
1247           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, albsol, rhcl, &           yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
1248           cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, &           tr_seri, zmasse, ncid_startphy)
1249           mp, upwd, dnwd, tr_seri, zmasse)  
1250        IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1251      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &           pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1252           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &           frac_impa, frac_nucl, pphis, airephy, dtphys)
          pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)  
1253    
1254      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1255      CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &      CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
          ue, uq)  
1256    
1257      ! diag. bilKP      ! diag. bilKP
1258    
1259      CALL transp_lay(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &      CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
1260           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1261    
1262      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
1263    
1264      ! conversion Ec -> E thermique      ! conversion Ec en énergie thermique
1265      DO k = 1, llm      DO k = 1, llm
1266         DO i = 1, klon         DO i = 1, klon
1267            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
# Line 1434  contains Line 1272  contains
1272         END DO         END DO
1273      END DO      END DO
1274    
1275      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1276         tit = 'after physic'         tit = 'after physic'
1277         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, &
1278              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1279         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1280         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1281         ! 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.
1282         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
# Line 1472  contains Line 1310  contains
1310      DO iq = 3, nqmx      DO iq = 3, nqmx
1311         DO k = 1, llm         DO k = 1, llm
1312            DO i = 1, klon            DO i = 1, klon
1313               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
1314            ENDDO            ENDDO
1315         ENDDO         ENDDO
1316      ENDDO      ENDDO
# Line 1485  contains Line 1323  contains
1323         ENDDO         ENDDO
1324      ENDDO      ENDDO
1325    
1326      ! Ecriture des sorties      CALL histwrite_phy("phis", pphis)
1327      call write_histins      CALL histwrite_phy("aire", airephy)
1328        CALL histwrite_phy("psol", paprs(:, 1))
1329        CALL histwrite_phy("precip", rain_fall + snow_fall)
1330        CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1331        CALL histwrite_phy("pluc", rain_con + snow_con)
1332        CALL histwrite_phy("tsol", zxtsol)
1333        CALL histwrite_phy("t2m", zt2m)
1334        CALL histwrite_phy("q2m", zq2m)
1335        CALL histwrite_phy("u10m", zu10m)
1336        CALL histwrite_phy("v10m", zv10m)
1337        CALL histwrite_phy("snow", snow_fall)
1338        CALL histwrite_phy("cdrm", cdragm)
1339        CALL histwrite_phy("cdrh", cdragh)
1340        CALL histwrite_phy("topl", toplw)
1341        CALL histwrite_phy("evap", evap)
1342        CALL histwrite_phy("sols", solsw)
1343        CALL histwrite_phy("soll", sollw)
1344        CALL histwrite_phy("solldown", sollwdown)
1345        CALL histwrite_phy("bils", bils)
1346        CALL histwrite_phy("sens", - sens)
1347        CALL histwrite_phy("fder", fder)
1348        CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1349        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1350        CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1351        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1352    
1353      ! Si c'est la fin, il faut conserver l'etat de redemarrage      DO nsrf = 1, nbsrf
1354      IF (lafin) THEN         CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)
1355         itau_phy = itau_phy + itap         CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1356         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &         CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1357              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &         CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1358              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &         CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1359              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &         CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1360              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)         CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1361      ENDIF         CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1362           CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1363      firstcal = .FALSE.      END DO
   
   contains  
   
     subroutine write_histins  
   
       ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09  
   
       use dimens_m, only: iim, jjm  
       USE histsync_m, ONLY: histsync  
       USE histwrite_m, ONLY: histwrite  
   
       real zout  
       integer itau_w ! pas de temps ecriture  
       REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)  
   
       !--------------------------------------------------  
   
       IF (ok_instan) THEN  
          ! Champs 2D:  
   
          zsto = dtphys * ecrit_ins  
          zout = dtphys * ecrit_ins  
          itau_w = itau_phy + itap  
   
          i = NINT(zout/zsto)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)  
          CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)  
   
          i = NINT(zout/zsto)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)  
          CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = paprs(i, 1)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)  
          !ccIM  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)  
          CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)  
          CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)  
          CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)  
          CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)  
   
          zx_tmp_fi2d(1:klon) = -1*sens(1:klon)  
          ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)  
          CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)  
   
          DO nsrf = 1, nbsrf  
             !XXX  
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
          END DO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)  
          CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)  
   
          !HBTM2  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)  
   
          ! Champs 3D:  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)  
          CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)  
          CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)  
1364    
1365           call histsync(nid_ins)      CALL histwrite_phy("albs", albsol)
1366        ENDIF      CALL histwrite_phy("rugs", zxrugs)
1367        CALL histwrite_phy("s_pblh", s_pblh)
1368        CALL histwrite_phy("s_pblt", s_pblt)
1369        CALL histwrite_phy("s_lcl", s_lcl)
1370        CALL histwrite_phy("s_capCL", s_capCL)
1371        CALL histwrite_phy("s_oliqCL", s_oliqCL)
1372        CALL histwrite_phy("s_cteiCL", s_cteiCL)
1373        CALL histwrite_phy("s_therm", s_therm)
1374        CALL histwrite_phy("s_trmb1", s_trmb1)
1375        CALL histwrite_phy("s_trmb2", s_trmb2)
1376        CALL histwrite_phy("s_trmb3", s_trmb3)
1377        if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1378        CALL histwrite_phy("temp", t_seri)
1379        CALL histwrite_phy("vitu", u_seri)
1380        CALL histwrite_phy("vitv", v_seri)
1381        CALL histwrite_phy("geop", zphi)
1382        CALL histwrite_phy("pres", play)
1383        CALL histwrite_phy("dtvdf", d_t_vdf)
1384        CALL histwrite_phy("dqvdf", d_q_vdf)
1385        CALL histwrite_phy("rhum", zx_rh)
1386    
1387        if (ok_instan) call histsync(nid_ins)
1388    
1389        IF (lafin) then
1390           call NF95_CLOSE(ncid_startphy)
1391           CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1392                fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1393                radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1394                t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1395                w01)
1396        end IF
1397    
1398      end subroutine write_histins      firstcal = .FALSE.
1399    
1400    END SUBROUTINE physiq    END SUBROUTINE physiq
1401    

Legend:
Removed from v.101  
changed lines
  Added in v.200

  ViewVC Help
Powered by ViewVC 1.1.21