/[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 13 by guez, Fri Jul 25 19:59:34 2008 UTC trunk/phylmd/physiq.f revision 90 by guez, Wed Mar 12 21:16:36 2014 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, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &
8         pplay, pphi, pphis, presnivs, u, v, t, qx, omega, d_u, d_v, &         u, v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps)
9         d_t, d_qx, d_ps, dudyn, PVteta)  
10        ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11      ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28      ! (subversion revision 678)
12    
13      ! Author : Z.X. Li (LMD/CNRS), date: 1993/08/18      ! Author: Z.X. Li (LMD/CNRS) 1993
14    
15      ! Objet: Moniteur general de la physique du modele      ! This is the main procedure for the "physics" part of the program.
16      !AA      Modifications quant aux traceurs :  
17      !AA                  -  uniformisation des parametrisations ds phytrac      use aaam_bud_m, only: aaam_bud
18      !AA                  -  stockage des moyennes des champs necessaires      USE abort_gcm_m, ONLY: abort_gcm
19      !AA                     en mode traceur off-line      use aeropt_m, only: aeropt
20        use ajsec_m, only: ajsec
21      USE ioipsl, only: ymds2ju, histwrite, histsync      USE calendar, ONLY: ymds2ju
22      use dimens_m, only: jjm, iim, llm      use calltherm_m, only: calltherm
23      use indicesol, only: nbsrf, is_ter, is_lic, is_sic, is_oce, &      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &
24           clnsurf, epsfra           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin
25      use dimphy, only: klon, nbtr      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
26      use conf_gcm_m, only: raz_date, offline, iphysiq           ok_orodr, ok_orolf, soil_model
27      use dimsoil, only: nsoilmx      USE clmain_m, ONLY: clmain
28      use temps, only: itau_phy, day_ref, annee_ref, itaufin      use clouds_gno_m, only: clouds_gno
29      use clesphys, only: ecrit_hf, ecrit_hf2mth, &      USE comgeomphy, ONLY: airephy, cuphy, cvphy
30           ecrit_ins, ecrit_mth, ecrit_day, &      USE concvl_m, ONLY: concvl
31           cdmmax, cdhmax, &      USE conf_gcm_m, ONLY: offline, raz_date
32           co2_ppm, ecrit_reg, ecrit_tra, ksta, ksta_ter, &      USE conf_phys_m, ONLY: conf_phys
33           ok_kzmin      use conflx_m, only: conflx
34      use clesphys2, only: iflag_con, ok_orolf, ok_orodr, nbapp_rad, &      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
35           cycle_diurne, new_oliq, soil_model      use diagcld2_m, only: diagcld2
36      use iniprint, only: prt_level      use diagetpq_m, only: diagetpq
37      use abort_gcm_m, only: abort_gcm      use diagphy_m, only: diagphy
38      use YOMCST, only: rcpd, rtt, rlvtt, rg, ra, rsigma, retv, romega      USE dimens_m, ONLY: llm, nqmx
39      use comgeomphy      USE dimphy, ONLY: klon, nbtr
40      use ctherm      USE dimsoil, ONLY: nsoilmx
41      use phytrac_m, only: phytrac      use drag_noro_m, only: drag_noro
42      use oasis_m      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
43      use radepsi      use fisrtilp_m, only: fisrtilp
44      use radopt      USE hgardfou_m, ONLY: hgardfou
45      use yoethf      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
46      use ini_hist, only: ini_histhf, ini_histday, ini_histins           nbsrf
47      use orbite_m, only: orbite, zenang      USE ini_histins_m, ONLY: ini_histins
48      use phyetat0_m, only: phyetat0, rlat, rlon      use newmicro_m, only: newmicro
49      use hgardfou_m, only: hgardfou      USE oasis_m, ONLY: ok_oasis
50      use conf_phys_m, only: conf_phys      USE orbite_m, ONLY: orbite, zenang
51        USE ozonecm_m, ONLY: ozonecm
52      ! Declaration des constantes et des fonctions thermodynamiques :      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
53      use fcttre, only: thermcep, foeew, qsats, qsatl      USE phyredem_m, ONLY: phyredem
54        USE phystokenc_m, ONLY: phystokenc
55      ! Variables argument:      USE phytrac_m, ONLY: phytrac
56        USE qcheck_m, ONLY: qcheck
57      INTEGER nq ! input nombre de traceurs (y compris vapeur d'eau)      use radlwsw_m, only: radlwsw
58      REAL, intent(in):: rdayvrai ! input numero du jour de l'experience      use readsulfate_m, only: readsulfate
59      REAL, intent(in):: gmtime ! heure de la journée en fraction de jour      use sugwd_m, only: sugwd
60      REAL, intent(in):: pdtphys ! pas d'integration pour la physique (seconde)      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
61      LOGICAL, intent(in):: firstcal ! first call to "calfis"      USE temps, ONLY: annee_ref, day_ref, itau_phy
62        use unit_nml_m, only: unit_nml
63        USE yoethf_m, ONLY: r2es, rvtmp2
64    
65        ! Arguments:
66    
67        REAL, intent(in):: rdayvrai
68        ! (elapsed time since January 1st 0h of the starting year, in days)
69    
70        REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
71        REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)
72      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
73    
74      REAL, intent(in):: paprs(klon, llm+1)      REAL, intent(in):: paprs(klon, llm + 1)
75      ! (pression pour chaque inter-couche, en Pa)      ! (pression pour chaque inter-couche, en Pa)
76        
77      REAL, intent(in):: pplay(klon, llm)      REAL, intent(in):: play(klon, llm)
78      ! (input pression pour le mileu de chaque couche (en Pa))      ! (input pression pour le mileu de chaque couche (en Pa))
79    
80      REAL pphi(klon, llm)        REAL, intent(in):: pphi(klon, llm)
81      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! (input geopotentiel de chaque couche (g z) (reference sol))
82    
83      REAL pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol
84    
85      REAL presnivs(llm)      REAL, intent(in):: u(klon, llm)
86      ! (input pressions approximat. des milieux couches ( en PA))      ! vitesse dans la direction X (de O a E) en m/s
87    
88      REAL u(klon, llm)  ! input vitesse dans la direction X (de O a E) en m/s      REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s
89      REAL v(klon, llm)  ! input vitesse Y (de S a N) en m/s      REAL, intent(in):: t(klon, llm) ! input temperature (K)
     REAL t(klon, llm)  ! input temperature (K)  
   
     REAL, intent(in):: qx(klon, llm, nq)  
     ! (humidite specifique (kg/kg) et fractions massiques des autres traceurs)  
   
     REAL omega(klon, llm)  ! input vitesse verticale en Pa/s  
     REAL d_u(klon, llm)  ! output tendance physique de "u" (m/s/s)  
     REAL d_v(klon, llm)  ! output tendance physique de "v" (m/s/s)  
     REAL d_t(klon, llm)  ! output tendance physique de "t" (K/s)  
     REAL d_qx(klon, llm, nq)  ! output tendance physique de "qx" (kg/kg/s)  
     REAL d_ps(klon)  ! output tendance physique de la pression au sol  
90    
91      INTEGER nbteta      REAL, intent(in):: qx(klon, llm, nqmx)
92      PARAMETER(nbteta=3)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
93    
94        REAL omega(klon, llm) ! input vitesse verticale en Pa/s
95        REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)
96        REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/s)
97        REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)
98        REAL d_qx(klon, llm, nqmx) ! output tendance physique de "qx" (kg/kg/s)
99        REAL d_ps(klon) ! output tendance physique de la pression au sol
100    
101      REAL PVteta(klon, nbteta)      LOGICAL:: firstcal = .true.
102      ! (output vorticite potentielle a des thetas constantes)  
103        INTEGER nbteta
104        PARAMETER(nbteta = 3)
105    
     LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE  
     PARAMETER (ok_cvl=.TRUE.)  
106      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
107      PARAMETER (ok_gust=.FALSE.)      PARAMETER (ok_gust = .FALSE.)
108    
109      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL check ! Verifier la conservation du modele en eau
110      PARAMETER (check=.FALSE.)      PARAMETER (check = .FALSE.)
111      LOGICAL ok_stratus ! Ajouter artificiellement les stratus  
112      PARAMETER (ok_stratus=.FALSE.)      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
113        ! Ajouter artificiellement les stratus
114    
115      ! Parametres lies au coupleur OASIS:      ! Parametres lies au coupleur OASIS:
116      INTEGER, SAVE :: npas, nexca      INTEGER, SAVE:: npas, nexca
117      logical rnpb      logical rnpb
118      parameter(rnpb=.true.)      parameter(rnpb = .true.)
119    
120      character(len=6), save:: ocean      character(len = 6):: ocean = 'force '
121      ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")      ! (type de mod\`ele oc\'ean \`a utiliser: "force" ou "slab" mais pas "couple")
122    
123      logical ok_ocean      ! "slab" ocean
124      SAVE ok_ocean      REAL, save:: tslab(klon) ! temperature of ocean slab
125        REAL, save:: seaice(klon) ! glace de mer (kg/m2)
126      !IM "slab" ocean      REAL fluxo(klon) ! flux turbulents ocean-glace de mer
127      REAL tslab(klon)    !Temperature du slab-ocean      REAL fluxg(klon) ! flux turbulents ocean-atmosphere
     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  
128    
129      ! Modele thermique du sol, a activer pour le cycle diurne:      ! Modele thermique du sol, a activer pour le cycle diurne:
130      logical ok_veget      logical:: ok_veget = .false. ! type de modele de vegetation utilise
     save ok_veget  
     LOGICAL ok_journe ! sortir le fichier journalier  
     save ok_journe  
   
     LOGICAL ok_mensuel ! sortir le fichier mensuel  
131    
132      LOGICAL ok_instan ! sortir le fichier instantane      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
133      save ok_instan      ! sorties journalieres, mensuelles et instantanees dans les
134        ! fichiers histday, histmth et histins
135    
136      LOGICAL ok_region ! sortir le fichier regional      LOGICAL ok_region ! sortir le fichier regional
137      PARAMETER (ok_region=.FALSE.)      PARAMETER (ok_region = .FALSE.)
138    
139      !     pour phsystoke avec thermiques      ! pour phsystoke avec thermiques
140      REAL fm_therm(klon, llm+1)      REAL fm_therm(klon, llm + 1)
141      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
142      real q2(klon, llm+1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
     save q2  
143    
144      INTEGER ivap          ! indice de traceurs pour vapeur d'eau      INTEGER ivap ! indice de traceurs pour vapeur d'eau
145      PARAMETER (ivap=1)      PARAMETER (ivap = 1)
146      INTEGER iliq          ! indice de traceurs pour eau liquide      INTEGER iliq ! indice de traceurs pour eau liquide
147      PARAMETER (iliq=2)      PARAMETER (iliq = 2)
148    
149      REAL t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
150      SAVE t_ancien, q_ancien      LOGICAL, save:: ancien_ok
     LOGICAL ancien_ok  
     SAVE ancien_ok  
151    
152      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
153      REAL d_q_dyn(klon, llm)  ! tendance dynamique pour "q" (kg/kg/s)      REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)
154    
155      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
156    
157      !IM Amip2 PV a theta constante      !IM Amip2 PV a theta constante
158    
159      CHARACTER(LEN=3) ctetaSTD(nbteta)      CHARACTER(LEN = 3) ctetaSTD(nbteta)
160      DATA ctetaSTD/'350', '380', '405'/      DATA ctetaSTD/'350', '380', '405'/
161      REAL rtetaSTD(nbteta)      REAL rtetaSTD(nbteta)
162      DATA rtetaSTD/350., 380., 405./      DATA rtetaSTD/350., 380., 405./
163    
164      !MI Amip2 PV a theta constante      !MI Amip2 PV a theta constante
165    
166      INTEGER klevp1      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
167      PARAMETER(klevp1=llm+1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
   
     REAL swdn0(klon, klevp1), swdn(klon, klevp1)  
     REAL swup0(klon, klevp1), swup(klon, klevp1)  
168      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
169    
170      REAL SWdn200clr(klon), SWdn200(klon)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
171      REAL SWup200clr(klon), SWup200(klon)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
     SAVE SWdn200clr, SWdn200, SWup200clr, SWup200  
   
     REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)  
     REAL lwup0(klon, klevp1), lwup(klon, klevp1)  
172      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
173    
     REAL LWdn200clr(klon), LWdn200(klon)  
     REAL LWup200clr(klon), LWup200(klon)  
     SAVE LWdn200clr, LWdn200, LWup200clr, LWup200  
   
174      !IM Amip2      !IM Amip2
175      ! variables a une pression donnee      ! variables a une pression donnee
176    
177      integer nlevSTD      integer nlevSTD
178      PARAMETER(nlevSTD=17)      PARAMETER(nlevSTD = 17)
179      real rlevSTD(nlevSTD)      real rlevSTD(nlevSTD)
180      DATA rlevSTD/100000., 92500., 85000., 70000., &      DATA rlevSTD/100000., 92500., 85000., 70000., &
181           60000., 50000., 40000., 30000., 25000., 20000., &           60000., 50000., 40000., 30000., 25000., 20000., &
182           15000., 10000., 7000., 5000., 3000., 2000., 1000./           15000., 10000., 7000., 5000., 3000., 2000., 1000./
183      CHARACTER(LEN=4) clevSTD(nlevSTD)      CHARACTER(LEN = 4) clevSTD(nlevSTD)
184      DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &      DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &
185           '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &           '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &
186           '70  ', '50  ', '30  ', '20  ', '10  '/           '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  
   
     REAL tsumSTD(klon, nlevSTD, nout)  
     REAL usumSTD(klon, nlevSTD, nout), vsumSTD(klon, nlevSTD, nout)  
     REAL wsumSTD(klon, nlevSTD, nout), phisumSTD(klon, nlevSTD, nout)  
     REAL qsumSTD(klon, nlevSTD, nout), rhsumSTD(klon, nlevSTD, nout)  
   
     SAVE tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD,  &  
          qsumSTD, rhsumSTD  
   
     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 uvsumSTD(klon, nlevSTD, nout)  
     real vqsumSTD(klon, nlevSTD, nout)  
     real vTsumSTD(klon, nlevSTD, nout)  
     real wqsumSTD(klon, nlevSTD, nout)  
   
     real vphiSTD(klon, nlevSTD)  
     real wTSTD(klon, nlevSTD)  
     real u2STD(klon, nlevSTD)  
     real v2STD(klon, nlevSTD)  
     real T2STD(klon, nlevSTD)  
   
     real vphisumSTD(klon, nlevSTD, nout)  
     real wTsumSTD(klon, nlevSTD, nout)  
     real u2sumSTD(klon, nlevSTD, nout)  
     real v2sumSTD(klon, nlevSTD, nout)  
     real T2sumSTD(klon, nlevSTD, nout)  
   
     SAVE uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD  
     SAVE vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD  
     !MI Amip2  
187    
188      ! prw: precipitable water      ! prw: precipitable water
189      real prw(klon)      real prw(klon)
# Line 265  contains Line 193  contains
193      REAL flwp(klon), fiwp(klon)      REAL flwp(klon), fiwp(klon)
194      REAL flwc(klon, llm), fiwc(klon, llm)      REAL flwc(klon, llm), fiwc(klon, llm)
195    
196      INTEGER l, kmax, lmax      INTEGER kmax, lmax
197      PARAMETER(kmax=8, lmax=8)      PARAMETER(kmax = 8, lmax = 8)
198      INTEGER kmaxm1, lmaxm1      INTEGER kmaxm1, lmaxm1
199      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)
200    
201      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)      REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
202      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./      DATA zx_tau/0., 0.3, 1.3, 3.6, 9.4, 23., 60./
203      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
204    
205      ! cldtopres pression au sommet des nuages      ! cldtopres pression au sommet des nuages
# Line 279  contains Line 207  contains
207      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
208    
209      ! taulev: numero du niveau de tau dans les sorties ISCCP      ! taulev: numero du niveau de tau dans les sorties ISCCP
210      CHARACTER(LEN=4) taulev(kmaxm1)      CHARACTER(LEN = 4) taulev(kmaxm1)
211    
212      DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/      DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/
213      CHARACTER(LEN=3) pclev(lmaxm1)      CHARACTER(LEN = 3) pclev(lmaxm1)
214      DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/      DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/
215    
216      CHARACTER(LEN=28) cnameisccp(lmaxm1, kmaxm1)      CHARACTER(LEN = 28) cnameisccp(lmaxm1, kmaxm1)
217      DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &      DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &
218           'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &           'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &
219           'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &           'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &
# Line 317  contains Line 245  contains
245      integer nid_hf, nid_hf3d      integer nid_hf, nid_hf3d
246      save nid_hf, nid_hf3d      save nid_hf, nid_hf3d
247    
     INTEGER        longcles  
     PARAMETER    ( longcles = 20 )  
   
248      ! Variables propres a la physique      ! Variables propres a la physique
249    
250      INTEGER, save:: radpas      INTEGER, save:: radpas
# Line 327  contains Line 252  contains
252      ! "physiq".)      ! "physiq".)
253    
254      REAL radsol(klon)      REAL radsol(klon)
255      SAVE radsol               ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
256    
257      INTEGER, SAVE:: itap ! number of calls to "physiq"      INTEGER, SAVE:: itap ! number of calls to "physiq"
258    
259      REAL ftsol(klon, nbsrf)      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
     SAVE ftsol                  ! temperature du sol  
260    
261      REAL ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
262      SAVE ftsoil                 ! temperature dans le sol      ! soil temperature of surface fraction
263    
264      REAL fevap(klon, nbsrf)      REAL, save:: fevap(klon, nbsrf) ! evaporation
     SAVE fevap                 ! evaporation  
265      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
266      SAVE fluxlat      SAVE fluxlat
267    
268      REAL fqsurf(klon, nbsrf)      REAL fqsurf(klon, nbsrf)
269      SAVE fqsurf                 ! humidite de l'air au contact de la surface      SAVE fqsurf ! humidite de l'air au contact de la surface
270    
271      REAL qsol(klon)      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol
     SAVE qsol                  ! hauteur d'eau dans le sol  
272    
273      REAL fsnow(klon, nbsrf)      REAL fsnow(klon, nbsrf)
274      SAVE fsnow                  ! epaisseur neigeuse      SAVE fsnow ! epaisseur neigeuse
275    
276      REAL falbe(klon, nbsrf)      REAL falbe(klon, nbsrf)
277      SAVE falbe                  ! albedo par type de surface      SAVE falbe ! albedo par type de surface
278      REAL falblw(klon, nbsrf)      REAL falblw(klon, nbsrf)
279      SAVE falblw                 ! albedo par type de surface      SAVE falblw ! albedo par type de surface
280    
281      ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
282      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
283      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
284      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 371  contains Line 293  contains
293      INTEGER igwd, idx(klon), itest(klon)      INTEGER igwd, idx(klon), itest(klon)
294    
295      REAL agesno(klon, nbsrf)      REAL agesno(klon, nbsrf)
296      SAVE agesno                 ! age de la neige      SAVE agesno ! age de la neige
297    
298      REAL run_off_lic_0(klon)      REAL run_off_lic_0(klon)
299      SAVE run_off_lic_0      SAVE run_off_lic_0
300      !KE43      !KE43
301      ! Variables liees a la convection de K. Emanuel (sb):      ! Variables liees a la convection de K. Emanuel (sb):
302    
303      REAL bas, top             ! cloud base and top levels      REAL bas, top ! cloud base and top levels
304      SAVE bas      SAVE bas
305      SAVE top      SAVE top
306    
307      REAL Ma(klon, llm)        ! undilute upward mass flux      REAL Ma(klon, llm) ! undilute upward mass flux
308      SAVE Ma      SAVE Ma
309      REAL qcondc(klon, llm)    ! in-cld water content from convect      REAL qcondc(klon, llm) ! in-cld water content from convect
310      SAVE qcondc      SAVE qcondc
311      REAL ema_work1(klon, llm), ema_work2(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
312      SAVE ema_work1, ema_work2      REAL, save:: wd(klon)
   
     REAL wd(klon) ! sb  
     SAVE wd       ! sb  
313    
314      ! Variables locales pour la couche limite (al1):      ! Variables locales pour la couche limite (al1):
315    
# Line 399  contains Line 318  contains
318      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
319      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
320    
321      !AA  Pour phytrac      ! Pour phytrac :
322      REAL ycoefh(klon, llm)    ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
323      REAL yu1(klon)            ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
324      REAL yv1(klon)            ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
325      REAL ffonte(klon, nbsrf)    !Flux thermique utilise pour fondre la neige      REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige
326      REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface      REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface
327      !                               !et necessaire pour limiter la      ! !et necessaire pour limiter la
328      !                               !hauteur de neige, en kg/m2/s      ! !hauteur de neige, en kg/m2/s
329      REAL zxffonte(klon), zxfqcalving(klon)      REAL zxffonte(klon), zxfqcalving(klon)
330    
331      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
# Line 418  contains Line 337  contains
337      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
338      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
339    
340      !AA      REAL, save:: rain_fall(klon) ! pluie
341      REAL rain_fall(klon) ! pluie      REAL, save:: snow_fall(klon) ! neige
     REAL snow_fall(klon) ! neige  
     save snow_fall, rain_fall  
     !IM cf FH pour Tiedtke 080604  
     REAL rain_tiedtke(klon), snow_tiedtke(klon)  
342    
343      REAL total_rain(klon), nday_rain(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
     save nday_rain  
344    
345      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon), devap(klon) ! evaporation and its derivative
346      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
347      REAL dlw(klon)    ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
348      SAVE dlw      SAVE dlw
349      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
350      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL fder(klon) ! Derive de flux (sensible et latente)
# Line 449  contains Line 363  contains
363      INTEGER julien      INTEGER julien
364    
365      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
366      REAL pctsrf(klon, nbsrf)      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
367      !IM      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE  
368    
     SAVE pctsrf                 ! sous-fraction du sol  
369      REAL albsol(klon)      REAL albsol(klon)
370      SAVE albsol                 ! albedo du sol total      SAVE albsol ! albedo du sol total
371      REAL albsollw(klon)      REAL albsollw(klon)
372      SAVE albsollw                 ! albedo du sol total      SAVE albsollw ! albedo du sol total
373    
374      REAL, SAVE:: wo(klon, llm) ! ozone      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
375    
376      ! Declaration des procedures appelees      ! Declaration des procedures appelees
377    
378      EXTERNAL alboc     ! calculer l'albedo sur ocean      EXTERNAL alboc ! calculer l'albedo sur ocean
     EXTERNAL ajsec     ! ajustement sec  
     EXTERNAL clmain    ! couche limite  
379      !KE43      !KE43
380      EXTERNAL conema3  ! convect4.3      EXTERNAL conema3 ! convect4.3
381      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)      EXTERNAL nuage ! calculer les proprietes radiatives
382      EXTERNAL nuage     ! calculer les proprietes radiatives      EXTERNAL transp ! transport total de l'eau et de l'energie
     EXTERNAL ozonecm   ! prescrire l'ozone  
     EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique  
     EXTERNAL radlwsw   ! rayonnements solaire et infrarouge  
     EXTERNAL transp    ! transport total de l'eau et de l'energie  
   
     EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression  
   
     EXTERNAL undefSTD  
     ! (somme les valeurs definies d'1 var a 1 niveau de pression)  
383    
384      ! Variables locales      ! Variables locales
385    
386      real clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
387      real clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
388    
389      save rnebcon, clwcon      REAL rhcl(klon, llm) ! humiditi relative ciel clair
390        REAL dialiq(klon, llm) ! eau liquide nuageuse
391      REAL rhcl(klon, llm)    ! humiditi relative ciel clair      REAL diafra(klon, llm) ! fraction nuageuse
392      REAL dialiq(klon, llm)  ! eau liquide nuageuse      REAL cldliq(klon, llm) ! eau liquide nuageuse
393      REAL diafra(klon, llm)  ! fraction nuageuse      REAL cldfra(klon, llm) ! fraction nuageuse
394      REAL cldliq(klon, llm)  ! eau liquide nuageuse      REAL cldtau(klon, llm) ! epaisseur optique
395      REAL cldfra(klon, llm)  ! fraction nuageuse      REAL cldemi(klon, llm) ! emissivite infrarouge
396      REAL cldtau(klon, llm)  ! epaisseur optique  
397      REAL cldemi(klon, llm)  ! emissivite infrarouge      REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
398        REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
399      REAL fluxq(klon, llm, nbsrf)   ! flux turbulent d'humidite      REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
400      REAL fluxt(klon, llm, nbsrf)   ! flux turbulent de chaleur      REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
     REAL fluxu(klon, llm, nbsrf)   ! flux turbulent de vitesse u  
     REAL fluxv(klon, llm, nbsrf)   ! flux turbulent de vitesse v  
401    
402      REAL zxfluxt(klon, llm)      REAL zxfluxt(klon, llm)
403      REAL zxfluxq(klon, llm)      REAL zxfluxq(klon, llm)
404      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
405      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
406    
407      REAL heat(klon, llm)    ! chauffage solaire      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
408      REAL heat0(klon, llm)   ! chauffage solaire ciel clair      ! les variables soient r\'emanentes.
409      REAL cool(klon, llm)    ! refroidissement infrarouge      REAL, save:: heat(klon, llm) ! chauffage solaire
410      REAL cool0(klon, llm)   ! refroidissement infrarouge ciel clair      REAL heat0(klon, llm) ! chauffage solaire ciel clair
411      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
412      real sollwdown(klon)    ! downward LW flux at surface      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
413      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
414        REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
415        real, save:: sollwdown(klon) ! downward LW flux at surface
416        REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
417      REAL albpla(klon)      REAL albpla(klon)
418      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
419      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
420      ! Le rayonnement n'est pas calcule tous les pas, il faut donc      SAVE albpla
421      !                      sauvegarder les sorties du rayonnement      SAVE heat0, cool0
     SAVE  heat, cool, albpla, topsw, toplw, solsw, sollw, sollwdown  
     SAVE  topsw0, toplw0, solsw0, sollw0, heat0, cool0  
422    
423      INTEGER itaprad      INTEGER itaprad
424      SAVE itaprad      SAVE itaprad
425    
426      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
427      REAL conv_t(klon, llm) ! convergence de la temperature(K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K/s)
428    
429      REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut      REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut
430      REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree      REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree
# Line 534  contains Line 434  contains
434      REAL dist, rmu0(klon), fract(klon)      REAL dist, rmu0(klon), fract(klon)
435      REAL zdtime ! pas de temps du rayonnement (s)      REAL zdtime ! pas de temps du rayonnement (s)
436      real zlongi      real zlongi
   
437      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
438      REAL za, zb      REAL za, zb
439      REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp      REAL zx_t, zx_qs, zdelta, zcor
440      real zqsat(klon, llm)      real zqsat(klon, llm)
441      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
442      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup=234.0)  
   
443      REAL zphi(klon, llm)      REAL zphi(klon, llm)
444    
445      !IM cf. AM Variables locales pour la CLA (hbtm2)      !IM cf. AM Variables locales pour la CLA (hbtm2)
446    
447      REAL pblh(klon, nbsrf)           ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
448      REAL plcl(klon, nbsrf)           ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
449      REAL capCL(klon, nbsrf)          ! CAPE de couche limite      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
450      REAL oliqCL(klon, nbsrf)          ! eau_liqu integree de couche limite      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
451      REAL cteiCL(klon, nbsrf)          ! cloud top instab. crit. couche limite      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
452      REAL pblt(klon, nbsrf)          ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
453      REAL therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
454      REAL trmb1(klon, nbsrf)          ! deep_cape      REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
455      REAL trmb2(klon, nbsrf)          ! inhibition      REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
456      REAL trmb3(klon, nbsrf)          ! Point Omega      REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
457      ! Grdeurs de sorties      ! Grdeurs de sorties
458      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
459      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
460      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
461      REAL s_trmb3(klon)      REAL s_trmb3(klon)
462    
463      ! Variables locales pour la convection de K. Emanuel (sb):      ! Variables locales pour la convection de K. Emanuel :
464    
465      REAL upwd(klon, llm)      ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
466      REAL dnwd(klon, llm)      ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
467      REAL dnwd0(klon, llm)     ! unsaturated downdraft mass flux      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
468      REAL tvp(klon, llm)       ! virtual temp of lifted parcel      REAL tvp(klon, llm) ! virtual temp of lifted parcel
469      REAL cape(klon)           ! CAPE      REAL cape(klon) ! CAPE
470      SAVE cape      SAVE cape
471    
472      REAL pbase(klon)          ! cloud base pressure      REAL pbase(klon) ! cloud base pressure
473      SAVE pbase      SAVE pbase
474      REAL bbase(klon)          ! cloud base buoyancy      REAL bbase(klon) ! cloud base buoyancy
475      SAVE bbase      SAVE bbase
476      REAL rflag(klon)          ! flag fonctionnement de convect      REAL rflag(klon) ! flag fonctionnement de convect
477      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
478      ! -- convect43:      ! -- convect43:
     INTEGER ntra              ! nb traceurs pour convect4.3  
479      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
480      REAL dplcldt(klon), dplcldr(klon)      REAL dplcldt(klon), dplcldr(klon)
481    
482      ! Variables du changement      ! Variables du changement
483    
484      ! con: convection      ! con: convection
485      ! lsc: condensation a grande echelle (Large-Scale-Condensation)      ! lsc: large scale condensation
486      ! ajs: ajustement sec      ! ajs: ajustement sec
487      ! eva: evaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
488      ! vdf: couche limite (Vertical DiFfusion)      ! vdf: vertical diffusion in boundary layer
489      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
490      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
491      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 599  contains Line 493  contains
493      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
494      REAL rneb(klon, llm)      REAL rneb(klon, llm)
495    
496      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
497      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
498      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
499      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
500      REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1)      REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
501      REAL prfl(klon, llm+1), psfl(klon, llm+1)      REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
502    
503      INTEGER ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
   
     SAVE ibas_con, itop_con  
504    
505      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon), rain_lsc(klon)
506      REAL snow_con(klon), snow_lsc(klon)      REAL snow_con(klon), snow_lsc(klon)
# Line 622  contains Line 514  contains
514      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)      REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
515      REAL d_t_lif(klon, llm)      REAL d_t_lif(klon, llm)
516    
517      REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)      REAL, save:: ratqs(klon, llm)
518      real ratqsbas, ratqshaut      real ratqss(klon, llm), ratqsc(klon, llm)
519      save ratqsbas, ratqshaut, ratqs      real:: ratqsbas = 0.01, ratqshaut = 0.3
520    
521      ! Parametres lies au nouveau schema de nuages (SB, PDF)      ! Parametres lies au nouveau schema de nuages (SB, PDF)
522      real fact_cldcon      real:: fact_cldcon = 0.375
523      real facttemps      real:: facttemps = 1.e-4
524      logical ok_newmicro      logical:: ok_newmicro = .true.
     save ok_newmicro  
     save fact_cldcon, facttemps  
525      real facteur      real facteur
526    
527      integer iflag_cldcon      integer:: iflag_cldcon = 1
     save iflag_cldcon  
   
528      logical ptconv(klon, llm)      logical ptconv(klon, llm)
529    
530      ! Variables liees a l'ecriture de la bande histoire physique      ! Variables locales pour effectuer les appels en s\'erie :
   
     integer itau_w   ! pas de temps ecriture = itap + itau_phy  
   
     ! Variables locales pour effectuer les appels en serie  
531    
532      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
533      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm), qs_seri(klon, llm)
# Line 653  contains Line 537  contains
537      REAL d_tr(klon, llm, nbtr)      REAL d_tr(klon, llm, nbtr)
538    
539      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
     INTEGER ndex2d(iim*(jjm + 1)), ndex3d(iim*(jjm + 1)*llm)  
540    
541      REAL zustrdr(klon), zvstrdr(klon)      REAL zustrdr(klon), zvstrdr(klon)
542      REAL zustrli(klon), zvstrli(klon)      REAL zustrli(klon), zvstrli(klon)
543      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
544      REAL aam, torsfc      REAL aam, torsfc
545    
546      REAL dudyn(iim+1, jjm + 1, llm)      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
   
     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)  
547    
548      INTEGER nid_day, nid_ins      INTEGER, SAVE:: nid_day, nid_ins
     SAVE nid_day, nid_ins  
549    
550      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.
551      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.
# Line 677  contains Line 554  contains
554    
555      REAL zsto      REAL zsto
556    
     character(len=20) modname  
     character(len=80) abort_message  
557      logical ok_sync      logical ok_sync
558      real date0      real date0
559    
560      !     Variables liees au bilan d'energie et d'enthalpi      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
561      REAL ztsol(klon)      REAL ztsol(klon)
562      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
563      REAL      d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
564      REAL      fs_bound, fq_bound      REAL fs_bound, fq_bound
565      SAVE      d_h_vcol_phy      REAL zero_v(klon)
566      REAL      zero_v(klon)      CHARACTER(LEN = 15) tit
567      CHARACTER(LEN=15) ztit      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
568      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
569      SAVE      ip_ebil  
570      DATA      ip_ebil/0/      REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
     INTEGER   if_ebil ! level for energy conserv. dignostics  
     SAVE      if_ebil  
     !+jld ec_conser  
     REAL d_t_ec(klon, llm)    ! tendance du a la conersion Ec -> E thermique  
571      REAL ZRCPD      REAL ZRCPD
572      !-jld ec_conser  
573      !IM: t2m, q2m, u10m, v10m      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
574      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)   !temperature, humidite a 2m      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
575      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m      REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
576      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
577      REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille  
578      !jq   Aerosol effects (Johannes Quaas, 27/11/2003)      ! Aerosol effects:
579      REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]  
580        REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
581      REAL sulfate_pi(klon, llm)  
582      ! (SO4 aerosol concentration [ug/m3] (pre-industrial value))      REAL, save:: sulfate_pi(klon, llm)
583      SAVE sulfate_pi      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
584    
585      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
586      ! (Cloud optical thickness for pre-industrial (pi) aerosols)      ! cloud optical thickness for pre-industrial (pi) aerosols
587    
588      REAL re(klon, llm)       ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
589      REAL fl(klon, llm)  ! denominator of re      REAL fl(klon, llm) ! denominator of re
590    
591      ! Aerosol optical properties      ! Aerosol optical properties
592      REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)      REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
593      REAL cg_ae(klon, llm, 2)      REAL, save:: cg_ae(klon, llm, 2)
594    
595      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.      REAL topswad(klon), solswad(klon) ! aerosol direct effect
596      ! ok_ade=T -ADE=topswad-topsw      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
597    
598      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.      REAL aerindex(klon) ! POLDER aerosol index
     ! ok_aie=T ->  
     !        ok_ade=T -AIE=topswai-topswad  
     !        ok_ade=F -AIE=topswai-topsw  
599    
600      REAL aerindex(klon)       ! POLDER aerosol index      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
601        LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
602    
603      ! Parameters      REAL:: bl95_b0 = 2., bl95_b1 = 0.2
604      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not      ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
605      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)      ! B). They link cloud droplet number concentration to aerosol mass
606        ! concentration.
607    
     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1  
608      SAVE u10m      SAVE u10m
609      SAVE v10m      SAVE v10m
610      SAVE t2m      SAVE t2m
611      SAVE q2m      SAVE q2m
612      SAVE ffonte      SAVE ffonte
613      SAVE fqcalving      SAVE fqcalving
     SAVE piz_ae  
     SAVE tau_ae  
     SAVE cg_ae  
614      SAVE rain_con      SAVE rain_con
615      SAVE snow_con      SAVE snow_con
616      SAVE topswai      SAVE topswai
# Line 753  contains Line 619  contains
619      SAVE solswad      SAVE solswad
620      SAVE d_u_con      SAVE d_u_con
621      SAVE d_v_con      SAVE d_v_con
622      SAVE rnebcon0  
623      SAVE clwcon0      real zmasse(klon, llm)
624      SAVE pblh      ! (column-density of mass of air in a cell, in kg m-2)
625      SAVE plcl  
626      SAVE capCL      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
627      SAVE oliqCL  
628      SAVE cteiCL      namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
629      SAVE pblt           fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &
630      SAVE therm           ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &
631      SAVE trmb1           nsplit_thermals
     SAVE trmb2  
     SAVE trmb3  
632    
633      !----------------------------------------------------------------      !----------------------------------------------------------------
634    
635      modname = 'physiq'      IF (if_ebil >= 1) zero_v = 0.
636      IF (if_ebil >= 1) THEN      ok_sync = .TRUE.
637         DO i=1, klon      IF (nqmx < 2) CALL abort_gcm('physiq', &
638            zero_v(i)=0.           'eaux vapeur et liquide sont indispensables', 1)
        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  
639    
640      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
641         !  initialiser         ! initialiser
642         u10m=0.         u10m = 0.
643         v10m=0.         v10m = 0.
644         t2m=0.         t2m = 0.
645         q2m=0.         q2m = 0.
646         ffonte=0.         ffonte = 0.
647         fqcalving=0.         fqcalving = 0.
648         piz_ae(:, :, :)=0.         piz_ae = 0.
649         tau_ae(:, :, :)=0.         tau_ae = 0.
650         cg_ae(:, :, :)=0.         cg_ae = 0.
651         rain_con(:)=0.         rain_con(:) = 0.
652         snow_con(:)=0.         snow_con(:) = 0.
653         bl95_b0=0.         topswai(:) = 0.
654         bl95_b1=0.         topswad(:) = 0.
655         topswai(:)=0.         solswai(:) = 0.
656         topswad(:)=0.         solswad(:) = 0.
657         solswai(:)=0.  
658         solswad(:)=0.         d_u_con = 0.
659           d_v_con = 0.
660         d_u_con = 0.0         rnebcon0 = 0.
661         d_v_con = 0.0         clwcon0 = 0.
662         rnebcon0 = 0.0         rnebcon = 0.
663         clwcon0 = 0.0         clwcon = 0.
664         rnebcon = 0.0  
665         clwcon = 0.0         pblh =0. ! Hauteur de couche limite
666           plcl =0. ! Niveau de condensation de la CLA
667         pblh   =0.        ! Hauteur de couche limite         capCL =0. ! CAPE de couche limite
668         plcl   =0.        ! Niveau de condensation de la CLA         oliqCL =0. ! eau_liqu integree de couche limite
669         capCL  =0.        ! CAPE de couche limite         cteiCL =0. ! cloud top instab. crit. couche limite
670         oliqCL =0.        ! eau_liqu integree de couche limite         pblt =0. ! T a la Hauteur de couche limite
671         cteiCL =0.        ! cloud top instab. crit. couche limite         therm =0.
672         pblt   =0.        ! T a la Hauteur de couche limite         trmb1 =0. ! deep_cape
673         therm  =0.         trmb2 =0. ! inhibition
674         trmb1  =0.        ! deep_cape         trmb3 =0. ! Point Omega
675         trmb2  =0.        ! inhibition  
676         trmb3  =0.        ! Point Omega         IF (if_ebil >= 1) d_h_vcol_phy = 0.
677    
678         IF (if_ebil >= 1) d_h_vcol_phy=0.         iflag_thermals = 0
679           nsplit_thermals = 1
680         ! appel a la lecture du run.def physique         print *, "Enter namelist 'physiq_nml'."
681           read(unit=*, nml=physiq_nml)
682         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &         write(unit_nml, nml=physiq_nml)
683              ok_instan, fact_cldcon, facttemps, ok_newmicro, &  
684              iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &         call conf_phys
             ok_ade, ok_aie,  &  
             bl95_b0, bl95_b1, &  
             iflag_thermals, nsplit_thermals)  
685    
686         ! Initialiser les compteurs:         ! Initialiser les compteurs:
687    
# Line 835  contains Line 689  contains
689         itap = 0         itap = 0
690         itaprad = 0         itaprad = 0
691         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &
692              seaice, fqsurf, qsol, fsnow, &              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &
693              falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollwdown, &              snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &
694              dlw, radsol, frugs, agesno, &              zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
695              zmea, zstd, zsig, zgam, zthe, zpic, zval, &              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
             t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon,  &  
             run_off_lic_0)  
696    
697         !   ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
698         q2(:, :, :)=1.e-8         q2 = 1e-8
699    
700         radpas = NINT( 86400. / pdtphys / nbapp_rad)         radpas = NINT(86400. / dtphys / nbapp_rad)
701    
702         ! on remet le calendrier a zero         ! on remet le calendrier a zero
703           IF (raz_date) itau_phy = 0
        IF (raz_date == 1) THEN  
           itau_phy = 0  
        ENDIF  
704    
705         PRINT *, 'cycle_diurne = ', cycle_diurne         PRINT *, 'cycle_diurne = ', cycle_diurne
706           CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &
707                ok_instan, ok_region)
708    
709         IF(ocean.NE.'force ') THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
710            ok_ocean=.TRUE.            print *, "Au minimum 4 appels par jour si cycle diurne"
711              call abort_gcm('physiq', &
712                   "Nombre d'appels au rayonnement insuffisant", 1)
713         ENDIF         ENDIF
714    
715         CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &         ! Initialisation pour le sch\'ema de convection d'Emanuel :
             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):  
716         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
717              ibas_con = 1
718            print *,"*** Convection de Kerry Emanuel 4.3  "            itop_con = 1
   
           !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  
   
719         ENDIF         ENDIF
720    
721         IF (ok_orodr) THEN         IF (ok_orodr) THEN
722            rugoro = MAX(1e-5, zstd * zsig / 2)            rugoro = MAX(1e-5, zstd * zsig / 2)
723            CALL SUGWD(klon, llm, paprs, pplay)            CALL SUGWD(paprs, play)
724         else         else
725            rugoro = 0.            rugoro = 0.
726         ENDIF         ENDIF
727    
728         lmt_pas = NINT(86400. / pdtphys)  ! tous les jours         lmt_pas = NINT(86400. / dtphys) ! tous les jours
729         print *, 'Number of time steps of "physics" per day: ', lmt_pas         print *, 'Number of time steps of "physics" per day: ', lmt_pas
730    
731         ecrit_ins = NINT(ecrit_ins/pdtphys)         ecrit_ins = NINT(ecrit_ins/dtphys)
732         ecrit_hf = NINT(ecrit_hf/pdtphys)         ecrit_hf = NINT(ecrit_hf/dtphys)
733         ecrit_day = NINT(ecrit_day/pdtphys)         ecrit_mth = NINT(ecrit_mth/dtphys)
734         ecrit_mth = NINT(ecrit_mth/pdtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
735         ecrit_tra = NINT(86400.*ecrit_tra/pdtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
        ecrit_reg = NINT(ecrit_reg/pdtphys)  
736    
737         ! Initialiser le couplage si necessaire         ! Initialiser le couplage si necessaire
738    
739         npas = 0         npas = 0
740         nexca = 0         nexca = 0
741    
742         print *,'AVANT HIST IFLAG_CON=', iflag_con         ! Initialisation des sorties
   
        !   Initialisation des sorties  
743    
744         call ini_histhf(pdtphys, presnivs, nid_hf, nid_hf3d)         call ini_histins(dtphys, ok_instan, nid_ins)
        call ini_histday(pdtphys, presnivs, ok_journe, nid_day)  
        call ini_histins(pdtphys, presnivs, ok_instan, nid_ins)  
745         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
746         !XXXPB Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
747         WRITE(*, *) 'physiq date0 : ', date0         print *, 'physiq date0: ', date0
748      ENDIF test_firstcal      ENDIF test_firstcal
749    
750      ! Mettre a zero des variables de sortie (pour securite)      ! Mettre a zero des variables de sortie (pour securite)
751    
752      DO i = 1, klon      DO i = 1, klon
753         d_ps(i) = 0.0         d_ps(i) = 0.
754      ENDDO      ENDDO
755      DO k = 1, llm      DO iq = 1, nqmx
        DO i = 1, klon  
           d_t(i, k) = 0.0  
           d_u(i, k) = 0.0  
           d_v(i, k) = 0.0  
        ENDDO  
     ENDDO  
     DO iq = 1, nq  
756         DO k = 1, llm         DO k = 1, llm
757            DO i = 1, klon            DO i = 1, klon
758               d_qx(i, k, iq) = 0.0               d_qx(i, k, iq) = 0.
759            ENDDO            ENDDO
760         ENDDO         ENDDO
761      ENDDO      ENDDO
762      da=0.      da = 0.
763      mp=0.      mp = 0.
764      phi(:, :, :)=0.      phi = 0.
765    
766      ! Ne pas affecter les valeurs entrees de u, v, h, et q      ! Ne pas affecter les valeurs entr\'ees de u, v, h, et q :
767    
768      DO k = 1, llm      DO k = 1, llm
769         DO i = 1, klon         DO i = 1, klon
770            t_seri(i, k)  = t(i, k)            t_seri(i, k) = t(i, k)
771            u_seri(i, k)  = u(i, k)            u_seri(i, k) = u(i, k)
772            v_seri(i, k)  = v(i, k)            v_seri(i, k) = v(i, k)
773            q_seri(i, k)  = qx(i, k, ivap)            q_seri(i, k) = qx(i, k, ivap)
774            ql_seri(i, k) = qx(i, k, iliq)            ql_seri(i, k) = qx(i, k, iliq)
775            qs_seri(i, k) = 0.            qs_seri(i, k) = 0.
776         ENDDO         ENDDO
777      ENDDO      ENDDO
778      IF (nq >= 3) THEN      IF (nqmx >= 3) THEN
779         tr_seri(:, :, :nq-2) = qx(:, :, 3:nq)         tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)
780      ELSE      ELSE
781         tr_seri(:, :, 1) = 0.         tr_seri(:, :, 1) = 0.
782      ENDIF      ENDIF
# Line 971  contains Line 791  contains
791      ENDDO      ENDDO
792    
793      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
794         ztit='after dynamic'         tit = 'after dynamics'
795         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
796              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
797              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
798         !     Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajout\'es dans la
799         !     on devrait avoir que la variation d'entalpie par la dynamique         !  dynamique, la variation d'enthalpie par la dynamique devrait
800         !     est egale a la variation de la physique au pas de temps precedent.         !  \^etre \'egale \`a la variation de la physique au pas de temps
801         !     Donc la somme de ces 2 variations devrait etre nulle.         !  pr\'ec\'edent.  Donc la somme de ces 2 variations devrait \^etre
802         call diagphy(airephy, ztit, ip_ebil &         !  nulle.
803              , zero_v, zero_v, zero_v, zero_v, zero_v &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
804              , zero_v, zero_v, zero_v, ztsol &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
805              , d_h_vcol+d_h_vcol_phy, d_qt, 0. &              d_qt, 0., fs_bound, fq_bound)
             , fs_bound, fq_bound )  
806      END IF      END IF
807    
808      ! Diagnostiquer la tendance dynamique      ! Diagnostic de la tendance dynamique :
   
809      IF (ancien_ok) THEN      IF (ancien_ok) THEN
810         DO k = 1, llm         DO k = 1, llm
811            DO i = 1, klon            DO i = 1, klon
812               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
813               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
814            ENDDO            ENDDO
815         ENDDO         ENDDO
816      ELSE      ELSE
817         DO k = 1, llm         DO k = 1, llm
818            DO i = 1, klon            DO i = 1, klon
819               d_t_dyn(i, k) = 0.0               d_t_dyn(i, k) = 0.
820               d_q_dyn(i, k) = 0.0               d_q_dyn(i, k) = 0.
821            ENDDO            ENDDO
822         ENDDO         ENDDO
823         ancien_ok = .TRUE.         ancien_ok = .TRUE.
824      ENDIF      ENDIF
825    
826      ! Ajouter le geopotentiel du sol:      ! Ajouter le geopotentiel du sol:
   
827      DO k = 1, llm      DO k = 1, llm
828         DO i = 1, klon         DO i = 1, klon
829            zphi(i, k) = pphi(i, k) + pphis(i)            zphi(i, k) = pphi(i, k) + pphis(i)
830         ENDDO         ENDDO
831      ENDDO      ENDDO
832    
833      ! Verifier les temperatures      ! Check temperatures:
   
834      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
835    
836      ! Incrementer le compteur de la physique      ! Incrementer le compteur de la physique
   
837      itap = itap + 1      itap = itap + 1
838      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(NINT(rdayvrai), 360)
839      if (julien == 0) julien = 360      if (julien == 0) julien = 360
840    
841      ! Mettre en action les conditions aux limites (albedo, sst, etc.).      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.  
842    
843      IF (MOD(itap - 1, lmt_pas) == 0) THEN      ! Mettre en action les conditions aux limites (albedo, sst etc.).
        CALL ozonecm(REAL(julien), rlat, paprs, wo)  
     ENDIF  
844    
845      ! Re-evaporer l'eau liquide nuageuse      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
846        wo = ozonecm(REAL(julien), paprs)
847    
848      DO k = 1, llm  ! re-evaporation de l'eau liquide nuageuse      ! \'Evaporation de l'eau liquide nuageuse :
849        DO k = 1, llm
850         DO i = 1, klon         DO i = 1, klon
851            zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))            zb = MAX(0., ql_seri(i, k))
852            zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))            t_seri(i, k) = t_seri(i, k) &
853            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  
854            q_seri(i, k) = q_seri(i, k) + zb            q_seri(i, k) = q_seri(i, k) + zb
           ql_seri(i, k) = 0.0  
855         ENDDO         ENDDO
856      ENDDO      ENDDO
857        ql_seri = 0.
858    
859      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
860         ztit='after reevap'         tit = 'after reevap'
861         CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
862              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
863              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
864         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
865              , zero_v, zero_v, zero_v, zero_v, zero_v &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
866              , zero_v, zero_v, zero_v, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
867    
868      END IF      END IF
869    
870      ! Appeler la diffusion verticale (programme de couche limite)      ! Appeler la diffusion verticale (programme de couche limite)
871    
872      DO i = 1, klon      DO i = 1, klon
873         zxrugs(i) = 0.0         zxrugs(i) = 0.
874      ENDDO      ENDDO
875      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
876         DO i = 1, klon         DO i = 1, klon
# Line 1079  contains Line 887  contains
887    
888      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), zlongi, dist)
889      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
890         zdtime = pdtphys * REAL(radpas)         zdtime = dtphys * REAL(radpas)
891         CALL zenang(zlongi, gmtime, zdtime, rmu0, fract)         CALL zenang(zlongi, time, zdtime, rmu0, fract)
892      ELSE      ELSE
893         rmu0 = -999.999         rmu0 = -999.999
894      ENDIF      ENDIF
895    
896      !     Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
897      albsol(:)=0.      albsol(:) = 0.
898      albsollw(:)=0.      albsollw(:) = 0.
899      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
900         DO i = 1, klon         DO i = 1, klon
901            albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)            albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)
# Line 1095  contains Line 903  contains
903         ENDDO         ENDDO
904      ENDDO      ENDDO
905    
906      !     Repartition sous maille des flux LW et SW      ! R\'epartition sous maille des flux longwave et shortwave
907      ! Repartition du longwave par sous-surface linearisee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
908    
909      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
910         DO i = 1, klon         DO i = 1, klon
911            fsollw(i, nsrf) = sollw(i) &            fsollw(i, nsrf) = sollw(i) &
912                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))
913            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))
914         ENDDO         ENDDO
915      ENDDO      ENDDO
916    
917      fder = dlw      fder = dlw
918    
919      CALL clmain(pdtphys, itap, date0, pctsrf, pctsrf_new, &      ! Couche limite:
920           t_seri, q_seri, u_seri, v_seri, &  
921           julien, rmu0, co2_ppm,  &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &
922           ok_veget, ocean, npas, nexca, ftsol, &           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &
923           soil_model, cdmmax, cdhmax, &           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
924           ksta, ksta_ter, ok_kzmin, ftsoil, qsol,  &           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &
925           paprs, pplay, fsnow, fqsurf, fevap, falbe, falblw, &           rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &
926           fluxlat, rain_fall, snow_fall, &           frugs, firstcal, agesno, rugoro, d_t_vdf, &
927           fsolsw, fsollw, sollwdown, fder, &           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &
928           rlon, rlat, cuphy, cvphy, frugs, &           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &
929           firstcal, lafin, agesno, rugoro, &           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
930           d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)
931           fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &  
932           q2, dsens, devap, &      ! Incr\'ementation des flux
933           ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &  
934           pblh, capCL, oliqCL, cteiCL, pblT, &      zxfluxt = 0.
935           therm, trmb1, trmb2, trmb3, plcl, &      zxfluxq = 0.
936           fqcalving, ffonte, run_off_lic_0, &      zxfluxu = 0.
937           fluxo, fluxg, tslab, seaice)      zxfluxv = 0.
   
     !XXX Incrementation des flux  
   
     zxfluxt=0.  
     zxfluxq=0.  
     zxfluxu=0.  
     zxfluxv=0.  
938      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
939         DO k = 1, llm         DO k = 1, llm
940            DO i = 1, klon            DO i = 1, klon
941               zxfluxt(i, k) = zxfluxt(i, k) +  &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
942                    fluxt(i, k, nsrf) * pctsrf( i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
943               zxfluxq(i, k) = zxfluxq(i, k) +  &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
944                    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)  
945            END DO            END DO
946         END DO         END DO
947      END DO      END DO
948      DO i = 1, klon      DO i = 1, klon
949         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
950         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
951         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
952      ENDDO      ENDDO
953    
# Line 1164  contains Line 961  contains
961      ENDDO      ENDDO
962    
963      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
964         ztit='after clmain'         tit = 'after clmain'
965         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
966              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
967              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
968         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
969              , zero_v, zero_v, zero_v, zero_v, sens &              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
970              , evap, zero_v, zero_v, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
971      END IF      END IF
972    
973      ! Incrementer la temperature du sol      ! Update surface temperature:
974    
975      DO i = 1, klon      DO i = 1, klon
976         zxtsol(i) = 0.0         zxtsol(i) = 0.
977         zxfluxlat(i) = 0.0         zxfluxlat(i) = 0.
978    
979         zt2m(i) = 0.0         zt2m(i) = 0.
980         zq2m(i) = 0.0         zq2m(i) = 0.
981         zu10m(i) = 0.0         zu10m(i) = 0.
982         zv10m(i) = 0.0         zv10m(i) = 0.
983         zxffonte(i) = 0.0         zxffonte(i) = 0.
984         zxfqcalving(i) = 0.0         zxfqcalving(i) = 0.
985    
986         s_pblh(i) = 0.0         s_pblh(i) = 0.
987         s_lcl(i) = 0.0         s_lcl(i) = 0.
988         s_capCL(i) = 0.0         s_capCL(i) = 0.
989         s_oliqCL(i) = 0.0         s_oliqCL(i) = 0.
990         s_cteiCL(i) = 0.0         s_cteiCL(i) = 0.
991         s_pblT(i) = 0.0         s_pblT(i) = 0.
992         s_therm(i) = 0.0         s_therm(i) = 0.
993         s_trmb1(i) = 0.0         s_trmb1(i) = 0.
994         s_trmb2(i) = 0.0         s_trmb2(i) = 0.
995         s_trmb3(i) = 0.0         s_trmb3(i) = 0.
996    
997         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +  &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
998              pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)  &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
999              THEN              'physiq : probl\`eme sous surface au point ', i, pctsrf(i, 1 : nbsrf)
           WRITE(*, *) 'physiq : pb sous surface au point ', i,  &  
                pctsrf(i, 1 : nbsrf)  
        ENDIF  
1000      ENDDO      ENDDO
1001      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1002         DO i = 1, klon         DO i = 1, klon
# Line 1217  contains Line 1009  contains
1009            zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)            zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
1010            zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)            zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
1011            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
1012            zxfqcalving(i) = zxfqcalving(i) +  &            zxfqcalving(i) = zxfqcalving(i) + &
1013                 fqcalving(i, nsrf)*pctsrf(i, nsrf)                 fqcalving(i, nsrf)*pctsrf(i, nsrf)
1014            s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)            s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
1015            s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)            s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
# Line 1236  contains Line 1028  contains
1028    
1029      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1030         DO i = 1, klon         DO i = 1, klon
1031            IF (pctsrf(i, nsrf)  <  epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
1032    
1033            IF (pctsrf(i, nsrf)  <  epsfra) t2m(i, nsrf) = zt2m(i)            IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
1034            IF (pctsrf(i, nsrf)  <  epsfra) q2m(i, nsrf) = zq2m(i)            IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
1035            IF (pctsrf(i, nsrf)  <  epsfra) u10m(i, nsrf) = zu10m(i)            IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
1036            IF (pctsrf(i, nsrf)  <  epsfra) v10m(i, nsrf) = zv10m(i)            IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
1037            IF (pctsrf(i, nsrf)  <  epsfra) ffonte(i, nsrf) = zxffonte(i)            IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
1038            IF (pctsrf(i, nsrf)  <  epsfra)  &            IF (pctsrf(i, nsrf) < epsfra) &
1039                 fqcalving(i, nsrf) = zxfqcalving(i)                 fqcalving(i, nsrf) = zxfqcalving(i)
1040            IF (pctsrf(i, nsrf)  <  epsfra) pblh(i, nsrf)=s_pblh(i)            IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
1041            IF (pctsrf(i, nsrf)  <  epsfra) plcl(i, nsrf)=s_lcl(i)            IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
1042            IF (pctsrf(i, nsrf)  <  epsfra) capCL(i, nsrf)=s_capCL(i)            IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
1043            IF (pctsrf(i, nsrf)  <  epsfra) oliqCL(i, nsrf)=s_oliqCL(i)            IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
1044            IF (pctsrf(i, nsrf)  <  epsfra) cteiCL(i, nsrf)=s_cteiCL(i)            IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
1045            IF (pctsrf(i, nsrf)  <  epsfra) pblT(i, nsrf)=s_pblT(i)            IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
1046            IF (pctsrf(i, nsrf)  <  epsfra) therm(i, nsrf)=s_therm(i)            IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
1047            IF (pctsrf(i, nsrf)  <  epsfra) trmb1(i, nsrf)=s_trmb1(i)            IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
1048            IF (pctsrf(i, nsrf)  <  epsfra) trmb2(i, nsrf)=s_trmb2(i)            IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
1049            IF (pctsrf(i, nsrf)  <  epsfra) trmb3(i, nsrf)=s_trmb3(i)            IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
1050         ENDDO         ENDDO
1051      ENDDO      ENDDO
1052    
1053      ! Calculer la derive du flux infrarouge      ! Calculer la derive du flux infrarouge
1054    
1055      DO i = 1, klon      DO i = 1, klon
1056         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
1057      ENDDO      ENDDO
1058    
1059      ! Appeler la convection (au choix)      ! Appeler la convection (au choix)
1060    
1061      DO k = 1, llm      DO k = 1, llm
1062         DO i = 1, klon         DO i = 1, klon
1063            conv_q(i, k) = d_q_dyn(i, k)  &            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys
1064                 + d_q_vdf(i, k)/pdtphys            conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys
           conv_t(i, k) = d_t_dyn(i, k)  &  
                + d_t_vdf(i, k)/pdtphys  
1065         ENDDO         ENDDO
1066      ENDDO      ENDDO
1067    
1068      IF (check) THEN      IF (check) THEN
1069         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1070         print *, "avantcon=", za         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)) &  
                   *(paprs(i, k)-paprs(i, k+1))/RG  
           ENDDO  
        ENDDO  
1071      ENDIF      ENDIF
1072      IF (iflag_con == 1) THEN  
1073         stop 'reactiver le call conlmd dans physiq.F'      if (iflag_con == 2) then
1074      ELSE IF (iflag_con == 2) THEN         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
1075         CALL conflx(pdtphys, paprs, pplay, t_seri, q_seri, &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
1076              conv_t, conv_q, zxfluxq(1, 1), omega, &              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
1077              d_t_con, d_q_con, rain_con, snow_con, &              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
1078              pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
1079              kcbot, kctop, kdtop, pmflxr, pmflxs)              kdtop, pmflxr, pmflxs)
1080         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
1081         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
1082         DO i = 1, klon         ibas_con = llm + 1 - kcbot
1083            ibas_con(i) = llm+1 - kcbot(i)         itop_con = llm + 1 - kctop
1084            itop_con(i) = llm+1 - kctop(i)      else
1085         ENDDO         ! iflag_con >= 3
     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  
1086    
1087         IF (.NOT. ok_gust) THEN         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &
1088            do i = 1, klon              v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &
1089               wd(i)=0.0              d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &
1090            enddo              itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &
1091         ENDIF              pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &
1092                wd, pmflxr, pmflxs, da, phi, mp, ntra=1)
1093           ! (number of tracers for the convection scheme of Kerry Emanuel:
1094           ! la partie traceurs est faite dans phytrac
1095           ! on met ntra = 1 pour limiter les appels mais on peut
1096           ! supprimer les calculs / ftra.)
1097    
1098           clwcon0 = qcondc
1099           mfu = upwd + dnwd
1100           IF (.NOT. ok_gust) wd = 0.
1101    
1102         ! Calcul des proprietes des nuages convectifs         ! Calcul des propri\'et\'es des nuages convectifs
1103    
1104         DO k = 1, llm         DO k = 1, llm
1105            DO i = 1, klon            DO i = 1, klon
1106               zx_t = t_seri(i, k)               zx_t = t_seri(i, k)
1107               IF (thermcep) THEN               IF (thermcep) THEN
1108                  zdelta = MAX(0., SIGN(1., rtt-zx_t))                  zdelta = MAX(0., SIGN(1., rtt-zx_t))
1109                  zx_qs  = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)                  zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)
1110                  zx_qs  = MIN(0.5, zx_qs)                  zx_qs = MIN(0.5, zx_qs)
1111                  zcor   = 1./(1.-retv*zx_qs)                  zcor = 1./(1.-retv*zx_qs)
1112                  zx_qs  = zx_qs*zcor                  zx_qs = zx_qs*zcor
1113               ELSE               ELSE
1114                  IF (zx_t < t_coup) THEN                  IF (zx_t < t_coup) THEN
1115                     zx_qs = qsats(zx_t)/pplay(i, k)                     zx_qs = qsats(zx_t)/play(i, k)
1116                  ELSE                  ELSE
1117                     zx_qs = qsatl(zx_t)/pplay(i, k)                     zx_qs = qsatl(zx_t)/play(i, k)
1118                  ENDIF                  ENDIF
1119               ENDIF               ENDIF
1120               zqsat(i, k)=zx_qs               zqsat(i, k) = zx_qs
1121            ENDDO            ENDDO
1122         ENDDO         ENDDO
1123    
1124         !   calcul des proprietes des nuages convectifs         ! calcul des proprietes des nuages convectifs
1125         clwcon0=fact_cldcon*clwcon0         clwcon0 = fact_cldcon * clwcon0
1126         call clouds_gno &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
1127              (klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, rnebcon0)              rnebcon0)
1128      ELSE  
1129         print *, "iflag_con non-prevu", iflag_con         mfd = 0.
1130         stop 1         pen_u = 0.
1131      ENDIF         pen_d = 0.
1132           pde_d = 0.
1133           pde_u = 0.
1134        END if
1135    
1136      DO k = 1, llm      DO k = 1, llm
1137         DO i = 1, klon         DO i = 1, klon
# Line 1391  contains Line 1143  contains
1143      ENDDO      ENDDO
1144    
1145      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1146         ztit='after convect'         tit = 'after convect'
1147         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1148              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1149              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1150         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1151              , zero_v, zero_v, zero_v, zero_v, zero_v &              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &
1152              , zero_v, rain_con, snow_con, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
1153      END IF      END IF
1154    
1155      IF (check) THEN      IF (check) THEN
1156         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1157         print *,"aprescon=", za         print *, "aprescon = ", za
1158         zx_t = 0.0         zx_t = 0.
1159         za = 0.0         za = 0.
1160         DO i = 1, klon         DO i = 1, klon
1161            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1162            zx_t = zx_t + (rain_con(i)+ &            zx_t = zx_t + (rain_con(i)+ &
1163                 snow_con(i))*airephy(i)/REAL(klon)                 snow_con(i))*airephy(i)/REAL(klon)
1164         ENDDO         ENDDO
1165         zx_t = zx_t/za*pdtphys         zx_t = zx_t/za*dtphys
1166         print *,"Precip=", zx_t         print *, "Precip = ", zx_t
1167      ENDIF      ENDIF
1168      IF (zx_ajustq) THEN  
1169         DO i = 1, klon      IF (iflag_con == 2) THEN
1170            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1171         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)) &  
                   *(paprs(i, k)-paprs(i, k+1))/RG  
           ENDDO  
        ENDDO  
        DO i = 1, klon  
           z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*pdtphys) &  
                /z_apres(i)  
        ENDDO  
1172         DO k = 1, llm         DO k = 1, llm
1173            DO i = 1, klon            DO i = 1, klon
1174               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  
1175                  q_seri(i, k) = q_seri(i, k) * z_factor(i)                  q_seri(i, k) = q_seri(i, k) * z_factor(i)
1176               ENDIF               ENDIF
1177            ENDDO            ENDDO
1178         ENDDO         ENDDO
1179      ENDIF      ENDIF
     zx_ajustq=.FALSE.  
1180    
1181      ! Convection seche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
1182    
1183      d_t_ajs=0.      d_t_ajs = 0.
1184      d_u_ajs=0.      d_u_ajs = 0.
1185      d_v_ajs=0.      d_v_ajs = 0.
1186      d_q_ajs=0.      d_q_ajs = 0.
1187      fm_therm=0.      fm_therm = 0.
1188      entr_therm=0.      entr_therm = 0.
1189    
1190      IF(prt_level>9)print *, &      if (iflag_thermals == 0) then
1191           'AVANT LA CONVECTION SECHE, iflag_thermals=' &         ! Ajustement sec
1192           , 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)  
1193         t_seri = t_seri + d_t_ajs         t_seri = t_seri + d_t_ajs
1194         q_seri = q_seri + d_q_ajs         q_seri = q_seri + d_q_ajs
1195      else      else
1196         !  Thermiques         ! Thermiques
1197         IF(prt_level>9)print *,'JUSTE AVANT, iflag_thermals=' &         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
1198              , iflag_thermals, '   nsplit_thermals=', nsplit_thermals              q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
        call calltherm(pdtphys &  
             , pplay, paprs, pphi &  
             , u_seri, v_seri, t_seri, q_seri &  
             , d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs &  
             , fm_therm, entr_therm)  
1199      endif      endif
1200    
1201      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1202         ztit='after dry_adjust'         tit = 'after dry_adjust'
1203         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1204              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1205              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1206      END IF      END IF
1207    
1208      !  Caclul des ratqs      ! Caclul des ratqs
1209    
1210      !   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
1211      !   on ecrase le tableau ratqsc calcule par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
1212      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
1213         do k=1, llm         do k = 1, llm
1214            do i=1, klon            do i = 1, klon
1215               if(ptconv(i, k)) then               if(ptconv(i, k)) then
1216                  ratqsc(i, k)=ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
1217                       +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)
1218               else               else
1219                  ratqsc(i, k)=0.                  ratqsc(i, k) = 0.
1220               endif               endif
1221            enddo            enddo
1222         enddo         enddo
1223      endif      endif
1224    
1225      !   ratqs stables      ! ratqs stables
1226      do k=1, llm      do k = 1, llm
1227         do i=1, klon         do i = 1, klon
1228            ratqss(i, k)=ratqsbas+(ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1229                 min((paprs(i, 1)-pplay(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1230         enddo         enddo
1231      enddo      enddo
1232    
1233      !  ratqs final      ! ratqs final
1234      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1235         !   les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
1236         !   ratqs final         ! ratqs final
1237         !   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
1238         !   relaxation des ratqs         ! relaxation des ratqs
1239         facteur=exp(-pdtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
1240         ratqs=max(ratqs*facteur, ratqss)         ratqs = max(ratqs, ratqsc)
        ratqs=max(ratqs, ratqsc)  
1241      else      else
1242         !   on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
1243         ratqs=ratqss         ratqs = ratqss
1244      endif      endif
1245    
1246      ! Appeler le processus de condensation a grande echelle      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1247      ! et le processus de precipitation           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
1248      CALL fisrtilp(pdtphys, paprs, pplay, &           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
1249           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)  
1250    
1251      WHERE (rain_lsc < 0) rain_lsc = 0.      WHERE (rain_lsc < 0) rain_lsc = 0.
1252      WHERE (snow_lsc < 0) snow_lsc = 0.      WHERE (snow_lsc < 0) snow_lsc = 0.
# Line 1541  contains Line 1261  contains
1261      ENDDO      ENDDO
1262      IF (check) THEN      IF (check) THEN
1263         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1264         print *,"apresilp=", za         print *, "apresilp = ", za
1265         zx_t = 0.0         zx_t = 0.
1266         za = 0.0         za = 0.
1267         DO i = 1, klon         DO i = 1, klon
1268            za = za + airephy(i)/REAL(klon)            za = za + airephy(i)/REAL(klon)
1269            zx_t = zx_t + (rain_lsc(i) &            zx_t = zx_t + (rain_lsc(i) &
1270                 + snow_lsc(i))*airephy(i)/REAL(klon)                 + snow_lsc(i))*airephy(i)/REAL(klon)
1271         ENDDO         ENDDO
1272         zx_t = zx_t/za*pdtphys         zx_t = zx_t/za*dtphys
1273         print *,"Precip=", zx_t         print *, "Precip = ", zx_t
1274      ENDIF      ENDIF
1275    
1276      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1277         ztit='after fisrt'         tit = 'after fisrt'
1278         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1279              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1280              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1281         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1282              , zero_v, zero_v, zero_v, zero_v, zero_v &              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &
1283              , zero_v, rain_lsc, snow_lsc, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
1284      END IF      END IF
1285    
1286      !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1287    
1288      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1289    
1290      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke      IF (iflag_cldcon <= -1) THEN
1291         snow_tiedtke=0.         ! seulement pour Tiedtke
1292           snow_tiedtke = 0.
1293         if (iflag_cldcon == -1) then         if (iflag_cldcon == -1) then
1294            rain_tiedtke=rain_con            rain_tiedtke = rain_con
1295         else         else
1296            rain_tiedtke=0.            rain_tiedtke = 0.
1297            do k=1, llm            do k = 1, llm
1298               do i=1, klon               do i = 1, klon
1299                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1300                     rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i, k)/pdtphys &                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &
1301                          *(paprs(i, k)-paprs(i, k+1))/rg                          *zmasse(i, k)
1302                  endif                  endif
1303               enddo               enddo
1304            enddo            enddo
1305         endif         endif
1306    
1307         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1308         CALL diagcld1(paprs, pplay, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1309              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1310         DO k = 1, llm         DO k = 1, llm
1311            DO i = 1, klon            DO i = 1, klon
1312               IF (diafra(i, k).GT.cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
1313                  cldliq(i, k) = dialiq(i, k)                  cldliq(i, k) = dialiq(i, k)
1314                  cldfra(i, k) = diafra(i, k)                  cldfra(i, k) = diafra(i, k)
1315               ENDIF               ENDIF
1316            ENDDO            ENDDO
1317         ENDDO         ENDDO
   
1318      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1319         ! On prend pour les nuages convectifs le max du calcul de la         ! On prend pour les nuages convectifs le maximum du calcul de
1320         ! 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
1321         ! facttemps         ! d'un facteur facttemps.
1322         facteur = pdtphys *facttemps         facteur = dtphys * facttemps
1323         do k=1, llm         do k = 1, llm
1324            do i=1, klon            do i = 1, klon
1325               rnebcon(i, k)=rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1326               if (rnebcon0(i, k)*clwcon0(i, k).gt.rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k) * clwcon0(i, k) &
1327                    then                    > rnebcon(i, k) * clwcon(i, k)) then
1328                  rnebcon(i, k)=rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
1329                  clwcon(i, k)=clwcon0(i, k)                  clwcon(i, k) = clwcon0(i, k)
1330               endif               endif
1331            enddo            enddo
1332         enddo         enddo
1333    
1334         !   On prend la somme des fractions nuageuses et des contenus en eau         ! On prend la somme des fractions nuageuses et des contenus en eau
1335         cldfra=min(max(cldfra, rnebcon), 1.)         cldfra = min(max(cldfra, rnebcon), 1.)
1336         cldliq=cldliq+rnebcon*clwcon         cldliq = cldliq + rnebcon*clwcon
   
1337      ENDIF      ENDIF
1338    
1339      ! 2. NUAGES STARTIFORMES      ! 2. Nuages stratiformes
1340    
1341      IF (ok_stratus) THEN      IF (ok_stratus) THEN
1342         CALL diagcld2(paprs, pplay, t_seri, q_seri, diafra, dialiq)         CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1343         DO k = 1, llm         DO k = 1, llm
1344            DO i = 1, klon            DO i = 1, klon
1345               IF (diafra(i, k).GT.cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
1346                  cldliq(i, k) = dialiq(i, k)                  cldliq(i, k) = dialiq(i, k)
1347                  cldfra(i, k) = diafra(i, k)                  cldfra(i, k) = diafra(i, k)
1348               ENDIF               ENDIF
# Line 1635  contains Line 1351  contains
1351      ENDIF      ENDIF
1352    
1353      ! Precipitation totale      ! Precipitation totale
   
1354      DO i = 1, klon      DO i = 1, klon
1355         rain_fall(i) = rain_con(i) + rain_lsc(i)         rain_fall(i) = rain_con(i) + rain_lsc(i)
1356         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
1357      ENDDO      ENDDO
1358    
1359      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1360         ztit="after diagcld"           dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1361         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &           d_h_vcol, d_qt, d_qw, d_ql, d_qs, 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  
1362    
1363        ! Humidit\'e relative pour diagnostic :
1364      DO k = 1, llm      DO k = 1, llm
1365         DO i = 1, klon         DO i = 1, klon
1366            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1367            IF (thermcep) THEN            IF (thermcep) THEN
1368               zdelta = MAX(0., SIGN(1., rtt-zx_t))               zdelta = MAX(0., SIGN(1., rtt-zx_t))
1369               zx_qs  = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)               zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)
1370               zx_qs  = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1371               zcor   = 1./(1.-retv*zx_qs)               zcor = 1./(1.-retv*zx_qs)
1372               zx_qs  = zx_qs*zcor               zx_qs = zx_qs*zcor
1373            ELSE            ELSE
1374               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
1375                  zx_qs = qsats(zx_t)/pplay(i, k)                  zx_qs = qsats(zx_t)/play(i, k)
1376               ELSE               ELSE
1377                  zx_qs = qsatl(zx_t)/pplay(i, k)                  zx_qs = qsatl(zx_t)/play(i, k)
1378               ENDIF               ENDIF
1379            ENDIF            ENDIF
1380            zx_rh(i, k) = q_seri(i, k)/zx_qs            zx_rh(i, k) = q_seri(i, k)/zx_qs
1381            zqsat(i, k)=zx_qs            zqsat(i, k) = zx_qs
1382         ENDDO         ENDDO
1383      ENDDO      ENDDO
1384      !jq - introduce the aerosol direct and first indirect radiative forcings  
1385      !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)      ! Introduce the aerosol direct and first indirect radiative forcings:
1386      IF (ok_ade.OR.ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1387         ! Get sulfate aerosol distribution         ! Get sulfate aerosol distribution :
1388         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(rdayvrai, firstcal, sulfate)
1389         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1390    
1391         ! Calculate aerosol optical properties (Olivier Boucher)         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1392         CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &              aerindex)
             tau_ae, piz_ae, cg_ae, aerindex)  
1393      ELSE      ELSE
1394         tau_ae(:, :, :)=0.0         tau_ae = 0.
1395         piz_ae(:, :, :)=0.0         piz_ae = 0.
1396         cg_ae(:, :, :)=0.0         cg_ae = 0.
1397      ENDIF      ENDIF
1398    
1399      ! Calculer les parametres optiques des nuages et quelques      ! Param\`etres optiques des nuages et quelques param\`etres pour diagnostics :
     ! parametres pour diagnostiques:  
   
1400      if (ok_newmicro) then      if (ok_newmicro) then
1401         CALL newmicro (paprs, pplay, ok_newmicro, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1402              t_seri, cldliq, cldfra, cldtau, cldemi, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1403              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)  
1404      else      else
1405         CALL nuage (paprs, pplay, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1406              t_seri, cldliq, cldfra, cldtau, cldemi, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1407              cldh, cldl, cldm, cldt, cldq, &              bl95_b1, cldtaupi, re, fl)
             ok_aie, &  
             sulfate, sulfate_pi, &  
             bl95_b0, bl95_b1, &  
             cldtaupi, re, fl)  
   
1408      endif      endif
1409    
1410      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
   
1411      IF (MOD(itaprad, radpas) == 0) THEN      IF (MOD(itaprad, radpas) == 0) THEN
1412         DO i = 1, klon         DO i = 1, klon
1413            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
# Line 1722  contains Line 1419  contains
1419                 + falblw(i, is_ter) * pctsrf(i, is_ter) &                 + falblw(i, is_ter) * pctsrf(i, is_ter) &
1420                 + falblw(i, is_sic) * pctsrf(i, is_sic)                 + falblw(i, is_sic) * pctsrf(i, is_sic)
1421         ENDDO         ENDDO
1422         ! nouveau rayonnement (compatible Arpege-IFS):         ! Rayonnement (compatible Arpege-IFS) :
1423         CALL radlwsw(dist, rmu0, fract,  &         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &
1424              paprs, pplay, zxtsol, albsol, albsollw, t_seri, q_seri, &              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &
1425              wo, &              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &
1426              cldfra, cldemi, cldtau, &              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &
1427              heat, heat0, cool, cool0, radsol, albpla, &              lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &
1428              topsw, toplw, solsw, sollw, &              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)
             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)  
1429         itaprad = 0         itaprad = 0
1430      ENDIF      ENDIF
1431      itaprad = itaprad + 1      itaprad = itaprad + 1
# Line 1746  contains Line 1434  contains
1434    
1435      DO k = 1, llm      DO k = 1, llm
1436         DO i = 1, klon         DO i = 1, klon
1437            t_seri(i, k) = t_seri(i, k) &            t_seri(i, k) = t_seri(i, k) + (heat(i, k)-cool(i, k)) * dtphys/86400.
                + (heat(i, k)-cool(i, k)) * pdtphys/86400.  
1438         ENDDO         ENDDO
1439      ENDDO      ENDDO
1440    
1441      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1442         ztit='after rad'         tit = 'after rad'
1443         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1444              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1445              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1446         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1447              , topsw, toplw, solsw, sollw, zero_v &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
1448              , zero_v, zero_v, zero_v, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
1449      END IF      END IF
1450    
1451      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
   
1452      DO i = 1, klon      DO i = 1, klon
1453         zxqsurf(i) = 0.0         zxqsurf(i) = 0.
1454         zxsnow(i) = 0.0         zxsnow(i) = 0.
1455      ENDDO      ENDDO
1456      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1457         DO i = 1, klon         DO i = 1, klon
# Line 1776  contains Line 1460  contains
1460         ENDDO         ENDDO
1461      ENDDO      ENDDO
1462    
1463      ! Calculer le bilan du sol et la derive de temperature (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1464    
1465      DO i = 1, klon      DO i = 1, klon
1466         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1467      ENDDO      ENDDO
1468    
1469      !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:  
1470    
1471      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1472         !  selection des points pour lesquels le shema est actif:         ! selection des points pour lesquels le shema est actif:
1473         igwd=0         igwd = 0
1474         DO i=1, klon         DO i = 1, klon
1475            itest(i)=0            itest(i) = 0
1476            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
1477               itest(i)=1               itest(i) = 1
1478               igwd=igwd+1               igwd = igwd + 1
1479               idx(igwd)=i               idx(igwd) = i
1480            ENDIF            ENDIF
1481         ENDDO         ENDDO
1482    
1483         CALL drag_noro(klon, llm, pdtphys, paprs, pplay, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1484              zmea, zstd, zsig, zgam, zthe, zpic, zval, &              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &
1485              igwd, idx, itest, &              zulow, zvlow, 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)  
1486    
1487         !  ajout des tendances         ! ajout des tendances
1488         DO k = 1, llm         DO k = 1, llm
1489            DO i = 1, klon            DO i = 1, klon
1490               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 1816  contains Line 1495  contains
1495      ENDIF      ENDIF
1496    
1497      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1498           ! S\'election des points pour lesquels le sch\'ema est actif :
1499         !  selection des points pour lesquels le shema est actif:         igwd = 0
1500         igwd=0         DO i = 1, klon
1501         DO i=1, klon            itest(i) = 0
1502            itest(i)=0            IF ((zpic(i) - zmea(i)) > 100.) THEN
1503            IF ((zpic(i)-zmea(i)).GT.100.) THEN               itest(i) = 1
1504               itest(i)=1               igwd = igwd + 1
1505               igwd=igwd+1               idx(igwd) = i
              idx(igwd)=i  
1506            ENDIF            ENDIF
1507         ENDDO         ENDDO
1508    
1509         CALL lift_noro(klon, llm, pdtphys, paprs, pplay, &         CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1510              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, &  
1511              d_t_lif, d_u_lif, d_v_lif)              d_t_lif, d_u_lif, d_v_lif)
1512    
1513         !  ajout des tendances         ! Ajout des tendances :
1514         DO k = 1, llm         DO k = 1, llm
1515            DO i = 1, klon            DO i = 1, klon
1516               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 1843  contains Line 1518  contains
1518               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)
1519            ENDDO            ENDDO
1520         ENDDO         ENDDO
1521        ENDIF
1522    
1523      ENDIF ! fin de test sur ok_orolf      ! Stress n\'ecessaires : toute la physique
   
     ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE  
1524    
1525      DO i = 1, klon      DO i = 1, klon
1526         zustrph(i)=0.         zustrph(i) = 0.
1527         zvstrph(i)=0.         zvstrph(i) = 0.
1528      ENDDO      ENDDO
1529      DO k = 1, llm      DO k = 1, llm
1530         DO i = 1, klon         DO i = 1, klon
1531            zustrph(i)=zustrph(i)+(u_seri(i, k)-u(i, k))/pdtphys* &            zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1532                 (paprs(i, k)-paprs(i, k+1))/rg                 * zmasse(i, k)
1533            zvstrph(i)=zvstrph(i)+(v_seri(i, k)-v(i, k))/pdtphys* &            zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1534                 (paprs(i, k)-paprs(i, k+1))/rg                 * zmasse(i, k)
1535         ENDDO         ENDDO
1536      ENDDO      ENDDO
1537    
1538      !IM calcul composantes axiales du moment angulaire et couple des montagnes      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &
1539             zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1540      CALL aaam_bud (27, klon, llm, gmtime, &  
1541           ra, rg, romega, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1542           rlat, rlon, pphis, &           2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1543           zustrdr, zustrli, zustrph, &           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1544           zvstrdr, zvstrli, zvstrph, &  
1545           paprs, u, v, &      ! Calcul des tendances traceurs
1546           aam, torsfc)      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &
1547             dtphys, u, t, paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, &
1548      IF (if_ebil >= 2) THEN           entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, &
1549         ztit='after orography'           albsol, rhcl, cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, &
1550         CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &           psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)
1551              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &  
1552              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1553      END IF           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1554             pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
     !AA Installation de l'interface online-offline pour traceurs  
   
     !   Calcul  des tendances traceurs  
   
     call phytrac(rnpb, itap, lmt_pas, julien,  gmtime, firstcal, lafin, nq-2, &  
          pdtphys, u, v, t, paprs, pplay, pmfu,  pmfd,  pen_u,  pde_u,  pen_d, &  
          pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &  
          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)  
   
     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  
1555    
1556      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1557        CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &
1558             ue, uq)
1559    
1560      CALL transp (paprs, zxtsol, &      ! diag. bilKP
          t_seri, q_seri, u_seri, v_seri, zphi, &  
          ve, vq, ue, uq)  
1561    
1562      !IM diag. bilKP      CALL transp_lay(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
   
     CALL transp_lay (paprs, zxtsol, &  
          t_seri, q_seri, u_seri, v_seri, zphi, &  
1563           ve_lay, vq_lay, ue_lay, uq_lay)           ve_lay, vq_lay, ue_lay, uq_lay)
1564    
1565      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
1566    
1567      !+jld ec_conser      ! conversion Ec -> E thermique
1568      DO k = 1, llm      DO k = 1, llm
1569         DO i = 1, klon         DO i = 1, klon
1570            ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i, k))            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1571            d_t_ec(i, k)=0.5/ZRCPD &            d_t_ec(i, k) = 0.5 / ZRCPD &
1572                 *(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)
1573            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)
1574            d_t_ec(i, k) = d_t_ec(i, k)/pdtphys            d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1575         END DO         END DO
1576      END DO      END DO
1577      !-jld ec_conser  
1578      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1579         ztit='after physic'         tit = 'after physic'
1580         CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1581              , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1582              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)              d_ql, d_qs, d_ec)
1583         !     Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1584         !     on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1585         !     est egale a la variation de la physique au pas de temps precedent.         ! est egale a la variation de la physique au pas de temps precedent.
1586         !     Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1587         call diagphy(airephy, ztit, ip_ebil &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1588              , topsw, toplw, solsw, sollw, sens &              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &
1589              , evap, rain_fall, snow_fall, ztsol &              fs_bound, fq_bound)
             , d_h_vcol, d_qt, d_ec &  
             , fs_bound, fq_bound )  
1590    
1591         d_h_vcol_phy=d_h_vcol         d_h_vcol_phy = d_h_vcol
1592    
1593      END IF      END IF
1594    
1595      !   SORTIES      ! SORTIES
   
     !IM Interpolation sur les niveaux de pression du NMC  
     call calcul_STDlev  
1596    
1597      !cc prw = eau precipitable      ! prw = eau precipitable
1598      DO i = 1, klon      DO i = 1, klon
1599         prw(i) = 0.         prw(i) = 0.
1600         DO k = 1, llm         DO k = 1, llm
1601            prw(i) = prw(i) + &            prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
                q_seri(i, k)*(paprs(i, k)-paprs(i, k+1))/RG  
1602         ENDDO         ENDDO
1603      ENDDO      ENDDO
1604    
     !IM initialisation + calculs divers diag AMIP2  
     call calcul_divers  
   
1605      ! Convertir les incrementations en tendances      ! Convertir les incrementations en tendances
1606    
1607      DO k = 1, llm      DO k = 1, llm
1608         DO i = 1, klon         DO i = 1, klon
1609            d_u(i, k) = ( u_seri(i, k) - u(i, k) ) / pdtphys            d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1610            d_v(i, k) = ( v_seri(i, k) - v(i, k) ) / pdtphys            d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1611            d_t(i, k) = ( t_seri(i, k)-t(i, k) ) / pdtphys            d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1612            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
1613            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
1614         ENDDO         ENDDO
1615      ENDDO      ENDDO
1616    
1617      IF (nq >= 3) THEN      IF (nqmx >= 3) THEN
1618         DO iq = 3, nq         DO iq = 3, nqmx
1619            DO  k = 1, llm            DO k = 1, llm
1620               DO  i = 1, klon               DO i = 1, klon
1621                  d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / pdtphys                  d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys
1622               ENDDO               ENDDO
1623            ENDDO            ENDDO
1624         ENDDO         ENDDO
1625      ENDIF      ENDIF
1626    
1627      ! 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:
   
1628      DO k = 1, llm      DO k = 1, llm
1629         DO i = 1, klon         DO i = 1, klon
1630            t_ancien(i, k) = t_seri(i, k)            t_ancien(i, k) = t_seri(i, k)
# Line 1993  contains Line 1632  contains
1632         ENDDO         ENDDO
1633      ENDDO      ENDDO
1634    
1635      !   Ecriture des sorties      ! Ecriture des sorties
   
     call write_histhf  
     call write_histday  
1636      call write_histins      call write_histins
1637    
1638      ! Si c'est la fin, il faut conserver l'etat de redemarrage      ! Si c'est la fin, il faut conserver l'etat de redemarrage
   
1639      IF (lafin) THEN      IF (lafin) THEN
1640         itau_phy = itau_phy + itap         itau_phy = itau_phy + itap
1641         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, &         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &
1642              ftsoil, tslab, seaice, fqsurf, qsol, &              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &
1643              fsnow, falbe, falblw, fevap, rain_fall, snow_fall, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
1644              solsw, sollwdown, dlw, &              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
1645              radsol, frugs, agesno, &              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
             zmea, zstd, zsig, zgam, zthe, zpic, zval, &  
             t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)  
1646      ENDIF      ENDIF
1647    
1648    contains      firstcal = .FALSE.
   
     subroutine calcul_STDlev  
   
       !     From phylmd/calcul_STDlev.h, v 1.1 2005/05/25 13:10:09  
   
       !IM on initialise les champs en debut du jour ou du mois  
   
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, tsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, usumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, vsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, wsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, phisumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, qsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, rhsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, uvsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, vqsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, vTsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, wqsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, vphisumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, wTsumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, u2sumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, v2sumSTD)  
       CALL ini_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, &  
            tnondef, T2sumSTD)  
   
       !IM on interpole sur les niveaux STD de pression a chaque pas de  
       !temps de la physique  
   
       DO k=1, nlevSTD  
   
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               t_seri, tlevSTD(:, k))  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               u_seri, ulevSTD(:, k))  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               v_seri, vlevSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=paprs(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., zx_tmp_fi3d, rlevSTD(k), &  
               omega, wlevSTD(:, k))  
   
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zphi/RG, philevSTD(:, k))  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               qx(:, :, ivap), qlevSTD(:, k))  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_rh*100., rhlevSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=u_seri(i, l)*v_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, uvSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=v_seri(i, l)*q_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, vqSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=v_seri(i, l)*t_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, vTSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=omega(i, l)*qx(i, l, ivap)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, wqSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=v_seri(i, l)*zphi(i, l)/RG  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, vphiSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=omega(i, l)*t_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, wTSTD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=u_seri(i, l)*u_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, u2STD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=v_seri(i, l)*v_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, v2STD(:, k))  
   
          DO l=1, llm  
             DO i=1, klon  
                zx_tmp_fi3d(i, l)=t_seri(i, l)*t_seri(i, l)  
             ENDDO !i  
          ENDDO !l  
          CALL plevel(klon, llm, .true., pplay, rlevSTD(k), &  
               zx_tmp_fi3d, T2STD(:, k))  
   
       ENDDO !k=1, nlevSTD  
   
       !IM on somme les valeurs definies a chaque pas de temps de la  
       ! physique ou toutes les 6 heures  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.TRUE.  
       CALL undefSTD(nlevSTD, itap, tlevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, tsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, ulevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, usumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, vlevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, vsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, wlevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, wsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, philevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, phisumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, qlevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, qsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, rhlevSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, rhsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, uvSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, uvsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, vqSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, vqsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, vTSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, vTsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, wqSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, wqsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, vphiSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, vphisumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, wTSTD, &  
            ecrit_hf, &  
            oknondef, tnondef, wTsumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, u2STD, &  
            ecrit_hf, &  
            oknondef, tnondef, u2sumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, v2STD, &  
            ecrit_hf, &  
            oknondef, tnondef, v2sumSTD)  
   
       oknondef(1:klon, 1:nlevSTD, 1:nout)=.FALSE.  
       CALL undefSTD(nlevSTD, itap, T2STD, &  
            ecrit_hf, &  
            oknondef, tnondef, T2sumSTD)  
   
       !IM on moyenne a la fin du jour ou du mois  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, tsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, usumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, vsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, wsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, phisumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, qsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, rhsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, uvsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, vqsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, vTsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, wqsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, vphisumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, wTsumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, u2sumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, v2sumSTD)  
   
       CALL moy_undefSTD(nlevSTD, itap, &  
            ecrit_day, ecrit_mth, ecrit_hf2mth, &  
            tnondef, T2sumSTD)  
   
       !IM interpolation a chaque pas de temps du SWup(clr) et  
       !SWdn(clr) a 200 hPa  
   
       CALL plevel(klon, klevp1, .true., paprs, 20000., &  
            swdn0, SWdn200clr)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            swdn, SWdn200)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            swup0, SWup200clr)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            swup, SWup200)  
   
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            lwdn0, LWdn200clr)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            lwdn, LWdn200)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            lwup0, LWup200clr)  
       CALL plevel(klon, klevp1, .false., paprs, 20000., &  
            lwup, LWup200)  
   
     end SUBROUTINE calcul_STDlev  
   
     !****************************************************  
   
     SUBROUTINE calcul_divers  
   
       ! From phylmd/calcul_divers.h, v 1.1 2005/05/25 13:10:09  
   
       ! initialisations diverses au "debut" du mois  
   
       IF(MOD(itap, ecrit_mth) == 1) THEN  
          DO i=1, klon  
             nday_rain(i)=0.  
          ENDDO  
       ENDIF  
   
       IF(MOD(itap, ecrit_day) == 0) THEN  
          !IM calcul total_rain, nday_rain  
          DO i = 1, klon  
             total_rain(i)=rain_fall(i)+snow_fall(i)    
             IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.  
          ENDDO  
       ENDIF  
   
     End SUBROUTINE calcul_divers  
   
     !***********************************************  
   
     subroutine write_histday  
   
       !     From phylmd/write_histday.h, v 1.3 2005/05/25 13:10:09  
   
       if (ok_journe) THEN  
   
          ndex2d = 0  
          ndex3d = 0  
   
          ! Champs 2D:  
   
          itau_w = itau_phy + itap  
   
          !   FIN ECRITURE DES CHAMPS 3D  
1649    
1650           if (ok_sync) then    contains
             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  
   
       ndex2d = 0  
       ndex3d = 0  
   
       itau_w = itau_phy + itap  
   
       call write_histhf3d  
   
       IF (ok_sync) THEN  
          call histsync(nid_hf)  
       ENDIF  
   
     end subroutine write_histhf  
   
     !***************************************************************  
1651    
1652      subroutine write_histins      subroutine write_histins
1653    
1654        ! From phylmd/write_histins.h, v 1.2 2005/05/25 13:10:09        ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09
1655    
1656          use dimens_m, only: iim, jjm
1657          USE histsync_m, ONLY: histsync
1658          USE histwrite_m, ONLY: histwrite
1659    
1660        real zout        real zout
1661          integer itau_w ! pas de temps ecriture
1662          REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1663    
1664        !--------------------------------------------------        !--------------------------------------------------
1665    
1666        IF (ok_instan) THEN        IF (ok_instan) THEN
   
          ndex2d = 0  
          ndex3d = 0  
   
1667           ! Champs 2D:           ! Champs 2D:
1668    
1669           zsto = pdtphys * ecrit_ins           zsto = dtphys * ecrit_ins
1670           zout = pdtphys * ecrit_ins           zout = dtphys * ecrit_ins
1671           itau_w = itau_phy + itap           itau_w = itau_phy + itap
1672    
1673           i = NINT(zout/zsto)           i = NINT(zout/zsto)
1674           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), pphis, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)
1675           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1676    
1677           i = NINT(zout/zsto)           i = NINT(zout/zsto)
1678           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), airephy, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)
1679           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1680    
1681           DO i = 1, klon           DO i = 1, klon
1682              zx_tmp_fi2d(i) = paprs(i, 1)              zx_tmp_fi2d(i) = paprs(i, 1)
1683           ENDDO           ENDDO
1684           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1685           CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)
1686    
1687           DO i = 1, klon           DO i = 1, klon
1688              zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)              zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
1689           ENDDO           ENDDO
1690           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1691           CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)
1692    
1693           DO i = 1, klon           DO i = 1, klon
1694              zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)              zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
1695           ENDDO           ENDDO
1696           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1697           CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)
1698    
1699           DO i = 1, klon           DO i = 1, klon
1700              zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)              zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
1701           ENDDO           ENDDO
1702           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1703           CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)
1704    
1705           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxtsol, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)
1706           CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)
1707           !ccIM           !ccIM
1708           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zt2m, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)
1709           CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)
1710    
1711           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zq2m, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)
1712           CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)
1713    
1714           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zu10m, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)
1715           CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)
1716    
1717           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zv10m, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)
1718           CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)
1719    
1720           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), snow_fall, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)
1721           CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)
1722    
1723           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragm, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)
1724           CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)
1725    
1726           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragh, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)
1727           CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)
1728    
1729           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), toplw, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)
1730           CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)
1731    
1732           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), evap, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)
1733           CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)
1734    
1735           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), solsw, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)
1736           CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)
1737    
1738           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollw, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)
1739           CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)
1740    
1741           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollwdown, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)
1742           CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)
               ndex2d)  
1743    
1744           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), bils, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)
1745           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1746    
1747           zx_tmp_fi2d(1:klon)=-1*sens(1:klon)           zx_tmp_fi2d(1:klon) = -1*sens(1:klon)
1748           !     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sens, zx_tmp_2d)           ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)
1749           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1750           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
1751    
1752           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), fder, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)
1753           CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)
1754    
1755           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_oce), zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)
1756           CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)
               ndex2d)  
1757    
1758           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_ter), zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)
1759           CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)
               ndex2d)  
1760    
1761           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_lic), zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)
1762           CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)
               ndex2d)  
1763    
1764           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_sic), zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)
1765           CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)
               ndex2d)  
1766    
1767           DO nsrf = 1, nbsrf           DO nsrf = 1, nbsrf
1768              !XXX              !XXX
1769              zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)*100.              zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.
1770              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1771              CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &
1772                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1773    
1774              zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)
1775              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1776              CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &
1777                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1778    
1779              zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)              zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)
1780              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1781              CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &
1782                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1783    
1784              zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)
1785              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1786              CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &
1787                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1788    
1789              zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)
1790              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1791              CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &
1792                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1793    
1794              zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)              zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)
1795              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1796              CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &
1797                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1798    
1799              zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)              zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)
1800              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1801              CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &
1802                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1803    
1804              zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)
1805              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1806              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1807                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1808    
1809              zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)
1810              CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1811              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1812                   zx_tmp_2d, iim*(jjm + 1), ndex2d)                   zx_tmp_2d)
1813    
1814           END DO           END DO
1815           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsol, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)
1816           CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
1817           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsollw, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)
1818           CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           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, iim*(jjm + 1), ndex2d)  
1819    
1820           !IM cf. AM 081204 BEG           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)
1821             CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
1822    
1823           !HBTM2           !HBTM2
1824    
1825           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblh, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)
1826           CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)
1827    
1828           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblt, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)
1829           CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)
1830    
1831           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_lcl, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)
1832           CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d, iim*(jjm + 1), ndex2d)           CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)
1833    
1834           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_capCL, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)
1835           CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)
               ndex2d)  
1836    
1837           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_oliqCL, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)
1838           CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)
               ndex2d)  
1839    
1840           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_cteiCL, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)
1841           CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)
               ndex2d)  
1842    
1843           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_therm, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)
1844           CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)
               ndex2d)  
1845    
1846           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb1, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)
1847           CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)
               ndex2d)  
1848    
1849           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb2, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)
1850           CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)
               ndex2d)  
1851    
1852           CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb3, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)
1853           CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d, iim*(jjm + 1), &           CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
               ndex2d)  
   
          !IM cf. AM 081204 END  
1854    
1855           ! Champs 3D:           ! Champs 3D:
1856    
1857           CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)
1858           CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d, &           CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)
               iim*(jjm + 1)*llm, ndex3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)  
          CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), zphi, zx_tmp_3d)  
          CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), pplay, zx_tmp_3d)  
          CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          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, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          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, &  
               iim*(jjm + 1)*llm, ndex3d)  
   
          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  
   
       ndex2d = 0  
       ndex3d = 0  
   
       itau_w = itau_phy + itap  
1859    
1860        ! Champs 3D:           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)
1861             CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)
1862    
1863        CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)
1864        CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d, &           CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)
            iim*(jjm + 1)*llm, ndex3d)  
1865    
1866        CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), qx(1, 1, ivap), zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)
1867        CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d, &           CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)
            iim*(jjm + 1)*llm, ndex3d)  
1868    
1869        CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)
1870        CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d, &           CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)
            iim*(jjm + 1)*llm, ndex3d)  
1871    
1872        CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)
1873        CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d, &           CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)
            iim*(jjm + 1)*llm, ndex3d)  
1874    
1875        if (nbtr >= 3) then           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)
1876           CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), tr_seri(1, 1, 3), &           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
               zx_tmp_3d)  
          CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d, iim*(jjm + 1)*llm, &  
               ndex3d)  
       end if  
1877    
1878        if (ok_sync) then           if (ok_sync) then
1879           call histsync(nid_hf3d)              call histsync(nid_ins)
1880        endif           endif
1881          ENDIF
1882    
1883      end subroutine write_histhf3d      end subroutine write_histins
1884    
1885    END SUBROUTINE physiq    END SUBROUTINE physiq
1886    
   !****************************************************  
   
   FUNCTION qcheck(klon, klev, paprs, q, ql, aire)  
   
     ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28  
   
     use YOMCST  
     IMPLICIT none  
   
     ! Calculer et imprimer l'eau totale. A utiliser pour verifier  
     ! la conservation de l'eau  
   
     INTEGER klon, klev  
     REAL, intent(in):: paprs(klon, klev+1)  
     real q(klon, klev), ql(klon, klev)  
     REAL aire(klon)  
     REAL qtotal, zx, qcheck  
     INTEGER i, k  
   
     zx = 0.0  
     DO i = 1, klon  
        zx = zx + aire(i)  
     ENDDO  
     qtotal = 0.0  
     DO k = 1, klev  
        DO i = 1, klon  
           qtotal = qtotal + (q(i, k)+ql(i, k)) * aire(i) &  
                *(paprs(i, k)-paprs(i, k+1))/RG  
        ENDDO  
     ENDDO  
   
     qcheck = qtotal/zx  
   
   END FUNCTION qcheck  
   
1887  end module physiq_m  end module physiq_m

Legend:
Removed from v.13  
changed lines
  Added in v.90

  ViewVC Help
Powered by ViewVC 1.1.21