/[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 186 by guez, Mon Mar 21 15:36: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, conv_emanuel, 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. Anne Mathieu variables pour la couche limite atmosphérique (hbtm)
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 455  contains Line 343  contains
343      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
344      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
345      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
346      ! Grdeurs de sorties      ! Grandeurs de sorties
347      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
348      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
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 502  contains Line 381  contains
381      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
382    
383      INTEGER, save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
384        real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
385    
386      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
387      REAL snow_con(klon), snow_lsc(klon)      REAL, save:: snow_con(klon) ! neige (mm / s)
388        real snow_lsc(klon)
389      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
390    
391      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
# Line 528  contains Line 409  contains
409      integer:: iflag_cldcon = 1      integer:: iflag_cldcon = 1
410      logical ptconv(klon, llm)      logical ptconv(klon, llm)
411    
412      ! Variables locales pour effectuer les appels en s\'erie :      ! Variables pour effectuer les appels en s\'erie :
413    
414      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
415      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
416      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
417        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
418    
419      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
420    
# Line 546  contains Line 425  contains
425    
426      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
427    
428      INTEGER, SAVE:: nid_day, nid_ins      INTEGER, SAVE:: nid_ins
429    
430      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.
431      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.
432      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.
433      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.
434    
     REAL zsto  
435      real date0      real date0
436    
437      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
438      REAL ztsol(klon)      REAL ztsol(klon)
439      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
440      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
441      REAL zero_v(klon)      REAL zero_v(klon)
442      CHARACTER(LEN = 15) tit      CHARACTER(LEN = 20) tit
443      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
444      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
445    
# Line 579  contains Line 456  contains
456      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
457    
458      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
459      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value      ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value
460    
461      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
462      ! cloud optical thickness for pre-industrial (pi) aerosols      ! cloud optical thickness for pre-industrial (pi) aerosols
# Line 611  contains Line 488  contains
488      SAVE ffonte      SAVE ffonte
489      SAVE fqcalving      SAVE fqcalving
490      SAVE rain_con      SAVE rain_con
     SAVE snow_con  
491      SAVE topswai      SAVE topswai
492      SAVE topswad      SAVE topswad
493      SAVE solswai      SAVE solswai
# 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         ! on remet le calendrier a zero         radpas = lmt_pas / nbapp_rad
576    
577           ! 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 (conv_emanuel) THEN
584            ibas_con = 1            ibas_con = 1
585            itop_con = 1            itop_con = 1
586         ENDIF         ENDIF
# 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
669    
670      ! Mettre en action les conditions aux limites (albedo, sst etc.).      ! Prescrire l'ozone :
   
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.  
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)
     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  
694    
695      ! calculs necessaires au calcul de l'albedo dans l'interface      ! Calculs nécessaires au calcul de l'albedo dans l'interface avec
696        ! la surface.
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, 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)
857    
858      DO k = 1, llm      ! Appeler la convection (au choix)
        DO i = 1, klon  
           conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys  
           conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys  
        ENDDO  
     ENDDO  
   
     IF (check) THEN  
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
   
     if (iflag_con == 2) then  
        z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)  
        CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &  
             q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &  
             d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &  
             mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &  
             kdtop, pmflxr, pmflxs)  
        WHERE (rain_con < 0.) rain_con = 0.  
        WHERE (snow_con < 0.) snow_con = 0.  
        ibas_con = llm + 1 - kcbot  
        itop_con = llm + 1 - kctop  
     else  
        ! iflag_con >= 3  
   
        CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &  
             v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &  
             d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &  
             itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &  
             pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &  
             wd, pmflxr, pmflxs, da, phi, mp, ntra=1)  
        ! (number of tracers for the convection scheme of Kerry Emanuel:  
        ! la partie traceurs est faite dans phytrac  
        ! on met ntra = 1 pour limiter les appels mais on peut  
        ! supprimer les calculs / ftra.)  
859    
860        if (conv_emanuel) then
861           da = 0.
862           mp = 0.
863           phi = 0.
864           CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
865                w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, &
866                itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, wd, &
867                pmflxr, da, phi, mp)
868           snow_con = 0.
869         clwcon0 = qcondc         clwcon0 = qcondc
870         mfu = upwd + dnwd         mfu = upwd + dnwd
871         IF (.NOT. ok_gust) wd = 0.         IF (.NOT. ok_gust) wd = 0.
872    
873         ! Calcul des propri\'et\'es des nuages convectifs         IF (thermcep) THEN
874              zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
875         DO k = 1, llm            zqsat = zqsat / (1. - retv * zqsat)
876            DO i = 1, klon         ELSE
877               zx_t = t_seri(i, k)            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
878               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  
879    
880         ! calcul des proprietes des nuages convectifs         ! Properties of convective clouds
881         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
882         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
883              rnebcon0)              rnebcon0)
884    
885           forall (i = 1:klon) ema_pct(i) = paprs(i,itop_con(i) + 1)
886         mfd = 0.         mfd = 0.
887         pen_u = 0.         pen_u = 0.
888         pen_d = 0.         pen_d = 0.
889         pde_d = 0.         pde_d = 0.
890         pde_u = 0.         pde_u = 0.
891        else
892           conv_q = d_q_dyn + d_q_vdf / dtphys
893           conv_t = d_t_dyn + d_t_vdf / dtphys
894           z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
895           CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
896                q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
897                d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
898                mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
899                kdtop, pmflxr, pmflxs)
900           WHERE (rain_con < 0.) rain_con = 0.
901           WHERE (snow_con < 0.) snow_con = 0.
902           ibas_con = llm + 1 - kcbot
903           itop_con = llm + 1 - kctop
904      END if      END if
905    
906      DO k = 1, llm      DO k = 1, llm
# 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 1153  contains Line 934  contains
934         print *, "Precip = ", zx_t         print *, "Precip = ", zx_t
935      ENDIF      ENDIF
936    
937      IF (iflag_con == 2) THEN      IF (.not. conv_emanuel) THEN
938         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
939         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
940         DO k = 1, llm         DO k = 1, llm
# 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
              idx(igwd) = i  
1233            ENDIF            ENDIF
1234         ENDDO         ENDDO
1235    
1236         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1237              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1238              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1239    
1240         ! ajout des tendances         ! ajout des tendances
1241         DO k = 1, llm         DO k = 1, llm
# Line 1486  contains Line 1252  contains
1252         igwd = 0         igwd = 0
1253         DO i = 1, klon         DO i = 1, klon
1254            itest(i) = 0            itest(i) = 0
1255            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
1256               itest(i) = 1               itest(i) = 1
1257               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1258            ENDIF            ENDIF
1259         ENDDO         ENDDO
1260    
# Line 1522  contains Line 1287  contains
1287         ENDDO         ENDDO
1288      ENDDO      ENDDO
1289    
1290      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1291           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1292    
1293      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1294           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, &
1295           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1296    
1297      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1298      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, t, &
1299           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, &
1300           entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, &
1301           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)  
1302    
1303      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1304           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, &
1305           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
1306    
1307      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1308      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)  
1309    
1310      ! diag. bilKP      ! diag. bilKP
1311    
1312      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, &
1313           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1314    
1315      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
# Line 1565  contains Line 1328  contains
1328      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1329         tit = 'after physic'         tit = 'after physic'
1330         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, &
1331              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)  
1332         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1333         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1334         ! 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.
1335         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1336         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1337              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)  
   
1338         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1339      END IF      END IF
1340    
1341      ! SORTIES      ! SORTIES
# Line 1601  contains Line 1360  contains
1360         ENDDO         ENDDO
1361      ENDDO      ENDDO
1362    
1363      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1364         DO iq = 3, nqmx         DO k = 1, llm
1365            DO k = 1, llm            DO i = 1, klon
1366               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  
1367            ENDDO            ENDDO
1368         ENDDO         ENDDO
1369      ENDIF      ENDDO
1370    
1371      ! 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:
1372      DO k = 1, llm      DO k = 1, llm
# Line 1619  contains Line 1376  contains
1376         ENDDO         ENDDO
1377      ENDDO      ENDDO
1378    
     ! Ecriture des sorties  
1379      call write_histins      call write_histins
1380    
1381      ! Si c'est la fin, il faut conserver l'etat de redemarrage      IF (lafin) then
1382      IF (lafin) THEN         call NF95_CLOSE(ncid_startphy)
1383         itau_phy = itau_phy + itap         CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1384         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1385              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1386              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &              t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1387              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              w01)
1388              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)      end IF
     ENDIF  
1389    
1390      firstcal = .FALSE.      firstcal = .FALSE.
1391    
# Line 1640  contains Line 1395  contains
1395    
1396        ! 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
1397    
1398          ! Ecriture des sorties
1399    
1400        use dimens_m, only: iim, jjm        use dimens_m, only: iim, jjm
1401          use gr_fi_ecrit_m, only: gr_fi_ecrit
1402        USE histsync_m, ONLY: histsync        USE histsync_m, ONLY: histsync
1403        USE histwrite_m, ONLY: histwrite        USE histwrite_m, ONLY: histwrite
1404    
1405        real zout        integer i, itau_w ! pas de temps ecriture
       integer itau_w ! pas de temps ecriture  
1406        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)
1407    
1408        !--------------------------------------------------        !--------------------------------------------------
# Line 1653  contains Line 1410  contains
1410        IF (ok_instan) THEN        IF (ok_instan) THEN
1411           ! Champs 2D:           ! Champs 2D:
1412    
          zsto = dtphys * ecrit_ins  
          zout = dtphys * ecrit_ins  
1413           itau_w = itau_phy + itap           itau_w = itau_phy + itap
1414    
          i = NINT(zout/zsto)  
1415           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)
1416           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1417    
          i = NINT(zout/zsto)  
1418           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)
1419           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1420    
# Line 1731  contains Line 1484  contains
1484           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)
1485           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1486    
1487           zx_tmp_fi2d(1:klon) = -1*sens(1:klon)           zx_tmp_fi2d(1:klon) = - sens(1:klon)
1488           ! 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)
1489           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)
1490           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
# Line 1793  contains Line 1546  contains
1546              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1547                   zx_tmp_2d)                   zx_tmp_2d)
1548    
1549              zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = falbe(:, nsrf)
1550              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)
1551              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1552                   zx_tmp_2d)                   zx_tmp_2d)
# Line 1801  contains Line 1554  contains
1554           END DO           END DO
1555           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)
1556           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)  
1557    
1558           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)
1559           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
# Line 1839  contains Line 1590  contains
1590           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)
1591           CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
1592    
1593             if (conv_emanuel) then
1594                CALL gr_fi_ecrit(1, klon, iim, jjm + 1, ema_pct, zx_tmp_2d)
1595                CALL histwrite(nid_ins, "ptop", itau_w, zx_tmp_2d)
1596             end if
1597    
1598           ! Champs 3D:           ! Champs 3D:
1599    
1600           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)
# Line 1862  contains Line 1618  contains
1618           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)
1619           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1620    
1621             CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zx_rh, zx_tmp_3d)
1622             CALL histwrite(nid_ins, "rhum", itau_w, zx_tmp_3d)
1623    
1624           call histsync(nid_ins)           call histsync(nid_ins)
1625        ENDIF        ENDIF
1626    
Legend:
Removed from v.92  
changed lines
  Added in v.186
  ViewVC Help
Powered by ViewVC 1.1.21