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

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

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

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

Legend:
Removed from v.10  
changed lines
  Added in v.180

  ViewVC Help
Powered by ViewVC 1.1.21