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

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

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

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

Legend:
Removed from v.35  
changed lines
  Added in v.195

  ViewVC Help
Powered by ViewVC 1.1.21