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

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

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

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

Legend:
Removed from v.53  
changed lines
  Added in v.200

  ViewVC Help
Powered by ViewVC 1.1.21