/[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/libf/phylmd/physiq.f90 revision 53 by guez, Fri Oct 7 13:11:58 2011 UTC trunk/Sources/phylmd/physiq.f revision 191 by guez, Mon May 9 19:56:28 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, d_ps, dudyn, PVteta)         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 (SVN revision 678)      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11      ! Author: Z.X. Li (LMD/CNRS) 1993      ! (subversion revision 678)
12    
13        ! Author: Z. X. Li (LMD/CNRS) 1993
14    
15      ! This is the main procedure for the "physics" part of the program.      ! This is the main procedure for the "physics" part of the program.
16    
17        use aaam_bud_m, only: aaam_bud
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
19        use aeropt_m, only: aeropt
20      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
     USE calendar, ONLY: ymds2ju  
21      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
22      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &      USE clesphys, ONLY: cdhmax, cdmmax, 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, ok_instan
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 comgeomphy, ONLY: airephy, cuphy, cvphy      use clouds_gno_m, only: clouds_gno
28        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, day_step, iphysiq
32      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
33        use conflx_m, only: conflx
34      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
35      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
36      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
37      USE dimens_m, ONLY: iim, jjm, llm, nqmx      use diagphy_m, only: diagphy
38      USE dimphy, ONLY: klon, nbtr      USE dimens_m, ONLY: llm, nqmx
39        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
45      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
46      USE histcom, ONLY: histsync      USE histsync_m, ONLY: histsync
47      USE histwrite_m, ONLY: histwrite      USE histwrite_phy_m, ONLY: histwrite_phy
48      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, &
49           nbsrf           nbsrf
50      USE ini_histhf_m, ONLY: ini_histhf      USE ini_histins_m, ONLY: ini_histins, nid_ins
51      USE ini_histday_m, ONLY: ini_histday      use netcdf95, only: NF95_CLOSE
52      USE ini_histins_m, ONLY: ini_histins      use newmicro_m, only: newmicro
53      USE oasis_m, ONLY: ok_oasis      use nr_util, only: assert
54      USE orbite_m, ONLY: orbite, zenang      use nuage_m, only: nuage
55        USE orbite_m, ONLY: orbite
56      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
57      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
58      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
59        USE phyredem0_m, ONLY: phyredem0
60      USE phystokenc_m, ONLY: phystokenc      USE phystokenc_m, ONLY: phystokenc
61      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
62      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
63      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
64      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      use readsulfate_m, only: readsulfate
65      USE temps, ONLY: annee_ref, day_ref, itau_phy      use readsulfate_preind_m, only: readsulfate_preind
66        use yoegwd, only: sugwd
67        USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt
68        use time_phylmdz, only: itap, increment_itap
69        use transp_m, only: transp
70        use transp_lay_m, only: transp_lay
71        use unit_nml_m, only: unit_nml
72        USE ymds2ju_m, ONLY: ymds2ju
73      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
74        use zenang_m, only: zenang
75    
76      ! Arguments:      logical, intent(in):: lafin ! dernier passage
77    
78      REAL, intent(in):: rdayvrai      integer, intent(in):: dayvrai
79      ! (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
80    
81      REAL, intent(in):: time ! heure de la journée en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
     REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)  
     logical, intent(in):: lafin ! dernier passage  
82    
83      REAL, intent(in):: paprs(klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
84      ! (pression pour chaque inter-couche, en Pa)      ! pression pour chaque inter-couche, en Pa
85    
86      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
87      ! (input pression pour le mileu de chaque couche (en Pa))      ! pression pour le mileu de chaque couche (en Pa)
88    
89      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
90      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! géopotentiel de chaque couche (référence sol)
91    
92      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
93    
94      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
95      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m/s
96    
97      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
98      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
99    
100      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
101      ! (humidité spécifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
102    
103      REAL omega(klon, llm) ! input vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
104      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
105      REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/s)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
106      REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)
     REAL d_qx(klon, llm, nqmx) ! output tendance physique de "qx" (kg/kg/s)  
     REAL d_ps(klon) ! output tendance physique de la pression au sol  
   
     LOGICAL:: firstcal = .true.  
107    
108      INTEGER nbteta      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
109      PARAMETER(nbteta = 3)      ! tendance physique de "qx" (s-1)
110    
111      REAL PVteta(klon, nbteta)      ! Local:
     ! (output vorticite potentielle a des thetas constantes)  
112    
113      LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE      LOGICAL:: firstcal = .true.
     PARAMETER (ok_cvl = .TRUE.)  
     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface  
     PARAMETER (ok_gust = .FALSE.)  
114    
115      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL, PARAMETER:: check = .FALSE.
116      PARAMETER (check = .FALSE.)      ! Verifier la conservation du modele en eau
117    
118      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
119      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
120    
     ! Parametres lies au coupleur OASIS:  
     INTEGER, SAVE:: npas, nexca  
     logical rnpb  
     parameter(rnpb = .true.)  
   
     character(len = 6), save:: ocean  
     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")  
   
     logical ok_ocean  
     SAVE ok_ocean  
   
     ! "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, save:: ok_veget  
     LOGICAL, save:: ok_journe ! sortir le fichier journalier  
   
     LOGICAL ok_mensuel ! sortir le fichier mensuel  
   
     LOGICAL ok_instan ! sortir le fichier instantane  
     save ok_instan  
   
     LOGICAL ok_region ! sortir le fichier regional  
     PARAMETER (ok_region = .FALSE.)  
   
121      ! pour phsystoke avec thermiques      ! pour phsystoke avec thermiques
122      REAL fm_therm(klon, llm + 1)      REAL fm_therm(klon, llm + 1)
123      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
124      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
125    
126      INTEGER ivap ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
127      PARAMETER (ivap = 1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     INTEGER iliq ! indice de traceurs pour eau liquide  
     PARAMETER (iliq = 2)  
128    
129      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
130      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
# Line 156  contains Line 134  contains
134    
135      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
136    
137      !IM Amip2 PV a theta constante      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
138        REAL swup0(klon, llm + 1), swup(klon, llm + 1)
     CHARACTER(LEN = 3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
   
     !MI Amip2 PV a theta constante  
   
     INTEGER klevp1  
     PARAMETER(klevp1 = llm + 1)  
   
     REAL swdn0(klon, klevp1), swdn(klon, klevp1)  
     REAL swup0(klon, klevp1), swup(klon, klevp1)  
139      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
140    
141      REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
142      REAL lwup0(klon, klevp1), lwup(klon, klevp1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
143      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
   
     !IM Amip2  
     ! variables a une pression donnee  
   
     integer nlevSTD  
     PARAMETER(nlevSTD = 17)  
     real rlevSTD(nlevSTD)  
     DATA rlevSTD/100000., 92500., 85000., 70000., &  
          60000., 50000., 40000., 30000., 25000., 20000., &  
          15000., 10000., 7000., 5000., 3000., 2000., 1000./  
     CHARACTER(LEN = 4) clevSTD(nlevSTD)  
     DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &  
          '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &  
          '70 ', '50 ', '30 ', '20 ', '10 '/  
144    
145      ! prw: precipitable water      ! prw: precipitable water
146      real prw(klon)      real prw(klon)
# Line 198  contains Line 150  contains
150      REAL flwp(klon), fiwp(klon)      REAL flwp(klon), fiwp(klon)
151      REAL flwc(klon, llm), fiwc(klon, llm)      REAL flwc(klon, llm), fiwc(klon, llm)
152    
     INTEGER kmax, lmax  
     PARAMETER(kmax = 8, lmax = 8)  
     INTEGER kmaxm1, lmaxm1  
     PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)  
   
     REAL zx_tau(kmaxm1), zx_pc(lmaxm1)  
     DATA zx_tau/0.0, 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.'/  
   
     !IM ISCCP simulator v3.4  
   
     integer nid_hf, nid_hf3d  
     save nid_hf, nid_hf3d  
   
153      ! Variables propres a la physique      ! Variables propres a la physique
154    
155      INTEGER, save:: radpas      INTEGER, save:: radpas
156      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
157      ! "physiq".)      ! "physiq".
158    
159      REAL radsol(klon)      REAL radsol(klon)
160      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
161    
     INTEGER, SAVE:: itap ! number of calls to "physiq"  
   
162      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
163    
164      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
165      ! soil temperature of surface fraction      ! soil temperature of surface fraction
166    
167      REAL fevap(klon, nbsrf)      REAL, save:: fevap(klon, nbsrf) ! evaporation
     SAVE fevap ! evaporation  
168      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
169      SAVE fluxlat      SAVE fluxlat
170    
171      REAL fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
172      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  
173    
174      REAL fsnow(klon, nbsrf)      REAL, save:: qsol(klon)
175      SAVE fsnow ! epaisseur neigeuse      ! column-density of water in soil, in kg m-2
176    
177      REAL falbe(klon, nbsrf)      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
178      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  
179    
180      ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
181      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
182      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
183      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 293  contains Line 186  contains
186      REAL, save:: zpic(klon) ! Maximum de l'OESM      REAL, save:: zpic(klon) ! Maximum de l'OESM
187      REAL, save:: zval(klon) ! Minimum de l'OESM      REAL, save:: zval(klon) ! Minimum de l'OESM
188      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
   
189      REAL zulow(klon), zvlow(klon)      REAL zulow(klon), zvlow(klon)
190        INTEGER igwd, itest(klon)
191    
192      INTEGER igwd, idx(klon), itest(klon)      REAL, save:: agesno(klon, nbsrf) ! age de la neige
193        REAL, save:: run_off_lic_0(klon)
     REAL agesno(klon, nbsrf)  
     SAVE agesno ! age de la neige  
   
     REAL run_off_lic_0(klon)  
     SAVE run_off_lic_0  
     !KE43  
     ! Variables liees a la convection de K. Emanuel (sb):  
   
     REAL bas, top ! cloud base and top levels  
     SAVE bas  
     SAVE top  
   
     REAL Ma(klon, llm) ! undilute upward mass flux  
     SAVE Ma  
     REAL qcondc(klon, llm) ! in-cld water content from convect  
     SAVE qcondc  
     REAL ema_work1(klon, llm), ema_work2(klon, llm)  
     SAVE ema_work1, ema_work2  
   
     REAL wd(klon) ! sb  
     SAVE wd ! sb  
   
     ! Variables locales pour la couche limite (al1):  
194    
195      ! Variables locales:      ! Variables li\'ees \`a la convection d'Emanuel :
196        REAL, save:: Ma(klon, llm) ! undilute upward mass flux
197        REAL, save:: qcondc(klon, llm) ! in-cld water content from convect
198        REAL, save:: sig1(klon, llm), w01(klon, llm)
199    
200        ! Variables pour la couche limite (Alain Lahellec) :
201      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
202      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
203    
204      !AA Pour phytrac      ! Pour phytrac :
205      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
206      REAL yu1(klon) ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
207      REAL yv1(klon) ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
208      REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige      REAL ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
209      REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface  
210      ! !et necessaire pour limiter la      REAL fqcalving(klon, nbsrf)
211      ! !hauteur de neige, en kg/m2/s      ! flux d'eau "perdue" par la surface et necessaire pour limiter la
212        ! hauteur de neige, en kg/m2/s
213    
214      REAL zxffonte(klon), zxfqcalving(klon)      REAL zxffonte(klon), zxfqcalving(klon)
215    
216      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
# Line 346  contains Line 222  contains
222      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
223      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
224    
225      !AA      REAL, save:: rain_fall(klon)
226      REAL rain_fall(klon) ! pluie      ! liquid water mass flux (kg/m2/s), positive down
227      REAL snow_fall(klon) ! neige  
228      save snow_fall, rain_fall      REAL, save:: snow_fall(klon)
229      !IM cf FH pour Tiedtke 080604      ! solid water mass flux (kg/m2/s), positive down
230    
231      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
232    
233      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon), devap(klon) ! evaporation and its derivative
234      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
235      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
236      SAVE dlw      SAVE dlw
237      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
238      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
     save fder  
239      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
240      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
241      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
242      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
243    
244      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
245      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
246    
247      ! Conditions aux limites      ! Conditions aux limites
248    
249      INTEGER julien      INTEGER julien
   
250      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
251      REAL pctsrf(klon, nbsrf)      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
252      !IM      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
253      REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE      REAL, save:: albsol(klon) ! albedo du sol total visible
   
     SAVE pctsrf ! sous-fraction du sol  
     REAL albsol(klon)  
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
254      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
255    
256      ! Declaration des procedures appelees      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
257        real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
     EXTERNAL alboc ! calculer l'albedo sur ocean  
     !KE43  
     EXTERNAL conema3 ! convect4.3  
     EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)  
     EXTERNAL nuage ! calculer les proprietes radiatives  
     EXTERNAL transp ! transport total de l'eau et de l'energie  
   
     ! Variables locales  
   
     real clwcon(klon, llm), rnebcon(klon, llm)  
     real clwcon0(klon, llm), rnebcon0(klon, llm)  
   
     save rnebcon, clwcon  
258    
259      REAL rhcl(klon, llm) ! humiditi relative ciel clair      REAL rhcl(klon, llm) ! humiditi relative ciel clair
260      REAL dialiq(klon, llm) ! eau liquide nuageuse      REAL dialiq(klon, llm) ! eau liquide nuageuse
# Line 420  contains Line 274  contains
274      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
275      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
276    
277      ! Le rayonnement n'est pas calcule tous les pas, il faut donc      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
278      ! que les variables soient rémanentes      ! les variables soient r\'emanentes.
279      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
280      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
281      REAL cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
282      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
283      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
284      real sollwdown(klon) ! downward LW flux at surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
285      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      real, save:: sollwdown(klon) ! downward LW flux at surface
286      REAL albpla(klon)      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
287      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL, save:: albpla(klon)
288      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
289      SAVE cool, albpla, topsw, toplw, solsw, sollw, sollwdown      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
     SAVE topsw0, toplw0, solsw0, sollw0, heat0, cool0  
   
     INTEGER itaprad  
     SAVE itaprad  
290    
291      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
292      REAL conv_t(klon, llm) ! convergence of temperature (K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K/s)
293    
294      REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut      REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
295      REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree      REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
296    
297      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
298    
299      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
300      REAL zdtime ! pas de temps du rayonnement (s)      real longi
     real zlongi  
   
301      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
302      REAL za, zb      REAL za, zb
303      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
304      real zqsat(klon, llm)      real zqsat(klon, llm)
305      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
306      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup = 234.0)  
   
307      REAL zphi(klon, llm)      REAL zphi(klon, llm)
308    
309      !IM cf. AM Variables locales pour la CLA (hbtm2)      ! cf. Anne Mathieu variables pour la couche limite atmosphérique (hbtm)
310    
311      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
312      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
# Line 472  contains Line 316  contains
316      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
317      REAL, SAVE:: therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
318      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
319      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
320      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
321      ! Grdeurs de sorties      ! Grandeurs de sorties
322      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
323      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
324      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
325      REAL s_trmb3(klon)      REAL s_trmb3(klon)
326    
327      ! Variables locales pour la convection de K. Emanuel (sb):      ! Variables pour la convection de K. Emanuel :
328    
329      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
330      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
331      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  
332      REAL cape(klon) ! CAPE      REAL cape(klon) ! CAPE
333      SAVE cape      SAVE cape
334    
     REAL pbase(klon) ! cloud base pressure  
     SAVE pbase  
     REAL bbase(klon) ! cloud base buoyancy  
     SAVE bbase  
     REAL rflag(klon) ! flag fonctionnement de convect  
335      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
     ! -- convect43:  
     INTEGER ntra ! nb traceurs pour convect4.3  
     REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)  
     REAL dplcldt(klon), dplcldr(klon)  
336    
337      ! Variables du changement      ! Variables du changement
338    
339      ! con: convection      ! con: convection
340      ! lsc: large scale condensation      ! lsc: large scale condensation
341      ! ajs: ajustement sec      ! ajs: ajustement sec
342      ! eva: évaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
343      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
344      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
345      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
# Line 514  contains Line 348  contains
348      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
349      REAL rneb(klon, llm)      REAL rneb(klon, llm)
350    
351      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
352      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
353      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
354      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
355      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
356      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
357    
358      INTEGER,save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
359        real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
360    
361      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
362      REAL snow_con(klon), snow_lsc(klon)      REAL, save:: snow_con(klon) ! neige (mm / s)
363        real snow_lsc(klon)
364      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
365    
366      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
# Line 535  contains Line 371  contains
371      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
372      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
373    
374      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
375      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
376      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
377    
378      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
379      real, save:: fact_cldcon      real:: fact_cldcon = 0.375
380      real, save:: facttemps      real:: facttemps = 1.e-4
381      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
382      real facteur      real facteur
383    
384      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
385      logical ptconv(klon, llm)      logical ptconv(klon, llm)
386    
387      ! Variables locales pour effectuer les appels en série :      ! Variables pour effectuer les appels en s\'erie :
388    
389      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
390      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
391      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
392        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
393    
394      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
395    
# Line 567  contains Line 398  contains
398      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
399      REAL aam, torsfc      REAL aam, torsfc
400    
     REAL dudyn(iim + 1, jjm + 1, llm)  
   
     REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique  
     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)  
   
     INTEGER, SAVE:: nid_day, nid_ins  
   
401      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.
402      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.
403      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.
404      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.
405    
     REAL zsto  
   
     character(len = 20) modname  
     character(len = 80) abort_message  
     logical ok_sync  
406      real date0      real date0
407    
408      ! Variables liées au bilan d'énergie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
409      REAL ztsol(klon)      REAL ztsol(klon)
410      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
411      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
412      REAL zero_v(klon)      REAL zero_v(klon)
413      CHARACTER(LEN = 15) ztit      CHARACTER(LEN = 20) tit
414      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
415      INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
416    
417      REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique      REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
418      REAL ZRCPD      REAL ZRCPD
419    
420      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
421      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
422      REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille      REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
423      REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille      REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
424      !jq Aerosol effects (Johannes Quaas, 27/11/2003)  
425      REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]      ! Aerosol effects:
426    
427        REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
428    
429      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
430      ! (SO4 aerosol concentration, in ug/m3, pre-industrial value)      ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value
431    
432      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
433      ! (Cloud optical thickness for pre-industrial (pi) aerosols)      ! cloud optical thickness for pre-industrial aerosols
434    
435      REAL re(klon, llm) ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
436      REAL fl(klon, llm) ! denominator of re      REAL fl(klon, llm) ! denominator of re
437    
438      ! Aerosol optical properties      ! Aerosol optical properties
439      REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)      REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
440      REAL cg_ae(klon, llm, 2)      REAL, save:: cg_ae(klon, llm, 2)
   
     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.  
     ! ok_ade = True -ADE = topswad-topsw  
441    
442      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.      REAL topswad(klon), solswad(klon) ! aerosol direct effect
443      ! ok_aie = True ->      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
     ! ok_ade = True -AIE = topswai-topswad  
     ! ok_ade = F -AIE = topswai-topsw  
444    
445      REAL aerindex(klon) ! POLDER aerosol index      REAL aerindex(klon) ! POLDER aerosol index
446    
447      ! Parameters      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
448      LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
449      REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)  
450        REAL:: bl95_b0 = 2., bl95_b1 = 0.2
451        ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
452        ! B). They link cloud droplet number concentration to aerosol mass
453        ! concentration.
454    
     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1  
455      SAVE u10m      SAVE u10m
456      SAVE v10m      SAVE v10m
457      SAVE t2m      SAVE t2m
458      SAVE q2m      SAVE q2m
459      SAVE ffonte      SAVE ffonte
460      SAVE fqcalving      SAVE fqcalving
     SAVE piz_ae  
     SAVE tau_ae  
     SAVE cg_ae  
461      SAVE rain_con      SAVE rain_con
     SAVE snow_con  
462      SAVE topswai      SAVE topswai
463      SAVE topswad      SAVE topswad
464      SAVE solswai      SAVE solswai
465      SAVE solswad      SAVE solswad
466      SAVE d_u_con      SAVE d_u_con
467      SAVE d_v_con      SAVE d_v_con
     SAVE rnebcon0  
     SAVE clwcon0  
468    
469      real zmasse(klon, llm)      real zmasse(klon, llm)
470      ! (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)
471    
472      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      integer, save:: ncid_startphy
473    
474        namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, &
475             iflag_cldcon, ratqsbas, ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, &
476             bl95_b1, iflag_thermals, nsplit_thermals
477    
478      !----------------------------------------------------------------      !----------------------------------------------------------------
479    
480      modname = 'physiq'      IF (if_ebil >= 1) zero_v = 0.
481      IF (if_ebil >= 1) THEN      IF (nqmx < 2) CALL abort_gcm('physiq', &
482         DO i = 1, klon           'eaux vapeur et liquide sont indispensables')
           zero_v(i) = 0.  
        END DO  
     END IF  
     ok_sync = .TRUE.  
     IF (nqmx < 2) THEN  
        abort_message = 'eaux vapeur et liquide sont indispensables'  
        CALL abort_gcm(modname, abort_message, 1)  
     ENDIF  
483    
484      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
485         ! initialiser         ! initialiser
# Line 684  contains Line 492  contains
492         piz_ae = 0.         piz_ae = 0.
493         tau_ae = 0.         tau_ae = 0.
494         cg_ae = 0.         cg_ae = 0.
495         rain_con(:) = 0.         rain_con = 0.
496         snow_con(:) = 0.         snow_con = 0.
497         bl95_b0 = 0.         topswai = 0.
498         bl95_b1 = 0.         topswad = 0.
499         topswai(:) = 0.         solswai = 0.
500         topswad(:) = 0.         solswad = 0.
501         solswai(:) = 0.  
502         solswad(:) = 0.         d_u_con = 0.
503           d_v_con = 0.
504         d_u_con = 0.0         rnebcon0 = 0.
505         d_v_con = 0.0         clwcon0 = 0.
506         rnebcon0 = 0.0         rnebcon = 0.
507         clwcon0 = 0.0         clwcon = 0.
        rnebcon = 0.0  
        clwcon = 0.0  
508    
509         pblh =0. ! Hauteur de couche limite         pblh =0. ! Hauteur de couche limite
510         plcl =0. ! Niveau de condensation de la CLA         plcl =0. ! Niveau de condensation de la CLA
# Line 708  contains Line 514  contains
514         pblt =0. ! T a la Hauteur de couche limite         pblt =0. ! T a la Hauteur de couche limite
515         therm =0.         therm =0.
516         trmb1 =0. ! deep_cape         trmb1 =0. ! deep_cape
517         trmb2 =0. ! inhibition         trmb2 =0. ! inhibition
518         trmb3 =0. ! Point Omega         trmb3 =0. ! Point Omega
519    
520         IF (if_ebil >= 1) d_h_vcol_phy = 0.         IF (if_ebil >= 1) d_h_vcol_phy = 0.
521    
522         ! appel a la lecture du run.def physique         iflag_thermals = 0
523           nsplit_thermals = 1
524           print *, "Enter namelist 'physiq_nml'."
525           read(unit=*, nml=physiq_nml)
526           write(unit_nml, nml=physiq_nml)
527    
528         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &         call conf_phys
             ok_instan, fact_cldcon, facttemps, ok_newmicro, &  
             iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &  
             ok_ade, ok_aie, &  
             bl95_b0, bl95_b1, &  
             iflag_thermals, nsplit_thermals)  
529    
530         ! Initialiser les compteurs:         ! Initialiser les compteurs:
531    
532         frugs = 0.         frugs = 0.
533         itap = 0         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
534         itaprad = 0              fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
535         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
536              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
537              snow_fall, solsw, sollwdown, dlw, radsol, frugs, agesno, zmea, &              w01, ncid_startphy)
             zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &  
             ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)  
538    
539         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
540         q2 = 1.e-8         q2 = 1e-8
   
        radpas = NINT(86400. / dtphys / nbapp_rad)  
   
        ! on remet le calendrier a zero  
        IF (raz_date) itau_phy = 0  
   
        PRINT *, 'cycle_diurne = ', cycle_diurne  
541    
542         IF(ocean.NE.'force ') THEN         lmt_pas = day_step / iphysiq
543            ok_ocean = .TRUE.         print *, 'Number of time steps of "physics" per day: ', lmt_pas
        ENDIF  
   
        CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &  
             ok_region)  
   
        IF (dtphys*REAL(radpas) > 21600..AND.cycle_diurne) THEN  
           print *,'Nbre d appels au rayonnement insuffisant'  
           print *,"Au minimum 4 appels par jour si cycle diurne"  
           abort_message = 'Nbre d appels au rayonnement insuffisant'  
           call abort_gcm(modname, abort_message, 1)  
        ENDIF  
        print *,"Clef pour la convection, iflag_con = ", iflag_con  
        print *,"Clef pour le driver de la convection, ok_cvl = ", &  
             ok_cvl  
   
        ! Initialisation pour la convection de K.E. (sb):  
        IF (iflag_con >= 3) THEN  
   
           print *,"*** Convection de Kerry Emanuel 4.3 "  
544    
545            !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG         radpas = lmt_pas / nbapp_rad
546            DO i = 1, klon         print *, "radpas = ", radpas
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END  
547    
548           ! Initialisation pour le sch\'ema de convection d'Emanuel :
549           IF (conv_emanuel) THEN
550              ibas_con = 1
551              itop_con = 1
552         ENDIF         ENDIF
553    
554         IF (ok_orodr) THEN         IF (ok_orodr) THEN
555            rugoro = MAX(1e-5, zstd * zsig / 2)            rugoro = MAX(1e-5, zstd * zsig / 2)
556            CALL SUGWD(klon, llm, paprs, play)            CALL SUGWD(paprs, play)
557         else         else
558            rugoro = 0.            rugoro = 0.
559         ENDIF         ENDIF
560    
        lmt_pas = NINT(86400. / dtphys) ! tous les jours  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
561         ecrit_ins = NINT(ecrit_ins/dtphys)         ecrit_ins = NINT(ecrit_ins/dtphys)
562         ecrit_hf = NINT(ecrit_hf/dtphys)         ecrit_hf = NINT(ecrit_hf/dtphys)
563         ecrit_mth = NINT(ecrit_mth/dtphys)         ecrit_mth = NINT(ecrit_mth/dtphys)
564         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
565         ecrit_reg = NINT(ecrit_reg/dtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
566    
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
        print *,'AVANT HIST IFLAG_CON = ', iflag_con  
   
567         ! Initialisation des sorties         ! Initialisation des sorties
568    
569         call ini_histhf(dtphys, nid_hf, nid_hf3d)         call ini_histins(dtphys)
570         call ini_histday(dtphys, ok_journe, nid_day, nqmx)         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
571         call ini_histins(dtphys, ok_instan, nid_ins)         ! Positionner date0 pour initialisation de ORCHIDEE
572         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         print *, 'physiq date0: ', date0
573         !XXXPB Positionner date0 pour initialisation de ORCHIDEE         CALL phyredem0(lmt_pas)
        WRITE(*, *) 'physiq date0: ', date0  
574      ENDIF test_firstcal      ENDIF test_firstcal
575    
576      ! Mettre a zero des variables de sortie (pour securite)      ! We will modify variables *_seri and we will not touch variables
577        ! u, v, t, qx:
578      DO i = 1, klon      t_seri = t
579         d_ps(i) = 0.0      u_seri = u
580      ENDDO      v_seri = v
581      DO iq = 1, nqmx      q_seri = qx(:, :, ivap)
582         DO k = 1, llm      ql_seri = qx(:, :, iliq)
583            DO i = 1, klon      tr_seri = qx(:, :, 3:nqmx)
              d_qx(i, k, iq) = 0.0  
           ENDDO  
        ENDDO  
     ENDDO  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
     ! Ne pas affecter les valeurs entrées de u, v, h, et q :  
   
     DO k = 1, llm  
        DO i = 1, klon  
           t_seri(i, k) = t(i, k)  
           u_seri(i, k) = u(i, k)  
           v_seri(i, k) = v(i, k)  
           q_seri(i, k) = qx(i, k, ivap)  
           ql_seri(i, k) = qx(i, k, iliq)  
           qs_seri(i, k) = 0.  
        ENDDO  
     ENDDO  
     IF (nqmx >= 3) THEN  
        tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)  
     ELSE  
        tr_seri(:, :, 1) = 0.  
     ENDIF  
584    
585      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  
586    
587      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
588         ztit = 'after dynamics'         tit = 'after dynamics'
589         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
590              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)
591              d_ql, d_qs, d_ec)         ! Comme les tendances de la physique sont ajout\'es dans la
592         ! Comme les tendances de la physique sont ajoutés dans la         ! dynamique, la variation d'enthalpie par la dynamique devrait
593         !  dynamique, la variation d'enthalpie par la dynamique devrait         ! \^etre \'egale \`a la variation de la physique au pas de temps
594         !  être égale à la variation de la physique au pas de temps         ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
595         !  précédent.  Donc la somme de ces 2 variations devrait être         ! nulle.
596         !  nulle.         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
        call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
597              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, &
598              d_qt, 0., fs_bound, fq_bound)              d_qt, 0.)
599      END IF      END IF
600    
601      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
# Line 876  contains Line 609  contains
609      ELSE      ELSE
610         DO k = 1, llm         DO k = 1, llm
611            DO i = 1, klon            DO i = 1, klon
612               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
613               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
614            ENDDO            ENDDO
615         ENDDO         ENDDO
616         ancien_ok = .TRUE.         ancien_ok = .TRUE.
# Line 893  contains Line 626  contains
626      ! Check temperatures:      ! Check temperatures:
627      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
628    
629      ! Incrementer le compteur de la physique      call increment_itap
630      itap = itap + 1      julien = MOD(dayvrai, 360)
     julien = MOD(NINT(rdayvrai), 360)  
631      if (julien == 0) julien = 360      if (julien == 0) julien = 360
632    
633      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.).  
634    
635      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone :
636      if (nqmx >= 5) then      wo = ozonecm(REAL(julien), paprs)
        wo = qx(:, :, 5) * zmasse / dobson_u / 1e3  
     else IF (MOD(itap - 1, lmt_pas) == 0) THEN  
        wo = ozonecm(REAL(julien), paprs)  
     ENDIF  
637    
638      ! Évaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
639      DO k = 1, llm      DO k = 1, llm
640         DO i = 1, klon         DO i = 1, klon
641            zb = MAX(0., ql_seri(i, k))            zb = MAX(0., ql_seri(i, k))
# Line 920  contains Line 646  contains
646      ENDDO      ENDDO
647      ql_seri = 0.      ql_seri = 0.
648    
649      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
650         ztit = 'after reevap'         tit = 'after reevap'
651         CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
652              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)
653              d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
654         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &              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)  
   
655      END IF      END IF
656    
657      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
658        zxrugs = sum(frugs * pctsrf, dim = 2)
     DO i = 1, klon  
        zxrugs(i) = 0.0  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)  
        ENDDO  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
659    
660      ! calculs necessaires au calcul de l'albedo dans l'interface      ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
661        ! la surface.
662    
663      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), longi, dist)
664      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
665         zdtime = dtphys * REAL(radpas)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
666      ELSE      ELSE
667         rmu0 = -999.999         mu0 = - 999.999
668      ENDIF      ENDIF
669    
670      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
671      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  
672    
673      ! Repartition sous maille des flux LW et SW      ! R\'epartition sous maille des flux longwave et shortwave
674      ! Repartition du longwave par sous-surface linearisee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
675    
676      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
677         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
678            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
679                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
680            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))      END forall
        ENDDO  
     ENDDO  
681    
682      fder = dlw      fder = dlw
683    
684      ! Couche limite:      ! Couche limite:
685    
686      CALL clmain(dtphys, itap, date0, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, v_seri, &
687           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, &           julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, &
688           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, &
689           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, firstcal, &
690           rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder, rlon, rlat, &           agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, &
691           cuphy, cvphy, frugs, firstcal, lafin, agesno, rugoro, d_t_vdf, &           fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, &
692           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, &
693           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
          pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &  
          fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)  
694    
695      ! Incrémentation des flux      ! Incr\'ementation des flux
696    
697      zxfluxt = 0.      zxfluxt = 0.
698      zxfluxq = 0.      zxfluxq = 0.
# Line 1002  contains Line 701  contains
701      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
702         DO k = 1, llm         DO k = 1, llm
703            DO i = 1, klon            DO i = 1, klon
704               zxfluxt(i, k) = zxfluxt(i, k) + &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
705                    fluxt(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
706               zxfluxq(i, k) = zxfluxq(i, k) + &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
707                    fluxq(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
              zxfluxu(i, k) = zxfluxu(i, k) + &  
                   fluxu(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxv(i, k) = zxfluxv(i, k) + &  
                   fluxv(i, k, nsrf) * pctsrf(i, nsrf)  
708            END DO            END DO
709         END DO         END DO
710      END DO      END DO
711      DO i = 1, klon      DO i = 1, klon
712         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
713         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
714         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
715      ENDDO      ENDDO
716    
# Line 1028  contains Line 723  contains
723         ENDDO         ENDDO
724      ENDDO      ENDDO
725    
726      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
727         ztit = 'after clmain'         tit = 'after clmain'
728         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
729              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)
730              d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
731         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &              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)  
732      END IF      END IF
733    
734      ! Update surface temperature:      ! Update surface temperature:
735    
736      DO i = 1, klon      DO i = 1, klon
737         zxtsol(i) = 0.0         zxtsol(i) = 0.
738         zxfluxlat(i) = 0.0         zxfluxlat(i) = 0.
739    
740         zt2m(i) = 0.0         zt2m(i) = 0.
741         zq2m(i) = 0.0         zq2m(i) = 0.
742         zu10m(i) = 0.0         zu10m(i) = 0.
743         zv10m(i) = 0.0         zv10m(i) = 0.
744         zxffonte(i) = 0.0         zxffonte(i) = 0.
745         zxfqcalving(i) = 0.0         zxfqcalving(i) = 0.
746    
747         s_pblh(i) = 0.0         s_pblh(i) = 0.
748         s_lcl(i) = 0.0         s_lcl(i) = 0.
749         s_capCL(i) = 0.0         s_capCL(i) = 0.
750         s_oliqCL(i) = 0.0         s_oliqCL(i) = 0.
751         s_cteiCL(i) = 0.0         s_cteiCL(i) = 0.
752         s_pblT(i) = 0.0         s_pblT(i) = 0.
753         s_therm(i) = 0.0         s_therm(i) = 0.
754         s_trmb1(i) = 0.0         s_trmb1(i) = 0.
755         s_trmb2(i) = 0.0         s_trmb2(i) = 0.
756         s_trmb3(i) = 0.0         s_trmb3(i) = 0.
   
        IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + &  
             pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.)  >  EPSFRA) &  
             THEN  
           WRITE(*, *) 'physiq : pb sous surface au point ', i, &  
                pctsrf(i, 1 : nbsrf)  
        ENDIF  
757      ENDDO      ENDDO
758    
759        call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
760    
761      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
762         DO i = 1, klon         DO i = 1, klon
763            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
# Line 1095  contains Line 784  contains
784         ENDDO         ENDDO
785      ENDDO      ENDDO
786    
787      ! 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 :
   
788      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
789         DO i = 1, klon         DO i = 1, klon
790            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
# Line 1121  contains Line 809  contains
809         ENDDO         ENDDO
810      ENDDO      ENDDO
811    
812      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
813    
814      DO i = 1, klon      DO i = 1, klon
815         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
816      ENDDO      ENDDO
817    
818      ! Appeler la convection (au choix)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
819    
820      DO k = 1, llm      ! Appeler la convection
        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  
     zx_ajustq = .FALSE.  
     IF (iflag_con == 2) zx_ajustq = .TRUE.  
     IF (zx_ajustq) THEN  
        DO i = 1, klon  
           z_avant(i) = 0.0  
        ENDDO  
        DO k = 1, llm  
           DO i = 1, klon  
              z_avant(i) = z_avant(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
     ENDIF  
821    
822      select case (iflag_con)      if (conv_emanuel) then
823      case (1)         da = 0.
824         print *, 'Réactiver l''appel à "conlmd" dans "physiq.F".'         mp = 0.
825         stop 1         phi = 0.
826      case (2)         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
827         CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, &
828              zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &              itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, &
829              pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &              da, phi, mp)
830              pmflxs)         snow_con = 0.
831         WHERE (rain_con < 0.) rain_con = 0.         clwcon0 = qcondc
832         WHERE (snow_con < 0.) snow_con = 0.         mfu = upwd + dnwd
833         DO i = 1, klon  
834            ibas_con(i) = llm + 1 - kcbot(i)         IF (thermcep) THEN
835            itop_con(i) = llm + 1 - kctop(i)            zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
836         ENDDO            zqsat = zqsat / (1. - retv * zqsat)
     case (3:)  
        ! 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.  
        ntra = 1  
        ! Schéma de convection modularisé et vectorisé :  
        ! (driver commun aux versions 3 et 4)  
   
        IF (ok_cvl) THEN  
           ! new driver for convectL  
           CALL concvl(iflag_con, dtphys, paprs, play, t_seri, q_seri, &  
                u_seri, v_seri, tr_seri, ntra, ema_work1, ema_work2, d_t_con, &  
                d_q_con, d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &  
                itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, pbase, &  
                bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, pmflxr, &  
                pmflxs, da, phi, mp)  
           clwcon0 = qcondc  
           pmfu = upwd + dnwd  
837         ELSE         ELSE
838            ! conema3 ne contient pas les traceurs            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
           CALL conema3(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, &  
                tr_seri, ntra, ema_work1, ema_work2, d_t_con, d_q_con, &  
                d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &  
                itop_con, upwd, dnwd, dnwd0, bas, top, Ma, cape, tvp, rflag, &  
                pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, clwcon0)  
839         ENDIF         ENDIF
840    
841         IF (.NOT. ok_gust) THEN         ! Properties of convective clouds
842            do i = 1, klon         clwcon0 = fact_cldcon * clwcon0
843               wd(i) = 0.0         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
844            enddo              rnebcon0)
845         ENDIF  
846           forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
847         ! Calcul des propriétés des nuages convectifs         mfd = 0.
848           pen_u = 0.
849         DO k = 1, llm         pen_d = 0.
850            DO i = 1, klon         pde_d = 0.
851               zx_t = t_seri(i, k)         pde_u = 0.
852               IF (thermcep) THEN      else
853                  zdelta = MAX(0., SIGN(1., rtt-zx_t))         conv_q = d_q_dyn + d_q_vdf / dtphys
854                  zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)         conv_t = d_t_dyn + d_t_vdf / dtphys
855                  zx_qs = MIN(0.5, zx_qs)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
856                  zcor = 1./(1.-retv*zx_qs)         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
857                  zx_qs = zx_qs*zcor              q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
858               ELSE              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
859                  IF (zx_t < t_coup) THEN              mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
860                     zx_qs = qsats(zx_t)/play(i, k)              kdtop, pmflxr, pmflxs)
861                  ELSE         WHERE (rain_con < 0.) rain_con = 0.
862                     zx_qs = qsatl(zx_t)/play(i, k)         WHERE (snow_con < 0.) snow_con = 0.
863                  ENDIF         ibas_con = llm + 1 - kcbot
864               ENDIF         itop_con = llm + 1 - kctop
865               zqsat(i, k) = zx_qs      END if
           ENDDO  
        ENDDO  
   
        ! calcul des proprietes des nuages convectifs  
        clwcon0 = fact_cldcon*clwcon0  
        call clouds_gno &  
             (klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, rnebcon0)  
     case default  
        print *, "iflag_con non-prevu", iflag_con  
        stop 1  
     END select  
866    
867      DO k = 1, llm      DO k = 1, llm
868         DO i = 1, klon         DO i = 1, klon
# Line 1244  contains Line 873  contains
873         ENDDO         ENDDO
874      ENDDO      ENDDO
875    
876      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
877         ztit = 'after convect'         tit = 'after convect'
878         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
879              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)
880              d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
881         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &              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)  
882      END IF      END IF
883    
884      IF (check) THEN      IF (check) THEN
885         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
886         print *,"aprescon = ", za         print *, "aprescon = ", za
887         zx_t = 0.0         zx_t = 0.
888         za = 0.0         za = 0.
889         DO i = 1, klon         DO i = 1, klon
890            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
891            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
892                 snow_con(i))*airephy(i)/REAL(klon)                 snow_con(i))*airephy(i)/REAL(klon)
893         ENDDO         ENDDO
894         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
895         print *,"Precip = ", zx_t         print *, "Precip = ", zx_t
896      ENDIF      ENDIF
897      IF (zx_ajustq) THEN  
898         DO i = 1, klon      IF (.not. conv_emanuel) THEN
899            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
900         ENDDO         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
        DO k = 1, llm  
           DO i = 1, klon  
              z_apres(i) = z_apres(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
        DO i = 1, klon  
           z_factor(i) = (z_avant(i)-(rain_con(i) + snow_con(i))*dtphys) &  
                /z_apres(i)  
        ENDDO  
901         DO k = 1, llm         DO k = 1, llm
902            DO i = 1, klon            DO i = 1, klon
903               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
# Line 1289  contains Line 906  contains
906            ENDDO            ENDDO
907         ENDDO         ENDDO
908      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
909    
910      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
911    
912      d_t_ajs = 0.      d_t_ajs = 0.
913      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1311  contains Line 927  contains
927              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
928      endif      endif
929    
930      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
931         ztit = 'after dry_adjust'         tit = 'after dry_adjust'
932         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
933              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)  
934      END IF      END IF
935    
936      ! Caclul des ratqs      ! Caclul des ratqs
937    
938      ! ratqs convectifs a l'ancienne en fonction de q(z = 0)-q / q      ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
939      ! on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
940      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
941         do k = 1, llm         do k = 1, llm
942            do i = 1, klon            do i = 1, klon
943               if(ptconv(i, k)) then               if(ptconv(i, k)) then
944                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
945                       +fact_cldcon*(q_seri(i, 1)-q_seri(i, k))/q_seri(i, k)                       * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
946               else               else
947                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
948               endif               endif
# Line 1338  contains Line 953  contains
953      ! ratqs stables      ! ratqs stables
954      do k = 1, llm      do k = 1, llm
955         do i = 1, klon         do i = 1, klon
956            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
957                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
958         enddo         enddo
959      enddo      enddo
960    
961      ! ratqs final      ! ratqs final
962      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
963         ! les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
964         ! ratqs final         ! ratqs final
965         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
966         ! relaxation des ratqs         ! relaxation des ratqs
967         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
968         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
969      else      else
970         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
971         ratqs = ratqss         ratqs = ratqss
972      endif      endif
973    
     ! Processus de condensation à grande echelle et processus de  
     ! précipitation :  
974      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
975           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
976           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
# Line 1376  contains Line 988  contains
988         ENDDO         ENDDO
989      ENDDO      ENDDO
990      IF (check) THEN      IF (check) THEN
991         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
992         print *,"apresilp = ", za         print *, "apresilp = ", za
993         zx_t = 0.0         zx_t = 0.
994         za = 0.0         za = 0.
995         DO i = 1, klon         DO i = 1, klon
996            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
997            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
998                 + snow_lsc(i))*airephy(i)/REAL(klon)                 + snow_lsc(i))*airephy(i)/REAL(klon)
999         ENDDO         ENDDO
1000         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
1001         print *,"Precip = ", zx_t         print *, "Precip = ", zx_t
1002      ENDIF      ENDIF
1003    
1004      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1005         ztit = 'after fisrt'         tit = 'after fisrt'
1006         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1007              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)
1008              d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1009         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &              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)  
1010      END IF      END IF
1011    
1012      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1013    
1014      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1015    
1016      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke      IF (iflag_cldcon <= - 1) THEN
1017           ! seulement pour Tiedtke
1018         snow_tiedtke = 0.         snow_tiedtke = 0.
1019         if (iflag_cldcon == -1) then         if (iflag_cldcon == - 1) then
1020            rain_tiedtke = rain_con            rain_tiedtke = rain_con
1021         else         else
1022            rain_tiedtke = 0.            rain_tiedtke = 0.
1023            do k = 1, llm            do k = 1, llm
1024               do i = 1, klon               do i = 1, klon
1025                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1026                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1027                          *zmasse(i, k)                          *zmasse(i, k)
1028                  endif                  endif
1029               enddo               enddo
# Line 1420  contains Line 1031  contains
1031         endif         endif
1032    
1033         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1034         CALL diagcld1(paprs, play, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1035              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1036         DO k = 1, llm         DO k = 1, llm
1037            DO i = 1, klon            DO i = 1, klon
1038               IF (diafra(i, k) > cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
# Line 1432  contains Line 1042  contains
1042            ENDDO            ENDDO
1043         ENDDO         ENDDO
1044      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1045         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1046         ! convection et du calcul du pas de temps précédent diminué d'un facteur         ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1047         ! facttemps         ! d'un facteur facttemps.
1048         facteur = dtphys *facttemps         facteur = dtphys * facttemps
1049         do k = 1, llm         do k = 1, llm
1050            do i = 1, klon            do i = 1, klon
1051               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1052               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1053                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1054                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1055                  clwcon(i, k) = clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1056               endif               endif
# Line 1467  contains Line 1077  contains
1077      ENDIF      ENDIF
1078    
1079      ! Precipitation totale      ! Precipitation totale
   
1080      DO i = 1, klon      DO i = 1, klon
1081         rain_fall(i) = rain_con(i) + rain_lsc(i)         rain_fall(i) = rain_con(i) + rain_lsc(i)
1082         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
1083      ENDDO      ENDDO
1084    
1085      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1086         ztit = "after diagcld"           dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1087         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &           d_qt, d_ec)
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
1088    
1089      ! Humidité relative pour diagnostic:      ! Humidit\'e relative pour diagnostic :
1090      DO k = 1, llm      DO k = 1, llm
1091         DO i = 1, klon         DO i = 1, klon
1092            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1093            IF (thermcep) THEN            IF (thermcep) THEN
1094               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)  
1095               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1096               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1. - retv*zx_qs)
1097               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
1098            ELSE            ELSE
1099               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
# Line 1503  contains Line 1108  contains
1108      ENDDO      ENDDO
1109    
1110      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
     ! Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)  
1111      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1112         ! Get sulfate aerosol distribution         ! Get sulfate aerosol distribution :
1113         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(dayvrai, time, firstcal, sulfate)
1114         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1115    
        ! Calculate aerosol optical properties (Olivier Boucher)  
1116         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, &
1117              aerindex)              aerindex)
1118      ELSE      ELSE
# Line 1518  contains Line 1121  contains
1121         cg_ae = 0.         cg_ae = 0.
1122      ENDIF      ENDIF
1123    
1124      ! Paramètres optiques des nuages et quelques paramètres pour      ! Param\`etres optiques des nuages et quelques param\`etres pour
1125      ! diagnostics :      ! diagnostics :
1126      if (ok_newmicro) then      if (ok_newmicro) then
1127         CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1128              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1129              fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, &              sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
             re, fl)  
1130      else      else
1131         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1132              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1133              bl95_b1, cldtaupi, re, fl)              bl95_b1, cldtaupi, re, fl)
1134      endif      endif
1135    
1136      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
1137      IF (MOD(itaprad, radpas) == 0) THEN         ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1138         DO i = 1, klon         ! Calcul de l'abedo moyen par maille
1139            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &         albsol = sum(falbe * pctsrf, dim = 2)
1140                 + falbe(i, is_lic) * pctsrf(i, is_lic) &  
1141                 + falbe(i, is_ter) * pctsrf(i, is_ter) &         ! Rayonnement (compatible Arpege-IFS) :
1142                 + falbe(i, is_sic) * pctsrf(i, is_sic)         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1143            albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1144                 + falblw(i, is_lic) * pctsrf(i, is_lic) &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1145                 + falblw(i, is_ter) * pctsrf(i, is_ter) &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1146                 + falblw(i, is_sic) * pctsrf(i, is_sic)              swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1147         ENDDO              solswad, cldtaupi, topswai, solswai)
        ! nouveau rayonnement (compatible Arpege-IFS):  
        CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &  
             albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &  
             heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &  
             sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &  
             lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &  
             cg_ae, topswad, solswad, cldtaupi, topswai, solswai)  
        itaprad = 0  
1148      ENDIF      ENDIF
     itaprad = itaprad + 1  
1149    
1150      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
1151    
1152      DO k = 1, llm      DO k = 1, llm
1153         DO i = 1, klon         DO i = 1, klon
1154            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.
1155         ENDDO         ENDDO
1156      ENDDO      ENDDO
1157    
1158      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1159         ztit = 'after rad'         tit = 'after rad'
1160         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1161              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)
1162              d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1163         call diagphy(airephy, ztit, ip_ebil, topsw, toplw, solsw, sollw, &              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)  
1164      END IF      END IF
1165    
1166      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
1167      DO i = 1, klon      DO i = 1, klon
1168         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1169         zxsnow(i) = 0.0         zxsnow(i) = 0.
1170      ENDDO      ENDDO
1171      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1172         DO i = 1, klon         DO i = 1, klon
# Line 1584  contains Line 1175  contains
1175         ENDDO         ENDDO
1176      ENDDO      ENDDO
1177    
1178      ! Calculer le bilan du sol et la dérive de température (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1179    
1180      DO i = 1, klon      DO i = 1, klon
1181         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1182      ENDDO      ENDDO
1183    
1184      ! Paramétrisation de l'orographie à l'échelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1185    
1186      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1187         ! selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
1188         igwd = 0         igwd = 0
1189         DO i = 1, klon         DO i = 1, klon
1190            itest(i) = 0            itest(i) = 0
1191            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.0)) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1192               itest(i) = 1               itest(i) = 1
1193               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1194            ENDIF            ENDIF
1195         ENDDO         ENDDO
1196    
1197         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1198              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1199              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1200    
1201         ! ajout des tendances         ! ajout des tendances
1202         DO k = 1, llm         DO k = 1, llm
# Line 1619  contains Line 1209  contains
1209      ENDIF      ENDIF
1210    
1211      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1212         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
1213         igwd = 0         igwd = 0
1214         DO i = 1, klon         DO i = 1, klon
1215            itest(i) = 0            itest(i) = 0
1216            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
1217               itest(i) = 1               itest(i) = 1
1218               igwd = igwd + 1               igwd = igwd + 1
              idx(igwd) = i  
1219            ENDIF            ENDIF
1220         ENDDO         ENDDO
1221    
# Line 1644  contains Line 1233  contains
1233         ENDDO         ENDDO
1234      ENDIF      ENDIF
1235    
1236      ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE      ! Stress n\'ecessaires : toute la physique
1237    
1238      DO i = 1, klon      DO i = 1, klon
1239         zustrph(i) = 0.         zustrph(i) = 0.
# Line 1652  contains Line 1241  contains
1241      ENDDO      ENDDO
1242      DO k = 1, llm      DO k = 1, llm
1243         DO i = 1, klon         DO i = 1, klon
1244            zustrph(i) = zustrph(i) + (u_seri(i, k)-u(i, k))/dtphys* zmasse(i, k)            zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1245            zvstrph(i) = zvstrph(i) + (v_seri(i, k)-v(i, k))/dtphys* zmasse(i, k)                 * zmasse(i, k)
1246              zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1247                   * zmasse(i, k)
1248         ENDDO         ENDDO
1249      ENDDO      ENDDO
1250    
1251      !IM calcul composantes axiales du moment angulaire et couple des montagnes      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1252             zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
     CALL aaam_bud(27, klon, llm, time, ra, rg, romega, rlat, rlon, pphis, &  
          zustrdr, zustrli, zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, &  
          aam, torsfc)  
1253    
1254      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1255         ztit = 'after orography'           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1256         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &           d_qt, d_ec)
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
1257    
1258      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1259      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, &      call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &
1260           nqmx-2, dtphys, u, t, paprs, play, pmfu, pmfd, pen_u, pde_u, &           play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &
1261           pen_d, pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &           yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
1262           frac_impa, frac_nucl, pphis, albsol, rhcl, cldfra, rneb, &           tr_seri, zmasse, ncid_startphy)
1263           diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, &  
1264           tr_seri, zmasse)      IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1265             pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1266      IF (offline) THEN           frac_impa, frac_nucl, pphis, airephy, dtphys)
        call phystokenc(dtphys, rlon, rlat, t, pmfu, pmfd, pen_u, pde_u, &  
             pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &  
             pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)  
     ENDIF  
1267    
1268      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1269      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)  
1270    
1271      ! diag. bilKP      ! diag. bilKP
1272    
1273      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, &
1274           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1275    
1276      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
# Line 1706  contains Line 1286  contains
1286         END DO         END DO
1287      END DO      END DO
1288    
1289      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1290         ztit = 'after physic'         tit = 'after physic'
1291         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1292              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)
1293              d_ql, d_qs, d_ec)         ! Comme les tendances de la physique sont ajoute dans la dynamique,
        ! Comme les tendances de la physique sont ajoute dans la dynamique,  
1294         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1295         ! 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.
1296         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1297         call diagphy(airephy, ztit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1298              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)  
   
1299         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1300      END IF      END IF
1301    
1302      ! SORTIES      ! SORTIES
1303    
1304      !cc prw = eau precipitable      ! prw = eau precipitable
1305      DO i = 1, klon      DO i = 1, klon
1306         prw(i) = 0.         prw(i) = 0.
1307         DO k = 1, llm         DO k = 1, llm
# Line 1745  contains Line 1321  contains
1321         ENDDO         ENDDO
1322      ENDDO      ENDDO
1323    
1324      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1325         DO iq = 3, nqmx         DO k = 1, llm
1326            DO k = 1, llm            DO i = 1, klon
1327               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  
1328            ENDDO            ENDDO
1329         ENDDO         ENDDO
1330      ENDIF      ENDDO
1331    
1332      ! 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:
1333      DO k = 1, llm      DO k = 1, llm
# Line 1763  contains Line 1337  contains
1337         ENDDO         ENDDO
1338      ENDDO      ENDDO
1339    
1340      ! Ecriture des sorties      CALL histwrite_phy("phis", pphis)
1341      call write_histhf      CALL histwrite_phy("aire", airephy)
1342      call write_histday      CALL histwrite_phy("psol", paprs(:, 1))
1343      call write_histins      CALL histwrite_phy("precip", rain_fall + snow_fall)
1344        CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1345      ! Si c'est la fin, il faut conserver l'etat de redemarrage      CALL histwrite_phy("pluc", rain_con + snow_con)
1346      IF (lafin) THEN      CALL histwrite_phy("tsol", zxtsol)
1347         itau_phy = itau_phy + itap      CALL histwrite_phy("t2m", zt2m)
1348         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &      CALL histwrite_phy("q2m", zq2m)
1349              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &      CALL histwrite_phy("u10m", zu10m)
1350              rain_fall, snow_fall, solsw, sollwdown, dlw, radsol, frugs, &      CALL histwrite_phy("v10m", zv10m)
1351              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &      CALL histwrite_phy("snow", snow_fall)
1352              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)      CALL histwrite_phy("cdrm", cdragm)
1353      ENDIF      CALL histwrite_phy("cdrh", cdragh)
1354        CALL histwrite_phy("topl", toplw)
1355      firstcal = .FALSE.      CALL histwrite_phy("evap", evap)
1356        CALL histwrite_phy("sols", solsw)
1357    contains      CALL histwrite_phy("soll", sollw)
1358        CALL histwrite_phy("solldown", sollwdown)
1359      subroutine write_histday      CALL histwrite_phy("bils", bils)
1360        CALL histwrite_phy("sens", - sens)
1361        use gr_phy_write_3d_m, only: gr_phy_write_3d      CALL histwrite_phy("fder", fder)
1362        integer itau_w ! pas de temps ecriture      CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1363        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1364        !------------------------------------------------      CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1365        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
       if (ok_journe) THEN  
          itau_w = itau_phy + itap  
          if (nqmx <= 4) then  
             call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &  
                  gr_phy_write_3d(wo) * 1e3)  
             ! (convert "wo" from kDU to DU)  
          end if  
          if (ok_sync) then  
             call histsync(nid_day)  
          endif  
       ENDIF  
   
     End subroutine write_histday  
   
     !****************************  
   
     subroutine write_histhf  
   
       ! From phylmd/write_histhf.h, version 1.5 2005/05/25 13:10:09  
   
       !------------------------------------------------  
   
       call write_histhf3d  
   
       IF (ok_sync) THEN  
          call histsync(nid_hf)  
       ENDIF  
   
     end subroutine write_histhf  
   
     !***************************************************************  
   
     subroutine write_histins  
   
       ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09  
   
       real zout  
       integer itau_w ! pas de temps ecriture  
   
       !--------------------------------------------------  
   
       IF (ok_instan) THEN  
          ! Champs 2D:  
   
          zsto = dtphys * ecrit_ins  
          zout = dtphys * ecrit_ins  
          itau_w = itau_phy + itap  
   
          i = NINT(zout/zsto)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)  
          CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)  
   
          i = NINT(zout/zsto)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)  
          CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = paprs(i, 1)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)  
1366    
1367           DO i = 1, klon      DO nsrf = 1, nbsrf
1368              zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)         CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)
1369           ENDDO         CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1370           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)         CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1371           CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)         CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1372           CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1373           DO i = 1, klon         CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1374              zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)         CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1375           ENDDO         CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1376           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)         CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1377           CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)      END DO
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)  
          !ccIM  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)  
          CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)  
          CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)  
          CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)  
          CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)  
   
          zx_tmp_fi2d(1:klon) = -1*sens(1:klon)  
          ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)  
          CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)  
   
          DO nsrf = 1, nbsrf  
             !XXX  
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
          END DO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)  
          CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)  
   
          !HBTM2  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)  
   
          ! Champs 3D:  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)  
          CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)  
          CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)  
   
          if (ok_sync) then  
             call histsync(nid_ins)  
          endif  
       ENDIF  
   
     end subroutine write_histins  
   
     !****************************************************  
   
     subroutine write_histhf3d  
   
       ! From phylmd/write_histhf3d.h, version 1.2 2005/05/25 13:10:09  
   
       integer itau_w ! pas de temps ecriture  
   
       !-------------------------------------------------------  
   
       itau_w = itau_phy + itap  
   
       ! Champs 3D:  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, qx(1, 1, ivap), zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)  
   
       if (nbtr >= 3) then  
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, tr_seri(1, 1, 3), &  
               zx_tmp_3d)  
          CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)  
       end if  
1378    
1379        if (ok_sync) then      CALL histwrite_phy("albs", albsol)
1380           call histsync(nid_hf3d)      CALL histwrite_phy("rugs", zxrugs)
1381        endif      CALL histwrite_phy("s_pblh", s_pblh)
1382        CALL histwrite_phy("s_pblt", s_pblt)
1383        CALL histwrite_phy("s_lcl", s_lcl)
1384        CALL histwrite_phy("s_capCL", s_capCL)
1385        CALL histwrite_phy("s_oliqCL", s_oliqCL)
1386        CALL histwrite_phy("s_cteiCL", s_cteiCL)
1387        CALL histwrite_phy("s_therm", s_therm)
1388        CALL histwrite_phy("s_trmb1", s_trmb1)
1389        CALL histwrite_phy("s_trmb2", s_trmb2)
1390        CALL histwrite_phy("s_trmb3", s_trmb3)
1391        if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1392        CALL histwrite_phy("temp", t_seri)
1393        CALL histwrite_phy("vitu", u_seri)
1394        CALL histwrite_phy("vitv", v_seri)
1395        CALL histwrite_phy("geop", zphi)
1396        CALL histwrite_phy("pres", play)
1397        CALL histwrite_phy("dtvdf", d_t_vdf)
1398        CALL histwrite_phy("dqvdf", d_q_vdf)
1399        CALL histwrite_phy("rhum", zx_rh)
1400    
1401        if (ok_instan) call histsync(nid_ins)
1402    
1403        IF (lafin) then
1404           call NF95_CLOSE(ncid_startphy)
1405           CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1406                fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1407                radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1408                t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1409                w01)
1410        end IF
1411    
1412      end subroutine write_histhf3d      firstcal = .FALSE.
1413    
1414    END SUBROUTINE physiq    END SUBROUTINE physiq
1415    

Legend:
Removed from v.53  
changed lines
  Added in v.191

  ViewVC Help
Powered by ViewVC 1.1.21