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