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

Diff of /trunk/phylmd/physiq.f

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21