/[lmdze]/trunk/phylmd/physiq.f90
ViewVC logotype

Diff of /trunk/phylmd/physiq.f90

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

trunk/libf/phylmd/physiq.f90 revision 56 by guez, Tue Jan 10 19:02:02 2012 UTC trunk/phylmd/physiq.f revision 322 by guez, Thu Jan 24 16:17:39 2019 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      use aaam_bud_m, only: aaam_bud
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
19      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
     USE calendar, ONLY: ymds2ju  
20      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
21      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &      USE clesphys, ONLY: cdhmax, cdmmax, ecrit_ins, ok_instan
22           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin      USE clesphys2, ONLY: conv_emanuel, nbapp_rad, new_oliq, ok_orodr, ok_orolf
23      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &      USE conf_interface_m, ONLY: conf_interface
24           ok_orodr, ok_orolf, soil_model      USE pbl_surface_m, ONLY: pbl_surface
25      USE clmain_m, ONLY: clmain      use clouds_gno_m, only: clouds_gno
26      USE comgeomphy, ONLY: airephy, cuphy, cvphy      use comconst, only: dtphys
27        USE comgeomphy, ONLY: airephy
28      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
29      USE conf_gcm_m, ONLY: offline, raz_date      USE conf_gcm_m, ONLY: lmt_pas
30      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
31        use conflx_m, only: conflx
32      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
33      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
34      use diagetpq_m, only: diagetpq      USE dimensions, ONLY: llm, nqmx
35      USE dimens_m, ONLY: iim, jjm, llm, nqmx      USE dimphy, ONLY: klon
     USE dimphy, ONLY: klon, nbtr  
36      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
37      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
38      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      use dynetat0_chosen_m, only: day_ref, annee_ref
39        USE fcttre, ONLY: foeew
40        use fisrtilp_m, only: fisrtilp
41      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
42      USE histcom, ONLY: histsync      USE histsync_m, ONLY: histsync
43      USE histwrite_m, ONLY: histwrite      USE histwrite_phy_m, ONLY: histwrite_phy
44      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, &
45           nbsrf           nbsrf
46      USE ini_histhf_m, ONLY: ini_histhf      USE ini_histins_m, ONLY: ini_histins, nid_ins
47      USE ini_histday_m, ONLY: ini_histday      use lift_noro_m, only: lift_noro
48      USE ini_histins_m, ONLY: ini_histins      use netcdf95, only: NF95_CLOSE
49      USE oasis_m, ONLY: ok_oasis      use newmicro_m, only: newmicro
50      USE orbite_m, ONLY: orbite, zenang      use nr_util, only: assert
51        use nuage_m, only: nuage
52        USE orbite_m, ONLY: orbite
53      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
54      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0
55      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
56      USE phystokenc_m, ONLY: phystokenc      USE phyredem0_m, ONLY: phyredem0
57      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
     USE qcheck_m, ONLY: qcheck  
58      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
59      use sugwd_m, only: sugwd      use yoegwd, only: sugwd
60      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt, rmo3, md
61      USE temps, ONLY: annee_ref, day_ref, itau_phy      use time_phylmdz, only: itap, increment_itap
62        use transp_m, only: transp
63        use transp_lay_m, only: transp_lay
64        use unit_nml_m, only: unit_nml
65        USE ymds2ju_m, ONLY: ymds2ju
66      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
67        use zenang_m, only: zenang
68    
     ! Arguments:  
   
     REAL, intent(in):: rdayvrai  
     ! (elapsed time since January 1st 0h of the starting year, in days)  
   
     REAL, intent(in):: time ! heure de la journée en fraction de jour  
     REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)  
69      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
70    
71      REAL, intent(in):: paprs(klon, llm + 1)      integer, intent(in):: dayvrai
72      ! (pression pour chaque inter-couche, en Pa)      ! current day number, based at value 1 on January 1st of annee_ref
73    
74      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
     ! (input pression pour le mileu de chaque couche (en Pa))  
75    
76      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
77      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! pression pour chaque inter-couche, en Pa
78    
79      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: play(:, :) ! (klon, llm)
80        ! pression pour le mileu de chaque couche (en Pa)
81    
82      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
83      ! vitesse dans la direction X (de O a E) en m/s      ! géopotentiel de chaque couche (référence sol)
84    
85      REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
     REAL, intent(in):: t(klon, llm) ! input temperature (K)  
86    
87      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: u(:, :) ! (klon, llm)
88      ! (humidité spécifique et fractions massiques des autres traceurs)      ! vitesse dans la direction X (de O a E) en m / s
89    
90      REAL omega(klon, llm) ! input vitesse verticale en Pa/s      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m / s
91      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
     REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/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  
92    
93      LOGICAL:: firstcal = .true.      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
94        ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
95    
96      INTEGER nbteta      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa / s
97      PARAMETER(nbteta = 3)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
98        REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
99        REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K / s)
100    
101      REAL PVteta(klon, nbteta)      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
102      ! (output vorticite potentielle a des thetas constantes)      ! tendance physique de "qx" (s-1)
103    
104      LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE      ! Local:
     PARAMETER (ok_cvl = .TRUE.)  
     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface  
     PARAMETER (ok_gust = .FALSE.)  
105    
106      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL:: firstcal = .true.
     PARAMETER (check = .FALSE.)  
107    
108      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
109      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
110    
111      ! Parametres lies au coupleur OASIS:      ! pour phystoke avec thermiques
     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.)  
   
     ! pour phsystoke avec thermiques  
112      REAL fm_therm(klon, llm + 1)      REAL fm_therm(klon, llm + 1)
113      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
114      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
115    
116      INTEGER ivap ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
117      PARAMETER (ivap = 1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     INTEGER iliq ! indice de traceurs pour eau liquide  
     PARAMETER (iliq = 2)  
118    
119      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
120      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
121    
122      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K / s)
123      REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)      REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg / kg / s)
124    
125      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
126    
127      !IM Amip2 PV a theta constante      REAL, save:: swdn0(klon, llm + 1), swdn(klon, llm + 1)
128        REAL, save:: swup0(klon, llm + 1), swup(klon, llm + 1)
129    
130      CHARACTER(LEN = 3) ctetaSTD(nbteta)      REAL, save:: lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
131      DATA ctetaSTD/'350', '380', '405'/      REAL, save:: lwup0(klon, llm + 1), lwup(klon, llm + 1)
     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)  
     SAVE swdn0, swdn, swup0, swup  
   
     REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)  
     REAL lwup0(klon, klevp1), lwup(klon, klevp1)  
     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 '/  
132    
133      ! prw: precipitable water      ! prw: precipitable water
134      real prw(klon)      real prw(klon)
135    
136      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg / m2)
137      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg / kg)
138      REAL flwp(klon), fiwp(klon)      REAL flwp(klon), fiwp(klon)
139      REAL flwc(klon, llm), fiwc(klon, llm)      REAL flwc(klon, llm), fiwc(klon, llm)
140    
     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  
   
141      ! Variables propres a la physique      ! Variables propres a la physique
142    
143      INTEGER, save:: radpas      INTEGER, save:: radpas
144      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
145      ! "physiq".)      ! "physiq".
146    
147      REAL radsol(klon)      REAL, save:: radsol(klon)
148      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      ! bilan radiatif net au sol (W/m2), positif vers le bas
149        
150      INTEGER, SAVE:: itap ! number of calls to "physiq"      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction, in K
   
     REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction  
151    
152      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
153      ! soil temperature of surface fraction      ! soil temperature of surface fraction
154    
155      REAL fevap(klon, nbsrf)      REAL fluxlat(klon, nbsrf) ! flux de chaleur latente, en W m-2
     SAVE fevap ! evaporation  
     REAL fluxlat(klon, nbsrf)  
     SAVE fluxlat  
   
     REAL fqsurf(klon, nbsrf)  
     SAVE fqsurf ! humidite de l'air au contact de la surface  
156    
157      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol      REAL, save:: fqsurf(klon, nbsrf)
158        ! humidite de l'air au contact de la surface
159    
160      REAL fsnow(klon, nbsrf)      REAL, save:: qsol(klon) ! column-density of water in soil, in kg m-2
161      SAVE fsnow ! epaisseur neigeuse      REAL, save:: fsnow(klon, nbsrf) ! \'epaisseur neigeuse
162        REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
163    
164      REAL falbe(klon, nbsrf)      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
     SAVE falbe ! albedo par type de surface  
     REAL falblw(klon, nbsrf)  
     SAVE falblw ! albedo par type de surface  
   
     ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :  
165      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
166      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
167      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 295  contains Line 170  contains
170      REAL, save:: zpic(klon) ! Maximum de l'OESM      REAL, save:: zpic(klon) ! Maximum de l'OESM
171      REAL, save:: zval(klon) ! Minimum de l'OESM      REAL, save:: zval(klon) ! Minimum de l'OESM
172      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
   
173      REAL zulow(klon), zvlow(klon)      REAL zulow(klon), zvlow(klon)
174        INTEGER ktest(klon)
175    
176      INTEGER igwd, idx(klon), itest(klon)      REAL, save:: agesno(klon, nbsrf) ! age de la neige
177        REAL, save:: run_off_lic_0(klon)
178    
179      REAL agesno(klon, nbsrf)      ! Variables li\'ees \`a la convection d'Emanuel :
180      SAVE agesno ! age de la neige      REAL, save:: Ma(klon, llm) ! undilute upward mass flux
181        REAL, save:: sig1(klon, llm), w01(klon, llm)
182    
183      REAL run_off_lic_0(klon)      ! Variables pour la couche limite (Alain Lahellec) :
184      SAVE run_off_lic_0      REAL cdragh(klon) ! drag coefficient pour T and Q
185      !KE43      REAL cdragm(klon) ! drag coefficient pour vent
     ! Variables liees a la convection de K. Emanuel (sb):  
186    
187      REAL bas, top ! cloud base and top levels      REAL coefh(klon, 2:llm) ! coef d'echange pour phytrac
     SAVE bas  
     SAVE top  
188    
189      REAL Ma(klon, llm) ! undilute upward mass flux      REAL, save:: ffonte(klon, nbsrf)
190      SAVE Ma      ! flux thermique utilise pour fondre la neige
     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  
191    
192      REAL wd(klon) ! sb      REAL fqcalving(klon, nbsrf)
193      SAVE wd ! sb      ! flux d'eau "perdue" par la surface et n\'ecessaire pour limiter
194        ! la hauteur de neige, en kg / m2 / s
195    
196      ! Variables locales pour la couche limite (al1):      REAL zxffonte(klon)
197    
198      ! Variables locales:      REAL, save:: pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
199        REAL, save:: pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
200    
201      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL, save:: pfrac_1nucl(klon, llm)
202      REAL cdragm(klon) ! drag coefficient pour vent      ! Produits des coefs lessi nucl (alpha = 1)
203    
204      !AA Pour phytrac      REAL frac_impa(klon, llm) ! fraction d'a\'erosols lessiv\'es (impaction)
     REAL ycoefh(klon, llm) ! coef d'echange pour phytrac  
     REAL yu1(klon) ! vents dans la premiere couche U  
     REAL yv1(klon) ! vents dans la premiere couche V  
     REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige  
     REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface  
     ! !et necessaire pour limiter la  
     ! !hauteur de neige, en kg/m2/s  
     REAL zxffonte(klon), zxfqcalving(klon)  
   
     REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction  
     save pfrac_impa  
     REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation  
     save pfrac_nucl  
     REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)  
     save pfrac_1nucl  
     REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)  
205      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
206    
207      !AA      REAL, save:: rain_fall(klon)
208      REAL rain_fall(klon) ! pluie      ! liquid water mass flux (kg / m2 / s), positive down
209      REAL snow_fall(klon) ! neige  
210      save snow_fall, rain_fall      REAL, save:: snow_fall(klon)
211      !IM cf FH pour Tiedtke 080604      ! solid water mass flux (kg / m2 / s), positive down
212    
213      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
214    
215      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon) ! flux d'\'evaporation au sol
216      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      real dflux_q(klon) ! derivative of the evaporation flux at the surface
217      REAL dlw(klon) ! derivee infra rouge      REAL sens(klon) ! flux de chaleur sensible au sol
218      SAVE dlw      real dflux_t(klon) ! derivee du flux de chaleur sensible au sol
219        REAL, save:: dlw(klon) ! derivative of infra-red flux
220      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
221      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL fder(klon) ! Derive de flux (sensible et latente)
     save fder  
222      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
223      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
224      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
225      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
226    
227      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
228      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
229    
230      ! Conditions aux limites      ! Conditions aux limites
231    
232      INTEGER julien      INTEGER julien
233        REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
234      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      REAL, save:: albsol(klon) ! albedo du sol total, visible, moyen par maille
     REAL pctsrf(klon, nbsrf)  
     !IM  
     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE  
   
     SAVE pctsrf ! sous-fraction du sol  
     REAL albsol(klon)  
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
235      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
236        real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
237    
238      ! Declaration des procedures appelees      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
239        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  
240    
241      REAL rhcl(klon, llm) ! humiditi relative ciel clair      REAL rhcl(klon, llm) ! humidit\'e relative ciel clair
242      REAL dialiq(klon, llm) ! eau liquide nuageuse      REAL dialiq(klon, llm) ! eau liquide nuageuse
243      REAL diafra(klon, llm) ! fraction nuageuse      REAL diafra(klon, llm) ! fraction nuageuse
244      REAL cldliq(klon, llm) ! eau liquide nuageuse      REAL cldliq(klon, llm) ! eau liquide nuageuse
# Line 412  contains Line 246  contains
246      REAL cldtau(klon, llm) ! epaisseur optique      REAL cldtau(klon, llm) ! epaisseur optique
247      REAL cldemi(klon, llm) ! emissivite infrarouge      REAL cldemi(klon, llm) ! emissivite infrarouge
248    
249      REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite      REAL flux_q(klon, nbsrf) ! flux turbulent d'humidite à la surface
     REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur  
     REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u  
     REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v  
   
     REAL zxfluxt(klon, llm)  
     REAL zxfluxq(klon, llm)  
     REAL zxfluxu(klon, llm)  
     REAL zxfluxv(klon, llm)  
250    
251      ! Le rayonnement n'est pas calcule tous les pas, il faut donc      REAL flux_t(klon, nbsrf)
252      ! que les variables soient rémanentes      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
253      REAL, save:: heat(klon, llm) ! chauffage solaire      ! vers le bas) à la surface
     REAL heat0(klon, llm) ! chauffage solaire ciel clair  
     REAL cool(klon, llm) ! refroidissement infrarouge  
     REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair  
     REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)  
     real sollwdown(klon) ! downward LW flux at surface  
     REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)  
     REAL albpla(klon)  
     REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface  
     REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface  
     SAVE cool, albpla, topsw, toplw, solsw, sollw, sollwdown  
     SAVE topsw0, toplw0, solsw0, sollw0, heat0, cool0  
   
     INTEGER itaprad  
     SAVE itaprad  
   
     REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)  
     REAL conv_t(klon, llm) ! convergence of temperature (K/s)  
   
     REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut  
     REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree  
   
     REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)  
   
     REAL dist, rmu0(klon), fract(klon)  
     REAL zdtime ! pas de temps du rayonnement (s)  
     real zlongi  
254    
255      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL flux_u(klon, nbsrf), flux_v(klon, nbsrf)
256      LOGICAL zx_ajustq      ! tension du vent (flux turbulent de vent) à la surface, en Pa
257    
258      REAL za, zb      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
259      REAL zx_t, zx_qs, zdelta, zcor      ! les variables soient r\'emanentes.
260        REAL, save:: heat(klon, llm) ! chauffage solaire
261        REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
262        REAL, save:: cool(klon, llm) ! refroidissement infrarouge
263        REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
264        REAL, save:: topsw(klon), toplw(klon), solsw(klon)
265    
266        REAL, save:: sollw(klon) ! surface net downward longwave flux, in W m-2
267        real, save:: sollwdown(klon) ! downward LW flux at surface
268        REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
269        REAL, save:: albpla(klon)
270    
271        REAL conv_q(klon, llm) ! convergence de l'humidite (kg / kg / s)
272        REAL conv_t(klon, llm) ! convergence of temperature (K / s)
273    
274        REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
275        REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
276    
277        REAL zxfluxlat(klon)
278        REAL dist, mu0(klon), fract(klon)
279        real longi
280        REAL z_avant(klon), z_apres(klon), z_factor(klon)
281        REAL zb
282        REAL zx_t, zx_qs, zcor
283      real zqsat(klon, llm)      real zqsat(klon, llm)
284      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
     REAL t_coup  
     PARAMETER (t_coup = 234.0)  
   
285      REAL zphi(klon, llm)      REAL zphi(klon, llm)
286    
287      !IM cf. AM Variables locales pour la CLA (hbtm2)      ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
288    
289      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
290      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
291      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
292      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
293      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
294      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T \`a la hauteur de couche limite
295      REAL, SAVE:: therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
296      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape      ! Grandeurs de sorties
     REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition  
     REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega  
     ! Grdeurs de sorties  
297      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
298      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
299      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon)
     REAL s_trmb3(klon)  
300    
301      ! Variables locales pour la convection de K. Emanuel (sb):      ! Variables pour la convection de K. Emanuel :
302    
303      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
304      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
305      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux      REAL, save:: cape(klon)
306      REAL tvp(klon, llm) ! virtual temp of lifted parcel  
     REAL cape(klon) ! CAPE  
     SAVE cape  
   
     REAL pbase(klon) ! cloud base pressure  
     SAVE pbase  
     REAL bbase(klon) ! cloud base buoyancy  
     SAVE bbase  
     REAL rflag(klon) ! flag fonctionnement de convect  
307      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)  
308    
309      ! Variables du changement      ! Variables du changement
310    
311      ! con: convection      ! con: convection
312      ! lsc: large scale condensation      ! lsc: large scale condensation
313      ! ajs: ajustement sec      ! ajs: ajustement sec
314      ! eva: évaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
315      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
316      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
317      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL, save:: d_u_con(klon, llm), d_v_con(klon, llm)
318      REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)      REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
319      REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)      REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
320      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
321      REAL rneb(klon, llm)      REAL rneb(klon, llm)
322    
323      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
324      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
325      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
326      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
327      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
328      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
329    
330      INTEGER,save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
331        real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
332    
333      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon)
334      REAL snow_con(klon), snow_lsc(klon)      real rain_lsc(klon)
335      REAL d_ts(klon, nbsrf)      REAL snow_con(klon) ! neige (mm / s)
336        real snow_lsc(klon)
337        REAL d_ts(klon, nbsrf) ! variation of ftsol
338    
339      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
340      REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)      REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
# Line 537  contains Line 344  contains
344      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
345      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
346    
347      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
348      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
349      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
350    
351      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
352      real, save:: fact_cldcon      real:: fact_cldcon = 0.375
353      real, save:: facttemps      real:: facttemps = 1.e-4
354      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
355      real facteur      real facteur
356    
357      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
358      logical ptconv(klon, llm)      logical ptconv(klon, llm)
359    
360      ! Variables locales pour effectuer les appels en série :      ! Variables pour effectuer les appels en s\'erie :
361    
362      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
363      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
364      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
365        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
366    
367      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
368    
369      REAL zustrdr(klon), zvstrdr(klon)      REAL zustrdr(klon), zvstrdr(klon)
370      REAL zustrli(klon), zvstrli(klon)      REAL zustrli(klon), zvstrli(klon)
     REAL zustrph(klon), zvstrph(klon)  
371      REAL aam, torsfc      REAL aam, torsfc
372    
     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  
   
373      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.
374      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.
375      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.
376      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.
377    
378      REAL zsto      REAL tsol(klon)
379    
380        REAL d_t_ec(klon, llm)
381        ! tendance due \`a la conversion d'\'energie cin\'etique en
382        ! énergie thermique
383    
384        REAL, save:: t2m(klon, nbsrf), q2m(klon, nbsrf)
385        ! temperature and humidity at 2 m
386    
387        REAL, save:: u10m_srf(klon, nbsrf), v10m_srf(klon, nbsrf)
388        ! composantes du vent \`a 10 m
389        
390        REAL zt2m(klon), zq2m(klon) ! température, humidité 2 m moyenne sur 1 maille
391        REAL u10m(klon), v10m(klon) ! vent \`a 10 m moyenn\' sur les sous-surfaces
392    
393      character(len = 20) modname      ! Aerosol effects:
     character(len = 80) abort_message  
     logical ok_sync  
     real date0  
   
     ! Variables liées au bilan d'énergie et d'enthalpie :  
     REAL ztsol(klon)  
     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec  
     REAL, SAVE:: d_h_vcol_phy  
     REAL fs_bound, fq_bound  
     REAL zero_v(klon)  
     CHARACTER(LEN = 15) ztit  
     INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics  
     INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics  
   
     REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique  
     REAL ZRCPD  
   
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m  
     REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille  
     REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille  
     !jq Aerosol effects (Johannes Quaas, 27/11/2003)  
     REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]  
   
     REAL, save:: sulfate_pi(klon, llm)  
     ! (SO4 aerosol concentration, in ug/m3, pre-industrial value)  
   
     REAL cldtaupi(klon, llm)  
     ! (Cloud optical thickness for pre-industrial (pi) aerosols)  
   
     REAL re(klon, llm) ! Cloud droplet effective radius  
     REAL fl(klon, llm) ! denominator of re  
   
     ! Aerosol optical properties  
     REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)  
     REAL cg_ae(klon, llm, 2)  
   
     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.  
     ! ok_ade = True -ADE = topswad-topsw  
   
     REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.  
     ! ok_aie = True ->  
     ! ok_ade = True -AIE = topswai-topswad  
     ! ok_ade = F -AIE = topswai-topsw  
   
     REAL aerindex(klon) ! POLDER aerosol index  
   
     ! Parameters  
     LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not  
     REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)  
   
     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1  
     SAVE u10m  
     SAVE v10m  
     SAVE t2m  
     SAVE q2m  
     SAVE ffonte  
     SAVE fqcalving  
     SAVE piz_ae  
     SAVE tau_ae  
     SAVE cg_ae  
     SAVE rain_con  
     SAVE snow_con  
     SAVE topswai  
     SAVE topswad  
     SAVE solswai  
     SAVE solswad  
     SAVE d_u_con  
     SAVE d_v_con  
     SAVE rnebcon0  
     SAVE clwcon0  
394    
395      real zmasse(klon, llm)      REAL, save:: topswad(klon), solswad(klon) ! aerosol direct effect
396        LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
397    
398        REAL:: bl95_b0 = 2., bl95_b1 = 0.2
399        ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
400        ! B). They link cloud droplet number concentration to aerosol mass
401        ! concentration.
402    
403        real zmasse(klon, llm)
404      ! (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)
405    
406      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      integer, save:: ncid_startphy
407    
408        namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, &
409             ratqsbas, ratqshaut, ok_ade, bl95_b0, bl95_b1, iflag_thermals, &
410             nsplit_thermals
411    
412      !----------------------------------------------------------------      !----------------------------------------------------------------
413    
414      modname = 'physiq'      IF (nqmx < 2) CALL abort_gcm('physiq', &
415      IF (if_ebil >= 1) THEN           'eaux vapeur et liquide sont indispensables')
        DO i = 1, klon  
           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  
416    
417      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
418         ! initialiser         ! initialiser
419         u10m = 0.         u10m_srf = 0.
420         v10m = 0.         v10m_srf = 0.
421         t2m = 0.         t2m = 0.
422         q2m = 0.         q2m = 0.
423         ffonte = 0.         ffonte = 0.
424         fqcalving = 0.         d_u_con = 0.
425         piz_ae = 0.         d_v_con = 0.
426         tau_ae = 0.         rnebcon0 = 0.
427         cg_ae = 0.         clwcon0 = 0.
428         rain_con(:) = 0.         rnebcon = 0.
429         snow_con(:) = 0.         clwcon = 0.
        bl95_b0 = 0.  
        bl95_b1 = 0.  
        topswai(:) = 0.  
        topswad(:) = 0.  
        solswai(:) = 0.  
        solswad(:) = 0.  
   
        d_u_con = 0.0  
        d_v_con = 0.0  
        rnebcon0 = 0.0  
        clwcon0 = 0.0  
        rnebcon = 0.0  
        clwcon = 0.0  
   
430         pblh =0. ! Hauteur de couche limite         pblh =0. ! Hauteur de couche limite
431         plcl =0. ! Niveau de condensation de la CLA         plcl =0. ! Niveau de condensation de la CLA
432         capCL =0. ! CAPE de couche limite         capCL =0. ! CAPE de couche limite
433         oliqCL =0. ! eau_liqu integree de couche limite         oliqCL =0. ! eau_liqu integree de couche limite
434         cteiCL =0. ! cloud top instab. crit. couche limite         cteiCL =0. ! cloud top instab. crit. couche limite
435         pblt =0. ! T a la Hauteur de couche limite         pblt =0.
436         therm =0.         therm =0.
437         trmb1 =0. ! deep_cape  
438         trmb2 =0. ! inhibition         iflag_thermals = 0
439         trmb3 =0. ! Point Omega         nsplit_thermals = 1
440           print *, "Enter namelist 'physiq_nml'."
441         IF (if_ebil >= 1) d_h_vcol_phy = 0.         read(unit=*, nml=physiq_nml)
442           write(unit_nml, nml=physiq_nml)
443         ! appel a la lecture du run.def physique  
444           call conf_phys
        call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &  
             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)  
445    
446         ! Initialiser les compteurs:         ! Initialiser les compteurs:
447    
448         frugs = 0.         frugs = 0.
449         itap = 0         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
450         itaprad = 0              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, &
451         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &              zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
452              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01, &
453              snow_fall, solsw, sollwdown, dlw, radsol, frugs, agesno, zmea, &              ncid_startphy)
             zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &  
             ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)  
454    
455         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
456         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  
   
        IF(ocean.NE.'force ') THEN  
           ok_ocean = .TRUE.  
        ENDIF  
457    
458         CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &         radpas = lmt_pas / nbapp_rad
459              ok_region)         print *, "radpas = ", radpas
   
        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 "  
   
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG  
           DO i = 1, klon  
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END  
460    
461           ! Initialisation pour le sch\'ema de convection d'Emanuel :
462           IF (conv_emanuel) THEN
463              ibas_con = 1
464              itop_con = 1
465         ENDIF         ENDIF
466    
467         IF (ok_orodr) THEN         IF (ok_orodr) THEN
# Line 783  contains Line 471  contains
471            rugoro = 0.            rugoro = 0.
472         ENDIF         ENDIF
473    
        lmt_pas = NINT(86400. / dtphys) ! tous les jours  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
        ecrit_ins = NINT(ecrit_ins/dtphys)  
        ecrit_hf = NINT(ecrit_hf/dtphys)  
        ecrit_mth = NINT(ecrit_mth/dtphys)  
        ecrit_tra = NINT(86400.*ecrit_tra/dtphys)  
        ecrit_reg = NINT(ecrit_reg/dtphys)  
   
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
        print *,'AVANT HIST IFLAG_CON = ', iflag_con  
   
474         ! Initialisation des sorties         ! Initialisation des sorties
475           call ini_histins(ok_newmicro)
476         call ini_histhf(dtphys, nid_hf, nid_hf3d)         CALL phyredem0
477         call ini_histday(dtphys, ok_journe, nid_day, nqmx)         call conf_interface
        call ini_histins(dtphys, ok_instan, nid_ins)  
        CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)  
        !XXXPB Positionner date0 pour initialisation de ORCHIDEE  
        WRITE(*, *) 'physiq date0: ', date0  
478      ENDIF test_firstcal      ENDIF test_firstcal
479    
480      ! Mettre a zero des variables de sortie (pour securite)      ! We will modify variables *_seri and we will not touch variables
481        ! u, v, t, qx:
482        t_seri = t
483        u_seri = u
484        v_seri = v
485        q_seri = qx(:, :, ivap)
486        ql_seri = qx(:, :, iliq)
487        tr_seri = qx(:, :, 3:nqmx)
488    
489      DO i = 1, klon      tsol = sum(ftsol * pctsrf, dim = 2)
        d_ps(i) = 0.0  
     ENDDO  
     DO iq = 1, nqmx  
        DO k = 1, llm  
           DO i = 1, klon  
              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  
   
     DO i = 1, klon  
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
   
     IF (if_ebil >= 1) THEN  
        ztit = 'after dynamics'  
        CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        ! Comme les tendances de la physique sont ajoutés dans la  
        !  dynamique, la variation d'enthalpie par la dynamique devrait  
        !  être égale à la variation de la physique au pas de temps  
        !  précédent.  Donc la somme de ces 2 variations devrait être  
        !  nulle.  
        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_h_vcol_phy, &  
             d_qt, 0., fs_bound, fq_bound)  
     END IF  
490    
491      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
492      IF (ancien_ok) THEN      IF (ancien_ok) THEN
# Line 878  contains Line 499  contains
499      ELSE      ELSE
500         DO k = 1, llm         DO k = 1, llm
501            DO i = 1, klon            DO i = 1, klon
502               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
503               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
504            ENDDO            ENDDO
505         ENDDO         ENDDO
506         ancien_ok = .TRUE.         ancien_ok = .TRUE.
# Line 895  contains Line 516  contains
516      ! Check temperatures:      ! Check temperatures:
517      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
518    
519      ! Incrementer le compteur de la physique      call increment_itap
520      itap = itap + 1      julien = MOD(dayvrai, 360)
     julien = MOD(NINT(rdayvrai), 360)  
521      if (julien == 0) julien = 360      if (julien == 0) julien = 360
522    
523      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
524    
525      ! Mettre en action les conditions aux limites (albedo, sst, etc.).      ! \'Evaporation de l'eau liquide nuageuse :
   
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.  
     if (nqmx >= 5) then  
        wo = qx(:, :, 5) * zmasse / dobson_u / 1e3  
     else IF (MOD(itap - 1, lmt_pas) == 0) THEN  
        wo = ozonecm(REAL(julien), paprs)  
     ENDIF  
   
     ! Évaporation de l'eau liquide nuageuse :  
526      DO k = 1, llm      DO k = 1, llm
527         DO i = 1, klon         DO i = 1, klon
528            zb = MAX(0., ql_seri(i, k))            zb = MAX(0., ql_seri(i, k))
# Line 922  contains Line 533  contains
533      ENDDO      ENDDO
534      ql_seri = 0.      ql_seri = 0.
535    
536      IF (if_ebil >= 2) THEN      frugs = MAX(frugs, 0.000015)
537         ztit = 'after reevap'      zxrugs = sum(frugs * pctsrf, dim = 2)
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        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, &  
             fs_bound, fq_bound)  
   
     END IF  
   
     ! Appeler la diffusion verticale (programme de couche limite)  
   
     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  
   
     ! calculs necessaires au calcul de l'albedo dans l'interface  
   
     CALL orbite(REAL(julien), zlongi, dist)  
     IF (cycle_diurne) THEN  
        zdtime = dtphys * REAL(radpas)  
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
     ELSE  
        rmu0 = -999.999  
     ENDIF  
   
     ! Calcul de l'abedo moyen par maille  
     albsol(:) = 0.  
     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  
   
     ! Repartition sous maille des flux LW et SW  
     ! Repartition du longwave par sous-surface linearisee  
538    
539      DO nsrf = 1, nbsrf      ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
540         DO i = 1, klon      ! la surface.
           fsollw(i, nsrf) = sollw(i) &  
                + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))  
           fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))  
        ENDDO  
     ENDDO  
541    
542      fder = dlw      CALL orbite(REAL(julien), longi, dist)
543        CALL zenang(longi, time, dtphys * radpas, mu0, fract)
544    
545      ! Couche limite:      CALL pbl_surface(pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
546             ftsol, cdmmax, cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
547      CALL clmain(dtphys, itap, date0, pctsrf, pctsrf_new, t_seri, q_seri, &           falbe, fluxlat, rain_fall, snow_fall, frugs, agesno, rugoro, d_t_vdf, &
548           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, &           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, flux_t, flux_q, flux_u, flux_v, &
549           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           cdragh, cdragm, q2, dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, &
550           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           v10m_srf, pblh, capCL, oliqCL, cteiCL, pblT, therm, plcl, fqcalving, &
551           rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder, rlon, rlat, &           ffonte, run_off_lic_0, albsol, sollw, solsw, tsol)
552           cuphy, cvphy, frugs, firstcal, lafin, agesno, rugoro, d_t_vdf, &  
553           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &      ! Incr\'ementation des flux
554           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &  
555           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &      sens = sum(flux_t * pctsrf, dim = 2)
556           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)      evap = - sum(flux_q * pctsrf, dim = 2)
557        fder = dlw + dflux_t + dflux_q
     ! Incrémentation des flux  
   
     zxfluxt = 0.  
     zxfluxq = 0.  
     zxfluxu = 0.  
     zxfluxv = 0.  
     DO nsrf = 1, nbsrf  
        DO k = 1, llm  
           DO i = 1, klon  
              zxfluxt(i, k) = zxfluxt(i, k) + &  
                   fluxt(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxq(i, k) = zxfluxq(i, k) + &  
                   fluxq(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)  
           END DO  
        END DO  
     END DO  
     DO i = 1, klon  
        sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol  
        evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol  
        fder(i) = dlw(i) + dsens(i) + devap(i)  
     ENDDO  
558    
559      DO k = 1, llm      DO k = 1, llm
560         DO i = 1, klon         DO i = 1, klon
# Line 1030  contains Line 565  contains
565         ENDDO         ENDDO
566      ENDDO      ENDDO
567    
568      IF (if_ebil >= 2) THEN      call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
569         ztit = 'after clmain'      ftsol = ftsol + d_ts ! update surface temperature
570         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &      tsol = sum(ftsol * pctsrf, dim = 2)
571              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &      zxfluxlat = sum(fluxlat * pctsrf, dim = 2)
572              d_ql, d_qs, d_ec)      zt2m = sum(t2m * pctsrf, dim = 2)
573         call diagphy(airephy, ztit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &      zq2m = sum(q2m * pctsrf, dim = 2)
574              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &      u10m = sum(u10m_srf * pctsrf, dim = 2)
575              fs_bound, fq_bound)      v10m = sum(v10m_srf * pctsrf, dim = 2)
576      END IF      zxffonte = sum(ffonte * pctsrf, dim = 2)
577        s_pblh = sum(pblh * pctsrf, dim = 2)
578      ! Update surface temperature:      s_lcl = sum(plcl * pctsrf, dim = 2)
579        s_capCL = sum(capCL * pctsrf, dim = 2)
580      DO i = 1, klon      s_oliqCL = sum(oliqCL * pctsrf, dim = 2)
581         zxtsol(i) = 0.0      s_cteiCL = sum(cteiCL * pctsrf, dim = 2)
582         zxfluxlat(i) = 0.0      s_pblT = sum(pblT * pctsrf, dim = 2)
583        s_therm = sum(therm * pctsrf, dim = 2)
584    
585         zt2m(i) = 0.0      ! Si une sous-fraction n'existe pas, elle prend la valeur moyenne :
        zq2m(i) = 0.0  
        zu10m(i) = 0.0  
        zv10m(i) = 0.0  
        zxffonte(i) = 0.0  
        zxfqcalving(i) = 0.0  
   
        s_pblh(i) = 0.0  
        s_lcl(i) = 0.0  
        s_capCL(i) = 0.0  
        s_oliqCL(i) = 0.0  
        s_cteiCL(i) = 0.0  
        s_pblT(i) = 0.0  
        s_therm(i) = 0.0  
        s_trmb1(i) = 0.0  
        s_trmb2(i) = 0.0  
        s_trmb3(i) = 0.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  
     ENDDO  
586      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
587         DO i = 1, klon         DO i = 1, klon
588            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)            IF (pctsrf(i, nsrf) < epsfra) then
589            zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)               ftsol(i, nsrf) = tsol(i)
590            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)               t2m(i, nsrf) = zt2m(i)
591                 q2m(i, nsrf) = zq2m(i)
592            zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)               u10m_srf(i, nsrf) = u10m(i)
593            zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)               v10m_srf(i, nsrf) = v10m(i)
594            zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)               ffonte(i, nsrf) = zxffonte(i)
595            zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)               pblh(i, nsrf) = s_pblh(i)
596            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)               plcl(i, nsrf) = s_lcl(i)
597            zxfqcalving(i) = zxfqcalving(i) + &               capCL(i, nsrf) = s_capCL(i)
598                 fqcalving(i, nsrf)*pctsrf(i, nsrf)               oliqCL(i, nsrf) = s_oliqCL(i)
599            s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)               cteiCL(i, nsrf) = s_cteiCL(i)
600            s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)               pblT(i, nsrf) = s_pblT(i)
601            s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)               therm(i, nsrf) = s_therm(i)
602            s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)            end IF
           s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)  
           s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)  
           s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)  
603         ENDDO         ENDDO
604      ENDDO      ENDDO
605    
606      ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne      dlw = - 4. * RSIGMA * tsol**3
607    
608      DO nsrf = 1, nbsrf      ! Appeler la convection
609         DO i = 1, klon  
610            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)      if (conv_emanuel) then
611           CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
612            IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)              d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
613            IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)              upwd, dnwd, Ma, cape, iflagctrl, clwcon0, pmflxr, da, phi, mp)
614            IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)         snow_con = 0.
615            IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)         mfu = upwd + dnwd
616            IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)  
617            IF (pctsrf(i, nsrf) < epsfra) &         zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
618                 fqcalving(i, nsrf) = zxfqcalving(i)         zqsat = zqsat / (1. - retv * zqsat)
619            IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)  
620            IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)         ! Properties of convective clouds
621            IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)         clwcon0 = fact_cldcon * clwcon0
622            IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
623            IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)              rnebcon0)
624            IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)  
625            IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)         forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
626            IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)         mfd = 0.
627            IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)         pen_u = 0.
628            IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)         pen_d = 0.
629         ENDDO         pde_d = 0.
630      ENDDO         pde_u = 0.
631        else
632      ! Calculer la derive du flux infrarouge         conv_q = d_q_dyn + d_q_vdf / dtphys
633           conv_t = d_t_dyn + d_t_vdf / dtphys
634      DO i = 1, klon         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
635         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         CALL conflx(paprs, play, t_seri(:, llm:1:- 1), q_seri(:, llm:1:- 1), &
636      ENDDO              conv_t, conv_q, - evap, omega, d_t_con, d_q_con, rain_con, &
637                snow_con, mfu(:, llm:1:- 1), mfd(:, llm:1:- 1), pen_u, pde_u, &
638      ! Appeler la convection (au choix)              pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, pmflxs)
   
     DO k = 1, llm  
        DO i = 1, klon  
           conv_q(i, k) = d_q_dyn(i, k) &  
                + d_q_vdf(i, k)/dtphys  
           conv_t(i, k) = d_t_dyn(i, k) &  
                + d_t_vdf(i, k)/dtphys  
        ENDDO  
     ENDDO  
     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  
   
     select case (iflag_con)  
     case (1)  
        print *, 'Réactiver l''appel à "conlmd" dans "physiq.F".'  
        stop 1  
     case (2)  
        CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &  
             zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &  
             pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &  
             pmflxs)  
639         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
640         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
641         DO i = 1, klon         ibas_con = llm + 1 - kcbot
642            ibas_con(i) = llm + 1 - kcbot(i)         itop_con = llm + 1 - kctop
643            itop_con(i) = llm + 1 - kctop(i)      END if
        ENDDO  
     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  
        ELSE  
           ! conema3 ne contient pas les traceurs  
           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)  
        ENDIF  
   
        IF (.NOT. ok_gust) THEN  
           do i = 1, klon  
              wd(i) = 0.0  
           enddo  
        ENDIF  
   
        ! Calcul des propriétés des nuages convectifs  
   
        DO k = 1, llm  
           DO i = 1, klon  
              zx_t = t_seri(i, k)  
              IF (thermcep) THEN  
                 zdelta = MAX(0., SIGN(1., rtt-zx_t))  
                 zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)  
                 zx_qs = MIN(0.5, zx_qs)  
                 zcor = 1./(1.-retv*zx_qs)  
                 zx_qs = zx_qs*zcor  
              ELSE  
                 IF (zx_t < t_coup) THEN  
                    zx_qs = qsats(zx_t)/play(i, k)  
                 ELSE  
                    zx_qs = qsatl(zx_t)/play(i, k)  
                 ENDIF  
              ENDIF  
              zqsat(i, k) = zx_qs  
           ENDDO  
        ENDDO  
   
        ! 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  
644    
645      DO k = 1, llm      DO k = 1, llm
646         DO i = 1, klon         DO i = 1, klon
# Line 1246  contains Line 651  contains
651         ENDDO         ENDDO
652      ENDDO      ENDDO
653    
654      IF (if_ebil >= 2) THEN      IF (.not. conv_emanuel) THEN
655         ztit = 'after convect'         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
656         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        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, &  
             fs_bound, fq_bound)  
     END IF  
   
     IF (check) THEN  
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *,"aprescon = ", za  
        zx_t = 0.0  
        za = 0.0  
        DO i = 1, klon  
           za = za + airephy(i)/REAL(klon)  
           zx_t = zx_t + (rain_con(i)+ &  
                snow_con(i))*airephy(i)/REAL(klon)  
        ENDDO  
        zx_t = zx_t/za*dtphys  
        print *,"Precip = ", zx_t  
     ENDIF  
     IF (zx_ajustq) THEN  
        DO i = 1, klon  
           z_apres(i) = 0.0  
        ENDDO  
        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  
657         DO k = 1, llm         DO k = 1, llm
658            DO i = 1, klon            DO i = 1, klon
659               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 1291  contains Line 662  contains
662            ENDDO            ENDDO
663         ENDDO         ENDDO
664      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
665    
666      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
667    
668      d_t_ajs = 0.      d_t_ajs = 0.
669      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1303  contains Line 673  contains
673      entr_therm = 0.      entr_therm = 0.
674    
675      if (iflag_thermals == 0) then      if (iflag_thermals == 0) then
        ! Ajustement sec  
676         CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)         CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
677         t_seri = t_seri + d_t_ajs         t_seri = t_seri + d_t_ajs
678         q_seri = q_seri + d_q_ajs         q_seri = q_seri + d_q_ajs
679      else      else
680         ! Thermiques         call calltherm(play, paprs, pphi, u_seri, v_seri, t_seri, q_seri, &
681         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_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)  
682      endif      endif
683    
     IF (if_ebil >= 2) THEN  
        ztit = 'after dry_adjust'  
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
   
684      ! Caclul des ratqs      ! Caclul des ratqs
685    
     ! ratqs convectifs a l'ancienne en fonction de q(z = 0)-q / q  
     ! on ecrase le tableau ratqsc calcule par clouds_gno  
686      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
687           ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
688           ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
689         do k = 1, llm         do k = 1, llm
690            do i = 1, klon            do i = 1, klon
691               if(ptconv(i, k)) then               if(ptconv(i, k)) then
692                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
693                       +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)
694               else               else
695                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
696               endif               endif
# Line 1340  contains Line 701  contains
701      ! ratqs stables      ! ratqs stables
702      do k = 1, llm      do k = 1, llm
703         do i = 1, klon         do i = 1, klon
704            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
705                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
706         enddo         enddo
707      enddo      enddo
708    
709      ! ratqs final      ! ratqs final
710      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
711         ! les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
712         ! ratqs final         ! ratqs final
713         ! 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
714         ! relaxation des ratqs         ! relaxation des ratqs
715         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
716         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
717      else      else
718         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
719         ratqs = ratqss         ratqs = ratqss
720      endif      endif
721    
722      ! Processus de condensation à grande echelle et processus de      CALL fisrtilp(paprs, play, t_seri, q_seri, ptconv, ratqs, d_t_lsc, &
723      ! précipitation :           d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, pfrac_impa, &
724      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &           pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, psfl, rhcl)
          d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &  
          pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &  
          psfl, rhcl)  
725    
726      WHERE (rain_lsc < 0) rain_lsc = 0.      WHERE (rain_lsc < 0) rain_lsc = 0.
727      WHERE (snow_lsc < 0) snow_lsc = 0.      WHERE (snow_lsc < 0) snow_lsc = 0.
# Line 1377  contains Line 734  contains
734            IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)            IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
735         ENDDO         ENDDO
736      ENDDO      ENDDO
     IF (check) THEN  
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *,"apresilp = ", za  
        zx_t = 0.0  
        za = 0.0  
        DO i = 1, klon  
           za = za + airephy(i)/REAL(klon)  
           zx_t = zx_t + (rain_lsc(i) &  
                + snow_lsc(i))*airephy(i)/REAL(klon)  
        ENDDO  
        zx_t = zx_t/za*dtphys  
        print *,"Precip = ", zx_t  
     ENDIF  
   
     IF (if_ebil >= 2) THEN  
        ztit = 'after fisrt'  
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        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, &  
             fs_bound, fq_bound)  
     END IF  
737    
738      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
739    
740      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
741    
742      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke      IF (iflag_cldcon <= - 1) THEN
743           ! seulement pour Tiedtke
744         snow_tiedtke = 0.         snow_tiedtke = 0.
745         if (iflag_cldcon == -1) then         if (iflag_cldcon == - 1) then
746            rain_tiedtke = rain_con            rain_tiedtke = rain_con
747         else         else
748            rain_tiedtke = 0.            rain_tiedtke = 0.
749            do k = 1, llm            do k = 1, llm
750               do i = 1, klon               do i = 1, klon
751                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
752                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k) / dtphys &
753                          *zmasse(i, k)                          * zmasse(i, k)
754                  endif                  endif
755               enddo               enddo
756            enddo            enddo
757         endif         endif
758    
759         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
760         CALL diagcld1(paprs, play, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
761              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
762         DO k = 1, llm         DO k = 1, llm
763            DO i = 1, klon            DO i = 1, klon
764               IF (diafra(i, k) > cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
# Line 1434  contains Line 768  contains
768            ENDDO            ENDDO
769         ENDDO         ENDDO
770      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
771         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
772         ! 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
773         ! facttemps         ! d'un facteur facttemps.
774         facteur = dtphys *facttemps         facteur = dtphys * facttemps
775         do k = 1, llm         do k = 1, llm
776            do i = 1, klon            do i = 1, klon
777               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
778               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
779                    then                    > rnebcon(i, k) * clwcon(i, k)) then
780                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
781                  clwcon(i, k) = clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
782               endif               endif
# Line 1451  contains Line 785  contains
785    
786         ! On prend la somme des fractions nuageuses et des contenus en eau         ! On prend la somme des fractions nuageuses et des contenus en eau
787         cldfra = min(max(cldfra, rnebcon), 1.)         cldfra = min(max(cldfra, rnebcon), 1.)
788         cldliq = cldliq + rnebcon*clwcon         cldliq = cldliq + rnebcon * clwcon
789      ENDIF      ENDIF
790    
791      ! 2. Nuages stratiformes      ! 2. Nuages stratiformes
# Line 1469  contains Line 803  contains
803      ENDIF      ENDIF
804    
805      ! Precipitation totale      ! Precipitation totale
   
806      DO i = 1, klon      DO i = 1, klon
807         rain_fall(i) = rain_con(i) + rain_lsc(i)         rain_fall(i) = rain_con(i) + rain_lsc(i)
808         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
809      ENDDO      ENDDO
810    
811      IF (if_ebil >= 2) THEN      ! Humidit\'e relative pour diagnostic :
        ztit = "after diagcld"  
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
   
     ! Humidité relative pour diagnostic:  
812      DO k = 1, llm      DO k = 1, llm
813         DO i = 1, klon         DO i = 1, klon
814            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
815            IF (thermcep) THEN            zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t) / play(i, k)
816               zdelta = MAX(0., SIGN(1., rtt-zx_t))            zx_qs = MIN(0.5, zx_qs)
817               zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)            zcor = 1. / (1. - retv * zx_qs)
818               zx_qs = MIN(0.5, zx_qs)            zx_qs = zx_qs * zcor
819               zcor = 1./(1.-retv*zx_qs)            zx_rh(i, k) = q_seri(i, k) / zx_qs
              zx_qs = zx_qs*zcor  
           ELSE  
              IF (zx_t < t_coup) THEN  
                 zx_qs = qsats(zx_t)/play(i, k)  
              ELSE  
                 zx_qs = qsatl(zx_t)/play(i, k)  
              ENDIF  
           ENDIF  
           zx_rh(i, k) = q_seri(i, k)/zx_qs  
820            zqsat(i, k) = zx_qs            zqsat(i, k) = zx_qs
821         ENDDO         ENDDO
822      ENDDO      ENDDO
823    
824      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Param\`etres optiques des nuages et quelques param\`etres pour
     ! Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)  
     IF (ok_ade .OR. ok_aie) THEN  
        ! Get sulfate aerosol distribution  
        CALL readsulfate(rdayvrai, firstcal, sulfate)  
        CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)  
   
        ! Calculate aerosol optical properties (Olivier Boucher)  
        CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &  
             aerindex)  
     ELSE  
        tau_ae = 0.  
        piz_ae = 0.  
        cg_ae = 0.  
     ENDIF  
   
     ! Paramètres optiques des nuages et quelques paramètres pour  
825      ! diagnostics :      ! diagnostics :
826      if (ok_newmicro) then      if (ok_newmicro) then
827         CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
828              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc)
             fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, &  
             re, fl)  
829      else      else
830         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
831              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq)
             bl95_b1, cldtaupi, re, fl)  
832      endif      endif
833    
834      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
835      IF (MOD(itaprad, radpas) == 0) THEN         wo = ozonecm(REAL(julien), paprs)
836         DO i = 1, klon         albsol = sum(falbe * pctsrf, dim = 2)
837            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &         CALL radlwsw(dist, mu0, fract, paprs, play, tsol, albsol, t_seri, &
838                 + falbe(i, is_lic) * pctsrf(i, is_lic) &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
839                 + falbe(i, is_ter) * pctsrf(i, is_ter) &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
840                 + falbe(i, is_sic) * pctsrf(i, is_sic)              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
841            albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &              swup0, swup, ok_ade, topswad, solswad)
                + falblw(i, is_lic) * pctsrf(i, is_lic) &  
                + falblw(i, is_ter) * pctsrf(i, is_ter) &  
                + falblw(i, is_sic) * pctsrf(i, is_sic)  
        ENDDO  
        ! 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  
842      ENDIF      ENDIF
     itaprad = itaprad + 1  
843    
844      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
   
845      DO k = 1, llm      DO k = 1, llm
846         DO i = 1, klon         DO i = 1, klon
847            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 &
848                   / 86400.
849         ENDDO         ENDDO
850      ENDDO      ENDDO
851    
852      IF (if_ebil >= 2) THEN      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
        ztit = 'after rad'  
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        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, &  
             fs_bound, fq_bound)  
     END IF  
   
     ! Calculer l'hydrologie de la surface  
853      DO i = 1, klon      DO i = 1, klon
854         zxqsurf(i) = 0.0         bils(i) = radsol(i) + sens(i) + zxfluxlat(i)
        zxsnow(i) = 0.0  
855      ENDDO      ENDDO
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)  
           zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
   
     ! Calculer le bilan du sol et la dérive de température (couplage)  
856    
857      DO i = 1, klon      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
        bils(i) = radsol(i) - sens(i) + zxfluxlat(i)  
     ENDDO  
   
     ! Paramétrisation de l'orographie à l'échelle sous-maille :  
858    
859      IF (ok_orodr) THEN      IF (ok_orodr) THEN
860         ! selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
        igwd = 0  
861         DO i = 1, klon         DO i = 1, klon
862            itest(i) = 0            ktest(i) = 0
863            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.0)) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
864               itest(i) = 1               ktest(i) = 1
              igwd = igwd + 1  
              idx(igwd) = i  
865            ENDIF            ENDIF
866         ENDDO         ENDDO
867    
868         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(paprs, play, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
869              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              ktest, t_seri, u_seri, v_seri, zulow, zvlow, zustrdr, zvstrdr, &
870              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              d_t_oro, d_u_oro, d_v_oro)
871    
872         ! ajout des tendances         ! ajout des tendances
873         DO k = 1, llm         DO k = 1, llm
# Line 1621  contains Line 880  contains
880      ENDIF      ENDIF
881    
882      IF (ok_orolf) THEN      IF (ok_orolf) THEN
883         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
        igwd = 0  
884         DO i = 1, klon         DO i = 1, klon
885            itest(i) = 0            ktest(i) = 0
886            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
887               itest(i) = 1               ktest(i) = 1
              igwd = igwd + 1  
              idx(igwd) = i  
888            ENDIF            ENDIF
889         ENDDO         ENDDO
890    
891         CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &         CALL lift_noro(paprs, play, zmea, zstd, zpic, ktest, t_seri, u_seri, &
892              itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &              v_seri, zulow, zvlow, zustrli, zvstrli, d_t_lif, d_u_lif, d_v_lif)
             d_t_lif, d_u_lif, d_v_lif)  
893    
894         ! Ajout des tendances :         ! Ajout des tendances :
895         DO k = 1, llm         DO k = 1, llm
# Line 1646  contains Line 901  contains
901         ENDDO         ENDDO
902      ENDIF      ENDIF
903    
904      ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE      CALL aaam_bud(rg, romega, pphis, zustrdr, zustrli, &
905             sum((u_seri - u) / dtphys * zmasse, dim = 2), zvstrdr, &
906      DO i = 1, klon           zvstrli, sum((v_seri - v) / dtphys * zmasse, dim = 2), paprs, u, v, &
907         zustrph(i) = 0.           aam, torsfc)
        zvstrph(i) = 0.  
     ENDDO  
     DO k = 1, llm  
        DO i = 1, klon  
           zustrph(i) = zustrph(i) + (u_seri(i, k)-u(i, k))/dtphys* zmasse(i, k)  
           zvstrph(i) = zvstrph(i) + (v_seri(i, k)-v(i, k))/dtphys* zmasse(i, k)  
        ENDDO  
     ENDDO  
   
     CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &  
          zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)  
   
     IF (if_ebil >= 2) THEN  
        ztit = 'after orography'  
        CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
     END IF  
908    
909      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
910      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, &      call phytrac(julien, time, firstcal, lafin, t, paprs, play, mfu, mfd, &
911           nqmx-2, dtphys, u, t, paprs, play, pmfu, pmfd, pen_u, pde_u, &           pde_u, pen_d, coefh, cdragh, fm_therm, entr_therm, u(:, 1), v(:, 1), &
912           pen_d, pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &           ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
913           frac_impa, frac_nucl, pphis, albsol, rhcl, cldfra, rneb, &           tr_seri, zmasse, ncid_startphy)
          diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, &  
          tr_seri, zmasse)  
   
     IF (offline) THEN  
        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  
914    
915      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
916      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)  
917    
918      ! diag. bilKP      ! diag. bilKP
919    
920      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, &
921           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
922    
923      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
924    
925      ! conversion Ec -> E thermique      ! conversion Ec en énergie thermique
926      DO k = 1, llm      DO k = 1, llm
927         DO i = 1, klon         DO i = 1, klon
928            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))            d_t_ec(i, k) = 0.5 / (RCPD * (1. + RVTMP2 * q_seri(i, k))) &
           d_t_ec(i, k) = 0.5 / ZRCPD &  
929                 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)                 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
930            t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)            t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
931            d_t_ec(i, k) = d_t_ec(i, k) / dtphys            d_t_ec(i, k) = d_t_ec(i, k) / dtphys
932         END DO         END DO
933      END DO      END DO
934    
     IF (if_ebil >= 1) THEN  
        ztit = 'after physic'  
        CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &  
             ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &  
             d_ql, d_qs, d_ec)  
        ! Comme les tendances de la physique sont ajoute dans la dynamique,  
        ! on devrait avoir que la variation d'entalpie par la dynamique  
        ! est egale a la variation de la physique au pas de temps precedent.  
        ! Donc la somme de ces 2 variations devrait etre nulle.  
        call diagphy(airephy, ztit, ip_ebil, topsw, toplw, solsw, sollw, sens, &  
             evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &  
             fs_bound, fq_bound)  
   
        d_h_vcol_phy = d_h_vcol  
   
     END IF  
   
935      ! SORTIES      ! SORTIES
936    
937      !cc prw = eau precipitable      ! prw = eau precipitable
938      DO i = 1, klon      DO i = 1, klon
939         prw(i) = 0.         prw(i) = 0.
940         DO k = 1, llm         DO k = 1, llm
941            prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)            prw(i) = prw(i) + q_seri(i, k) * zmasse(i, k)
942         ENDDO         ENDDO
943      ENDDO      ENDDO
944    
# Line 1744  contains Line 954  contains
954         ENDDO         ENDDO
955      ENDDO      ENDDO
956    
957      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
958         DO iq = 3, nqmx         DO k = 1, llm
959            DO k = 1, llm            DO i = 1, klon
960               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  
961            ENDDO            ENDDO
962         ENDDO         ENDDO
963      ENDIF      ENDDO
964    
965      ! 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:
966      DO k = 1, llm      DO k = 1, llm
# Line 1762  contains Line 970  contains
970         ENDDO         ENDDO
971      ENDDO      ENDDO
972    
973      ! Ecriture des sorties      CALL histwrite_phy("phis", pphis)
974      call write_histhf      CALL histwrite_phy("aire", airephy)
975      call write_histday      CALL histwrite_phy("psol", paprs(:, 1))
976      call write_histins      CALL histwrite_phy("precip", rain_fall + snow_fall)
977        CALL histwrite_phy("plul", rain_lsc + snow_lsc)
978      ! Si c'est la fin, il faut conserver l'etat de redemarrage      CALL histwrite_phy("pluc", rain_con + snow_con)
979      IF (lafin) THEN      CALL histwrite_phy("tsol", tsol)
980         itau_phy = itau_phy + itap      CALL histwrite_phy("t2m", zt2m)
981         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &      CALL histwrite_phy("q2m", zq2m)
982              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &      CALL histwrite_phy("u10m", u10m)
983              rain_fall, snow_fall, solsw, sollwdown, dlw, radsol, frugs, &      CALL histwrite_phy("v10m", v10m)
984              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &      CALL histwrite_phy("snow", snow_fall)
985              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)      CALL histwrite_phy("cdrm", cdragm)
986      ENDIF      CALL histwrite_phy("cdrh", cdragh)
987        CALL histwrite_phy("topl", toplw)
988      firstcal = .FALSE.      CALL histwrite_phy("evap", evap)
989        CALL histwrite_phy("sols", solsw)
990    contains      CALL histwrite_phy("rls", sollw)
991        CALL histwrite_phy("solldown", sollwdown)
992      subroutine write_histday      CALL histwrite_phy("bils", bils)
993        CALL histwrite_phy("sens", sens)
994        use gr_phy_write_3d_m, only: gr_phy_write_3d      CALL histwrite_phy("fder", fder)
995        integer itau_w ! pas de temps ecriture      CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
996        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
997        !------------------------------------------------      CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
998        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
999        if (ok_journe) THEN      CALL histwrite_phy("zxfqcalving", sum(fqcalving * pctsrf, dim = 2))
1000           itau_w = itau_phy + itap      CALL histwrite_phy("albs", albsol)
1001           if (nqmx <= 4) then      CALL histwrite_phy("tro3", wo * dobson_u * 1e3 / zmasse / rmo3 * md)
1002              call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &      CALL histwrite_phy("rugs", zxrugs)
1003                   gr_phy_write_3d(wo) * 1e3)      CALL histwrite_phy("s_pblh", s_pblh)
1004              ! (convert "wo" from kDU to DU)      CALL histwrite_phy("s_pblt", s_pblt)
1005           end if      CALL histwrite_phy("s_lcl", s_lcl)
1006           if (ok_sync) then      CALL histwrite_phy("s_capCL", s_capCL)
1007              call histsync(nid_day)      CALL histwrite_phy("s_oliqCL", s_oliqCL)
1008           endif      CALL histwrite_phy("s_cteiCL", s_cteiCL)
1009        ENDIF      CALL histwrite_phy("s_therm", s_therm)
1010        CALL histwrite_phy("temp", t_seri)
1011      End subroutine write_histday      CALL histwrite_phy("vitu", u_seri)
1012        CALL histwrite_phy("vitv", v_seri)
1013      !****************************      CALL histwrite_phy("geop", zphi)
1014        CALL histwrite_phy("pres", play)
1015      subroutine write_histhf      CALL histwrite_phy("dtvdf", d_t_vdf)
1016        CALL histwrite_phy("dqvdf", d_q_vdf)
1017        ! From phylmd/write_histhf.h, version 1.5 2005/05/25 13:10:09      CALL histwrite_phy("rhum", zx_rh)
1018        CALL histwrite_phy("d_t_ec", d_t_ec)
1019        CALL histwrite_phy("dtsw0", heat0 / 86400.)
1020        CALL histwrite_phy("dtlw0", - cool0 / 86400.)
1021        CALL histwrite_phy("msnow", sum(fsnow * pctsrf, dim = 2))
1022        call histwrite_phy("qsurf", sum(fqsurf * pctsrf, dim = 2))
1023        call histwrite_phy("flat", zxfluxlat)
1024    
1025        !------------------------------------------------      DO nsrf = 1, nbsrf
1026           CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1027        call write_histhf3d         CALL histwrite_phy("sens_"//clnsurf(nsrf), flux_t(:, nsrf))
1028           CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1029        IF (ok_sync) THEN         CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1030           call histsync(nid_hf)         CALL histwrite_phy("taux_"//clnsurf(nsrf), flux_u(:, nsrf))
1031        ENDIF         CALL histwrite_phy("tauy_"//clnsurf(nsrf), flux_v(:, nsrf))
1032           CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1033      end subroutine write_histhf         CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1034           CALL histwrite_phy("u10m_"//clnsurf(nsrf), u10m_srf(:, nsrf))
1035      !***************************************************************         CALL histwrite_phy("v10m_"//clnsurf(nsrf), v10m_srf(:, nsrf))
1036        END DO
     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)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)  
          !ccIM  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)  
          CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)  
          CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)  
          CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)  
          CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)  
   
          zx_tmp_fi2d(1:klon) = -1*sens(1:klon)  
          ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)  
          CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)  
   
          DO nsrf = 1, nbsrf  
             !XXX  
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
          END DO  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)  
          CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)  
   
          !HBTM2  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)  
          CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)  
   
          ! Champs 3D:  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)  
          CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)  
          CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)  
          CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)  
   
          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  
1037    
1038        if (ok_sync) then      if (conv_emanuel) then
1039           call histsync(nid_hf3d)         CALL histwrite_phy("ptop", ema_pct)
1040        endif         CALL histwrite_phy("dnwd0", - mp)
1041        end if
1042    
1043        if (ok_instan) call histsync(nid_ins)
1044    
1045        IF (lafin) then
1046           call NF95_CLOSE(ncid_startphy)
1047           CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
1048                rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, &
1049                zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
1050                rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
1051        end IF
1052    
1053      end subroutine write_histhf3d      firstcal = .FALSE.
1054    
1055    END SUBROUTINE physiq    END SUBROUTINE physiq
1056    

Legend:
Removed from v.56  
changed lines
  Added in v.322

  ViewVC Help
Powered by ViewVC 1.1.21