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

Diff of /trunk/phylmd/physiq.f

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21