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

Diff of /trunk/phylmd/physiq.f

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

trunk/libf/phylmd/physiq.f90 revision 18 by guez, Thu Aug 7 12:29:13 2008 UTC trunk/Sources/phylmd/physiq.f revision 202 by guez, Wed Jun 8 12:23:41 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(nq, firstcal, lafin, rdayvrai, gmtime, pdtphys, paprs, &    SUBROUTINE physiq(lafin, dayvrai, time, paprs, play, pphi, pphis, u, v, t, &
8         pplay, pphi, pphis, presnivs, 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)  
   
     ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28  
   
     ! Author : Z.X. Li (LMD/CNRS), date: 1993/08/18  
   
     ! 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 ioipsl, only: ymds2ju, histwrite, histsync  
     use dimens_m, only: jjm, iim, llm  
     use indicesol, only: nbsrf, is_ter, is_lic, is_sic, is_oce, &  
          clnsurf, epsfra  
     use dimphy, only: klon, nbtr  
     use conf_gcm_m, only: raz_date, offline, iphysiq  
     use dimsoil, only: nsoilmx  
     use temps, only: itau_phy, day_ref, annee_ref, itaufin  
     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 iniprint, only: prt_level  
     use abort_gcm_m, only: abort_gcm  
     use YOMCST, only: rcpd, rtt, rlvtt, rg, ra, rsigma, retv, romega  
     use comgeomphy  
     use ctherm  
     use phytrac_m, only: phytrac  
     use oasis_m  
     use radepsi  
     use radopt  
     use yoethf  
     use ini_hist, only: ini_histhf, ini_histday, ini_histins  
     use orbite_m, only: orbite, zenang  
     use phyetat0_m, only: phyetat0, rlat, rlon  
     use hgardfou_m, only: hgardfou  
     use conf_phys_m, only: conf_phys  
     use phyredem_m, only: phyredem  
     use qcheck_m, only: qcheck  
   
     ! Declaration des constantes et des fonctions thermodynamiques :  
     use fcttre, only: thermcep, foeew, qsats, qsatl  
   
     ! Variables argument:  
   
     INTEGER, intent(in):: nq ! nombre de traceurs (y compris vapeur d'eau)  
   
     REAL, intent(in):: rdayvrai  
     ! (elapsed time since January 1st 0h of the starting year, in days)  
   
     REAL, intent(in):: gmtime ! heure de la journée en fraction de jour  
     REAL, intent(in):: pdtphys ! pas d'integration pour la physique (seconde)  
     LOGICAL, intent(in):: firstcal ! first call to "calfis"  
     logical, intent(in):: lafin ! dernier passage  
9    
10      REAL, intent(in):: paprs(klon, llm+1)      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11      ! (pression pour chaque inter-couche, en Pa)      ! (subversion revision 678)
12    
13      REAL, intent(in):: pplay(klon, llm)      ! Author: Z. X. Li (LMD/CNRS) 1993
     ! (input pression pour le mileu de chaque couche (en Pa))  
14    
15      REAL pphi(klon, llm)        ! This is the main procedure for the "physics" part of the program.
16      ! (input geopotentiel de chaque couche (g z) (reference sol))  
17        use aaam_bud_m, only: aaam_bud
18        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_ins, ksta, ksta_ter, ok_kzmin, &
22             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, lmt_pas
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    
73      REAL pphis(klon) ! input geopotentiel du sol      logical, intent(in):: lafin ! dernier passage
   
     REAL presnivs(llm)  
     ! (input pressions approximat. des milieux couches ( en PA))  
74    
75      REAL u(klon, llm)  ! input vitesse dans la direction X (de O a E) en m/s      integer, intent(in):: dayvrai
76      REAL v(klon, llm)  ! input vitesse Y (de S a N) en m/s      ! current day number, based at value 1 on January 1st of annee_ref
     REAL t(klon, llm)  ! input temperature (K)  
77    
78      REAL, intent(in):: qx(klon, llm, nq)      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
     ! (humidite specifique (kg/kg) et fractions massiques des autres traceurs)  
79    
80      REAL omega(klon, llm)  ! input vitesse verticale en Pa/s      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
81      REAL d_u(klon, llm)  ! output tendance physique de "u" (m/s/s)      ! pression pour chaque inter-couche, en Pa
     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, nq)  ! output tendance physique de "qx" (kg/kg/s)  
     REAL d_ps(klon)  ! output tendance physique de la pression au sol  
82    
83      INTEGER nbteta      REAL, intent(in):: play(:, :) ! (klon, llm)
84      PARAMETER(nbteta=3)      ! pression pour le mileu de chaque couche (en Pa)
85    
86      REAL PVteta(klon, nbteta)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
87      ! (output vorticite potentielle a des thetas constantes)      ! géopotentiel de chaque couche (référence sol)
88    
89      LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
     PARAMETER (ok_cvl=.TRUE.)  
     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface  
     PARAMETER (ok_gust=.FALSE.)  
90    
91      LOGICAL check ! Verifier la conservation du modele en eau      REAL, intent(in):: u(:, :) ! (klon, llm)
92      PARAMETER (check=.FALSE.)      ! vitesse dans la direction X (de O a E) en m / s
     LOGICAL ok_stratus ! Ajouter artificiellement les stratus  
     PARAMETER (ok_stratus=.FALSE.)  
93    
94      ! Parametres lies au coupleur OASIS:      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m / s
95      INTEGER, SAVE :: npas, nexca      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
     logical rnpb  
     parameter(rnpb=.true.)  
96    
97      character(len=6), save:: ocean      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
98      ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
99    
100      logical ok_ocean      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa / s
101      SAVE ok_ocean      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      !IM "slab" ocean      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
106      REAL tslab(klon)    !Temperature du slab-ocean      ! tendance physique de "qx" (s-1)
     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  
107    
108      ! Modele thermique du sol, a activer pour le cycle diurne:      ! Local:
     logical, save:: ok_veget  
     LOGICAL, save:: ok_journe ! sortir le fichier journalier  
109    
110      LOGICAL ok_mensuel ! sortir le fichier mensuel      LOGICAL:: firstcal = .true.
111    
112      LOGICAL ok_instan ! sortir le fichier instantane      LOGICAL, PARAMETER:: check = .FALSE.
113      save ok_instan      ! Verifier la conservation du modele en eau
114    
115      LOGICAL ok_region ! sortir le fichier regional      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
116      PARAMETER (ok_region=.FALSE.)      ! Ajouter artificiellement les stratus
117    
118      !     pour phsystoke avec thermiques      ! pour phystoke avec thermiques
119      REAL fm_therm(klon, llm+1)      REAL fm_therm(klon, llm + 1)
120      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
121      real q2(klon, llm+1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
     save q2  
122    
123      INTEGER ivap          ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
124      PARAMETER (ivap=1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     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  
125    
126      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
127      REAL d_q_dyn(klon, llm)  ! tendance dynamique pour "q" (kg/kg/s)      LOGICAL, save:: ancien_ok
128    
129      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K / s)
130        REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg / kg / s)
     !IM Amip2 PV a theta constante  
   
     CHARACTER(LEN=3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
131    
132      !MI Amip2 PV a theta constante      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
   
     INTEGER klevp1  
     PARAMETER(klevp1=llm+1)  
133    
134      REAL swdn0(klon, klevp1), swdn(klon, klevp1)      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
135      REAL swup0(klon, klevp1), swup(klon, klevp1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
136      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
137    
138      REAL SWdn200clr(klon), SWdn200(klon)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
139      REAL SWup200clr(klon), SWup200(klon)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
140      SAVE SWdn200clr, SWdn200, SWup200clr, SWup200      SAVE lwdn0, lwdn, lwup0, lwup
   
     REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)  
     REAL lwup0(klon, klevp1), lwup(klon, klevp1)  
     SAVE lwdn0, lwdn, lwup0, lwup  
   
     REAL LWdn200clr(klon), LWdn200(klon)  
     REAL LWup200clr(klon), LWup200(klon)  
     SAVE LWdn200clr, LWdn200, LWup200clr, LWup200  
   
     !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  '/  
   
     real tlevSTD(klon, nlevSTD), qlevSTD(klon, nlevSTD)  
     real rhlevSTD(klon, nlevSTD), philevSTD(klon, nlevSTD)  
     real ulevSTD(klon, nlevSTD), vlevSTD(klon, nlevSTD)  
     real wlevSTD(klon, nlevSTD)  
   
     ! nout : niveau de output des variables a une pression donnee  
     INTEGER nout  
     PARAMETER(nout=3) !nout=1 : day; =2 : mth; =3 : NMC  
   
     logical oknondef(klon, nlevSTD, nout)  
     real tnondef(klon, nlevSTD, nout)  
     save tnondef  
   
     ! les produits uvSTD, vqSTD, .., T2STD sont calcules  
     ! a partir des valeurs instantannees toutes les 6 h  
     ! qui sont moyennees sur le mois  
   
     real uvSTD(klon, nlevSTD)  
     real vqSTD(klon, nlevSTD)  
     real vTSTD(klon, nlevSTD)  
     real wqSTD(klon, nlevSTD)  
   
     real vphiSTD(klon, nlevSTD)  
     real wTSTD(klon, nlevSTD)  
     real u2STD(klon, nlevSTD)  
     real v2STD(klon, nlevSTD)  
     real T2STD(klon, nlevSTD)  
141    
142      ! prw: precipitable water      ! prw: precipitable water
143      real prw(klon)      real prw(klon)
144    
145      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg / m2)
146      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg / kg)
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 l, 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  
   
     INTEGER        longcles  
     PARAMETER    ( longcles = 20 )  
   
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 344  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)
     REAL agesno(klon, nbsrf)  
     SAVE agesno                 ! age de la neige  
   
     REAL run_off_lic_0(klon)  
     SAVE run_off_lic_0  
     !KE43  
     ! Variables liees a la convection de K. Emanuel (sb):  
191    
192      REAL bas, top             ! cloud base and top levels      ! Variables li\'ees \`a la convection d'Emanuel :
193      SAVE bas      REAL, save:: Ma(klon, llm) ! undilute upward mass flux
194      SAVE top      REAL, save:: qcondc(klon, llm) ! in-cld water content from convect
195        REAL, save:: sig1(klon, llm), w01(klon, llm)
     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 397  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        REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
248      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      REAL, save:: albsol(klon) ! albedo du sol total visible
     REAL pctsrf(klon, nbsrf)  
     !IM  
     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE  
   
     SAVE pctsrf                 ! sous-fraction du sol  
     REAL albsol(klon)  
     SAVE albsol                 ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw                 ! albedo du sol total  
   
249      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
250    
251      ! Declaration des procedures appelees      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
252        real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
253    
254      EXTERNAL alboc     ! calculer l'albedo sur ocean      REAL rhcl(klon, llm) ! humiditi relative ciel clair
255      EXTERNAL ajsec     ! ajustement sec      REAL dialiq(klon, llm) ! eau liquide nuageuse
256      EXTERNAL clmain    ! couche limite      REAL diafra(klon, llm) ! fraction nuageuse
257      !KE43      REAL cldliq(klon, llm) ! eau liquide nuageuse
258      EXTERNAL conema3  ! convect4.3      REAL cldfra(klon, llm) ! fraction nuageuse
259      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)      REAL cldtau(klon, llm) ! epaisseur optique
260      EXTERNAL nuage     ! calculer les proprietes radiatives      REAL cldemi(klon, llm) ! emissivite infrarouge
261      EXTERNAL ozonecm   ! prescrire l'ozone  
262      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge      REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
263      EXTERNAL transp    ! transport total de l'eau et de l'energie      REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
264        REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
265      ! Variables locales      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  
266    
267      REAL zxfluxt(klon, llm)      REAL zxfluxt(klon, llm)
268      REAL zxfluxq(klon, llm)      REAL zxfluxq(klon, llm)
269      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
270      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
271    
272      REAL heat(klon, llm)    ! chauffage solaire      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
273      REAL heat0(klon, llm)   ! chauffage solaire ciel clair      ! les variables soient r\'emanentes.
274      REAL cool(klon, llm)    ! refroidissement infrarouge      REAL, save:: heat(klon, llm) ! chauffage solaire
275      REAL cool0(klon, llm)   ! refroidissement infrarouge ciel clair      REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
276      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
277      real sollwdown(klon)    ! downward LW flux at surface      REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
278      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
279      REAL albpla(klon)      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
280      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface      real, save:: sollwdown(klon) ! downward LW flux at surface
281      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
282      ! Le rayonnement n'est pas calcule tous les pas, il faut donc      REAL, save:: albpla(klon)
283      !                      sauvegarder les sorties du rayonnement      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
284      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  
285    
286      INTEGER itaprad      REAL conv_q(klon, llm) ! convergence de l'humidite (kg / kg / s)
287      SAVE itaprad      REAL conv_t(klon, llm) ! convergence of temperature (K / s)
288    
289      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
290      REAL conv_t(klon, llm) ! convergence de la temperature(K/s)      REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
   
     REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut  
     REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree  
291    
292      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
293    
294      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
295      REAL zdtime ! pas de temps du rayonnement (s)      real longi
     real zlongi  
   
296      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
297      REAL za, zb      REAL za, zb
298      REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp      REAL zx_t, zx_qs, zcor
299      real zqsat(klon, llm)      real zqsat(klon, llm)
300      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
301      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup=234.0)  
   
302      REAL zphi(klon, llm)      REAL zphi(klon, llm)
303    
304      !IM cf. AM Variables locales pour la CLA (hbtm2)      ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
305    
306      REAL pblh(klon, nbsrf)           ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
307      REAL plcl(klon, nbsrf)           ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
308      REAL capCL(klon, nbsrf)          ! CAPE de couche limite      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
309      REAL oliqCL(klon, nbsrf)          ! eau_liqu integree de couche limite      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
310      REAL cteiCL(klon, nbsrf)          ! cloud top instab. crit. couche limite      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
311      REAL pblt(klon, nbsrf)          ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
312      REAL therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
313      REAL trmb1(klon, nbsrf)          ! deep_cape      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
314      REAL trmb2(klon, nbsrf)          ! inhibition      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
315      REAL trmb3(klon, nbsrf)          ! Point Omega      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
316      ! Grdeurs de sorties      ! Grandeurs de sorties
317      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
318      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
319      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
320      REAL s_trmb3(klon)      REAL s_trmb3(klon)
321    
322      ! Variables locales pour la convection de K. Emanuel (sb):      ! Variables pour la convection de K. Emanuel :
323    
324      REAL upwd(klon, llm)      ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
325      REAL dnwd(klon, llm)      ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
326      REAL dnwd0(klon, llm)     ! unsaturated downdraft mass flux      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
327      REAL tvp(klon, llm)       ! virtual temp of lifted parcel      REAL cape(klon) ! CAPE
     REAL cape(klon)           ! CAPE  
328      SAVE cape      SAVE cape
329    
330      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)  
331    
332      ! Variables du changement      ! Variables du changement
333    
334      ! con: convection      ! con: convection
335      ! lsc: condensation a grande echelle (Large-Scale-Condensation)      ! lsc: large scale condensation
336      ! ajs: ajustement sec      ! ajs: ajustement sec
337      ! eva: evaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
338      ! vdf: couche limite (Vertical DiFfusion)      ! vdf: vertical diffusion in boundary layer
339      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
340      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
341      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 569  contains Line 343  contains
343      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
344      REAL rneb(klon, llm)      REAL rneb(klon, llm)
345    
346      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
347      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
348      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
349      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
350      REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1)      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
351      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)  
352    
353      SAVE ibas_con, itop_con      INTEGER, save:: ibas_con(klon), itop_con(klon)
354        real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
355    
356      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
357      REAL snow_con(klon), snow_lsc(klon)      REAL, save:: snow_con(klon) ! neige (mm / s)
358        real snow_lsc(klon)
359      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
360    
361      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
# Line 592  contains Line 366  contains
366      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
367      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
368    
369      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
370      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
371      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
372    
373      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
374      real, save:: fact_cldcon      real:: fact_cldcon = 0.375
375      real, save:: facttemps      real:: facttemps = 1.e-4
376      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
377      real facteur      real facteur
378    
379      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
380      logical ptconv(klon, llm)      logical ptconv(klon, llm)
381    
382      ! Variables locales pour effectuer les appels en serie      ! Variables pour effectuer les appels en s\'erie :
383    
384      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
385      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
386      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
387        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
388    
389      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
390    
# Line 624  contains Line 393  contains
393      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
394      REAL aam, torsfc      REAL aam, torsfc
395    
     REAL dudyn(iim+1, jjm + 1, llm)  
   
     REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique  
     REAL zx_tmp_fi3d(klon, llm) ! variable temporaire pour champs 3D  
   
     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)  
   
     INTEGER, SAVE:: nid_day, nid_ins  
   
396      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.
397      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.
398      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.
399      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.
400    
     REAL zsto  
   
     character(len=20) modname  
     character(len=80) abort_message  
     logical ok_sync  
401      real date0      real date0
402    
403      !     Variables liees au bilan d'energie et d'enthalpi      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
404      REAL ztsol(klon)      REAL ztsol(klon)
405      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
406      REAL      d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
407      REAL      fs_bound, fq_bound      REAL zero_v(klon)
408      SAVE      d_h_vcol_phy      CHARACTER(LEN = 20) tit
409      REAL      zero_v(klon)      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
410      CHARACTER(LEN=15) ztit      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
411      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.  
412      SAVE      ip_ebil      REAL d_t_ec(klon, llm)
413      DATA      ip_ebil/0/      ! tendance due \`a la conversion Ec en énergie thermique
414      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)
   
     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.  
     ! ok_ade=T -ADE=topswad-topsw  
438    
439      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.      REAL topswad(klon), solswad(klon) ! aerosol direct effect
440      ! ok_aie=T ->      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
     !        ok_ade=T -AIE=topswai-topswad  
     !        ok_ade=F -AIE=topswai-topsw  
441    
442      REAL aerindex(klon)       ! POLDER aerosol index      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
443        LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
444    
445      ! Parameters      REAL:: bl95_b0 = 2., bl95_b1 = 0.2
446      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not      ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
447      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)      ! B). They link cloud droplet number concentration to aerosol mass
448        ! concentration.
449    
     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 (nq  <  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)
533              dlw, radsol, frugs, agesno, &  
534              zmea, zstd, zsig, zgam, zthe, zpic, zval, &         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
535              t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon,  &         q2 = 1e-8
536              run_off_lic_0)  
537           radpas = lmt_pas / nbapp_rad
538         !   ATTENTION : il faudra a terme relire q2 dans l'etat initial         print *, "radpas = ", radpas
539         q2(:, :, :)=1.e-8  
540           ! Initialisation pour le sch\'ema de convection d'Emanuel :
541         radpas = NINT( 86400. / pdtphys / nbapp_rad)         IF (conv_emanuel) THEN
542              ibas_con = 1
543         ! on remet le calendrier a zero            itop_con = 1
        IF (raz_date) itau_phy = 0  
   
        PRINT *, 'cycle_diurne = ', cycle_diurne  
   
        IF(ocean.NE.'force ') THEN  
           ok_ocean=.TRUE.  
        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  
   
           print *,"*** Convection de Kerry Emanuel 4.3  "  
   
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG  
           DO i = 1, klon  
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
           !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END  
   
544         ENDIF         ENDIF
545    
546         IF (ok_orodr) THEN         IF (ok_orodr) THEN
547            rugoro = MAX(1e-5, zstd * zsig / 2)            rugoro = MAX(1e-5, zstd * zsig / 2)
548            CALL SUGWD(klon, llm, paprs, pplay)            CALL SUGWD(paprs, play)
549         else         else
550            rugoro = 0.            rugoro = 0.
551         ENDIF         ENDIF
552    
553         lmt_pas = NINT(86400. / pdtphys)  ! tous les jours         ecrit_ins = NINT(ecrit_ins / dtphys)
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
        ecrit_ins = NINT(ecrit_ins/pdtphys)  
        ecrit_hf = NINT(ecrit_hf/pdtphys)  
        ecrit_mth = NINT(ecrit_mth/pdtphys)  
        ecrit_tra = NINT(86400.*ecrit_tra/pdtphys)  
        ecrit_reg = NINT(ecrit_reg/pdtphys)  
   
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
        print *,'AVANT HIST IFLAG_CON=', iflag_con  
   
        !   Initialisation des sorties  
   
        call ini_histhf(pdtphys, presnivs, nid_hf, nid_hf3d)  
        call ini_histday(pdtphys, presnivs, ok_journe, nid_day, nq)  
        call ini_histins(pdtphys, presnivs, 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  
     ENDIF test_firstcal  
554    
555      ! Mettre a zero des variables de sortie (pour securite)         ! Initialisation des sorties
556    
557      DO i = 1, klon         call ini_histins(dtphys)
558         d_ps(i) = 0.0         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
559      ENDDO         ! Positionner date0 pour initialisation de ORCHIDEE
560      DO k = 1, llm         print *, 'physiq date0: ', date0
561         DO i = 1, klon         CALL phyredem0
562            d_t(i, k) = 0.0      ENDIF test_firstcal
           d_u(i, k) = 0.0  
           d_v(i, k) = 0.0  
        ENDDO  
     ENDDO  
     DO iq = 1, nq  
        DO k = 1, llm  
           DO i = 1, klon  
              d_qx(i, k, iq) = 0.0  
           ENDDO  
        ENDDO  
     ENDDO  
     da=0.  
     mp=0.  
     phi(:, :, :)=0.  
563    
564      ! Ne pas affecter les valeurs entrees de u, v, h, et q      ! We will modify variables *_seri and we will not touch variables
565        ! u, v, t, qx:
566        t_seri = t
567        u_seri = u
568        v_seri = v
569        q_seri = qx(:, :, ivap)
570        ql_seri = qx(:, :, iliq)
571        tr_seri = qx(:, :, 3:nqmx)
572    
573      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 (nq >= 3) THEN  
        tr_seri(:, :, :nq-2) = qx(:, :, 3:nq)  
     ELSE  
        tr_seri(:, :, 1) = 0.  
     ENDIF  
   
     DO i = 1, klon  
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
574    
575      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
576         ztit='after dynamic'         tit = 'after dynamics'
577         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
578              , 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)
579              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)         ! Comme les tendances de la physique sont ajout\'es dans la
580         !     Comme les tendances de la physique sont ajoute dans la dynamique,         ! dynamique, la variation d'enthalpie par la dynamique devrait
581         !     on devrait avoir que la variation d'entalpie par la dynamique         ! \^etre \'egale \`a la variation de la physique au pas de temps
582         !     est egale a la variation de la physique au pas de temps precedent.         ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
583         !     Donc la somme de ces 2 variations devrait etre nulle.         ! nulle.
584         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
585              , zero_v, zero_v, zero_v, zero_v, zero_v &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
586              , zero_v, zero_v, zero_v, ztsol &              d_qt, 0.)
             , d_h_vcol+d_h_vcol_phy, d_qt, 0. &  
             , fs_bound, fq_bound )  
587      END IF      END IF
588    
589      ! Diagnostiquer la tendance dynamique      ! Diagnostic de la tendance dynamique :
   
590      IF (ancien_ok) THEN      IF (ancien_ok) THEN
591         DO k = 1, llm         DO k = 1, llm
592            DO i = 1, klon            DO i = 1, klon
593               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
594               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
595            ENDDO            ENDDO
596         ENDDO         ENDDO
597      ELSE      ELSE
598         DO k = 1, llm         DO k = 1, llm
599            DO i = 1, klon            DO i = 1, klon
600               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
601               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
602            ENDDO            ENDDO
603         ENDDO         ENDDO
604         ancien_ok = .TRUE.         ancien_ok = .TRUE.
605      ENDIF      ENDIF
606    
607      ! Ajouter le geopotentiel du sol:      ! Ajouter le geopotentiel du sol:
   
608      DO k = 1, llm      DO k = 1, llm
609         DO i = 1, klon         DO i = 1, klon
610            zphi(i, k) = pphi(i, k) + pphis(i)            zphi(i, k) = pphi(i, k) + pphis(i)
611         ENDDO         ENDDO
612      ENDDO      ENDDO
613    
614      ! Verifier les temperatures      ! Check temperatures:
   
615      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
616    
617      ! Incrementer le compteur de la physique      call increment_itap
618        julien = MOD(dayvrai, 360)
     itap = itap + 1  
     julien = MOD(NINT(rdayvrai), 360)  
619      if (julien == 0) julien = 360      if (julien == 0) julien = 360
620    
621      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
622    
623      ! Mettre en action les conditions aux limites (albedo, sst, etc.).      ! Prescrire l'ozone :
624      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      wo = ozonecm(REAL(julien), paprs)
625    
626      if (nq >= 5) then      ! \'Evaporation de l'eau liquide nuageuse :
627         wo = qx(:, :, 5) * zmasse / dobson_u / 1e3      DO k = 1, llm
     else IF (MOD(itap - 1, lmt_pas) == 0) THEN  
        CALL ozonecm(REAL(julien), rlat, paprs, wo)  
     ENDIF  
   
     ! Re-evaporer l'eau liquide nuageuse  
   
     DO k = 1, llm  ! re-evaporation de l'eau liquide nuageuse  
628         DO i = 1, klon         DO i = 1, klon
629            zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))            zb = MAX(0., ql_seri(i, k))
630            zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))            t_seri(i, k) = t_seri(i, k) &
631            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  
632            q_seri(i, k) = q_seri(i, k) + zb            q_seri(i, k) = q_seri(i, k) + zb
           ql_seri(i, k) = 0.0  
633         ENDDO         ENDDO
634      ENDDO      ENDDO
635        ql_seri = 0.
636    
637      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
638         ztit='after reevap'         tit = 'after reevap'
639         CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
640              , 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)
641              , 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, &
642         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 )  
   
643      END IF      END IF
644    
645      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
646        zxrugs = sum(frugs * pctsrf, dim = 2)
647    
648      DO i = 1, klon      ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
649         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  
650    
651      ! calculs necessaires au calcul de l'albedo dans l'interface      CALL orbite(REAL(julien), longi, dist)
   
     CALL orbite(REAL(julien), zlongi, dist)  
652      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
653         zdtime = pdtphys * REAL(radpas)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
        CALL zenang(zlongi, gmtime, zdtime, rmu0, fract)  
654      ELSE      ELSE
655         rmu0 = -999.999         mu0 = - 999.999
656      ENDIF      ENDIF
657    
658      !     Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
659      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  
660    
661      !     Repartition sous maille des flux LW et SW      ! R\'epartition sous maille des flux longwave et shortwave
662      ! Repartition du longwave par sous-surface linearisee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
663    
664      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
665         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
666            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
667                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
668            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))      END forall
        ENDDO  
     ENDDO  
669    
670      fder = dlw      fder = dlw
671    
672      CALL clmain(pdtphys, itap, date0, pctsrf, pctsrf_new, &      ! Couche limite:
673           t_seri, q_seri, u_seri, v_seri, &  
674           julien, rmu0, co2_ppm,  &      CALL clmain(dtphys, pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
675           ok_veget, ocean, npas, nexca, ftsol, &           ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
676           soil_model, cdmmax, cdhmax, &           paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, rain_fall, &
677           ksta, ksta_ter, ok_kzmin, ftsoil, qsol,  &           snow_fall, fsolsw, fsollw, fder, rlat, frugs, agesno, rugoro, &
678           paprs, pplay, fsnow, fqsurf, fevap, falbe, falblw, &           d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, &
679           fluxlat, rain_fall, snow_fall, &           fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, &
680           fsolsw, fsollw, sollwdown, fder, &           u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, &
681           rlon, rlat, cuphy, cvphy, frugs, &           trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
682           firstcal, lafin, agesno, rugoro, &  
683           d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &      ! Incr\'ementation des flux
684           fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &  
685           q2, dsens, devap, &      zxfluxt = 0.
686           ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &      zxfluxq = 0.
687           pblh, capCL, oliqCL, cteiCL, pblT, &      zxfluxu = 0.
688           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.  
689      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
690         DO k = 1, llm         DO k = 1, llm
691            DO i = 1, klon            DO i = 1, klon
692               zxfluxt(i, k) = zxfluxt(i, k) +  &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
693                    fluxt(i, k, nsrf) * pctsrf( i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
694               zxfluxq(i, k) = zxfluxq(i, k) +  &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
695                    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)  
696            END DO            END DO
697         END DO         END DO
698      END DO      END DO
699      DO i = 1, klon      DO i = 1, klon
700         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
701         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
702         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
703      ENDDO      ENDDO
704    
# Line 1130  contains Line 711  contains
711         ENDDO         ENDDO
712      ENDDO      ENDDO
713    
714      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
715         ztit='after clmain'         tit = 'after clmain'
716         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
717              , 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)
718              , 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, &
719         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 )  
720      END IF      END IF
721    
722      ! Incrementer la temperature du sol      ! Update surface temperature:
723    
724      DO i = 1, klon      DO i = 1, klon
725         zxtsol(i) = 0.0         zxfluxlat(i) = 0.
        zxfluxlat(i) = 0.0  
726    
727         zt2m(i) = 0.0         zt2m(i) = 0.
728         zq2m(i) = 0.0         zq2m(i) = 0.
729         zu10m(i) = 0.0         zu10m(i) = 0.
730         zv10m(i) = 0.0         zv10m(i) = 0.
731         zxffonte(i) = 0.0         zxffonte(i) = 0.
732         zxfqcalving(i) = 0.0         zxfqcalving(i) = 0.
733    
734         s_pblh(i) = 0.0         s_pblh(i) = 0.
735         s_lcl(i) = 0.0         s_lcl(i) = 0.
736         s_capCL(i) = 0.0         s_capCL(i) = 0.
737         s_oliqCL(i) = 0.0         s_oliqCL(i) = 0.
738         s_cteiCL(i) = 0.0         s_cteiCL(i) = 0.
739         s_pblT(i) = 0.0         s_pblT(i) = 0.
740         s_therm(i) = 0.0         s_therm(i) = 0.
741         s_trmb1(i) = 0.0         s_trmb1(i) = 0.
742         s_trmb2(i) = 0.0         s_trmb2(i) = 0.
743         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  
744      ENDDO      ENDDO
745    
746        call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
747    
748        ftsol = ftsol + d_ts
749        zxtsol = sum(ftsol * pctsrf, dim = 2)
750      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
751         DO i = 1, klon         DO i = 1, klon
752            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf) * pctsrf(i, nsrf)
753            zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
754            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)            zt2m(i) = zt2m(i) + t2m(i, nsrf) * pctsrf(i, nsrf)
755              zq2m(i) = zq2m(i) + q2m(i, nsrf) * pctsrf(i, nsrf)
756            zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)            zu10m(i) = zu10m(i) + u10m(i, nsrf) * pctsrf(i, nsrf)
757            zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)            zv10m(i) = zv10m(i) + v10m(i, nsrf) * pctsrf(i, nsrf)
758            zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf) * pctsrf(i, nsrf)
759            zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)            zxfqcalving(i) = zxfqcalving(i) + &
760            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)                 fqcalving(i, nsrf) * pctsrf(i, nsrf)
761            zxfqcalving(i) = zxfqcalving(i) +  &            s_pblh(i) = s_pblh(i) + pblh(i, nsrf) * pctsrf(i, nsrf)
762                 fqcalving(i, nsrf)*pctsrf(i, nsrf)            s_lcl(i) = s_lcl(i) + plcl(i, nsrf) * pctsrf(i, nsrf)
763            s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)            s_capCL(i) = s_capCL(i) + capCL(i, nsrf) * pctsrf(i, nsrf)
764            s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)            s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) * pctsrf(i, nsrf)
765            s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)            s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) * pctsrf(i, nsrf)
766            s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)            s_pblT(i) = s_pblT(i) + pblT(i, nsrf) * pctsrf(i, nsrf)
767            s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)            s_therm(i) = s_therm(i) + therm(i, nsrf) * pctsrf(i, nsrf)
768            s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)            s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) * pctsrf(i, nsrf)
769            s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)            s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) * pctsrf(i, nsrf)
770            s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)            s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) * pctsrf(i, nsrf)
           s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)  
771         ENDDO         ENDDO
772      ENDDO      ENDDO
773    
774      ! 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 :
   
775      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
776         DO i = 1, klon         DO i = 1, klon
777            IF (pctsrf(i, nsrf)  <  epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
778    
779            IF (pctsrf(i, nsrf)  <  epsfra) t2m(i, nsrf) = zt2m(i)            IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
780            IF (pctsrf(i, nsrf)  <  epsfra) q2m(i, nsrf) = zq2m(i)            IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
781            IF (pctsrf(i, nsrf)  <  epsfra) u10m(i, nsrf) = zu10m(i)            IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
782            IF (pctsrf(i, nsrf)  <  epsfra) v10m(i, nsrf) = zv10m(i)            IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
783            IF (pctsrf(i, nsrf)  <  epsfra) ffonte(i, nsrf) = zxffonte(i)            IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
784            IF (pctsrf(i, nsrf)  <  epsfra)  &            IF (pctsrf(i, nsrf) < epsfra) &
785                 fqcalving(i, nsrf) = zxfqcalving(i)                 fqcalving(i, nsrf) = zxfqcalving(i)
786            IF (pctsrf(i, nsrf)  <  epsfra) pblh(i, nsrf)=s_pblh(i)            IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
787            IF (pctsrf(i, nsrf)  <  epsfra) plcl(i, nsrf)=s_lcl(i)            IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
788            IF (pctsrf(i, nsrf)  <  epsfra) capCL(i, nsrf)=s_capCL(i)            IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
789            IF (pctsrf(i, nsrf)  <  epsfra) oliqCL(i, nsrf)=s_oliqCL(i)            IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
790            IF (pctsrf(i, nsrf)  <  epsfra) cteiCL(i, nsrf)=s_cteiCL(i)            IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
791            IF (pctsrf(i, nsrf)  <  epsfra) pblT(i, nsrf)=s_pblT(i)            IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
792            IF (pctsrf(i, nsrf)  <  epsfra) therm(i, nsrf)=s_therm(i)            IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
793            IF (pctsrf(i, nsrf)  <  epsfra) trmb1(i, nsrf)=s_trmb1(i)            IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
794            IF (pctsrf(i, nsrf)  <  epsfra) trmb2(i, nsrf)=s_trmb2(i)            IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
795            IF (pctsrf(i, nsrf)  <  epsfra) trmb3(i, nsrf)=s_trmb3(i)            IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
796         ENDDO         ENDDO
797      ENDDO      ENDDO
798    
799      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
800    
801      DO i = 1, klon      DO i = 1, klon
802         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
803      ENDDO      ENDDO
804    
805      ! Appeler la convection (au choix)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
806    
807      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  
808    
809         IF (.NOT. ok_gust) THEN      if (conv_emanuel) then
810            do i = 1, klon         da = 0.
811               wd(i)=0.0         mp = 0.
812            enddo         phi = 0.
813           CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
814                d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
815                upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, da, phi, mp)
816           snow_con = 0.
817           clwcon0 = qcondc
818           mfu = upwd + dnwd
819    
820           IF (thermcep) THEN
821              zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
822              zqsat = zqsat / (1. - retv * zqsat)
823           ELSE
824              zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
825         ENDIF         ENDIF
826    
827         ! Calcul des proprietes des nuages convectifs         ! Properties of convective clouds
828           clwcon0 = fact_cldcon * clwcon0
829         DO k = 1, llm         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
830            DO i = 1, klon              rnebcon0)
831               zx_t = t_seri(i, k)  
832               IF (thermcep) THEN         forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
833                  zdelta = MAX(0., SIGN(1., rtt-zx_t))         mfd = 0.
834                  zx_qs  = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)         pen_u = 0.
835                  zx_qs  = MIN(0.5, zx_qs)         pen_d = 0.
836                  zcor   = 1./(1.-retv*zx_qs)         pde_d = 0.
837                  zx_qs  = zx_qs*zcor         pde_u = 0.
838               ELSE      else
839                  IF (zx_t < t_coup) THEN         conv_q = d_q_dyn + d_q_vdf / dtphys
840                     zx_qs = qsats(zx_t)/pplay(i, k)         conv_t = d_t_dyn + d_t_vdf / dtphys
841                  ELSE         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
842                     zx_qs = qsatl(zx_t)/pplay(i, k)         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
843                  ENDIF              q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
844               ENDIF              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
845               zqsat(i, k)=zx_qs              mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
846            ENDDO              kdtop, pmflxr, pmflxs)
847         ENDDO         WHERE (rain_con < 0.) rain_con = 0.
848           WHERE (snow_con < 0.) snow_con = 0.
849         !   calcul des proprietes des nuages convectifs         ibas_con = llm + 1 - kcbot
850         clwcon0=fact_cldcon*clwcon0         itop_con = llm + 1 - kctop
851         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  
852    
853      DO k = 1, llm      DO k = 1, llm
854         DO i = 1, klon         DO i = 1, klon
# Line 1355  contains Line 859  contains
859         ENDDO         ENDDO
860      ENDDO      ENDDO
861    
862      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
863         ztit='after convect'         tit = 'after convect'
864         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
865              , 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)
866              , 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, &
867         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 )  
868      END IF      END IF
869    
870      IF (check) THEN      IF (check) THEN
871         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
872         print *,"aprescon=", za         print *, "aprescon = ", za
873         zx_t = 0.0         zx_t = 0.
874         za = 0.0         za = 0.
875         DO i = 1, klon         DO i = 1, klon
876            za = za + airephy(i)/REAL(klon)            za = za + airephy(i) / REAL(klon)
877            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
878                 snow_con(i))*airephy(i)/REAL(klon)                 snow_con(i)) * airephy(i) / REAL(klon)
879         ENDDO         ENDDO
880         zx_t = zx_t/za*pdtphys         zx_t = zx_t / za * dtphys
881         print *,"Precip=", zx_t         print *, "Precip = ", zx_t
882      ENDIF      ENDIF
883      IF (zx_ajustq) THEN  
884         DO i = 1, klon      IF (.not. conv_emanuel) THEN
885            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
886         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  
887         DO k = 1, llm         DO k = 1, llm
888            DO i = 1, klon            DO i = 1, klon
889               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  
890                  q_seri(i, k) = q_seri(i, k) * z_factor(i)                  q_seri(i, k) = q_seri(i, k) * z_factor(i)
891               ENDIF               ENDIF
892            ENDDO            ENDDO
893         ENDDO         ENDDO
894      ENDIF      ENDIF
     zx_ajustq=.FALSE.  
895    
896      ! Convection seche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
897    
898      d_t_ajs=0.      d_t_ajs = 0.
899      d_u_ajs=0.      d_u_ajs = 0.
900      d_v_ajs=0.      d_v_ajs = 0.
901      d_q_ajs=0.      d_q_ajs = 0.
902      fm_therm=0.      fm_therm = 0.
903      entr_therm=0.      entr_therm = 0.
904    
905      IF(prt_level>9)print *, &      if (iflag_thermals == 0) then
906           'AVANT LA CONVECTION SECHE, iflag_thermals=' &         ! Ajustement sec
907           , 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)  
908         t_seri = t_seri + d_t_ajs         t_seri = t_seri + d_t_ajs
909         q_seri = q_seri + d_q_ajs         q_seri = q_seri + d_q_ajs
910      else      else
911         !  Thermiques         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
912         IF(prt_level>9)print *,'JUSTE AVANT, iflag_thermals=' &              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
             , iflag_thermals, '   nsplit_thermals=', nsplit_thermals  
        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)  
913      endif      endif
914    
915      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
916         ztit='after dry_adjust'         tit = 'after dry_adjust'
917         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
918              , 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)  
919      END IF      END IF
920    
921      !  Caclul des ratqs      ! Caclul des ratqs
922    
923      !   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
924      !   on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
925      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
926         do k=1, llm         do k = 1, llm
927            do i=1, klon            do i = 1, klon
928               if(ptconv(i, k)) then               if(ptconv(i, k)) then
929                  ratqsc(i, k)=ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
930                       +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)
931               else               else
932                  ratqsc(i, k)=0.                  ratqsc(i, k) = 0.
933               endif               endif
934            enddo            enddo
935         enddo         enddo
936      endif      endif
937    
938      !   ratqs stables      ! ratqs stables
939      do k=1, llm      do k = 1, llm
940         do i=1, klon         do i = 1, klon
941            ratqss(i, k)=ratqsbas+(ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
942                 min((paprs(i, 1)-pplay(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
943         enddo         enddo
944      enddo      enddo
945    
946      !  ratqs final      ! ratqs final
947      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
948         !   les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
949         !   ratqs final         ! ratqs final
950         !   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
951         !   relaxation des ratqs         ! relaxation des ratqs
952         facteur=exp(-pdtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
953         ratqs=max(ratqs*facteur, ratqss)         ratqs = max(ratqs, ratqsc)
        ratqs=max(ratqs, ratqsc)  
954      else      else
955         !   on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
956         ratqs=ratqss         ratqs = ratqss
957      endif      endif
958    
959      ! Appeler le processus de condensation a grande echelle      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
960      ! et le processus de precipitation           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
961      CALL fisrtilp(pdtphys, paprs, pplay, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
962           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)  
963    
964      WHERE (rain_lsc < 0) rain_lsc = 0.      WHERE (rain_lsc < 0) rain_lsc = 0.
965      WHERE (snow_lsc < 0) snow_lsc = 0.      WHERE (snow_lsc < 0) snow_lsc = 0.
# Line 1505  contains Line 973  contains
973         ENDDO         ENDDO
974      ENDDO      ENDDO
975      IF (check) THEN      IF (check) THEN
976         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
977         print *,"apresilp=", za         print *, "apresilp = ", za
978         zx_t = 0.0         zx_t = 0.
979         za = 0.0         za = 0.
980         DO i = 1, klon         DO i = 1, klon
981            za = za + airephy(i)/REAL(klon)            za = za + airephy(i) / REAL(klon)
982            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
983                 + snow_lsc(i))*airephy(i)/REAL(klon)                 + snow_lsc(i)) * airephy(i) / REAL(klon)
984         ENDDO         ENDDO
985         zx_t = zx_t/za*pdtphys         zx_t = zx_t / za * dtphys
986         print *,"Precip=", zx_t         print *, "Precip = ", zx_t
987      ENDIF      ENDIF
988    
989      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
990         ztit='after fisrt'         tit = 'after fisrt'
991         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
992              , 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)
993              , 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, &
994         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 )  
995      END IF      END IF
996    
997      !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
998    
999      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1000    
1001      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke      IF (iflag_cldcon <= - 1) THEN
1002         snow_tiedtke=0.         ! seulement pour Tiedtke
1003         if (iflag_cldcon == -1) then         snow_tiedtke = 0.
1004            rain_tiedtke=rain_con         if (iflag_cldcon == - 1) then
1005              rain_tiedtke = rain_con
1006         else         else
1007            rain_tiedtke=0.            rain_tiedtke = 0.
1008            do k=1, llm            do k = 1, llm
1009               do i=1, klon               do i = 1, klon
1010                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1011                     rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i, k)/pdtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k) / dtphys &
1012                          *zmasse(i, k)                          * zmasse(i, k)
1013                  endif                  endif
1014               enddo               enddo
1015            enddo            enddo
1016         endif         endif
1017    
1018         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1019         CALL diagcld1(paprs, pplay, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1020              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1021         DO k = 1, llm         DO k = 1, llm
1022            DO i = 1, klon            DO i = 1, klon
1023               IF (diafra(i, k).GT.cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
1024                  cldliq(i, k) = dialiq(i, k)                  cldliq(i, k) = dialiq(i, k)
1025                  cldfra(i, k) = diafra(i, k)                  cldfra(i, k) = diafra(i, k)
1026               ENDIF               ENDIF
1027            ENDDO            ENDDO
1028         ENDDO         ENDDO
   
1029      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1030         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1031         ! 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
1032         ! facttemps         ! d'un facteur facttemps.
1033         facteur = pdtphys *facttemps         facteur = dtphys * facttemps
1034         do k=1, llm         do k = 1, llm
1035            do i=1, klon            do i = 1, klon
1036               rnebcon(i, k)=rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1037               if (rnebcon0(i, k)*clwcon0(i, k).gt.rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1038                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1039                  rnebcon(i, k)=rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1040                  clwcon(i, k)=clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1041               endif               endif
1042            enddo            enddo
1043         enddo         enddo
1044    
1045         !   On prend la somme des fractions nuageuses et des contenus en eau         ! On prend la somme des fractions nuageuses et des contenus en eau
1046         cldfra=min(max(cldfra, rnebcon), 1.)         cldfra = min(max(cldfra, rnebcon), 1.)
1047         cldliq=cldliq+rnebcon*clwcon         cldliq = cldliq + rnebcon * clwcon
   
1048      ENDIF      ENDIF
1049    
1050      ! 2. NUAGES STARTIFORMES      ! 2. Nuages stratiformes
1051    
1052      IF (ok_stratus) THEN      IF (ok_stratus) THEN
1053         CALL diagcld2(paprs, pplay, t_seri, q_seri, diafra, dialiq)         CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1054         DO k = 1, llm         DO k = 1, llm
1055            DO i = 1, klon            DO i = 1, klon
1056               IF (diafra(i, k).GT.cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
1057                  cldliq(i, k) = dialiq(i, k)                  cldliq(i, k) = dialiq(i, k)
1058                  cldfra(i, k) = diafra(i, k)                  cldfra(i, k) = diafra(i, k)
1059               ENDIF               ENDIF
# Line 1600  contains Line 1062  contains
1062      ENDIF      ENDIF
1063    
1064      ! Precipitation totale      ! Precipitation totale
   
1065      DO i = 1, klon      DO i = 1, klon
1066         rain_fall(i) = rain_con(i) + rain_lsc(i)         rain_fall(i) = rain_con(i) + rain_lsc(i)
1067         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
1068      ENDDO      ENDDO
1069    
1070      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1071         ztit="after diagcld"           dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1072         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  
1073    
1074        ! Humidit\'e relative pour diagnostic :
1075      DO k = 1, llm      DO k = 1, llm
1076         DO i = 1, klon         DO i = 1, klon
1077            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1078            IF (thermcep) THEN            IF (thermcep) THEN
1079               zdelta = MAX(0., SIGN(1., rtt-zx_t))               zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t) / play(i, k)
1080               zx_qs  = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)               zx_qs = MIN(0.5, zx_qs)
1081               zx_qs  = MIN(0.5, zx_qs)               zcor = 1. / (1. - retv * zx_qs)
1082               zcor   = 1./(1.-retv*zx_qs)               zx_qs = zx_qs * zcor
              zx_qs  = zx_qs*zcor  
1083            ELSE            ELSE
1084               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
1085                  zx_qs = qsats(zx_t)/pplay(i, k)                  zx_qs = qsats(zx_t) / play(i, k)
1086               ELSE               ELSE
1087                  zx_qs = qsatl(zx_t)/pplay(i, k)                  zx_qs = qsatl(zx_t) / play(i, k)
1088               ENDIF               ENDIF
1089            ENDIF            ENDIF
1090            zx_rh(i, k) = q_seri(i, k)/zx_qs            zx_rh(i, k) = q_seri(i, k) / zx_qs
1091            zqsat(i, k)=zx_qs            zqsat(i, k) = zx_qs
1092         ENDDO         ENDDO
1093      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  
1094    
1095      ! Calculer les parametres optiques des nuages et quelques      ! Introduce the aerosol direct and first indirect radiative forcings:
1096      ! parametres pour diagnostiques:      tau_ae = 0.
1097        piz_ae = 0.
1098        cg_ae = 0.
1099    
1100        ! Param\`etres optiques des nuages et quelques param\`etres pour
1101        ! diagnostics :
1102      if (ok_newmicro) then      if (ok_newmicro) then
1103         CALL newmicro (paprs, pplay, ok_newmicro, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1104              t_seri, cldliq, cldfra, cldtau, cldemi, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1105              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)  
1106      else      else
1107         CALL nuage (paprs, pplay, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1108              t_seri, cldliq, cldfra, cldtau, cldemi, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1109              cldh, cldl, cldm, cldt, cldq, &              bl95_b1, cldtaupi, re, fl)
             ok_aie, &  
             sulfate, sulfate_pi, &  
             bl95_b0, bl95_b1, &  
             cldtaupi, re, fl)  
   
1110      endif      endif
1111    
1112      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
1113           ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1114      IF (MOD(itaprad, radpas) == 0) THEN         ! Calcul de l'abedo moyen par maille
1115         DO i = 1, klon         albsol = sum(falbe * pctsrf, dim = 2)
1116            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &  
1117                 + falbe(i, is_lic) * pctsrf(i, is_lic) &         ! Rayonnement (compatible Arpege-IFS) :
1118                 + falbe(i, is_ter) * pctsrf(i, is_ter) &         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1119                 + falbe(i, is_sic) * pctsrf(i, is_sic)              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1120            albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1121                 + falblw(i, is_lic) * pctsrf(i, is_lic) &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1122                 + falblw(i, is_ter) * pctsrf(i, is_ter) &              swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1123                 + 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  
1124      ENDIF      ENDIF
     itaprad = itaprad + 1  
1125    
1126      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
1127    
1128      DO k = 1, llm      DO k = 1, llm
1129         DO i = 1, klon         DO i = 1, klon
1130            t_seri(i, k) = t_seri(i, k) &            t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys &
1131                 + (heat(i, k)-cool(i, k)) * pdtphys/86400.                 / 86400.
1132         ENDDO         ENDDO
1133      ENDDO      ENDDO
1134    
1135      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1136         ztit='after rad'         tit = 'after rad'
1137         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1138              , 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)
1139              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1140         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 )  
1141      END IF      END IF
1142    
1143      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
   
1144      DO i = 1, klon      DO i = 1, klon
1145         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1146         zxsnow(i) = 0.0         zxsnow(i) = 0.
1147      ENDDO      ENDDO
1148      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1149         DO i = 1, klon         DO i = 1, klon
1150            zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)            zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf) * pctsrf(i, nsrf)
1151            zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)            zxsnow(i) = zxsnow(i) + fsnow(i, nsrf) * pctsrf(i, nsrf)
1152         ENDDO         ENDDO
1153      ENDDO      ENDDO
1154    
1155      ! Calculer le bilan du sol et la derive de temperature (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1156    
1157      DO i = 1, klon      DO i = 1, klon
1158         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1159      ENDDO      ENDDO
1160    
1161      !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:  
1162    
1163      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1164         !  selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
1165         igwd=0         igwd = 0
1166         DO i=1, klon         DO i = 1, klon
1167            itest(i)=0            itest(i) = 0
1168            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
1169               itest(i)=1               itest(i) = 1
1170               igwd=igwd+1               igwd = igwd + 1
              idx(igwd)=i  
1171            ENDIF            ENDIF
1172         ENDDO         ENDDO
1173    
1174         CALL drag_noro(klon, llm, pdtphys, paprs, pplay, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1175              zmea, zstd, zsig, zgam, zthe, zpic, zval, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1176              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)  
1177    
1178         !  ajout des tendances         ! ajout des tendances
1179         DO k = 1, llm         DO k = 1, llm
1180            DO i = 1, klon            DO i = 1, klon
1181               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 1781  contains Line 1186  contains
1186      ENDIF      ENDIF
1187    
1188      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1189           ! S\'election des points pour lesquels le sch\'ema est actif :
1190         !  selection des points pour lesquels le shema est actif:         igwd = 0
1191         igwd=0         DO i = 1, klon
1192         DO i=1, klon            itest(i) = 0
1193            itest(i)=0            IF (zpic(i) - zmea(i) > 100.) THEN
1194            IF ((zpic(i)-zmea(i)).GT.100.) THEN               itest(i) = 1
1195               itest(i)=1               igwd = igwd + 1
              igwd=igwd+1  
              idx(igwd)=i  
1196            ENDIF            ENDIF
1197         ENDDO         ENDDO
1198    
1199         CALL lift_noro(klon, llm, pdtphys, paprs, pplay, &         CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1200              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, &  
1201              d_t_lif, d_u_lif, d_v_lif)              d_t_lif, d_u_lif, d_v_lif)
1202    
1203         !  ajout des tendances         ! Ajout des tendances :
1204         DO k = 1, llm         DO k = 1, llm
1205            DO i = 1, klon            DO i = 1, klon
1206               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 1808  contains Line 1208  contains
1208               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)
1209            ENDDO            ENDDO
1210         ENDDO         ENDDO
1211        ENDIF
1212    
1213      ENDIF ! fin de test sur ok_orolf      ! Stress n\'ecessaires : toute la physique
   
     ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE  
1214    
1215      DO i = 1, klon      DO i = 1, klon
1216         zustrph(i)=0.         zustrph(i) = 0.
1217         zvstrph(i)=0.         zvstrph(i) = 0.
1218      ENDDO      ENDDO
1219      DO k = 1, llm      DO k = 1, llm
1220         DO i = 1, klon         DO i = 1, klon
1221            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 &
1222            zvstrph(i)=zvstrph(i)+(v_seri(i, k)-v(i, k))/pdtphys* zmasse(i, k)                 * zmasse(i, k)
1223              zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1224                   * zmasse(i, k)
1225         ENDDO         ENDDO
1226      ENDDO      ENDDO
1227    
1228      !IM calcul composantes axiales du moment angulaire et couple des montagnes      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1229             zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1230    
1231      CALL aaam_bud(27, klon, llm, gmtime, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1232           ra, rg, romega, &           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1233           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  
   
     !AA Installation de l'interface online-offline pour traceurs  
1234    
1235      !   Calcul  des tendances traceurs      ! Calcul des tendances traceurs
1236        call phytrac(julien, time, firstcal, lafin, dtphys, t, paprs, play, mfu, &
1237             mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, &
1238             pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, tr_seri, &
1239             zmasse, ncid_startphy)
1240    
1241      call phytrac(rnpb, itap, lmt_pas, julien,  gmtime, firstcal, lafin, nq-2, &      IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1242           pdtphys, u, v, t, paprs, pplay, pmfu,  pmfd,  pen_u,  pde_u,  pen_d, &           pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1243           pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &           frac_impa, frac_nucl, pphis, airephy, dtphys)
          frac_impa,  frac_nucl, presnivs, pphis, pphi, albsol, rhcl, cldfra, &  
          rneb,  diafra,  cldliq, itop_con, ibas_con, pmflxr, pmflxs, prfl, &  
          psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)  
   
     IF (offline) THEN  
   
        print*, 'Attention on met a 0 les thermiques pour phystoke'  
        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  
1244    
1245      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1246        CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
1247    
1248      CALL transp (paprs, zxtsol, &      ! diag. bilKP
          t_seri, q_seri, u_seri, v_seri, zphi, &  
          ve, vq, ue, uq)  
1249    
1250      !IM diag. bilKP      CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
   
     CALL transp_lay (paprs, zxtsol, &  
          t_seri, q_seri, u_seri, v_seri, zphi, &  
1251           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1252    
1253      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
1254    
1255      !+jld ec_conser      ! conversion Ec en énergie thermique
1256      DO k = 1, llm      DO k = 1, llm
1257         DO i = 1, klon         DO i = 1, klon
1258            ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i, k))            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1259            d_t_ec(i, k)=0.5/ZRCPD &            d_t_ec(i, k) = 0.5 / ZRCPD &
1260                 *(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)
1261            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)
1262            d_t_ec(i, k) = d_t_ec(i, k)/pdtphys            d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1263         END DO         END DO
1264      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  
1265    
1266        IF (if_ebil >= 1) THEN
1267           tit = 'after physic'
1268           CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1269                ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1270           ! Comme les tendances de la physique sont ajoute dans la dynamique,
1271           ! on devrait avoir que la variation d'entalpie par la dynamique
1272           ! est egale a la variation de la physique au pas de temps precedent.
1273           ! Donc la somme de ces 2 variations devrait etre nulle.
1274           call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1275                evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
1276           d_h_vcol_phy = d_h_vcol
1277      END IF      END IF
1278    
1279      !   SORTIES      ! SORTIES
1280    
1281      !cc prw = eau precipitable      ! prw = eau precipitable
1282      DO i = 1, klon      DO i = 1, klon
1283         prw(i) = 0.         prw(i) = 0.
1284         DO k = 1, llm         DO k = 1, llm
1285            prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)            prw(i) = prw(i) + q_seri(i, k) * zmasse(i, k)
1286         ENDDO         ENDDO
1287      ENDDO      ENDDO
1288    
# Line 1922  contains Line 1290  contains
1290    
1291      DO k = 1, llm      DO k = 1, llm
1292         DO i = 1, klon         DO i = 1, klon
1293            d_u(i, k) = ( u_seri(i, k) - u(i, k) ) / pdtphys            d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1294            d_v(i, k) = ( v_seri(i, k) - v(i, k) ) / pdtphys            d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1295            d_t(i, k) = ( t_seri(i, k)-t(i, k) ) / pdtphys            d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1296            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
1297            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
1298         ENDDO         ENDDO
1299      ENDDO      ENDDO
1300    
1301      IF (nq >= 3) THEN      DO iq = 3, nqmx
1302         DO iq = 3, nq         DO k = 1, llm
1303            DO  k = 1, llm            DO i = 1, klon
1304               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  
1305            ENDDO            ENDDO
1306         ENDDO         ENDDO
1307      ENDIF      ENDDO
1308    
1309      ! 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:
1310      DO k = 1, llm      DO k = 1, llm
# Line 1948  contains Line 1314  contains
1314         ENDDO         ENDDO
1315      ENDDO      ENDDO
1316    
1317      !   Ecriture des sorties      CALL histwrite_phy("phis", pphis)
1318      call write_histhf      CALL histwrite_phy("aire", airephy)
1319      call write_histday      CALL histwrite_phy("psol", paprs(:, 1))
1320      call write_histins      CALL histwrite_phy("precip", rain_fall + snow_fall)
1321        CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1322      ! Si c'est la fin, il faut conserver l'etat de redemarrage      CALL histwrite_phy("pluc", rain_con + snow_con)
1323      IF (lafin) THEN      CALL histwrite_phy("tsol", zxtsol)
1324         itau_phy = itau_phy + itap      CALL histwrite_phy("t2m", zt2m)
1325         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, &      CALL histwrite_phy("q2m", zq2m)
1326              ftsoil, tslab, seaice, fqsurf, qsol, &      CALL histwrite_phy("u10m", zu10m)
1327              fsnow, falbe, falblw, fevap, rain_fall, snow_fall, &      CALL histwrite_phy("v10m", zv10m)
1328              solsw, sollwdown, dlw, &      CALL histwrite_phy("snow", snow_fall)
1329              radsol, frugs, agesno, &      CALL histwrite_phy("cdrm", cdragm)
1330              zmea, zstd, zsig, zgam, zthe, zpic, zval, &      CALL histwrite_phy("cdrh", cdragh)
1331              t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)      CALL histwrite_phy("topl", toplw)
1332      ENDIF      CALL histwrite_phy("evap", evap)
1333        CALL histwrite_phy("sols", solsw)
1334    contains      CALL histwrite_phy("soll", sollw)
1335        CALL histwrite_phy("solldown", sollwdown)
1336      subroutine write_histday      CALL histwrite_phy("bils", bils)
1337        CALL histwrite_phy("sens", - sens)
1338        use grid_change, only: gr_phy_write_3d      CALL histwrite_phy("fder", fder)
1339        integer itau_w  ! pas de temps ecriture      CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1340        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1341        !------------------------------------------------      CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1342        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
       if (ok_journe) THEN  
          itau_w = itau_phy + itap  
          if (nq <= 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:  
1343    
1344           zsto = pdtphys * ecrit_ins      DO nsrf = 1, nbsrf
1345           zout = pdtphys * ecrit_ins         CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf) * 100.)
1346           itau_w = itau_phy + itap         CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1347           CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1348           i = NINT(zout/zsto)         CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1349           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), pphis, zx_tmp_2d)         CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1350           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)         CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1351           CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1352           i = NINT(zout/zsto)         CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1353           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), airephy, zx_tmp_2d)         CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1354           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)      END DO
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = paprs(i, 1)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)  
   
          DO i = 1, klon  
             zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)  
          ENDDO  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxtsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)  
          !ccIM  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zt2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zq2m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zu10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zv10m, zx_tmp_2d)  
          CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), snow_fall, zx_tmp_2d)  
          CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragm, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragh, zx_tmp_2d)  
          CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), toplw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), evap, zx_tmp_2d)  
          CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), solsw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollwdown, zx_tmp_2d)  
          CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), bils, zx_tmp_2d)  
          CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)  
   
          zx_tmp_fi2d(1:klon)=-1*sens(1:klon)  
          !     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sens, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
          CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), fder, zx_tmp_2d)  
          CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_oce), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_ter), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_lic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_sic), zx_tmp_2d)  
          CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)  
   
          DO nsrf = 1, nbsrf  
             !XXX  
             zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)*100.  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
             zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)  
             CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)  
             CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &  
                  zx_tmp_2d)  
   
          END DO  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsol, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)  
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)  
   
          CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxrugs, zx_tmp_2d)  
          CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)  
   
          !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  
1355    
1356        if (ok_sync) then      CALL histwrite_phy("albs", albsol)
1357           call histsync(nid_hf3d)      CALL histwrite_phy("rugs", zxrugs)
1358        endif      CALL histwrite_phy("s_pblh", s_pblh)
1359        CALL histwrite_phy("s_pblt", s_pblt)
1360        CALL histwrite_phy("s_lcl", s_lcl)
1361        CALL histwrite_phy("s_capCL", s_capCL)
1362        CALL histwrite_phy("s_oliqCL", s_oliqCL)
1363        CALL histwrite_phy("s_cteiCL", s_cteiCL)
1364        CALL histwrite_phy("s_therm", s_therm)
1365        CALL histwrite_phy("s_trmb1", s_trmb1)
1366        CALL histwrite_phy("s_trmb2", s_trmb2)
1367        CALL histwrite_phy("s_trmb3", s_trmb3)
1368        if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1369        CALL histwrite_phy("temp", t_seri)
1370        CALL histwrite_phy("vitu", u_seri)
1371        CALL histwrite_phy("vitv", v_seri)
1372        CALL histwrite_phy("geop", zphi)
1373        CALL histwrite_phy("pres", play)
1374        CALL histwrite_phy("dtvdf", d_t_vdf)
1375        CALL histwrite_phy("dqvdf", d_q_vdf)
1376        CALL histwrite_phy("rhum", zx_rh)
1377    
1378        if (ok_instan) call histsync(nid_ins)
1379    
1380        IF (lafin) then
1381           call NF95_CLOSE(ncid_startphy)
1382           CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1383                fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1384                radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1385                t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1386                w01)
1387        end IF
1388    
1389      end subroutine write_histhf3d      firstcal = .FALSE.
1390    
1391    END SUBROUTINE physiq    END SUBROUTINE physiq
1392    

Legend:
Removed from v.18  
changed lines
  Added in v.202

  ViewVC Help
Powered by ViewVC 1.1.21