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

Legend:
Removed from v.12  
changed lines
  Added in v.79

  ViewVC Help
Powered by ViewVC 1.1.21