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

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

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

trunk/phylmd/physiq.f revision 92 by guez, Wed Mar 26 18:16:05 2014 UTC trunk/Sources/phylmd/physiq.f revision 178 by guez, Fri Mar 11 18:47:26 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    
# Line 19  contains Line 19  contains
19      use aeropt_m, only: aeropt      use aeropt_m, only: aeropt
20      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
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, ecrit_hf, ecrit_ins, ecrit_mth, &
23           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin           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      use clouds_gno_m, only: clouds_gno
28      USE comgeomphy, ONLY: airephy, cuphy, cvphy      use comconst, only: dtphys
29        USE comgeomphy, ONLY: airephy
30      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
31      USE conf_gcm_m, ONLY: offline, raz_date      USE conf_gcm_m, ONLY: offline, raz_date, day_step, iphysiq
32      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
33      use conflx_m, only: conflx      use conflx_m, only: conflx
34      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
# Line 35  contains Line 36  contains
36      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
37      use diagphy_m, only: diagphy      use diagphy_m, only: diagphy
38      USE dimens_m, ONLY: llm, nqmx      USE dimens_m, ONLY: llm, nqmx
39      USE dimphy, ONLY: klon, nbtr      USE dimphy, ONLY: klon
40      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
41      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
42        use dynetat0_m, only: day_ref, annee_ref
43      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
44      use fisrtilp_m, only: fisrtilp      use fisrtilp_m, only: fisrtilp
45      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
46      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, &
47           nbsrf           nbsrf
48      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins
49        use netcdf95, only: NF95_CLOSE
50      use newmicro_m, only: newmicro      use newmicro_m, only: newmicro
51      USE oasis_m, ONLY: ok_oasis      use nuage_m, only: nuage
52      USE orbite_m, ONLY: orbite, zenang      USE orbite_m, ONLY: orbite
53      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
54      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
55      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
56        USE phyredem0_m, ONLY: phyredem0
57      USE phystokenc_m, ONLY: phystokenc      USE phystokenc_m, ONLY: phystokenc
58      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
59      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
60      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
61      use readsulfate_m, only: readsulfate      use readsulfate_m, only: readsulfate
62      use sugwd_m, only: sugwd      use readsulfate_preind_m, only: readsulfate_preind
63      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      use yoegwd, only: sugwd
64      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt
65        use transp_m, only: transp
66        use transp_lay_m, only: transp_lay
67      use unit_nml_m, only: unit_nml      use unit_nml_m, only: unit_nml
68      USE ymds2ju_m, ONLY: ymds2ju      USE ymds2ju_m, ONLY: ymds2ju
69      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
70        use zenang_m, only: zenang
71    
72      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
73    
74      REAL, intent(in):: rdayvrai      integer, intent(in):: dayvrai
75      ! (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
76    
77      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)  
78    
79      REAL, intent(in):: paprs(klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
80      ! (pression pour chaque inter-couche, en Pa)      ! pression pour chaque inter-couche, en Pa
81    
82      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
83      ! (input pression pour le mileu de chaque couche (en Pa))      ! pression pour le mileu de chaque couche (en Pa)
84    
85      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
86      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! géopotentiel de chaque couche (référence sol)
87    
88      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
89    
90      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
91      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m/s
92    
93      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
94      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
95    
96      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
97      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
98    
99      REAL, intent(in):: omega(klon, llm) ! vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
100      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m s-2)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
101      REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m s-2)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
102      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)
103      REAL, intent(out):: d_qx(klon, llm, nqmx) ! tendance physique de "qx" (s-1)  
104        REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
105        ! tendance physique de "qx" (s-1)
106    
107      ! Local:      ! Local:
108    
109      LOGICAL:: firstcal = .true.      LOGICAL:: firstcal = .true.
110    
     INTEGER nbteta  
     PARAMETER(nbteta = 3)  
   
111      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
112      PARAMETER (ok_gust = .FALSE.)      PARAMETER (ok_gust = .FALSE.)
113    
114      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL, PARAMETER:: check = .FALSE.
115      PARAMETER (check = .FALSE.)      ! Verifier la conservation du modele en eau
116    
117      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
118      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
119    
     ! Parametres lies au coupleur OASIS:  
     INTEGER, SAVE:: npas, nexca  
     logical rnpb  
     parameter(rnpb = .true.)  
   
     character(len = 6):: ocean = 'force '  
     ! (type de mod\`ele oc\'ean \`a utiliser: "force" ou "slab" mais  
     ! pas "couple")  
   
     ! "slab" ocean  
     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  
   
     ! 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 142  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 155  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    
     ! Amip2 PV a theta constante  
   
     CHARACTER(LEN = 3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
   
     ! Amip2 PV a theta constante  
   
143      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
144      REAL swup0(klon, llm + 1), swup(klon, llm + 1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
145      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
# Line 177  contains Line 153  contains
153    
154      integer nlevSTD      integer nlevSTD
155      PARAMETER(nlevSTD = 17)      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 '/  
156    
157      ! prw: precipitable water      ! prw: precipitable water
158      real prw(klon)      real prw(klon)
# Line 197  contains Line 165  contains
165      INTEGER kmax, lmax      INTEGER kmax, lmax
166      PARAMETER(kmax = 8, lmax = 8)      PARAMETER(kmax = 8, lmax = 8)
167      INTEGER kmaxm1, lmaxm1      INTEGER kmaxm1, lmaxm1
168      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)      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  
   
     integer nid_hf, nid_hf3d  
     save nid_hf, nid_hf3d  
169    
170      ! Variables propres a la physique      ! Variables propres a la physique
171    
172      INTEGER, save:: radpas      INTEGER, save:: radpas
173      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
174      ! "physiq".)      ! "physiq".
175    
176      REAL radsol(klon)      REAL radsol(klon)
177      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
178    
179      INTEGER, SAVE:: itap ! number of calls to "physiq"      INTEGER:: itap = 0 ! number of calls to "physiq"
180    
181      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
182    
# Line 266  contains Line 187  contains
187      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
188      SAVE fluxlat      SAVE fluxlat
189    
190      REAL fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
191      SAVE fqsurf ! humidite de l'air au contact de la surface      ! humidite de l'air au contact de la surface
192    
193      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol      REAL, save:: qsol(klon)
194        ! column-density of water in soil, in kg m-2
195    
196      REAL fsnow(klon, nbsrf)      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
197      SAVE fsnow ! epaisseur neigeuse      REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
   
     REAL falbe(klon, nbsrf)  
     SAVE falbe ! albedo par type de surface  
     REAL falblw(klon, nbsrf)  
     SAVE falblw ! albedo par type de surface  
198    
199      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
200      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
# Line 288  contains Line 205  contains
205      REAL, save:: zpic(klon) ! Maximum de l'OESM      REAL, save:: zpic(klon) ! Maximum de l'OESM
206      REAL, save:: zval(klon) ! Minimum de l'OESM      REAL, save:: zval(klon) ! Minimum de l'OESM
207      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
   
208      REAL zulow(klon), zvlow(klon)      REAL zulow(klon), zvlow(klon)
209        INTEGER igwd, itest(klon)
     INTEGER igwd, idx(klon), itest(klon)  
210    
211      REAL agesno(klon, nbsrf)      REAL agesno(klon, nbsrf)
212      SAVE agesno ! age de la neige      SAVE agesno ! age de la neige
# Line 301  contains Line 216  contains
216      !KE43      !KE43
217      ! Variables liees a la convection de K. Emanuel (sb):      ! Variables liees a la convection de K. Emanuel (sb):
218    
     REAL bas, top ! cloud base and top levels  
     SAVE bas  
     SAVE top  
   
219      REAL Ma(klon, llm) ! undilute upward mass flux      REAL Ma(klon, llm) ! undilute upward mass flux
220      SAVE Ma      SAVE Ma
221      REAL qcondc(klon, llm) ! in-cld water content from convect      REAL qcondc(klon, llm) ! in-cld water content from convect
# Line 312  contains Line 223  contains
223      REAL, save:: sig1(klon, llm), w01(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
224      REAL, save:: wd(klon)      REAL, save:: wd(klon)
225    
226      ! Variables locales pour la couche limite (al1):      ! Variables pour la couche limite (al1):
   
     ! Variables locales:  
227    
228      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
229      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
# Line 338  contains Line 247  contains
247      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
248      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
249    
250      REAL, save:: rain_fall(klon) ! pluie      REAL, save:: rain_fall(klon)
251      REAL, save:: snow_fall(klon) ! neige      ! liquid water mass flux (kg/m2/s), positive down
252    
253        REAL, save:: snow_fall(klon)
254        ! solid water mass flux (kg/m2/s), positive down
255    
256      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
257    
# Line 348  contains Line 260  contains
260      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
261      SAVE dlw      SAVE dlw
262      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
263      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
     save fder  
264      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
265      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
266      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
267      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
268    
269      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
270      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
271    
272      ! Conditions aux limites      ! Conditions aux limites
273    
274      INTEGER julien      INTEGER julien
   
275      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
276      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
277      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
278        REAL, save:: albsol(klon) ! albedo du sol total visible
     REAL albsol(klon)  
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
279      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
280    
     ! Declaration des procedures appelees  
   
     EXTERNAL alboc ! calculer l'albedo sur ocean  
     !KE43  
     EXTERNAL conema3 ! convect4.3  
     EXTERNAL nuage ! calculer les proprietes radiatives  
     EXTERNAL transp ! transport total de l'eau et de l'energie  
   
     ! Variables locales  
   
281      real, save:: clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
282      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
283    
# Line 408  contains Line 302  contains
302      ! 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
303      ! les variables soient r\'emanentes.      ! les variables soient r\'emanentes.
304      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
305      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
306      REAL, save:: cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
307      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
308      REAL, save:: topsw(klon), toplw(klon), solsw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
309      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
310      real, save:: sollwdown(klon) ! downward LW flux at surface      real, save:: sollwdown(klon) ! downward LW flux at surface
311      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
312      REAL albpla(klon)      REAL, save:: albpla(klon)
313      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
314      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
     SAVE albpla  
     SAVE heat0, cool0  
   
     INTEGER itaprad  
     SAVE itaprad  
315    
316      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
317      REAL conv_t(klon, llm) ! convergence of temperature (K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K/s)
# Line 432  contains Line 321  contains
321    
322      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
323    
324      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
325      REAL zdtime ! pas de temps du rayonnement (s)      real longi
     real zlongi  
326      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
327      REAL za, zb      REAL za, zb
328      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
329      real zqsat(klon, llm)      real zqsat(klon, llm)
330      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
331      REAL, PARAMETER:: t_coup = 234.      REAL, PARAMETER:: t_coup = 234.
332      REAL zphi(klon, llm)      REAL zphi(klon, llm)
333    
334      ! cf. AM Variables locales pour la CLA (hbtm2)      ! cf. AM Variables pour la CLA (hbtm2)
335    
336      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
337      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
# Line 461  contains Line 349  contains
349      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
350      REAL s_trmb3(klon)      REAL s_trmb3(klon)
351    
352      ! Variables locales pour la convection de K. Emanuel :      ! Variables pour la convection de K. Emanuel :
353    
354      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
355      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
356      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  
357      REAL cape(klon) ! CAPE      REAL cape(klon) ! CAPE
358      SAVE cape      SAVE cape
359    
     REAL pbase(klon) ! cloud base pressure  
     SAVE pbase  
     REAL bbase(klon) ! cloud base buoyancy  
     SAVE bbase  
     REAL rflag(klon) ! flag fonctionnement de convect  
360      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)  
361    
362      ! Variables du changement      ! Variables du changement
363    
# Line 528  contains Line 407  contains
407      integer:: iflag_cldcon = 1      integer:: iflag_cldcon = 1
408      logical ptconv(klon, llm)      logical ptconv(klon, llm)
409    
410      ! Variables locales pour effectuer les appels en s\'erie :      ! Variables pour effectuer les appels en s\'erie :
411    
412      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
413      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
414      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
415        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
416    
417      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
418    
# Line 546  contains Line 423  contains
423    
424      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
425    
426      INTEGER, SAVE:: nid_day, nid_ins      INTEGER, SAVE:: nid_ins
427    
428      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.
429      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.
430      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.
431      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.
432    
     REAL zsto  
433      real date0      real date0
434    
435      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
436      REAL ztsol(klon)      REAL ztsol(klon)
437      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
438      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
439      REAL zero_v(klon)      REAL zero_v(klon)
440      CHARACTER(LEN = 15) tit      CHARACTER(LEN = 20) tit
441      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
442      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
443    
# Line 579  contains Line 454  contains
454      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
455    
456      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
457      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value      ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value
458    
459      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
460      ! cloud optical thickness for pre-industrial (pi) aerosols      ! cloud optical thickness for pre-industrial (pi) aerosols
# Line 623  contains Line 498  contains
498      ! (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)
499    
500      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
501        integer, save:: ncid_startphy, itau_phy
502    
503      namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &      namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
504           fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &           facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
505           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  
506    
507      !----------------------------------------------------------------      !----------------------------------------------------------------
508    
509      IF (if_ebil >= 1) zero_v = 0.      IF (if_ebil >= 1) zero_v = 0.
510      IF (nqmx < 2) CALL abort_gcm('physiq', &      IF (nqmx < 2) CALL abort_gcm('physiq', &
511           'eaux vapeur et liquide sont indispensables', 1)           'eaux vapeur et liquide sont indispensables')
512    
513      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
514         ! initialiser         ! initialiser
# Line 646  contains Line 521  contains
521         piz_ae = 0.         piz_ae = 0.
522         tau_ae = 0.         tau_ae = 0.
523         cg_ae = 0.         cg_ae = 0.
524         rain_con(:) = 0.         rain_con = 0.
525         snow_con(:) = 0.         snow_con = 0.
526         topswai(:) = 0.         topswai = 0.
527         topswad(:) = 0.         topswad = 0.
528         solswai(:) = 0.         solswai = 0.
529         solswad(:) = 0.         solswad = 0.
530    
531         d_u_con = 0.         d_u_con = 0.
532         d_v_con = 0.         d_v_con = 0.
# Line 684  contains Line 559  contains
559         ! Initialiser les compteurs:         ! Initialiser les compteurs:
560    
561         frugs = 0.         frugs = 0.
562         itap = 0         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
563         itaprad = 0              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
564         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
565              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
566              snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &              run_off_lic_0, sig1, w01, ncid_startphy, itau_phy)
             zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &  
             ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)  
567    
568         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
569         q2 = 1e-8         q2 = 1e-8
570    
571         radpas = NINT(86400. / dtphys / nbapp_rad)         lmt_pas = day_step / iphysiq
572           print *, 'Number of time steps of "physics" per day: ', lmt_pas
573    
574           radpas = lmt_pas / nbapp_rad
575    
576         ! on remet le calendrier a zero         ! On remet le calendrier a zero
577         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
578    
579         PRINT *, 'cycle_diurne = ', cycle_diurne         CALL printflag(radpas, ok_journe, ok_instan, ok_region)
        CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &  
             ok_instan, ok_region)  
   
        IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN  
           print *, "Au minimum 4 appels par jour si cycle diurne"  
           call abort_gcm('physiq', &  
                "Nombre d'appels au rayonnement insuffisant", 1)  
        ENDIF  
580    
581         ! Initialisation pour le sch\'ema de convection d'Emanuel :         ! Initialisation pour le sch\'ema de convection d'Emanuel :
582         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
# Line 723  contains Line 591  contains
591            rugoro = 0.            rugoro = 0.
592         ENDIF         ENDIF
593    
        lmt_pas = NINT(86400. / dtphys) ! tous les jours  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
594         ecrit_ins = NINT(ecrit_ins/dtphys)         ecrit_ins = NINT(ecrit_ins/dtphys)
595         ecrit_hf = NINT(ecrit_hf/dtphys)         ecrit_hf = NINT(ecrit_hf/dtphys)
596         ecrit_mth = NINT(ecrit_mth/dtphys)         ecrit_mth = NINT(ecrit_mth/dtphys)
597         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
598         ecrit_reg = NINT(ecrit_reg/dtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
599    
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
600         ! Initialisation des sorties         ! Initialisation des sorties
601    
602         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins, itau_phy)
603         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
604         ! Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
605         print *, 'physiq date0: ', date0         print *, 'physiq date0: ', date0
606           CALL phyredem0(lmt_pas, itau_phy)
607      ENDIF test_firstcal      ENDIF test_firstcal
608    
     ! Mettre a zero des variables de sortie (pour securite)  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
609      ! We will modify variables *_seri and we will not touch variables      ! We will modify variables *_seri and we will not touch variables
610      ! u, v, h, q:      ! u, v, t, qx:
611      DO k = 1, llm      t_seri = t
612         DO i = 1, klon      u_seri = u
613            t_seri(i, k) = t(i, k)      v_seri = v
614            u_seri(i, k) = u(i, k)      q_seri = qx(:, :, ivap)
615            v_seri(i, k) = v(i, k)      ql_seri = qx(:, :, iliq)
616            q_seri(i, k) = qx(i, k, ivap)      tr_seri = qx(:, :, 3:nqmx)
           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  
617    
618      DO i = 1, klon      ztsol = sum(ftsol * pctsrf, dim = 2)
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
619    
620      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
621         tit = 'after dynamics'         tit = 'after dynamics'
622         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, &
623              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)  
624         ! Comme les tendances de la physique sont ajout\'es dans la         ! Comme les tendances de la physique sont ajout\'es dans la
625         !  dynamique, la variation d'enthalpie par la dynamique devrait         !  dynamique, la variation d'enthalpie par la dynamique devrait
626         !  \^etre \'egale \`a la variation de la physique au pas de temps         !  \^etre \'egale \`a la variation de la physique au pas de temps
# Line 789  contains Line 628  contains
628         !  nulle.         !  nulle.
629         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, &
630              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, &
631              d_qt, 0., fs_bound, fq_bound)              d_qt, 0.)
632      END IF      END IF
633    
634      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
# Line 820  contains Line 659  contains
659      ! Check temperatures:      ! Check temperatures:
660      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
661    
662      ! Incrementer le compteur de la physique      ! Incrémenter le compteur de la physique
663      itap = itap + 1      itap = itap + 1
664      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(dayvrai, 360)
665      if (julien == 0) julien = 360      if (julien == 0) julien = 360
666    
667      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
   
     ! Mettre en action les conditions aux limites (albedo, sst etc.).  
668    
669      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone :
670      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
671    
672      ! \'Evaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
# Line 846  contains Line 683  contains
683      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
684         tit = 'after reevap'         tit = 'after reevap'
685         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, &
686              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)  
687         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, &
688              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)  
   
689      END IF      END IF
690    
691      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
692        zxrugs = sum(frugs * pctsrf, dim = 2)
     DO i = 1, klon  
        zxrugs(i) = 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  
693    
694      ! calculs necessaires au calcul de l'albedo dans l'interface      ! Calculs nécessaires au calcul de l'albedo dans l'interface avec
695        ! la surface.
696    
697      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), longi, dist)
698      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
699         zdtime = dtphys * REAL(radpas)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
700      ELSE      ELSE
701         rmu0 = -999.999         mu0 = - 999.999
702      ENDIF      ENDIF
703    
704      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
705      albsol(:) = 0.      albsol = sum(falbe * pctsrf, dim = 2)
     albsollw(:) = 0.  
     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  
706    
707      ! R\'epartition sous maille des flux longwave et shortwave      ! R\'epartition sous maille des flux longwave et shortwave
708      ! R\'epartition du longwave par sous-surface lin\'earis\'ee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
709    
710      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
711         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
712            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
713                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
714            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))      END forall
        ENDDO  
     ENDDO  
715    
716      fder = dlw      fder = dlw
717    
718      ! Couche limite:      ! Couche limite:
719    
720      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
721           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &           v_seri, julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, &
722           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, &
723           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           fluxlat, rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, &
724           rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &           firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &
725           frugs, firstcal, agesno, rugoro, d_t_vdf, &           fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, &
726           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, &
727           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           pblT, therm, trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, &
728           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           run_off_lic_0)
          fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)  
729    
730      ! Incr\'ementation des flux      ! Incr\'ementation des flux
731    
# Line 950  contains Line 761  contains
761      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
762         tit = 'after clmain'         tit = 'after clmain'
763         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, &
764              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)  
765         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, &
766              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)  
767      END IF      END IF
768    
769      ! Update surface temperature:      ! Update surface temperature:
# Line 983  contains Line 792  contains
792    
793         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) &
794              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
795              'physiq : probl\`eme sous surface au point ', i, pctsrf(i, 1 : nbsrf)              'physiq : probl\`eme sous surface au point ', i, &
796                pctsrf(i, 1 : nbsrf)
797      ENDDO      ENDDO
798      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
799         DO i = 1, klon         DO i = 1, klon
# Line 1011  contains Line 821  contains
821         ENDDO         ENDDO
822      ENDDO      ENDDO
823    
824      ! 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 :
   
825      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
826         DO i = 1, klon         DO i = 1, klon
827            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
# Line 1037  contains Line 846  contains
846         ENDDO         ENDDO
847      ENDDO      ENDDO
848    
849      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
850    
851      DO i = 1, klon      DO i = 1, klon
852         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
853      ENDDO      ENDDO
854    
855      ! Appeler la convection (au choix)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
   
     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  
     ENDDO  
856    
857      IF (check) THEN      ! Appeler la convection (au choix)
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
858    
859      if (iflag_con == 2) then      if (iflag_con == 2) then
860           conv_q = d_q_dyn + d_q_vdf / dtphys
861           conv_t = d_t_dyn + d_t_vdf / dtphys
862         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
863         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
864              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &              q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
865              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
866              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &              mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
867              kdtop, pmflxr, pmflxs)              kdtop, pmflxr, pmflxs)
868         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
869         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
# Line 1071  contains Line 872  contains
872      else      else
873         ! iflag_con >= 3         ! iflag_con >= 3
874    
875         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &         da = 0.
876              v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &         mp = 0.
877              d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &         phi = 0.
878              itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
879              pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &
880              wd, pmflxr, pmflxs, da, phi, mp, ntra=1)              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &
881         ! (number of tracers for the convection scheme of Kerry Emanuel:              qcondc, wd, pmflxr, pmflxs, da, phi, mp)
        ! la partie traceurs est faite dans phytrac  
        ! on met ntra = 1 pour limiter les appels mais on peut  
        ! supprimer les calculs / ftra.)  
   
882         clwcon0 = qcondc         clwcon0 = qcondc
883         mfu = upwd + dnwd         mfu = upwd + dnwd
884         IF (.NOT. ok_gust) wd = 0.         IF (.NOT. ok_gust) wd = 0.
885    
886         ! Calcul des propri\'et\'es des nuages convectifs         IF (thermcep) THEN
887              zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
888         DO k = 1, llm            zqsat = zqsat / (1. - retv * zqsat)
889            DO i = 1, klon         ELSE
890               zx_t = t_seri(i, k)            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
891               IF (thermcep) THEN         ENDIF
                 zdelta = MAX(0., SIGN(1., rtt-zx_t))  
                 zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)  
                 zx_qs = MIN(0.5, zx_qs)  
                 zcor = 1./(1.-retv*zx_qs)  
                 zx_qs = zx_qs*zcor  
              ELSE  
                 IF (zx_t < t_coup) THEN  
                    zx_qs = qsats(zx_t)/play(i, k)  
                 ELSE  
                    zx_qs = qsatl(zx_t)/play(i, k)  
                 ENDIF  
              ENDIF  
              zqsat(i, k) = zx_qs  
           ENDDO  
        ENDDO  
892    
893         ! calcul des proprietes des nuages convectifs         ! Properties of convective clouds
894         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
895         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
896              rnebcon0)              rnebcon0)
# Line 1132  contains Line 914  contains
914      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
915         tit = 'after convect'         tit = 'after convect'
916         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, &
917              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)  
918         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, &
919              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)  
920      END IF      END IF
921    
922      IF (check) THEN      IF (check) THEN
923         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
924         print *, "aprescon = ", za         print *, "aprescon = ", za
925         zx_t = 0.         zx_t = 0.
926         za = 0.         za = 0.
# Line 1188  contains Line 968  contains
968      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
969         tit = 'after dry_adjust'         tit = 'after dry_adjust'
970         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, &
971              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)  
972      END IF      END IF
973    
974      ! Caclul des ratqs      ! Caclul des ratqs
# Line 1247  contains Line 1026  contains
1026         ENDDO         ENDDO
1027      ENDDO      ENDDO
1028      IF (check) THEN      IF (check) THEN
1029         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1030         print *, "apresilp = ", za         print *, "apresilp = ", za
1031         zx_t = 0.         zx_t = 0.
1032         za = 0.         za = 0.
# Line 1263  contains Line 1042  contains
1042      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1043         tit = 'after fisrt'         tit = 'after fisrt'
1044         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, &
1045              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)  
1046         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, &
1047              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)  
1048      END IF      END IF
1049    
1050      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1051    
1052      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1053    
1054      IF (iflag_cldcon <= -1) THEN      IF (iflag_cldcon <= - 1) THEN
1055         ! seulement pour Tiedtke         ! seulement pour Tiedtke
1056         snow_tiedtke = 0.         snow_tiedtke = 0.
1057         if (iflag_cldcon == -1) then         if (iflag_cldcon == - 1) then
1058            rain_tiedtke = rain_con            rain_tiedtke = rain_con
1059         else         else
1060            rain_tiedtke = 0.            rain_tiedtke = 0.
1061            do k = 1, llm            do k = 1, llm
1062               do i = 1, klon               do i = 1, klon
1063                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1064                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1065                          *zmasse(i, k)                          *zmasse(i, k)
1066                  endif                  endif
1067               enddo               enddo
# Line 1344  contains Line 1121  contains
1121      ENDDO      ENDDO
1122    
1123      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, &
1124           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, &
1125           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1126    
1127      ! Humidit\'e relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
1128      DO k = 1, llm      DO k = 1, llm
1129         DO i = 1, klon         DO i = 1, klon
1130            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1131            IF (thermcep) THEN            IF (thermcep) THEN
1132               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)  
1133               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1134               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1. - retv*zx_qs)
1135               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
1136            ELSE            ELSE
1137               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
# Line 1372  contains Line 1148  contains
1148      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
1149      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1150         ! Get sulfate aerosol distribution :         ! Get sulfate aerosol distribution :
1151         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(dayvrai, time, firstcal, sulfate)
1152         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1153    
1154         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1155              aerindex)              aerindex)
# Line 1383  contains Line 1159  contains
1159         cg_ae = 0.         cg_ae = 0.
1160      ENDIF      ENDIF
1161    
1162      ! Param\`etres optiques des nuages et quelques param\`etres pour diagnostics :      ! Param\`etres optiques des nuages et quelques param\`etres pour
1163        ! diagnostics :
1164      if (ok_newmicro) then      if (ok_newmicro) then
1165         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1166              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
# Line 1394  contains Line 1171  contains
1171              bl95_b1, cldtaupi, re, fl)              bl95_b1, cldtaupi, re, fl)
1172      endif      endif
1173    
1174      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
1175      IF (MOD(itaprad, radpas) == 0) THEN         ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1176         DO i = 1, klon         ! Calcul de l'abedo moyen par maille
1177            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &         albsol = sum(falbe * pctsrf, dim = 2)
1178                 + 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  
1179         ! Rayonnement (compatible Arpege-IFS) :         ! Rayonnement (compatible Arpege-IFS) :
1180         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1181              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1182              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1183              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1184              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, &
1185              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)              solswad, cldtaupi, topswai, solswai)
        itaprad = 0  
1186      ENDIF      ENDIF
     itaprad = itaprad + 1  
1187    
1188      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
1189    
1190      DO k = 1, llm      DO k = 1, llm
1191         DO i = 1, klon         DO i = 1, klon
1192            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.
1193         ENDDO         ENDDO
1194      ENDDO      ENDDO
1195    
1196      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1197         tit = 'after rad'         tit = 'after rad'
1198         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, &
1199              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)  
1200         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1201              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)  
1202      END IF      END IF
1203    
1204      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
# Line 1456  contains Line 1222  contains
1222      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1223    
1224      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1225         ! selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
1226         igwd = 0         igwd = 0
1227         DO i = 1, klon         DO i = 1, klon
1228            itest(i) = 0            itest(i) = 0
1229            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1230               itest(i) = 1               itest(i) = 1
1231               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1232            ENDIF            ENDIF
1233         ENDDO         ENDDO
1234    
1235         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1236              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1237              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1238    
1239         ! ajout des tendances         ! ajout des tendances
1240         DO k = 1, llm         DO k = 1, llm
# Line 1486  contains Line 1251  contains
1251         igwd = 0         igwd = 0
1252         DO i = 1, klon         DO i = 1, klon
1253            itest(i) = 0            itest(i) = 0
1254            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
1255               itest(i) = 1               itest(i) = 1
1256               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1257            ENDIF            ENDIF
1258         ENDDO         ENDDO
1259    
# Line 1522  contains Line 1286  contains
1286         ENDDO         ENDDO
1287      ENDDO      ENDDO
1288    
1289      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1290           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1291    
1292      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1293           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, &
1294           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1295    
1296      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1297      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, t, &
1298           dtphys, u, t, paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, &           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1299           entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, &
1300           albsol, rhcl, cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, &           dnwd, tr_seri, zmasse, ncid_startphy, nid_ins, itau_phy)
          psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)  
1301    
1302      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1303           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1304           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
1305    
1306      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1307      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)  
1308    
1309      ! diag. bilKP      ! diag. bilKP
1310    
1311      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, &
1312           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1313    
1314      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
# Line 1565  contains Line 1327  contains
1327      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1328         tit = 'after physic'         tit = 'after physic'
1329         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, &
1330              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)  
1331         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1332         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1333         ! 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.
1334         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1335         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1336              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)  
   
1337         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1338      END IF      END IF
1339    
1340      ! SORTIES      ! SORTIES
# Line 1601  contains Line 1359  contains
1359         ENDDO         ENDDO
1360      ENDDO      ENDDO
1361    
1362      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1363         DO iq = 3, nqmx         DO k = 1, llm
1364            DO k = 1, llm            DO i = 1, klon
1365               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  
1366            ENDDO            ENDDO
1367         ENDDO         ENDDO
1368      ENDIF      ENDDO
1369    
1370      ! 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:
1371      DO k = 1, llm      DO k = 1, llm
# Line 1619  contains Line 1375  contains
1375         ENDDO         ENDDO
1376      ENDDO      ENDDO
1377    
     ! Ecriture des sorties  
1378      call write_histins      call write_histins
1379    
1380      ! Si c'est la fin, il faut conserver l'etat de redemarrage      IF (lafin) then
1381      IF (lafin) THEN         call NF95_CLOSE(ncid_startphy)
1382         itau_phy = itau_phy + itap         CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1383         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1384              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1385              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &              t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1386              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              w01)
1387              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)      end IF
     ENDIF  
1388    
1389      firstcal = .FALSE.      firstcal = .FALSE.
1390    
# Line 1640  contains Line 1394  contains
1394    
1395        ! 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
1396    
1397          ! Ecriture des sorties
1398    
1399        use dimens_m, only: iim, jjm        use dimens_m, only: iim, jjm
1400        USE histsync_m, ONLY: histsync        USE histsync_m, ONLY: histsync
1401        USE histwrite_m, ONLY: histwrite        USE histwrite_m, ONLY: histwrite
1402    
1403        real zout        integer i, itau_w ! pas de temps ecriture
       integer itau_w ! pas de temps ecriture  
1404        REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)        REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1405    
1406        !--------------------------------------------------        !--------------------------------------------------
# Line 1653  contains Line 1408  contains
1408        IF (ok_instan) THEN        IF (ok_instan) THEN
1409           ! Champs 2D:           ! Champs 2D:
1410    
          zsto = dtphys * ecrit_ins  
          zout = dtphys * ecrit_ins  
1411           itau_w = itau_phy + itap           itau_w = itau_phy + itap
1412    
          i = NINT(zout/zsto)  
1413           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)
1414           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1415    
          i = NINT(zout/zsto)  
1416           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)
1417           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1418    
# Line 1731  contains Line 1482  contains
1482           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)
1483           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1484    
1485           zx_tmp_fi2d(1:klon) = -1*sens(1:klon)           zx_tmp_fi2d(1:klon) = - sens(1:klon)
1486           ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)           ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)
1487           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1488           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
# Line 1793  contains Line 1544  contains
1544              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1545                   zx_tmp_2d)                   zx_tmp_2d)
1546    
1547              zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = falbe(:, nsrf)
1548              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1549              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1550                   zx_tmp_2d)                   zx_tmp_2d)
# Line 1801  contains Line 1552  contains
1552           END DO           END DO
1553           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)
1554           CALL histwrite(nid_ins, "albs", itau_w, 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)  
1555    
1556           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)
1557           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
# Line 1862  contains Line 1611  contains
1611           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)
1612           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1613    
1614             CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zx_rh, zx_tmp_3d)
1615             CALL histwrite(nid_ins, "rhum", itau_w, zx_tmp_3d)
1616    
1617           call histsync(nid_ins)           call histsync(nid_ins)
1618        ENDIF        ENDIF
1619    
Legend:
Removed from v.92  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21