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

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

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

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

Legend:
Removed from v.7  
changed lines
  Added in v.191

  ViewVC Help
Powered by ViewVC 1.1.21