source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90 @ 3604

Last change on this file since 3604 was 3604, checked in by oboucher, 5 years ago

addition of flag_volc_surfstrat for VolMIP configurations

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 172.7 KB
Line 
1!
2! $Id: physiq_mod.F90 3604 2019-11-19 22:04:56Z oboucher $
3!
4!#define IO_DEBUG
5MODULE physiq_mod
6
7  IMPLICIT NONE
8
9CONTAINS
10
11  SUBROUTINE physiq (nlon,nlev, &
12       debut,lafin,pdtphys_, &
13       paprs,pplay,pphi,pphis,presnivs, &
14       u,v,rot,t,qx, &
15       flxmass_w, &
16       d_u, d_v, d_t, d_qx, d_ps)
17
18    use assert_m, only: assert
19    USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
20         histwrite, ju2ymds, ymds2ju, getin
21    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
22    USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
23         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
24    USE write_field_phy
25    USE dimphy
26    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
27    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
28    USE mod_phys_lmdz_para
29    USE iophy
30    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
31    USE phystokenc_mod, ONLY: offline, phystokenc
32    USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
33    USE vampir
34    USE pbl_surface_mod, ONLY : pbl_surface
35    USE change_srf_frac_mod
36    USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
37    USE tropopause_m,     ONLY: dyn_tropopause
38#ifdef CPP_Dust
39    USE phytracr_spl_mod, ONLY: phytracr_spl
40#endif
41#ifdef CPP_StratAer
42    USE strataer_mod, ONLY: strataer_init
43#endif
44    USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
45       ! [Variables internes non sauvegardees de la physique]
46       ! Variables locales pour effectuer les appels en serie
47       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
48       ! Dynamic tendencies (diagnostics)
49       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
50       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
51       ! Physic tendencies
52       d_t_con,d_q_con,d_u_con,d_v_con, &
53       d_tr, &                              !! to be removed?? (jyg)
54       d_t_wake,d_q_wake, &
55       d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
56       d_t_ajsb,d_q_ajsb, &
57       d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
58       d_t_ajs_w,d_q_ajs_w, &
59       d_t_ajs_x,d_q_ajs_x, &
60       !
61       d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
62       d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
63       d_t_lscst,d_q_lscst, &
64       d_t_lscth,d_q_lscth, &
65       plul_st,plul_th, &
66       !
67       d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
68       d_t_vdf_w,d_q_vdf_w, &
69       d_t_vdf_x,d_q_vdf_x, &
70       d_ts, &
71       !
72       d_t_oli,d_u_oli,d_v_oli, &
73       d_t_oro,d_u_oro,d_v_oro, &
74       d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
75       d_t_lif,d_u_lif,d_v_lif, &
76       d_t_ec, &
77       !
78       du_gwd_hines,dv_gwd_hines,d_t_hin, &
79       dv_gwd_rando,dv_gwd_front, &
80       east_gwstress,west_gwstress, &
81       d_q_ch4, &
82       !  Special RRTM
83       ZLWFT0_i,ZSWFT0_i,ZFLDN0,  &
84       ZFLUP0,ZFSDN0,ZFSUP0,      &
85       !
86       topswad_aero,solswad_aero,   &
87       topswai_aero,solswai_aero,   &
88       topswad0_aero,solswad0_aero, &
89       !LW additional
90       toplwad_aero,sollwad_aero,   &
91       toplwai_aero,sollwai_aero,   &
92       toplwad0_aero,sollwad0_aero, &
93       !
94       topsw_aero,solsw_aero,       &
95       topsw0_aero,solsw0_aero,     &
96       topswcf_aero,solswcf_aero,   &
97       tausum_aero,tau3d_aero,      &
98       drytausum_aero,              &
99       !
100       !variables CFMIP2/CMIP5
101       topswad_aerop, solswad_aerop,   &
102       topswai_aerop, solswai_aerop,   &
103       topswad0_aerop, solswad0_aerop, &
104       topsw_aerop, topsw0_aerop,      &
105       solsw_aerop, solsw0_aerop,      &
106       topswcf_aerop, solswcf_aerop,   &
107       !LW diagnostics
108       toplwad_aerop, sollwad_aerop,   &
109       toplwai_aerop, sollwai_aerop,   &
110       toplwad0_aerop, sollwad0_aerop, &
111       !
112       ptstar, pt0, slp, &
113       !
114       bils, &
115       !
116       cldh, cldl,cldm, cldq, cldt,      &
117       JrNt,                             &
118       dthmin, evap, fder, plcl, plfc,   &
119       prw, prlw, prsw,                  &
120       s_lcl, s_pblh, s_pblt, s_therm,   &
121       cdragm, cdragh,                   &
122       zustar, zu10m, zv10m, rh2m, qsat2m, &
123       zq2m, zt2m, weak_inversion, &
124       zq2m_cor,zt2m_cor,zu10m_cor,zv10m_cor, & ! pour corriger d'un bug
125       zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
126       t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
127       !
128       s_pblh_x, s_pblh_w, &
129       s_lcl_x, s_lcl_w,   &
130       !
131       slab_wfbils, tpot, tpote,               &
132       ue, uq, ve, vq, zxffonte,               &
133       uwat, vwat,                             &
134       zxfqcalving, zxfluxlat,                 &
135       zxrunofflic,                            &
136       zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
137       rain_lsc, rain_num,                     &
138       !
139       sens_x, sens_w, &
140       zxfluxlat_x, zxfluxlat_w, &
141       !
142       dtvdf_x, dtvdf_w, &
143       dqvdf_x, dqvdf_w, &
144       pbl_tke_input, &
145       t_therm, q_therm, u_therm, v_therm, &
146       cdragh_x, cdragh_w, &
147       cdragm_x, cdragm_w, &
148       kh, kh_x, kh_w, &
149       !
150       wake_k, &
151       alp_wake, &
152       wake_h, wake_omg, &
153                       ! tendencies of delta T and delta q:
154       d_deltat_wk, d_deltaq_wk, &         ! due to wakes
155       d_deltat_wk_gw, d_deltaq_wk_gw, &   ! due to wake induced gravity waves
156       d_deltat_vdf, d_deltaq_vdf, &       ! due to vertical diffusion
157       d_deltat_the, d_deltaq_the, &       ! due to thermals
158       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
159                       ! tendencies of wake fractional area and wake number per unit area:
160       d_s_wk,  d_dens_wk, &             ! due to wakes
161!!!       d_s_vdf, d_dens_vdf, &            ! due to vertical diffusion
162!!!       d_s_the, d_dens_the, &            ! due to thermals
163       !                                 
164       ptconv, ratqsc, &
165       wbeff, convoccur, zmax_th, &
166       sens, flwp, fiwp,  &
167       alp_bl_conv,alp_bl_det,  &
168       alp_bl_fluct_m,alp_bl_fluct_tke,  &
169       alp_bl_stat, n2, s2,  &
170       proba_notrig, random_notrig,  &
171       !
172       dnwd0,  &
173       omega,  &
174       epmax_diag,  &
175
176       !    Deep convective variables used in phytrac
177       pmflxr, pmflxs,  &
178       wdtrainA, wdtrainM,  &
179       upwd, dnwd, &
180       ep,  &
181       da, mp, &
182       phi, &
183       wght_cvfd, &
184       phi2, &
185       d1a, dam, &
186       ev, &
187       elij, &
188       clw, &
189       epmlmMm, eplaMm, &
190       sij, &
191
192       cldemi,  &
193       cldfra, cldtau, fiwc,  &
194       fl, re, flwc,  &
195       ref_liq, ref_ice, theta,  &
196       ref_liq_pi, ref_ice_pi,  &
197       zphi, zx_rh,  &
198       pmfd, pmfu,  &
199       !
200       t2m, fluxlat,  &
201       fsollw, evap_pot,  &
202       fsolsw, wfbils, wfbilo,  &
203       wfevap, wfrain, wfsnow,  & 
204       prfl, psfl, fraca, Vprecip,  &
205       zw2,  &
206       
207       fluxu, fluxv,  &
208       fluxt,  &
209
210       uwriteSTD, vwriteSTD, &                !pour calcul_STDlev.h
211       wwriteSTD, phiwriteSTD, &              !pour calcul_STDlev.h
212       qwriteSTD, twriteSTD, rhwriteSTD, &    !pour calcul_STDlev.h
213       
214       beta_prec,  &
215       rneb,  &
216       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
217       !
218    USE phys_state_var_mod ! Variables sauvegardees de la physique
219#ifdef CPP_Dust
220  USE phys_output_write_spl_mod
221#else
222    USE phys_output_var_mod ! Variables pour les ecritures des sorties
223#endif
224
225    USE phys_output_write_mod
226    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
227    USE phys_output_mod
228    USE phys_output_ctrlout_mod
229    use open_climoz_m, only: open_climoz ! ozone climatology from a file
230    use regr_pr_time_av_m, only: regr_pr_time_av
231    use netcdf95, only: nf95_close
232    !IM for NMC files
233    !     use netcdf, only: nf90_fill_real
234    use netcdf, only: nf90_fill_real
235    use mod_phys_lmdz_mpi_data, only: is_mpi_root
236    USE aero_mod
237    use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
238    use conf_phys_m, only: conf_phys
239    use radlwsw_m, only: radlwsw
240    use phyaqua_mod, only: zenang_an
241    USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
242         start_time, pdtphys, day_ini
243    USE tracinca_mod, ONLY: config_inca
244#ifdef CPP_XIOS
245    USE wxios, ONLY: missing_val, missing_val_omp
246    USE xios, ONLY: xios_get_field_attr, xios_field_is_active
247#endif
248#ifdef REPROBUS
249    USE CHEM_REP, ONLY : Init_chem_rep_xjour
250#endif
251    USE indice_sol_mod
252    USE phytrac_mod, ONLY : phytrac
253
254#ifdef CPP_RRTM
255    USE YOERAD, ONLY : NRADLP
256    USE YOESW, ONLY : RSUN
257#endif
258    USE ioipsl_getin_p_mod, ONLY : getin_p
259
260#ifndef CPP_XIOS
261    USE paramLMDZ_phy_mod
262#endif
263
264    USE cmp_seri_mod
265    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
266  &      fl_ebil, fl_cor_ebil
267
268    !IM stations CFMIP
269    USE CFMIP_point_locations
270    use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
271    use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
272    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps
273
274
275    IMPLICIT none
276    !>======================================================================
277    !!
278    !! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
279    !!
280    !! Objet: Moniteur general de la physique du modele
281    !!AA      Modifications quant aux traceurs :
282    !!AA                  -  uniformisation des parametrisations ds phytrac
283    !!AA                  -  stockage des moyennes des champs necessaires
284    !!AA                     en mode traceur off-line
285    !!======================================================================
286    !!   CLEFS CPP POUR LES IO
287    !!   =====================
288#define histNMC
289    !!======================================================================
290    !!    modif   ( P. Le Van ,  12/10/98 )
291    !!
292    !!  Arguments:
293    !!
294    !! nlon----input-I-nombre de points horizontaux
295    !! nlev----input-I-nombre de couches verticales, doit etre egale a klev
296    !! debut---input-L-variable logique indiquant le premier passage
297    !! lafin---input-L-variable logique indiquant le dernier passage
298    !! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
299    !! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
300    !! pdtphys-input-R-pas d'integration pour la physique (seconde)
301    !! paprs---input-R-pression pour chaque inter-couche (en Pa)
302    !! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
303    !! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
304    !! pphis---input-R-geopotentiel du sol
305    !! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
306    !! u-------input-R-vitesse dans la direction X (de O a E) en m/s
307    !! v-------input-R-vitesse Y (de S a N) en m/s
308    !! t-------input-R-temperature (K)
309    !! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
310    !! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
311    !! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
312    !! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
313    !! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
314    !! flxmass_w -input-R- flux de masse verticale
315    !! d_u-----output-R-tendance physique de "u" (m/s/s)
316    !! d_v-----output-R-tendance physique de "v" (m/s/s)
317    !! d_t-----output-R-tendance physique de "t" (K/s)
318    !! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
319    !! d_ps----output-R-tendance physique de la pression au sol
320    !!======================================================================
321    integer jjmp1
322    !  parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
323    !  integer iip1
324    !  parameter (iip1=iim+1)
325
326    include "regdim.h"
327    include "dimsoil.h"
328    include "clesphys.h"
329    include "thermcell.h"
330    include "dimpft.h"
331    !======================================================================
332    LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
333    !$OMP THREADPRIVATE(ok_volcan)
334    INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato
335    !$OMP THREADPRIVATE(flag_volc_surfstrat)
336    LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
337    PARAMETER (ok_cvl=.TRUE.)
338    LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
339    PARAMETER (ok_gust=.FALSE.)
340    INTEGER iflag_radia     ! active ou non le rayonnement (MPL)
341    SAVE iflag_radia
342    !$OMP THREADPRIVATE(iflag_radia)
343    !======================================================================
344    LOGICAL check ! Verifier la conservation du modele en eau
345    PARAMETER (check=.FALSE.)
346    LOGICAL ok_stratus ! Ajouter artificiellement les stratus
347    PARAMETER (ok_stratus=.FALSE.)
348    !======================================================================
349    REAL amn, amx
350    INTEGER igout
351    !======================================================================
352    ! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
353    ! en attente du codage des cles par Fred
354    ! iflag_cycle_diurne est initialise par conf_phys et se trouve
355    ! dans clesphys.h (IM)
356    !======================================================================
357    ! Modele thermique du sol, a activer pour le cycle diurne:
358    !cc      LOGICAL soil_model
359    !cc      PARAMETER (soil_model=.FALSE.)
360    !======================================================================
361    ! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
362    ! le calcul du rayonnement est celle apres la precipitation des nuages.
363    ! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
364    ! la condensation et la precipitation. Cette cle augmente les impacts
365    ! radiatifs des nuages.
366    !cc      LOGICAL new_oliq
367    !cc      PARAMETER (new_oliq=.FALSE.)
368    !======================================================================
369    ! Clefs controlant deux parametrisations de l'orographie:
370    !c      LOGICAL ok_orodr
371    !cc      PARAMETER (ok_orodr=.FALSE.)
372    !cc      LOGICAL ok_orolf
373    !cc      PARAMETER (ok_orolf=.FALSE.)
374    !======================================================================
375    LOGICAL ok_journe ! sortir le fichier journalier
376    save ok_journe
377    !$OMP THREADPRIVATE(ok_journe)
378    !
379    LOGICAL ok_mensuel ! sortir le fichier mensuel
380    save ok_mensuel
381    !$OMP THREADPRIVATE(ok_mensuel)
382    !
383    LOGICAL ok_instan ! sortir le fichier instantane
384    save ok_instan
385    !$OMP THREADPRIVATE(ok_instan)
386    !
387    LOGICAL ok_LES ! sortir le fichier LES
388    save ok_LES                           
389    !$OMP THREADPRIVATE(ok_LES)                 
390    !
391    LOGICAL callstats ! sortir le fichier stats
392    save callstats                           
393    !$OMP THREADPRIVATE(callstats)                 
394    !
395    LOGICAL ok_region ! sortir le fichier regional
396    PARAMETER (ok_region=.FALSE.)
397    !======================================================================
398    real seuil_inversion
399    save seuil_inversion
400    !$OMP THREADPRIVATE(seuil_inversion)
401    integer iflag_ratqs
402    save iflag_ratqs
403    !$OMP THREADPRIVATE(iflag_ratqs)
404    real facteur
405
406    REAL wmax_th(klon)
407    REAL tau_overturning_th(klon)
408
409    integer lmax_th(klon)
410    integer limbas(klon)
411    real ratqscth(klon,klev)
412    real ratqsdiff(klon,klev)
413    real zqsatth(klon,klev)
414
415    !======================================================================
416    !
417    INTEGER ivap          ! indice de traceurs pour vapeur d'eau
418    PARAMETER (ivap=1)
419    INTEGER iliq          ! indice de traceurs pour eau liquide
420    PARAMETER (iliq=2)
421    !CR: on ajoute la phase glace
422    INTEGER isol          ! indice de traceurs pour eau glace
423    PARAMETER (isol=3)
424    !
425    !
426    ! Variables argument:
427    !
428    INTEGER nlon
429    INTEGER nlev
430    REAL,INTENT(IN) :: pdtphys_
431    ! NB: pdtphys to be used in physics is in time_phylmdz_mod
432    LOGICAL debut, lafin
433    REAL paprs(klon,klev+1)
434    REAL pplay(klon,klev)
435    REAL pphi(klon,klev)
436    REAL pphis(klon)
437    REAL presnivs(klev)
438!JLD    REAL znivsig(klev)
439!JLD    real pir
440
441    REAL u(klon,klev)
442    REAL v(klon,klev)
443
444    REAL, intent(in):: rot(klon, klev)
445    ! relative vorticity, in s-1, needed for frontal waves
446
447    REAL t(klon,klev),thetal(klon,klev)
448    ! thetal: ligne suivante a decommenter si vous avez les fichiers
449    !     MPL 20130625
450    ! fth_fonctions.F90 et parkind1.F90
451    ! sinon thetal=theta
452    !     REAL fth_thetae,fth_thetav,fth_thetal
453    REAL qx(klon,klev,nqtot)
454    REAL flxmass_w(klon,klev)
455    REAL d_u(klon,klev)
456    REAL d_v(klon,klev)
457    REAL d_t(klon,klev)
458    REAL d_qx(klon,klev,nqtot)
459    REAL d_ps(klon)
460  ! variables pour tend_to_tke
461    REAL duadd(klon,klev)
462    REAL dvadd(klon,klev)
463    REAL dtadd(klon,klev)
464
465#ifndef CPP_XIOS
466    REAL, SAVE :: missing_val=nf90_fill_real
467#endif
468!!   Variables moved to phys_local_var_mod
469!!    ! Variables pour le transport convectif
470!!    real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
471!!    real wght_cvfd(klon,klev)
472!!    ! Variables pour le lessivage convectif
473!!    ! RomP >>>
474!!    real phi2(klon,klev,klev)
475!!    real d1a(klon,klev),dam(klon,klev)
476!!    real ev(klon,klev)
477!!    real clw(klon,klev),elij(klon,klev,klev)
478!!    real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
479!!    ! RomP <<<
480    !IM definition dynamique o_trac dans phys_output_open
481    !      type(ctrl_out) :: o_trac(nqtot)
482
483    ! variables a une pression donnee
484    !
485    include "declare_STDlev.h"
486    !
487    !
488    include "radopt.h"
489    !
490    !
491    INTEGER debug
492    INTEGER n
493    !ym      INTEGER npoints
494    !ym      PARAMETER(npoints=klon)
495    !
496    INTEGER nregISCtot
497    PARAMETER(nregISCtot=1)
498    !
499    ! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
500    ! sur 1 region rectangulaire y compris pour 1 point
501    ! imin_debut : indice minimum de i; nbpti : nombre de points en
502    ! direction i (longitude)
503    ! jmin_debut : indice minimum de j; nbptj : nombre de points en
504    ! direction j (latitude)
505!JLD    INTEGER imin_debut, nbpti
506!JLD    INTEGER jmin_debut, nbptj
507    !IM: region='3d' <==> sorties en global
508    CHARACTER*3 region
509    PARAMETER(region='3d')
510    logical ok_hf
511    !
512    save ok_hf
513    !$OMP THREADPRIVATE(ok_hf)
514
515    INTEGER,PARAMETER :: longcles=20
516    REAL,SAVE :: clesphy0(longcles)
517    !$OMP THREADPRIVATE(clesphy0)
518    !
519    ! Variables propres a la physique
520    INTEGER itap
521    SAVE itap                   ! compteur pour la physique
522    !$OMP THREADPRIVATE(itap)
523
524    INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
525    !$OMP THREADPRIVATE(abortphy)
526    !
527    REAL,save ::  solarlong0
528    !$OMP THREADPRIVATE(solarlong0)
529
530    !
531    !  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
532    !
533    !IM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
534    REAL zulow(klon),zvlow(klon)
535    !
536    INTEGER igwd,idx(klon),itest(klon)
537    !
538    !      REAL,allocatable,save :: run_off_lic_0(:)
539    ! !$OMP THREADPRIVATE(run_off_lic_0)
540    !ym      SAVE run_off_lic_0
541    !KE43
542    ! Variables liees a la convection de K. Emanuel (sb):
543    !
544    REAL bas, top             ! cloud base and top levels
545    SAVE bas
546    SAVE top
547    !$OMP THREADPRIVATE(bas, top)
548    !------------------------------------------------------------------
549    ! Upmost level reached by deep convection and related variable (jyg)
550    !
551    INTEGER izero
552    INTEGER k_upper_cv
553    !------------------------------------------------------------------
554    ! Compteur de l'occurence de cvpas=1
555    INTEGER Ncvpaseq1
556    SAVE Ncvpaseq1
557    !$OMP THREADPRIVATE(Ncvpaseq1)
558    !
559    !==========================================================================
560    !CR04.12.07: on ajoute les nouvelles variables du nouveau schema
561    !de convection avec poches froides
562    ! Variables li\'ees \`a la poche froide (jyg)
563
564    REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
565    !
566    REAL wape_prescr, fip_prescr
567    INTEGER it_wape_prescr
568    SAVE wape_prescr, fip_prescr, it_wape_prescr
569    !$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
570    !
571    ! variables supplementaires de concvl
572    REAL Tconv(klon,klev)
573!!    variable moved to phys_local_var_mod
574!!    REAL sij(klon,klev,klev)
575!!    !
576!!    ! variables pour tester la conservation de l'energie dans concvl
577!!    REAL, DIMENSION(klon,klev)     :: d_t_con_sat
578!!    REAL, DIMENSION(klon,klev)     :: d_q_con_sat
579!!    REAL, DIMENSION(klon,klev)     :: dql_sat
580
581    real, save :: alp_bl_prescr=0.
582    real, save :: ale_bl_prescr=0.
583
584    real, save :: wake_s_min_lsp=0.1
585
586    !$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
587    !$OMP THREADPRIVATE(wake_s_min_lsp)
588
589
590    real ok_wk_lsp(klon)
591
592    !RC
593    ! Variables li\'ees \`a la poche froide (jyg et rr)
594
595    INTEGER,  SAVE               :: iflag_wake_tend  ! wake: if =0, then wake state variables are
596                                                     ! updated within calwake
597    !$OMP THREADPRIVATE(iflag_wake_tend)
598    INTEGER,  SAVE               :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
599                                                        ! power provided by the wakes; else, Alp_wk is the
600                                                        ! lifting power conditionned on the presence of a
601                                                        ! gust-front in the grid cell.
602    !$OMP THREADPRIVATE(iflag_alp_wk_cond)
603    INTEGER,  SAVE               :: iflag_bug_t2m_ipslcm61=0 !
604    !$OMP THREADPRIVATE(iflag_bug_t2m_ipslcm61)
605
606    REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
607    REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
608
609    REAL wake_dth(klon,klev)        ! wake : temp pot difference
610
611    REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta
612    ! transported by LS omega
613    REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of
614    ! large scale omega
615    REAL wake_dtKE(klon,klev)       ! Wake : differential heating
616    ! (wake - unpertubed) CONV
617    REAL wake_dqKE(klon,klev)       ! Wake : differential moistening
618    ! (wake - unpertubed) CONV
619    REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
620    REAL wake_spread(klon,klev)     ! spreading term in wake_delt
621    !
622    !pourquoi y'a pas de save??
623    !
624!!!    INTEGER, SAVE, DIMENSION(klon)   :: wake_k
625!!!    !$OMP THREADPRIVATE(wake_k)
626    !
627    !jyg<
628    !cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
629    !>jyg
630
631    REAL wake_fip_0(klon)           ! Average Front Incoming Power (unconditionned)
632    REAL wake_gfl(klon)             ! Gust Front Length
633!!!    REAL wake_dens(klon)         ! moved to phys_state_var_mod
634    !
635    !
636    REAL dt_dwn(klon,klev)
637    REAL dq_dwn(klon,klev)
638    REAL M_dwn(klon,klev)
639    REAL M_up(klon,klev)
640    REAL dt_a(klon,klev)
641    REAL dq_a(klon,klev)
642    REAL d_t_adjwk(klon,klev)                !jyg
643    REAL d_q_adjwk(klon,klev)                !jyg
644    LOGICAL,SAVE :: ok_adjwk=.FALSE.
645    !$OMP THREADPRIVATE(ok_adjwk)
646    INTEGER,SAVE :: iflag_adjwk=0            !jyg
647    !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
648    REAL,SAVE :: oliqmax=999.,oicemax=999.
649    !$OMP THREADPRIVATE(oliqmax,oicemax)
650    REAL, SAVE :: alp_offset
651    !$OMP THREADPRIVATE(alp_offset)
652    REAL, SAVE :: dtcon_multistep_max=1.e6
653    !$OMP THREADPRIVATE(dtcon_multistep_max)
654    REAL, SAVE :: dqcon_multistep_max=1.e6
655    !$OMP THREADPRIVATE(dqcon_multistep_max)
656
657 
658    !
659    !RR:fin declarations poches froides
660    !==========================================================================
661
662    REAL ztv(klon,klev),ztva(klon,klev)
663    REAL zpspsk(klon,klev)
664    REAL ztla(klon,klev),zqla(klon,klev)
665    REAL zthl(klon,klev)
666
667    !cc nrlmd le 10/04/2012
668
669    !--------Stochastic Boundary Layer Triggering: ALE_BL--------
670    !---Propri\'et\'es du thermiques au LCL
671    real zlcl_th(klon)          ! Altitude du LCL calcul\'e
672    ! continument (pcon dans
673    ! thermcell_main.F90)
674    real fraca0(klon)           ! Fraction des thermiques au LCL
675    real w0(klon)               ! Vitesse des thermiques au LCL
676    real w_conv(klon)           ! Vitesse verticale de grande \'echelle au LCL
677    real tke0(klon,klev+1)      ! TKE au d\'ebut du pas de temps
678    real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL
679    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL
680
681!JLD    !---D\'eclenchement stochastique
682!JLD    integer :: tau_trig(klon)
683
684    REAL,SAVE :: random_notrig_max=1.
685    !$OMP THREADPRIVATE(random_notrig_max)
686
687    !--------Statistical Boundary Layer Closure: ALP_BL--------
688    !---Profils de TKE dans et hors du thermique
689    real therm_tke_max(klon,klev)   ! Profil de TKE dans les thermiques
690    real env_tke_max(klon,klev)     ! Profil de TKE dans l'environnement
691
692    !-------Activer les tendances de TKE due a l'orograp??ie---------
693     INTEGER, SAVE :: addtkeoro
694    !$OMP THREADPRIVATE(addtkeoro)
695     REAL, SAVE :: alphatkeoro
696    !$OMP THREADPRIVATE(alphatkeoro)
697     LOGICAL, SAVE :: smallscales_tkeoro
698    !$OMP THREADPRIVATE(smallscales_tkeoro)
699
700
701
702    !cc fin nrlmd le 10/04/2012
703
704    ! Variables locales pour la couche limite (al1):
705    !
706    !Al1      REAL pblh(klon)           ! Hauteur de couche limite
707    !Al1      SAVE pblh
708    !34EK
709    !
710    ! Variables locales:
711    !
712    !AA
713    !AA  Pour phytrac
714    REAL u1(klon)             ! vents dans la premiere couche U
715    REAL v1(klon)             ! vents dans la premiere couche V
716
717    !@$$      LOGICAL offline           ! Controle du stockage ds "physique"
718    !@$$      PARAMETER (offline=.false.)
719    !@$$      INTEGER physid
720    REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
721    REAL frac_nucl(klon,klev) ! idem (nucleation)
722    ! RomP >>>
723    REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
724    ! RomP <<<
725    REAL          :: calday
726
727    !IM cf FH pour Tiedtke 080604
728    REAL rain_tiedtke(klon),snow_tiedtke(klon)
729    !
730    !IM 050204 END
731    REAL devap(klon) ! evaporation et sa derivee
732    REAL dsens(klon) ! chaleur sensible et sa derivee
733
734    !
735    ! Conditions aux limites
736    !
737    !
738    REAL :: day_since_equinox
739    ! Date de l'equinoxe de printemps
740    INTEGER, parameter :: mth_eq=3, day_eq=21
741    REAL :: jD_eq
742
743    LOGICAL, parameter :: new_orbit = .true.
744
745    !
746    INTEGER lmt_pas
747    SAVE lmt_pas                ! frequence de mise a jour
748    !$OMP THREADPRIVATE(lmt_pas)
749    real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
750    !     (column-density of mass of air in a cell, in kg m-2)
751    real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
752
753    !IM sorties
754    REAL un_jour
755    PARAMETER(un_jour=86400.)
756    INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
757    SAVE itapm1    !mis a jour le dernier pas de temps du mois en cours
758    !$OMP THREADPRIVATE(itapm1)
759    !======================================================================
760    !
761    ! Declaration des procedures appelees
762    !
763    EXTERNAL angle     ! calculer angle zenithal du soleil
764    EXTERNAL alboc     ! calculer l'albedo sur ocean
765    EXTERNAL ajsec     ! ajustement sec
766    EXTERNAL conlmd    ! convection (schema LMD)
767    !KE43
768    EXTERNAL conema3  ! convect4.3
769    EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
770    !AA
771    ! JBM (3/14) fisrtilp_tr not loaded
772    ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
773    !                          ! stockage des coefficients necessaires au
774    !                          ! lessivage OFF-LINE et ON-LINE
775    EXTERNAL hgardfou  ! verifier les temperatures
776    EXTERNAL nuage     ! calculer les proprietes radiatives
777    !C      EXTERNAL o3cm      ! initialiser l'ozone
778    EXTERNAL orbite    ! calculer l'orbite terrestre
779    EXTERNAL phyetat0  ! lire l'etat initial de la physique
780    EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
781    EXTERNAL suphel    ! initialiser certaines constantes
782    EXTERNAL transp    ! transport total de l'eau et de l'energie
783    !IM
784    EXTERNAL haut2bas  !variables de haut en bas
785    EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
786    EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
787    !     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
788    ! EXTERNAL moyglo_aire
789    ! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
790    ! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
791    !
792    !
793    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
794    ! Local variables
795    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796    !
797    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
798    REAL dialiq(klon,klev)  ! eau liquide nuageuse
799    REAL diafra(klon,klev)  ! fraction nuageuse
800    REAL cldliq(klon,klev)  ! eau liquide nuageuse
801    !
802    !XXX PB
803    REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
804    !
805    REAL zxfluxt(klon, klev)
806    REAL zxfluxq(klon, klev)
807    REAL zxfluxu(klon, klev)
808    REAL zxfluxv(klon, klev)
809
810    ! Le rayonnement n'est pas calcule tous les pas, il faut donc
811    !                      sauvegarder les sorties du rayonnement
812    !ym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
813    !ym      SAVE  sollwdownclr, toplwdown, toplwdownclr
814    !ym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
815    !
816    INTEGER itaprad
817    SAVE itaprad
818    !$OMP THREADPRIVATE(itaprad)
819    !
820    REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
821    REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
822    !
823#ifdef INCA
824    REAL zxsnow_dummy(klon)
825#endif
826    REAL zsav_tsol(klon)
827    !
828    REAL dist, rmu0(klon), fract(klon)
829    REAL zrmu0(klon), zfract(klon)
830    REAL zdtime, zdtime1, zdtime2, zlongi
831    !
832    REAL qcheck
833    REAL z_avant(klon), z_apres(klon), z_factor(klon)
834    LOGICAL zx_ajustq
835    !
836    REAL za
837    REAL zx_t, zx_qs, zdelta, zcor
838    real zqsat(klon,klev)
839    !
840    INTEGER i, k, iq, j, nsrf, ll, l
841    !
842    REAL t_coup
843    PARAMETER (t_coup=234.0)
844
845    !ym A voir plus tard !!
846    !ym      REAL zx_relief(iim,jjmp1)
847    !ym      REAL zx_aire(iim,jjmp1)
848    !
849    ! Grandeurs de sorties
850    REAL s_capCL(klon)
851    REAL s_oliqCL(klon), s_cteiCL(klon)
852    REAL s_trmb1(klon), s_trmb2(klon)
853    REAL s_trmb3(klon)
854
855    ! La convection n'est pas calculee tous les pas, il faut donc
856    !                      sauvegarder les sorties de la convection
857    !ym      SAVE 
858    !ym      SAVE 
859    !ym      SAVE 
860    !
861    INTEGER itapcv, itapwk
862    SAVE itapcv, itapwk
863    !$OMP THREADPRIVATE(itapcv, itapwk)
864
865    !KE43
866    ! Variables locales pour la convection de K. Emanuel (sb):
867
868    REAL tvp(klon,klev)       ! virtual temp of lifted parcel
869    CHARACTER*40 capemaxcels  !max(CAPE)
870
871    REAL rflag(klon)          ! flag fonctionnement de convect
872    INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
873
874    ! -- convect43:
875    INTEGER ntra              ! nb traceurs pour convect4.3
876    REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
877    REAL dplcldt(klon), dplcldr(klon)
878    !?     .     condm_con(klon,klev),conda_con(klon,klev),
879    !?     .     mr_con(klon,klev),ep_con(klon,klev)
880    !?     .    ,sadiab(klon,klev),wadiab(klon,klev)
881    ! --
882    !34EK
883    !
884    ! Variables du changement
885    !
886    ! con: convection
887    ! lsc: condensation a grande echelle (Large-Scale-Condensation)
888    ! ajs: ajustement sec
889    ! eva: evaporation de l'eau liquide nuageuse
890    ! vdf: couche limite (Vertical DiFfusion)
891    !
892    ! tendance nulles
893    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
894    REAL, dimension(klon)     :: dsig0, ddens0
895    INTEGER, dimension(klon)  :: wkoccur1
896    ! tendance buffer pour appel de add_phys_tend
897    REAL, DIMENSION(klon,klev)  :: d_q_ch4_dtime
898    !
899    ! Flag pour pouvoir ne pas ajouter les tendances.
900    ! Par defaut, les tendances doivente etre ajoutees et
901    ! flag_inhib_tend = 0
902    ! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
903    ! croissant de print quand la valeur du flag augmente
904    !!! attention, ce flag doit etre change avec prudence !!!
905    INTEGER :: flag_inhib_tend = 0 !  0 is the default value
906!!    INTEGER :: flag_inhib_tend = 2
907    !
908    ! Logical switch to a bug : reseting to 0 convective variables at the
909    ! begining of physiq.
910    LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
911    !$OMP THREADPRIVATE(ok_bug_cv_trac)
912
913    !
914    !********************************************************
915    !     declarations
916
917    !********************************************************
918    !IM 081204 END
919    !
920    REAL pen_u(klon,klev), pen_d(klon,klev)
921    REAL pde_u(klon,klev), pde_d(klon,klev)
922    INTEGER kcbot(klon), kctop(klon), kdtop(klon)
923    !
924    real ratqsbas,ratqshaut,tau_ratqs
925    save ratqsbas,ratqshaut,tau_ratqs
926    !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
927    REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
928    !$OMP THREADPRIVATE(ratqsp0, ratqsdp)
929
930    ! Parametres lies au nouveau schema de nuages (SB, PDF)
931    real fact_cldcon
932    real facttemps
933    logical ok_newmicro
934    save ok_newmicro
935    !$OMP THREADPRIVATE(ok_newmicro)
936    !real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
937    save fact_cldcon,facttemps
938    !$OMP THREADPRIVATE(fact_cldcon,facttemps)
939
940    integer iflag_cld_th
941    save iflag_cld_th
942    !$OMP THREADPRIVATE(iflag_cld_th)
943!IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
944    !IM cf. AM 081204 BEG
945    logical ptconvth(klon,klev)
946    !IM cf. AM 081204 END
947    !
948    ! Variables liees a l'ecriture de la bande histoire physique
949    !
950    !======================================================================
951    !
952
953    !
954!JLD    integer itau_w   ! pas de temps ecriture = itap + itau_phy
955    !
956    !
957    ! Variables locales pour effectuer les appels en serie
958    !
959    !IM RH a 2m (la surface)
960    REAL Lheat
961
962    INTEGER        length
963    PARAMETER    ( length = 100 )
964    REAL tabcntr0( length       )
965    !
966!JLD    INTEGER ndex2d(nbp_lon*nbp_lat)
967    !IM
968    !
969    !IM AMIP2 BEG
970!JLD    REAL moyglo, mountor
971    !IM 141004 BEG
972    REAL zustrdr(klon), zvstrdr(klon)
973    REAL zustrli(klon), zvstrli(klon)
974    REAL zustrph(klon), zvstrph(klon)
975    REAL aam, torsfc
976    !IM 141004 END
977    !IM 190504 BEG
978    !  INTEGER imp1jmp1
979    !  PARAMETER(imp1jmp1=(iim+1)*jjmp1)
980    !ym A voir plus tard
981    !  REAL zx_tmp((nbp_lon+1)*nbp_lat)
982    !  REAL airedyn(nbp_lon+1,nbp_lat)
983    !IM 190504 END
984!JLD    LOGICAL ok_msk
985!JLD    REAL msk(klon)
986    !ym A voir plus tard
987    !ym      REAL zm_wo(jjmp1, klev)
988    !IM AMIP2 END
989    !
990    REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
991    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
992!JLD    REAL zx_tmp_2d(nbp_lon,nbp_lat)
993!JLD    REAL zx_lon(nbp_lon,nbp_lat)
994!JLD    REAL zx_lat(nbp_lon,nbp_lat)
995    !
996    INTEGER nid_ctesGCM
997    SAVE nid_ctesGCM
998    !$OMP THREADPRIVATE(nid_ctesGCM)
999    !
1000    !IM 280405 BEG
1001    !  INTEGER nid_bilKPins, nid_bilKPave
1002    !  SAVE nid_bilKPins, nid_bilKPave
1003    !  !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1004    !
1005    REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1006    REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1007    REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1008    REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1009    !
1010!JLD    REAL zjulian
1011!JLD    SAVE zjulian
1012!JLD!$OMP THREADPRIVATE(zjulian)
1013
1014!JLD    INTEGER nhori, nvert
1015!JLD    REAL zsto
1016!JLD    REAL zstophy, zout
1017
1018    character*20 modname
1019    character*80 abort_message
1020    logical, save ::  ok_sync, ok_sync_omp
1021    !$OMP THREADPRIVATE(ok_sync)
1022    real date0
1023
1024    ! essai writephys
1025    integer fid_day, fid_mth, fid_ins
1026    parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1027    integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1028    parameter (prof2d_on = 1, prof3d_on = 2, &
1029         prof2d_av = 3, prof3d_av = 4)
1030    REAL ztsol(klon)
1031    REAL q2m(klon,nbsrf)  ! humidite a 2m
1032
1033    !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1034    CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
1035    CHARACTER*40 tinst, tave
1036    REAL cldtaupi(klon,klev) ! Cloud optical thickness for
1037    ! pre-industrial (pi) aerosols
1038
1039    INTEGER :: naero
1040    ! Aerosol optical properties
1041    CHARACTER*4, DIMENSION(naero_grp) :: rfname
1042    REAL, DIMENSION(klon,klev)     :: mass_solu_aero ! total mass
1043    ! concentration
1044    ! for all soluble
1045    ! aerosols[ug/m3]
1046    REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi
1047    ! - " - (pre-industrial value)
1048
1049    ! Parameters
1050    LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
1051    LOGICAL ok_alw            ! Apply aerosol LW effect or not
1052    LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1053    REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
1054    SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1055    !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
1056    LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
1057    ! false : lecture des aerosol dans un fichier
1058    !$OMP THREADPRIVATE(aerosol_couple)   
1059    LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3
1060    ! false : use offline chemistry O3
1061    !$OMP THREADPRIVATE(chemistry_couple)   
1062    INTEGER, SAVE :: flag_aerosol
1063    !$OMP THREADPRIVATE(flag_aerosol)
1064    LOGICAL, SAVE :: flag_bc_internal_mixture
1065    !$OMP THREADPRIVATE(flag_bc_internal_mixture)
1066    LOGICAL, SAVE :: new_aod
1067    !$OMP THREADPRIVATE(new_aod)
1068    !
1069    !--STRAT AEROSOL
1070    INTEGER, SAVE :: flag_aerosol_strat
1071    !$OMP THREADPRIVATE(flag_aerosol_strat)
1072    !
1073    !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1074    LOGICAL, SAVE :: flag_aer_feedback
1075    !$OMP THREADPRIVATE(flag_aer_feedback)
1076
1077    !c-fin STRAT AEROSOL
1078    !
1079    ! Declaration des constantes et des fonctions thermodynamiques
1080    !
1081    LOGICAL,SAVE :: first=.true.
1082    !$OMP THREADPRIVATE(first)
1083
1084    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1085    ! Note that pressure vectors are in Pa and in stricly ascending order
1086    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
1087    !     (let it keep the default OpenMP shared attribute)
1088    !     Allowed values are 0, 1 and 2
1089    !     0: do not read an ozone climatology
1090    !     1: read a single ozone climatology that will be used day and night
1091    !     2: read two ozone climatologies, the average day and night
1092    !     climatology and the daylight climatology
1093    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
1094    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1095    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1096    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
1097    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1098                                  = ["tro3         ","tro3_daylight"]
1099    ! vars_climoz(1:read_climoz): variables names in climoz file.
1100    ! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1101    REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1102                 ! the ozone fields, old method.
1103
1104    include "YOMCST.h"
1105    include "YOETHF.h"
1106    include "FCTTRE.h"
1107    !IM 100106 BEG : pouvoir sortir les ctes de la physique
1108    include "conema3.h"
1109    include "fisrtilp.h"
1110    include "nuage.h"
1111    include "compbl.h"
1112    !IM 100106 END : pouvoir sortir les ctes de la physique
1113    !
1114    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1115    ! Declarations pour Simulateur COSP
1116    !============================================================
1117#ifdef CPP_COSP
1118    real :: mr_ozone(klon,klev)
1119#endif
1120    !IM stations CFMIP
1121    INTEGER, SAVE :: nCFMIP
1122    !$OMP THREADPRIVATE(nCFMIP)
1123    INTEGER, PARAMETER :: npCFMIP=120
1124    INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1125    REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1126    !$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1127    INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1128    REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1129    !$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1130    INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1131    !$OMP THREADPRIVATE(iGCM, jGCM)
1132    logical, dimension(nfiles)            :: phys_out_filestations
1133    logical, parameter :: lNMC=.FALSE.
1134
1135    !IM betaCRF
1136    REAL, SAVE :: pfree, beta_pbl, beta_free
1137    !$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1138    REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
1139    !$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
1140    LOGICAL, SAVE :: mskocean_beta
1141    !$OMP THREADPRIVATE(mskocean_beta)
1142    REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1143    ! cldemirad pour evaluer les
1144    ! retros liees aux CRF
1145    REAL, dimension(klon, klev) :: cldtaurad   ! epaisseur optique
1146    ! pour radlwsw pour
1147    ! tester "CRF off"
1148    REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1149    ! pour radlwsw pour
1150    ! tester "CRF off"
1151    REAL, dimension(klon, klev) :: cldemirad   ! emissivite pour
1152    ! radlwsw pour tester
1153    ! "CRF off"
1154    REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
1155
1156#ifdef INCA
1157    ! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1158    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1159    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1160    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1161    REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1162#endif
1163    REAL, DIMENSION(klon,nbtr) :: init_source
1164
1165
1166    !lwoff=y : offset LW CRE for radiation code and other schemes
1167    REAL, SAVE :: betalwoff
1168    !OMP THREADPRIVATE(betalwoff)
1169!
1170    INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1171    REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1172    REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1173    integer iostat
1174
1175    REAL zzz
1176    !albedo SB >>>
1177    real,dimension(6),save :: SFRWL
1178    !albedo SB <<<
1179
1180    !--OB variables for mass fixer (hard coded for now)
1181    logical, parameter :: mass_fixer=.false.
1182    real qql1(klon),qql2(klon),corrqql
1183
1184    REAL pi
1185
1186    pi = 4. * ATAN(1.)
1187
1188    ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1189    jjmp1=nbp_lat
1190
1191    !======================================================================
1192    ! Gestion calendrier : mise a jour du module phys_cal_mod
1193    !
1194    pdtphys=pdtphys_
1195    CALL update_time(pdtphys)
1196
1197    !======================================================================
1198    ! Ecriture eventuelle d'un profil verticale en entree de la physique.
1199    ! Utilise notamment en 1D mais peut etre active egalement en 3D
1200    ! en imposant la valeur de igout.
1201    !======================================================================d
1202    IF (prt_level.ge.1) THEN
1203       igout=klon/2+1/klon
1204       write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1205       write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1206            longitude_deg(igout)
1207       write(lunout,*) &
1208            'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1209       write(lunout,*) &
1210            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1211
1212       write(lunout,*) 'paprs, play, phi, u, v, t'
1213       DO k=1,klev
1214          write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1215               u(igout,k),v(igout,k),t(igout,k)
1216       ENDDO
1217       write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
1218       DO k=1,klev
1219          write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
1220       ENDDO
1221    ENDIF
1222
1223    ! Quick check on pressure levels:
1224    call assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1225            "physiq_mod paprs bad order")
1226
1227    IF (first) THEN
1228       !CR:nvelles variables convection/poches froides
1229
1230       print*, '================================================='
1231       print*, 'Allocation des variables locales et sauvegardees'
1232       CALL phys_local_var_init
1233       !
1234       pasphys=pdtphys
1235       !     appel a la lecture du run.def physique
1236       CALL conf_phys(ok_journe, ok_mensuel, &
1237            ok_instan, ok_hf, &
1238            ok_LES, &
1239            callstats, &
1240            solarlong0,seuil_inversion, &
1241            fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1242            iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1243            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
1244            chemistry_couple, flag_aerosol, flag_aerosol_strat,         &
1245            flag_aer_feedback, new_aod, &
1246            flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1247                                ! nv flags pour la convection et les
1248                                ! poches froides
1249            read_climoz, &
1250            alp_offset)
1251       CALL phys_state_var_init(read_climoz)
1252       CALL phys_output_var_init
1253
1254#ifdef CPP_StratAer
1255       CALL strataer_init
1256#endif
1257       print*, '================================================='
1258       !
1259       !CR: check sur le nb de traceurs de l eau
1260       IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1261          WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1262               '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
1263          STOP
1264       ENDIF
1265
1266       Ncvpaseq1 = 0
1267       dnwd0=0.0
1268       ftd=0.0
1269       fqd=0.0
1270       cin=0.
1271       !ym Attention pbase pas initialise dans concvl !!!!
1272       pbase=0
1273       !IM 180608
1274
1275       itau_con=0
1276       first=.false.
1277
1278    ENDIF  ! first
1279
1280    !ym => necessaire pour iflag_con != 2   
1281    pmfd(:,:) = 0.
1282    pen_u(:,:) = 0.
1283    pen_d(:,:) = 0.
1284    pde_d(:,:) = 0.
1285    pde_u(:,:) = 0.
1286    aam=0.
1287    d_t_adjwk(:,:)=0
1288    d_q_adjwk(:,:)=0
1289
1290    alp_bl_conv(:)=0.
1291
1292    torsfc=0.
1293    forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1294
1295    modname = 'physiq'
1296
1297    IF (debut) THEN
1298       CALL suphel ! initialiser constantes et parametres phys.
1299! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1300       tau_gl=5.
1301       CALL getin_p('tau_gl', tau_gl)
1302! tau_gl : constante de rappel de la temperature a la surface de la glace - en secondes
1303       tau_gl=86400.*tau_gl
1304       print*,'debut physiq_mod tau_gl=',tau_gl
1305       iflag_bug_t2m_ipslcm61 = 1
1306       CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
1307       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1308       CALL getin_p('random_notrig_max',random_notrig_max)
1309       CALL getin_p('ok_adjwk',ok_adjwk)
1310       IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
1311       ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1312                      ! 1 => convective adjustment but state variables are unchanged
1313                      ! 2 => convective adjustment and state variables are changed
1314       CALL getin_p('iflag_adjwk',iflag_adjwk)
1315       CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1316       CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
1317       CALL getin_p('oliqmax',oliqmax)
1318       CALL getin_p('oicemax',oicemax)
1319       CALL getin_p('ratqsp0',ratqsp0)
1320       CALL getin_p('ratqsdp',ratqsdp)
1321       iflag_wake_tend = 0
1322       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1323       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1324                                  ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1325       CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1326       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
1327       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1328       CALL getin_p('fl_ebil',fl_ebil)
1329       fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1330       CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1331       iflag_phytrac = 1 ! by default we do want to call phytrac
1332       CALL getin_p('iflag_phytrac',iflag_phytrac)
1333       nvm_lmdz = 13
1334       CALL getin_p('NVM',nvm_lmdz)
1335    ENDIF
1336
1337    IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1338
1339
1340    !======================================================================
1341    ! Gestion calendrier : mise a jour du module phys_cal_mod
1342    !
1343    !     CALL phys_cal_update(jD_cur,jH_cur)
1344
1345    !
1346    ! Si c'est le debut, il faut initialiser plusieurs choses
1347    !          ********
1348    !
1349    IF (debut) THEN
1350       !rv CRinitialisation de wght_th et lalim_conv pour la
1351       !definition de la couche alimentation de la convection a partir
1352       !des caracteristiques du thermique
1353       wght_th(:,:)=1.
1354       lalim_conv(:)=1
1355       !RC
1356       ustar(:,:)=0.
1357!       u10m(:,:)=0.
1358!       v10m(:,:)=0.
1359       rain_con(:)=0.
1360       snow_con(:)=0.
1361       topswai(:)=0.
1362       topswad(:)=0.
1363       solswai(:)=0.
1364       solswad(:)=0.
1365
1366       wmax_th(:)=0.
1367       tau_overturning_th(:)=0.
1368
1369       IF (type_trac == 'inca') THEN
1370          ! jg : initialisation jusqu'au ces variables sont dans restart
1371          ccm(:,:,:) = 0.
1372          tau_aero(:,:,:,:) = 0.
1373          piz_aero(:,:,:,:) = 0.
1374          cg_aero(:,:,:,:) = 0.
1375
1376          config_inca='none' ! default
1377          CALL getin_p('config_inca',config_inca)
1378
1379       ELSE
1380          config_inca='none' ! default
1381       ENDIF
1382
1383       IF (aerosol_couple .AND. (config_inca /= "aero" &
1384            .AND. config_inca /= "aeNP ")) THEN
1385          abort_message &
1386               = 'if aerosol_couple is activated, config_inca need to be ' &
1387               // 'aero or aeNP'
1388          CALL abort_physic (modname,abort_message,1)
1389       ENDIF
1390
1391
1392
1393       rnebcon0(:,:) = 0.0
1394       clwcon0(:,:) = 0.0
1395       rnebcon(:,:) = 0.0
1396       clwcon(:,:) = 0.0
1397
1398       !
1399       print*,'iflag_coupl,iflag_clos,iflag_wake', &
1400            iflag_coupl,iflag_clos,iflag_wake
1401       print*,'iflag_cycle_diurne', iflag_cycle_diurne
1402       !
1403       IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1404          abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1405          CALL abort_physic (modname,abort_message,1)
1406       ENDIF
1407       !
1408       !
1409       ! Initialiser les compteurs:
1410       !
1411       itap    = 0
1412       itaprad = 0
1413       itapcv = 0
1414       itapwk = 0
1415
1416       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1417       !! Un petit travail \`a faire ici.
1418       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1419
1420       IF (iflag_pbl>1) THEN
1421          PRINT*, "Using method MELLOR&YAMADA"
1422       ENDIF
1423
1424       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1425       ! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1426       ! phylmd plutot que dyn3d
1427       ! Attention : la version precedente n'etait pas tres propre.
1428       ! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1429       ! pour obtenir le meme resultat.
1430!jyg for fh<
1431!!       dtime=pdtphys
1432       dtime=NINT(pdtphys)
1433       WRITE(lunout,*) 'Pas de temps dtime pdtphys ',dtime,pdtphys
1434       IF (abs(dtime-pdtphys)>1.e-10) THEN
1435          abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1436          CALL abort_physic(modname,abort_message,1)
1437       ENDIF
1438!>jyg
1439       IF (MOD(NINT(86400./dtime),nbapp_rad).EQ.0) THEN
1440          radpas = NINT( 86400./dtime)/nbapp_rad
1441       ELSE
1442          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1443               'multiple de nbapp_rad'
1444          WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1445               'mais 1+1<>2'
1446          abort_message='nbre de pas de temps physique n est pas multiple ' &
1447               // 'de nbapp_rad'
1448          CALL abort_physic(modname,abort_message,1)
1449       ENDIF
1450       IF (nbapp_cv .EQ. 0) nbapp_cv=86400./dtime
1451       IF (nbapp_wk .EQ. 0) nbapp_wk=86400./dtime
1452       print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1453       IF (MOD(NINT(86400./dtime),nbapp_cv).EQ.0) THEN
1454          cvpas_0 = NINT( 86400./dtime)/nbapp_cv
1455          cvpas = cvpas_0
1456       print *,'physiq, cvpas ',cvpas
1457       ELSE
1458          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1459               'multiple de nbapp_cv'
1460          WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1461               'mais 1+1<>2'
1462          abort_message='nbre de pas de temps physique n est pas multiple ' &
1463               // 'de nbapp_cv'
1464          call abort_physic(modname,abort_message,1)
1465       ENDIF
1466       IF (MOD(NINT(86400./dtime),nbapp_wk).EQ.0) THEN
1467          wkpas = NINT( 86400./dtime)/nbapp_wk
1468       print *,'physiq, wkpas ',wkpas
1469       ELSE
1470          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1471               'multiple de nbapp_wk'
1472          WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1473               'mais 1+1<>2'
1474          abort_message='nbre de pas de temps physique n est pas multiple ' &
1475               // 'de nbapp_wk'
1476          call abort_physic(modname,abort_message,1)
1477       ENDIF
1478       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1479
1480       CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1481!jyg<
1482       IF (klon_glo==1) THEN
1483          IF (iflag_pbl > 1) THEN         
1484              pbl_tke(:,:,is_ave) = 0.
1485              DO nsrf=1,nbsrf
1486                DO k = 1,klev+1
1487                     pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1488                         +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1489                ENDDO
1490              ENDDO
1491          ELSE   ! (iflag_pbl > 1)
1492              pbl_tke(:,:,:) = 0.
1493          ENDIF  ! (iflag_pbl > 1)
1494!>jyg
1495       ENDIF
1496       !IM begin
1497       print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1498            ,ratqs(1,1)
1499       !IM end
1500
1501
1502       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1503       !
1504       ! on remet le calendrier a zero
1505       !
1506       IF (raz_date .eq. 1) THEN
1507          itau_phy = 0
1508       ENDIF
1509
1510       CALL printflag( tabcntr0,radpas,ok_journe, &
1511            ok_instan, ok_region )
1512       !
1513       IF (ABS(dtime-pdtphys).GT.0.001) THEN
1514          WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1515               pdtphys
1516          abort_message='Pas physique n est pas correct '
1517          !           call abort_physic(modname,abort_message,1)
1518          dtime=pdtphys
1519       ENDIF
1520       IF (nlon .NE. klon) THEN
1521          WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
1522               klon
1523          abort_message='nlon et klon ne sont pas coherents'
1524          CALL abort_physic(modname,abort_message,1)
1525       ENDIF
1526       IF (nlev .NE. klev) THEN
1527          WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1528               klev
1529          abort_message='nlev et klev ne sont pas coherents'
1530          CALL abort_physic(modname,abort_message,1)
1531       ENDIF
1532       !
1533       IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1534          WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1535          WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1536          abort_message='Nbre d appels au rayonnement insuffisant'
1537          CALL abort_physic(modname,abort_message,1)
1538       ENDIF
1539       WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1540       WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1541            ok_cvl
1542       !
1543       !KE43
1544       ! Initialisation pour la convection de K.E. (sb):
1545       IF (iflag_con.GE.3) THEN
1546
1547          WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1548          WRITE(lunout,*) &
1549               "On va utiliser le melange convectif des traceurs qui"
1550          WRITE(lunout,*)"est calcule dans convect4.3"
1551          WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1552
1553          DO i = 1, klon
1554             ema_cbmf(i) = 0.
1555             ema_pcb(i)  = 0.
1556             ema_pct(i)  = 0.
1557             !          ema_workcbmf(i) = 0.
1558          ENDDO
1559          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1560          DO i = 1, klon
1561             ibas_con(i) = 1
1562             itop_con(i) = 1
1563          ENDDO
1564          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1565          !================================================================
1566          !CR:04.12.07: initialisations poches froides
1567          ! Controle de ALE et ALP pour la fermeture convective (jyg)
1568          IF (iflag_wake>=1) THEN
1569             CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1570                  ,alp_bl_prescr, ale_bl_prescr)
1571             ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1572             !        print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1573             !
1574             ! Initialize tendencies of wake state variables (for some flag values
1575             ! they are not computed).
1576             d_deltat_wk(:,:) = 0.
1577             d_deltaq_wk(:,:) = 0.
1578             d_deltat_wk_gw(:,:) = 0.
1579             d_deltaq_wk_gw(:,:) = 0.
1580             d_deltat_vdf(:,:) = 0.
1581             d_deltaq_vdf(:,:) = 0.
1582             d_deltat_the(:,:) = 0.
1583             d_deltaq_the(:,:) = 0.
1584             d_deltat_ajs_cv(:,:) = 0.
1585             d_deltaq_ajs_cv(:,:) = 0.
1586             d_s_wk(:) = 0.
1587             d_dens_wk(:) = 0.
1588          ENDIF
1589
1590          !        do i = 1,klon
1591          !           Ale_bl(i)=0.
1592          !           Alp_bl(i)=0.
1593          !        enddo
1594
1595          !===================================================================
1596          !IM stations CFMIP
1597          nCFMIP=npCFMIP
1598          OPEN(98,file='npCFMIP_param.data',status='old', &
1599               form='formatted',iostat=iostat)
1600          IF (iostat == 0) THEN
1601             READ(98,*,end=998) nCFMIP
1602998          CONTINUE
1603             CLOSE(98)
1604             CONTINUE
1605             IF(nCFMIP.GT.npCFMIP) THEN
1606                print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1607                CALL abort_physic("physiq", "", 1)
1608             ELSE
1609                print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1610             ENDIF
1611
1612             !
1613             ALLOCATE(tabCFMIP(nCFMIP))
1614             ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1615             ALLOCATE(tabijGCM(nCFMIP))
1616             ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1617             ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1618             !
1619             ! lecture des nCFMIP stations CFMIP, de leur numero
1620             ! et des coordonnees geographiques lonCFMIP, latCFMIP
1621             !
1622             CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
1623                  lonCFMIP, latCFMIP)
1624             !
1625             ! identification des
1626             ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1627             ! grille de LMDZ
1628             ! 2) indices points tabijGCM de la grille physique 1d sur
1629             ! klon points
1630             ! 3) indices iGCM, jGCM de la grille physique 2d
1631             !
1632             CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1633                  tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1634             !
1635          ELSE
1636             ALLOCATE(tabijGCM(0))
1637             ALLOCATE(lonGCM(0), latGCM(0))
1638             ALLOCATE(iGCM(0), jGCM(0))
1639          ENDIF
1640       ELSE
1641          ALLOCATE(tabijGCM(0))
1642          ALLOCATE(lonGCM(0), latGCM(0))
1643          ALLOCATE(iGCM(0), jGCM(0))
1644       ENDIF
1645
1646       DO i=1,klon
1647          rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1648       ENDDO
1649
1650       !34EK
1651       IF (ok_orodr) THEN
1652
1653          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1654          ! FH sans doute a enlever de finitivement ou, si on le
1655          ! garde, l'activer justement quand ok_orodr = false.
1656          ! ce rugoro est utilise par la couche limite et fait double emploi
1657          ! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1658          !           DO i=1,klon
1659          !             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1660          !           ENDDO
1661          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1662          IF (ok_strato) THEN
1663             CALL SUGWD_strato(klon,klev,paprs,pplay)
1664          ELSE
1665             CALL SUGWD(klon,klev,paprs,pplay)
1666          ENDIF
1667
1668          DO i=1,klon
1669             zuthe(i)=0.
1670             zvthe(i)=0.
1671             IF (zstd(i).gt.10.) THEN
1672                zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1673                zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1674             ENDIF
1675          ENDDO
1676       ENDIF
1677       !
1678       !
1679       lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1680       WRITE(lunout,*)'La frequence de lecture surface est de ',  &
1681            lmt_pas
1682       !
1683       capemaxcels = 't_max(X)'
1684       t2mincels = 't_min(X)'
1685       t2maxcels = 't_max(X)'
1686       tinst = 'inst(X)'
1687       tave = 'ave(X)'
1688       !IM cf. AM 081204 BEG
1689       write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1690       !IM cf. AM 081204 END
1691       !
1692       !=============================================================
1693       !   Initialisation des sorties
1694       !=============================================================
1695
1696#ifdef CPP_XIOS
1697! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1698! initialised at that moment
1699       ! Get "missing_val" value from XML files (from temperature variable)
1700       !$OMP MASTER
1701       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1702       !$OMP END MASTER
1703       !$OMP BARRIER
1704       missing_val=missing_val_omp
1705#endif
1706
1707#ifdef CPP_IOIPSL
1708
1709       !$OMP MASTER
1710       ! FH : if ok_sync=.true. , the time axis is written at each time step
1711       ! in the output files. Only at the end in the opposite case
1712       ok_sync_omp=.false.
1713       CALL getin('ok_sync',ok_sync_omp)
1714       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1715            iGCM,jGCM,lonGCM,latGCM, &
1716            jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
1717            type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1718            ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1719            read_climoz, phys_out_filestations, &
1720            new_aod, aerosol_couple, &
1721            flag_aerosol_strat, pdtphys, paprs, pphis,  &
1722            pplay, lmax_th, ptconv, ptconvth, ivap,  &
1723            d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1724       !$OMP END MASTER
1725       !$OMP BARRIER
1726       ok_sync=ok_sync_omp
1727
1728       freq_outNMC(1) = ecrit_files(7)
1729       freq_outNMC(2) = ecrit_files(8)
1730       freq_outNMC(3) = ecrit_files(9)
1731       WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1732       WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1733       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1734
1735#ifndef CPP_XIOS
1736       CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
1737#endif
1738
1739#endif
1740       ecrit_reg = ecrit_reg * un_jour
1741       ecrit_tra = ecrit_tra * un_jour
1742
1743       !XXXPB Positionner date0 pour initialisation de ORCHIDEE
1744       date0 = jD_ref
1745       WRITE(*,*) 'physiq date0 : ',date0
1746       !
1747       !
1748       !
1749       ! Prescrire l'ozone dans l'atmosphere
1750       !
1751       !
1752       !c         DO i = 1, klon
1753       !c         DO k = 1, klev
1754       !c            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1755       !c         ENDDO
1756       !c         ENDDO
1757       !
1758       IF (type_trac == 'inca') THEN
1759#ifdef INCA
1760          CALL VTe(VTphysiq)
1761          CALL VTb(VTinca)
1762          calday = REAL(days_elapsed) + jH_cur
1763          WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1764
1765          CALL chemini(  &
1766               rg, &
1767               ra, &
1768               cell_area, &
1769               latitude_deg, &
1770               longitude_deg, &
1771               presnivs, &
1772               calday, &
1773               klon, &
1774               nqtot, &
1775               nqo, &
1776               pdtphys, &
1777               annee_ref, &
1778               year_cur, &
1779               day_ref,  &
1780               day_ini, &
1781               start_time, &
1782               itau_phy, &
1783               date0, &
1784               io_lon, &
1785               io_lat, &
1786               chemistry_couple, &
1787               init_source, &
1788               init_tauinca, &
1789               init_pizinca, &
1790               init_cginca, &
1791               init_ccminca)
1792
1793          ! initialisation des variables depuis le restart de inca
1794          ccm(:,:,:) = init_ccminca
1795          tau_aero(:,:,:,:) = init_tauinca
1796          piz_aero(:,:,:,:) = init_pizinca
1797          cg_aero(:,:,:,:) = init_cginca
1798
1799
1800          CALL VTe(VTinca)
1801          CALL VTb(VTphysiq)
1802#endif
1803       ENDIF
1804       !
1805       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1806       ! Nouvelle initialisation pour le rayonnement RRTM
1807       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1808
1809       CALL iniradia(klon,klev,paprs(1,1:klev+1))
1810
1811       !$omp single
1812       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
1813           press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1814       !$omp end single
1815       !
1816       !IM betaCRF
1817       pfree=70000. !Pa
1818       beta_pbl=1.
1819       beta_free=1.
1820       lon1_beta=-180.
1821       lon2_beta=+180.
1822       lat1_beta=90.
1823       lat2_beta=-90.
1824       mskocean_beta=.FALSE.
1825
1826       !albedo SB >>>
1827       select case(nsw)
1828       case(2)
1829          SFRWL(1)=0.45538747
1830          SFRWL(2)=0.54461211
1831       case(4)
1832          SFRWL(1)=0.45538747
1833          SFRWL(2)=0.32870591
1834          SFRWL(3)=0.18568763
1835          SFRWL(4)=3.02191470E-02
1836       case(6)
1837          SFRWL(1)=1.28432794E-03
1838          SFRWL(2)=0.12304168
1839          SFRWL(3)=0.33106142
1840          SFRWL(4)=0.32870591
1841          SFRWL(5)=0.18568763
1842          SFRWL(6)=3.02191470E-02
1843       end select
1844
1845
1846       !albedo SB <<<
1847
1848       OPEN(99,file='beta_crf.data',status='old', &
1849            form='formatted',err=9999)
1850       READ(99,*,end=9998) pfree
1851       READ(99,*,end=9998) beta_pbl
1852       READ(99,*,end=9998) beta_free
1853       READ(99,*,end=9998) lon1_beta
1854       READ(99,*,end=9998) lon2_beta
1855       READ(99,*,end=9998) lat1_beta
1856       READ(99,*,end=9998) lat2_beta
1857       READ(99,*,end=9998) mskocean_beta
18589998   Continue
1859       CLOSE(99)
18609999   Continue
1861       WRITE(*,*)'pfree=',pfree
1862       WRITE(*,*)'beta_pbl=',beta_pbl
1863       WRITE(*,*)'beta_free=',beta_free
1864       WRITE(*,*)'lon1_beta=',lon1_beta
1865       WRITE(*,*)'lon2_beta=',lon2_beta
1866       WRITE(*,*)'lat1_beta=',lat1_beta
1867       WRITE(*,*)'lat2_beta=',lat2_beta
1868       WRITE(*,*)'mskocean_beta=',mskocean_beta
1869
1870      !lwoff=y : offset LW CRE for radiation code and other schemes
1871      !lwoff=y : betalwoff=1.
1872      betalwoff=0.
1873      IF (ok_lwoff) THEN
1874         betalwoff=1.
1875      ENDIF
1876      WRITE(*,*)'ok_lwoff=',ok_lwoff
1877      !
1878      !lwoff=y to begin only sollw and sollwdown are set up to CS values
1879      sollw = sollw + betalwoff * (sollw0 - sollw)
1880      sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
1881                    sollwdown(:))
1882    ENDIF
1883    !
1884    !   ****************     Fin  de   IF ( debut  )   ***************
1885    !
1886    !
1887    ! Incrementer le compteur de la physique
1888    !
1889    itap   = itap + 1
1890    IF (is_master .OR. prt_level > 9) THEN
1891      IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
1892         WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
1893         WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
1894 100     FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
1895      ENDIF
1896    ENDIF
1897    !
1898    !
1899    ! Update fraction of the sub-surfaces (pctsrf) and
1900    ! initialize, where a new fraction has appeared, all variables depending
1901    ! on the surface fraction.
1902    !
1903    CALL change_srf_frac(itap, dtime, days_elapsed+1,  &
1904         pctsrf, fevap, z0m, z0h, agesno,              &
1905         falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1906
1907    ! Update time and other variables in Reprobus
1908    IF (type_trac == 'repr') THEN
1909#ifdef REPROBUS
1910       CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1911       print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1912       CALL Rtime(debut)
1913#endif
1914    ENDIF
1915
1916
1917    ! Tendances bidons pour les processus qui n'affectent pas certaines
1918    ! variables.
1919    du0(:,:)=0.
1920    dv0(:,:)=0.
1921    dt0 = 0.
1922    dq0(:,:)=0.
1923    dql0(:,:)=0.
1924    dqi0(:,:)=0.
1925    dsig0(:) = 0.
1926    ddens0(:) = 0.
1927    wkoccur1(:)=1
1928    !
1929    ! Mettre a zero des variables de sortie (pour securite)
1930    !
1931    DO i = 1, klon
1932       d_ps(i) = 0.0
1933    ENDDO
1934    DO k = 1, klev
1935       DO i = 1, klon
1936          d_t(i,k) = 0.0
1937          d_u(i,k) = 0.0
1938          d_v(i,k) = 0.0
1939       ENDDO
1940    ENDDO
1941    DO iq = 1, nqtot
1942       DO k = 1, klev
1943          DO i = 1, klon
1944             d_qx(i,k,iq) = 0.0
1945          ENDDO
1946       ENDDO
1947    ENDDO
1948    beta_prec_fisrt(:,:)=0.
1949    beta_prec(:,:)=0.
1950    !
1951    !   Output variables from the convective scheme should not be set to 0
1952    !   since convection is not always called at every time step.
1953    IF (ok_bug_cv_trac) THEN
1954      da(:,:)=0.
1955      mp(:,:)=0.
1956      phi(:,:,:)=0.
1957      ! RomP >>>
1958      phi2(:,:,:)=0.
1959      epmlmMm(:,:,:)=0.
1960      eplaMm(:,:)=0.
1961      d1a(:,:)=0.
1962      dam(:,:)=0.
1963      pmflxr(:,:)=0.
1964      pmflxs(:,:)=0.
1965      ! RomP <<<
1966    ENDIF
1967
1968    !
1969    ! Ne pas affecter les valeurs entrees de u, v, h, et q
1970    !
1971    DO k = 1, klev
1972       DO i = 1, klon
1973          t_seri(i,k)  = t(i,k)
1974          u_seri(i,k)  = u(i,k)
1975          v_seri(i,k)  = v(i,k)
1976          q_seri(i,k)  = qx(i,k,ivap)
1977          ql_seri(i,k) = qx(i,k,iliq)
1978          !CR: ATTENTION, on rajoute la variable glace
1979          IF (nqo.eq.2) THEN
1980             qs_seri(i,k) = 0.
1981          ELSE IF (nqo.eq.3) THEN
1982             qs_seri(i,k) = qx(i,k,isol)
1983          ENDIF
1984       ENDDO
1985    ENDDO
1986    !
1987    !--OB mass fixer
1988    IF (mass_fixer) THEN
1989    !--store initial water burden
1990    qql1(:)=0.0
1991    DO k = 1, klev
1992      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
1993    ENDDO
1994    ENDIF
1995    !--fin mass fixer
1996
1997    tke0(:,:)=pbl_tke(:,:,is_ave)
1998    !CR:Nombre de traceurs de l'eau: nqo
1999    !  IF (nqtot.GE.3) THEN
2000    IF (nqtot.GE.(nqo+1)) THEN
2001       !     DO iq = 3, nqtot       
2002       DO iq = nqo+1, nqtot 
2003          DO  k = 1, klev
2004             DO  i = 1, klon
2005                !              tr_seri(i,k,iq-2) = qx(i,k,iq)
2006                tr_seri(i,k,iq-nqo) = qx(i,k,iq)
2007             ENDDO
2008          ENDDO
2009       ENDDO
2010    ELSE
2011       DO k = 1, klev
2012          DO i = 1, klon
2013             tr_seri(i,k,1) = 0.0
2014          ENDDO
2015       ENDDO
2016    ENDIF
2017    !
2018    DO i = 1, klon
2019       ztsol(i) = 0.
2020    ENDDO
2021    DO nsrf = 1, nbsrf
2022       DO i = 1, klon
2023          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2024       ENDDO
2025    ENDDO
2026    ! Initialize variables used for diagnostic purpose
2027    IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2028
2029    ! Diagnostiquer la tendance dynamique
2030    !
2031    IF (ancien_ok) THEN
2032    !
2033       d_u_dyn(:,:)  = (u_seri(:,:)-u_ancien(:,:))/dtime
2034       d_v_dyn(:,:)  = (v_seri(:,:)-v_ancien(:,:))/dtime
2035       d_t_dyn(:,:)  = (t_seri(:,:)-t_ancien(:,:))/dtime
2036       d_q_dyn(:,:)  = (q_seri(:,:)-q_ancien(:,:))/dtime
2037       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/dtime
2038       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/dtime
2039       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2040       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/dtime
2041       CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2042       d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/dtime
2043       CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2044       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/dtime
2045       ! !! RomP >>>   td dyn traceur
2046       IF (nqtot.GT.nqo) THEN     ! jyg
2047          DO iq = nqo+1, nqtot      ! jyg
2048              d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/dtime ! jyg
2049          ENDDO
2050       ENDIF
2051       ! !! RomP <<<
2052    ELSE
2053       d_u_dyn(:,:)  = 0.0
2054       d_v_dyn(:,:)  = 0.0
2055       d_t_dyn(:,:)  = 0.0
2056       d_q_dyn(:,:)  = 0.0
2057       d_ql_dyn(:,:) = 0.0
2058       d_qs_dyn(:,:) = 0.0
2059       d_q_dyn2d(:)  = 0.0
2060       d_ql_dyn2d(:) = 0.0
2061       d_qs_dyn2d(:) = 0.0
2062       ! !! RomP >>>   td dyn traceur
2063       IF (nqtot.GT.nqo) THEN                                       ! jyg
2064          DO iq = nqo+1, nqtot                                      ! jyg
2065              d_tr_dyn(:,:,iq-nqo)= 0.0                             ! jyg
2066          ENDDO
2067       ENDIF
2068       ! !! RomP <<<
2069       ancien_ok = .TRUE.
2070    ENDIF
2071    !
2072    ! Ajouter le geopotentiel du sol:
2073    !
2074    DO k = 1, klev
2075       DO i = 1, klon
2076          zphi(i,k) = pphi(i,k) + pphis(i)
2077       ENDDO
2078    ENDDO
2079    !
2080    ! Verifier les temperatures
2081    !
2082    !IM BEG
2083    IF (check) THEN
2084       amn=MIN(ftsol(1,is_ter),1000.)
2085       amx=MAX(ftsol(1,is_ter),-1000.)
2086       DO i=2, klon
2087          amn=MIN(ftsol(i,is_ter),amn)
2088          amx=MAX(ftsol(i,is_ter),amx)
2089       ENDDO
2090       !
2091       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2092    ENDIF !(check) THEN
2093    !IM END
2094    !
2095    CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2096    IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2097
2098    !
2099    !IM BEG
2100    IF (check) THEN
2101       amn=MIN(ftsol(1,is_ter),1000.)
2102       amx=MAX(ftsol(1,is_ter),-1000.)
2103       DO i=2, klon
2104          amn=MIN(ftsol(i,is_ter),amn)
2105          amx=MAX(ftsol(i,is_ter),amx)
2106       ENDDO
2107       !
2108       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2109    ENDIF !(check) THEN
2110    !IM END
2111    !
2112    ! Mettre en action les conditions aux limites (albedo, sst, etc.).
2113    ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2114    !
2115    ! Update ozone if day change
2116    IF (MOD(itap-1,lmt_pas) == 0) THEN
2117       IF (read_climoz <= 0) THEN
2118          ! Once per day, update ozone from Royer:
2119          IF (solarlong0<-999.) then
2120             ! Generic case with evolvoing season
2121             zzz=real(days_elapsed+1)
2122          ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2123             ! Particular case with annual mean insolation
2124             zzz=real(90) ! could be revisited
2125             IF (read_climoz/=-1) THEN
2126                abort_message ='read_climoz=-1 is recommended when ' &
2127                     // 'solarlong0=1000.'
2128                CALL abort_physic (modname,abort_message,1)
2129             ENDIF
2130          ELSE
2131             ! Case where the season is imposed with solarlong0
2132             zzz=real(90) ! could be revisited
2133          ENDIF
2134
2135          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2136       ELSE
2137          !--- ro3i = elapsed days number since current year 1st january, 0h
2138          ro3i=days_elapsed+jh_cur-jh_1jan
2139          !--- scaling for old style files (360 records)
2140          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2141          IF(adjust_tropopause) THEN
2142             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
2143                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2144                      time_climoz ,  longitude_deg,   latitude_deg,          &
2145                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2146          ELSE
2147             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
2148                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2149                      time_climoz )
2150          END IF
2151          ! Convert from mole fraction of ozone to column density of ozone in a
2152          ! cell, in kDU:
2153          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2154               * zmasse / dobson_u / 1e3
2155          ! (By regridding ozone values for LMDZ only once a day, we
2156          ! have already neglected the variation of pressure in one
2157          ! day. So do not recompute "wo" at each time step even if
2158          ! "zmasse" changes a little.)
2159       ENDIF
2160    ENDIF
2161    !
2162    ! Re-evaporer l'eau liquide nuageuse
2163    !
2164     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2165   &         d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2166
2167     CALL add_phys_tend &
2168            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2169               'eva',abortphy,flag_inhib_tend,itap,0)
2170    call prt_enerbil('eva',itap)
2171
2172    !=========================================================================
2173    ! Calculs de l'orbite.
2174    ! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2175    ! doit donc etre plac\'e avant radlwsw et pbl_surface
2176
2177    ! !!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2178    CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2179    day_since_equinox = (jD_cur + jH_cur) - jD_eq
2180    !
2181    !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
2182    !   solarlong0
2183    IF (solarlong0<-999.) THEN
2184       IF (new_orbit) THEN
2185          ! calcul selon la routine utilisee pour les planetes
2186          CALL solarlong(day_since_equinox, zlongi, dist)
2187       ELSE
2188          ! calcul selon la routine utilisee pour l'AR4
2189          CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2190       ENDIF
2191    ELSE
2192       zlongi=solarlong0  ! longitude solaire vraie
2193       dist=1.            ! distance au soleil / moyenne
2194    ENDIF
2195
2196    IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2197
2198
2199    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2200    ! Calcul de l'ensoleillement :
2201    ! ============================
2202    ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2203    ! l'annee a partir d'une formule analytique.
2204    ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2205    ! non nul aux poles.
2206    IF (abs(solarlong0-1000.)<1.e-4) THEN
2207       CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2208            latitude_deg,longitude_deg,rmu0,fract)
2209       swradcorr(:) = 1.0
2210       JrNt(:) = 1.0
2211       zrmu0(:) = rmu0(:)
2212    ELSE
2213       ! recode par Olivier Boucher en sept 2015
2214       SELECT CASE (iflag_cycle_diurne)
2215       CASE(0) 
2216          !  Sans cycle diurne
2217          CALL angle(zlongi, latitude_deg, fract, rmu0)
2218          swradcorr = 1.0
2219          JrNt = 1.0
2220          zrmu0 = rmu0
2221       CASE(1) 
2222          !  Avec cycle diurne sans application des poids
2223          !  bit comparable a l ancienne formulation cycle_diurne=true
2224          !  on integre entre gmtime et gmtime+radpas
2225          zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2226          CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2227               latitude_deg,longitude_deg,rmu0,fract)
2228          zrmu0 = rmu0
2229          swradcorr = 1.0
2230          ! Calcul du flag jour-nuit
2231          JrNt = 0.0
2232          WHERE (fract.GT.0.0) JrNt = 1.0
2233       CASE(2) 
2234          !  Avec cycle diurne sans application des poids
2235          !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2236          !  Comme cette routine est appele a tous les pas de temps de
2237          !  la physique meme si le rayonnement n'est pas appele je
2238          !  remonte en arriere les radpas-1 pas de temps
2239          !  suivant. Petite ruse avec MOD pour prendre en compte le
2240          !  premier pas de temps de la physique pendant lequel
2241          !  itaprad=0
2242          zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)     
2243          zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
2244          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2245               latitude_deg,longitude_deg,rmu0,fract)
2246          !
2247          ! Calcul des poids
2248          !
2249          zdtime1=-dtime !--on corrige le rayonnement pour representer le
2250          zdtime2=0.0    !--pas de temps de la physique qui se termine
2251          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2252               latitude_deg,longitude_deg,zrmu0,zfract)
2253          swradcorr = 0.0
2254          WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2255               swradcorr=zfract/fract*zrmu0/rmu0
2256          ! Calcul du flag jour-nuit
2257          JrNt = 0.0
2258          WHERE (zfract.GT.0.0) JrNt = 1.0
2259       END SELECT
2260    ENDIF
2261    sza_o = ACOS (rmu0) *180./pi
2262
2263    IF (mydebug) THEN
2264       CALL writefield_phy('u_seri',u_seri,nbp_lev)
2265       CALL writefield_phy('v_seri',v_seri,nbp_lev)
2266       CALL writefield_phy('t_seri',t_seri,nbp_lev)
2267       CALL writefield_phy('q_seri',q_seri,nbp_lev)
2268    ENDIF
2269
2270    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2271    ! Appel au pbl_surface : Planetary Boudary Layer et Surface
2272    ! Cela implique tous les interactions des sous-surfaces et la
2273    ! partie diffusion turbulent du couche limit.
2274    !
2275    ! Certains varibales de sorties de pbl_surface sont utiliser que pour
2276    ! ecriture des fihiers hist_XXXX.nc, ces sont :
2277    !   qsol,      zq2m,      s_pblh,  s_lcl,
2278    !   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
2279    !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
2280    !   zu10m,     zv10m,   fder,
2281    !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
2282    !   frugs,     agesno,    fsollw,  fsolsw,
2283    !   d_ts,      fevap,     fluxlat, t2m,
2284    !   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
2285    !
2286    ! Certains ne sont pas utiliser du tout :
2287    !   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2288    !
2289
2290    ! Calcul de l'humidite de saturation au niveau du sol
2291
2292
2293
2294    IF (iflag_pbl/=0) THEN
2295
2296       !jyg+nrlmd<
2297!!jyg       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2298       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2299          print *,'debut du splitting de la PBL'
2300       ENDIF
2301       ! !!
2302       !>jyg+nrlmd
2303       !
2304       !-------gustiness calculation-------!
2305       IF (iflag_gusts==0) THEN
2306          gustiness(1:klon)=0
2307       ELSE IF (iflag_gusts==1) THEN
2308          gustiness(1:klon)=f_gust_bl*ale_bl(1:klon)+f_gust_wk*ale_wake(1:klon)
2309       ELSE IF (iflag_gusts==2) THEN
2310          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
2311          ! ELSE IF (iflag_gusts==2) THEN
2312          !    do i = 1, klon
2313          !       gustiness(i)=f_gust_bl*ale_bl(i)+sigma_wk(i)*f_gust_wk&
2314          !           *ale_wake(i) !! need to make sigma_wk accessible here
2315          !    enddo
2316          ! ELSE IF (iflag_gusts==3) THEN
2317          !    do i = 1, klon
2318          !       gustiness(i)=f_gust_bl*alp_bl(i)+f_gust_wk*alp_wake(i)
2319          !    enddo
2320       ENDIF
2321
2322
2323
2324       CALL pbl_surface(  &
2325            dtime,     date0,     itap,    days_elapsed+1, &
2326            debut,     lafin, &
2327            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
2328            zsig,      sollwdown, pphi,    cldt,      &
2329            rain_fall, snow_fall, solsw,   sollw,     &
2330            gustiness,                                &
2331            t_seri,    q_seri,    u_seri,  v_seri,    &
2332                                !nrlmd+jyg<
2333            wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2334                                !>nrlmd+jyg
2335            pplay,     paprs,     pctsrf,             &
2336            ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2337                                !albedo SB <<<
2338            cdragh,    cdragm,  u1,    v1,            &
2339                                !albedo SB >>>
2340                                ! albsol1,   albsol2,   sens,    evap,      &
2341            albsol_dir,   albsol_dif,   sens,    evap,   & 
2342                                !albedo SB <<<
2343            albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
2344            zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
2345            d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
2346                                !nrlmd<
2347                                !jyg<
2348            d_t_vdf_w, d_q_vdf_w, &
2349            d_t_vdf_x, d_q_vdf_x, &
2350            sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2351                                !>jyg
2352            delta_tsurf,wake_dens, &
2353            cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2354            kh,kh_x,kh_w, &
2355                                !>nrlmd
2356            coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2357            slab_wfbils,                 &
2358            qsol,      zq2m,      s_pblh,  s_lcl, &
2359                                !jyg<
2360            s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2361                                !>jyg
2362            s_capCL,   s_oliqCL,  s_cteiCL,s_pblT, &
2363            s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
2364            zustar, zu10m,     zv10m,   fder, &
2365            zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
2366            z0m, z0h,     agesno,    fsollw,  fsolsw, &
2367            d_ts,      fevap,     fluxlat, t2m, &
2368            wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2369            fluxt,   fluxu,  fluxv, &
2370            dsens,     devap,     zxsnow, &
2371            zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
2372                                !nrlmd+jyg<
2373            wake_delta_pbl_TKE, &
2374                                !>nrlmd+jyg
2375             treedrg )
2376!FC
2377       !
2378       !  Add turbulent diffusion tendency to the wake difference variables
2379!!jyg       IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2380       IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2381!jyg<
2382          d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2383          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2384          CALL add_wake_tend &
2385             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy)
2386       ELSE
2387          d_deltat_vdf(:,:) = 0.
2388          d_deltaq_vdf(:,:) = 0.
2389!>jyg
2390       ENDIF
2391        if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
2392          CALL borne_var_surf( klon,nbsrf,                     &
2393            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),    &
2394            ftsol,pctsrf,                                       &
2395            t2m, q2m, u10m, v10m,                               &
2396            zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor)
2397        ELSE
2398          zt2m_cor(:)=zt2m(:)
2399          zq2m_cor(:)=zq2m(:)
2400          zu10m_cor(:)=zu10m(:)
2401          zv10m_cor(:)=zv10m(:)
2402        ENDIF
2403
2404
2405
2406
2407
2408       !---------------------------------------------------------------------
2409       ! ajout des tendances de la diffusion turbulente
2410       IF (klon_glo==1) THEN
2411          CALL add_pbl_tend &
2412               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2413               'vdf',abortphy,flag_inhib_tend,itap)
2414       ELSE
2415          CALL add_phys_tend &
2416               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2417               'vdf',abortphy,flag_inhib_tend,itap,0)
2418       ENDIF
2419       call prt_enerbil('vdf',itap)
2420       !--------------------------------------------------------------------
2421
2422       IF (mydebug) THEN
2423          CALL writefield_phy('u_seri',u_seri,nbp_lev)
2424          CALL writefield_phy('v_seri',v_seri,nbp_lev)
2425          CALL writefield_phy('t_seri',t_seri,nbp_lev)
2426          CALL writefield_phy('q_seri',q_seri,nbp_lev)
2427       ENDIF
2428
2429       !albedo SB >>>
2430       albsol1=0.
2431       albsol2=0.
2432       falb1=0.
2433       falb2=0.
2434       SELECT CASE(nsw)
2435       CASE(2)
2436          albsol1=albsol_dir(:,1)
2437          albsol2=albsol_dir(:,2)
2438          falb1=falb_dir(:,1,:)
2439          falb2=falb_dir(:,2,:)
2440       CASE(4)
2441          albsol1=albsol_dir(:,1)
2442          albsol2=albsol_dir(:,2)*SFRWL(2)+albsol_dir(:,3)*SFRWL(3) &
2443               +albsol_dir(:,4)*SFRWL(4)
2444          albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2445          falb1=falb_dir(:,1,:)
2446          falb2=falb_dir(:,2,:)*SFRWL(2)+falb_dir(:,3,:)*SFRWL(3) &
2447               +falb_dir(:,4,:)*SFRWL(4)
2448          falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2449       CASE(6)
2450          albsol1=albsol_dir(:,1)*SFRWL(1)+albsol_dir(:,2)*SFRWL(2) &
2451               +albsol_dir(:,3)*SFRWL(3)
2452          albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2453          albsol2=albsol_dir(:,4)*SFRWL(4)+albsol_dir(:,5)*SFRWL(5) &
2454               +albsol_dir(:,6)*SFRWL(6)
2455          albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2456          falb1=falb_dir(:,1,:)*SFRWL(1)+falb_dir(:,2,:)*SFRWL(2) &
2457               +falb_dir(:,3,:)*SFRWL(3)
2458          falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2459          falb2=falb_dir(:,4,:)*SFRWL(4)+falb_dir(:,5,:)*SFRWL(5) &
2460               +falb_dir(:,6,:)*SFRWL(6)
2461          falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2462       END SELECt
2463       !albedo SB <<<
2464
2465
2466       CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2467            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2468
2469    ENDIF
2470    ! =================================================================== c
2471    !   Calcul de Qsat
2472
2473    DO k = 1, klev
2474       DO i = 1, klon
2475          zx_t = t_seri(i,k)
2476          IF (thermcep) THEN
2477             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2478             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2479             zx_qs  = MIN(0.5,zx_qs)
2480             zcor   = 1./(1.-retv*zx_qs)
2481             zx_qs  = zx_qs*zcor
2482          ELSE
2483             !!           IF (zx_t.LT.t_coup) THEN             !jyg
2484             IF (zx_t.LT.rtt) THEN                  !jyg
2485                zx_qs = qsats(zx_t)/pplay(i,k)
2486             ELSE
2487                zx_qs = qsatl(zx_t)/pplay(i,k)
2488             ENDIF
2489          ENDIF
2490          zqsat(i,k)=zx_qs
2491       ENDDO
2492    ENDDO
2493
2494    IF (prt_level.ge.1) THEN
2495       write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
2496       write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2497    ENDIF
2498    !
2499    ! Appeler la convection (au choix)
2500    !
2501    DO k = 1, klev
2502       DO i = 1, klon
2503          conv_q(i,k) = d_q_dyn(i,k)  &
2504               + d_q_vdf(i,k)/dtime
2505          conv_t(i,k) = d_t_dyn(i,k)  &
2506               + d_t_vdf(i,k)/dtime
2507       ENDDO
2508    ENDDO
2509    IF (check) THEN
2510       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2511       WRITE(lunout,*) "avantcon=", za
2512    ENDIF
2513    zx_ajustq = .FALSE.
2514    IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2515    IF (zx_ajustq) THEN
2516       DO i = 1, klon
2517          z_avant(i) = 0.0
2518       ENDDO
2519       DO k = 1, klev
2520          DO i = 1, klon
2521             z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2522                  *(paprs(i,k)-paprs(i,k+1))/RG
2523          ENDDO
2524       ENDDO
2525    ENDIF
2526
2527    ! Calcule de vitesse verticale a partir de flux de masse verticale
2528    DO k = 1, klev
2529       DO i = 1, klon
2530          omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2531       ENDDO
2532    ENDDO
2533
2534    IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2535         omega(igout, :)
2536    !
2537    ! Appel de la convection tous les "cvpas"
2538    !
2539!!jyg    IF (MOD(itapcv,cvpas).EQ.0) THEN
2540!!    print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2541!!                       itapcv, cvpas, itap-1, cvpas_0
2542    IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
2543
2544    !
2545    ! Mettre a zero des variables de sortie (pour securite)
2546    !
2547    pmflxr(:,:) = 0.
2548    pmflxs(:,:) = 0.
2549    wdtrainA(:,:) = 0.
2550    wdtrainM(:,:) = 0.
2551    upwd(:,:) = 0.
2552    dnwd(:,:) = 0.
2553    ep(:,:) = 0.
2554    da(:,:)=0.
2555    mp(:,:)=0.
2556    wght_cvfd(:,:)=0.
2557    phi(:,:,:)=0.
2558    phi2(:,:,:)=0.
2559    epmlmMm(:,:,:)=0.
2560    eplaMm(:,:)=0.
2561    d1a(:,:)=0.
2562    dam(:,:)=0.
2563    elij(:,:,:)=0.
2564    ev(:,:)=0.
2565    clw(:,:)=0.
2566    sij(:,:,:)=0.
2567    !
2568    IF (iflag_con.EQ.1) THEN
2569       abort_message ='reactiver le call conlmd dans physiq.F'
2570       CALL abort_physic (modname,abort_message,1)
2571       !     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2572       !    .             d_t_con, d_q_con,
2573       !    .             rain_con, snow_con, ibas_con, itop_con)
2574    ELSE IF (iflag_con.EQ.2) THEN
2575       CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2576            conv_t, conv_q, -evap, omega, &
2577            d_t_con, d_q_con, rain_con, snow_con, &
2578            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2579            kcbot, kctop, kdtop, pmflxr, pmflxs)
2580       d_u_con = 0.
2581       d_v_con = 0.
2582
2583       WHERE (rain_con < 0.) rain_con = 0.
2584       WHERE (snow_con < 0.) snow_con = 0.
2585       DO i = 1, klon
2586          ibas_con(i) = klev+1 - kcbot(i)
2587          itop_con(i) = klev+1 - kctop(i)
2588       ENDDO
2589    ELSE IF (iflag_con.GE.3) THEN
2590       ! nb of tracers for the KE convection:
2591       ! MAF la partie traceurs est faite dans phytrac
2592       ! on met ntra=1 pour limiter les appels mais on peut
2593       ! supprimer les calculs / ftra.
2594       ntra = 1
2595
2596       !=======================================================================
2597       !ajout pour la parametrisation des poches froides: calcul de
2598       !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2599       IF (iflag_wake>=1) THEN
2600         DO k=1,klev
2601            DO i=1,klon
2602                t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2603                q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2604                t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2605                q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2606            ENDDO
2607         ENDDO
2608       ELSE
2609               t_w(:,:) = t_seri(:,:)
2610                q_w(:,:) = q_seri(:,:)
2611                t_x(:,:) = t_seri(:,:)
2612                q_x(:,:) = q_seri(:,:)
2613       ENDIF
2614       !
2615       !jyg<
2616       ! Perform dry adiabatic adjustment on wake profile
2617       ! The corresponding tendencies are added to the convective tendencies
2618       ! after the call to the convective scheme.
2619       IF (iflag_wake>=1) then
2620          IF (iflag_adjwk >= 1) THEN
2621             limbas(:) = 1
2622             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2623                  d_t_adjwk, d_q_adjwk)
2624             !
2625             DO k=1,klev
2626                DO i=1,klon
2627                   IF (wake_s(i) .GT. 1.e-3) THEN
2628                      t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2629                      q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2630                      d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2631                      d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2632                   ELSE
2633                      d_deltat_ajs_cv(i,k) = 0.
2634                      d_deltaq_ajs_cv(i,k) = 0.
2635                   ENDIF
2636                ENDDO
2637             ENDDO
2638             IF (iflag_adjwk == 2) THEN
2639               CALL add_wake_tend &
2640                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2641             ENDIF  ! (iflag_adjwk == 2)
2642          ENDIF  ! (iflag_adjwk >= 1)
2643       ENDIF ! (iflag_wake>=1)
2644       !>jyg
2645       !
2646       
2647!!      print *,'physiq. q_w(1,k), q_x(1,k) ', &
2648!!             (k, q_w(1,k), q_x(1,k),k=1,25)
2649
2650!jyg<
2651       CALL alpale( debut, itap, dtime, paprs, omega, t_seri,   &
2652                    alp_offset, it_wape_prescr,  wape_prescr, fip_prescr, &
2653                    ale_bl_prescr, alp_bl_prescr, &
2654                    wake_pe, wake_fip,  &
2655                    Ale_bl, Ale_bl_trig, Alp_bl, &
2656                    Ale, Alp , Ale_wake, Alp_wake)
2657!>jyg
2658!
2659       ! sb, oct02:
2660       ! Schema de convection modularise et vectorise:
2661       ! (driver commun aux versions 3 et 4)
2662       !
2663       IF (ok_cvl) THEN ! new driver for convectL
2664          !
2665          !jyg<
2666          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2667          ! Calculate the upmost level of deep convection loops: k_upper_cv
2668          !  (near 22 km)
2669          k_upper_cv = klev
2670          !izero = klon/2+1/klon
2671          !DO k = klev,1,-1
2672          !   IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2673          !ENDDO
2674          ! FH : nouveau calcul base sur un profil global sans quoi
2675          ! le modele etait sensible au decoupage de domaines
2676          DO k = klev,1,-1
2677             IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
2678          ENDDO
2679          IF (prt_level .ge. 5) THEN
2680             Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2681                  k_upper_cv
2682          ENDIF
2683          !
2684          !>jyg
2685          IF (type_trac == 'repr') THEN
2686             nbtr_tmp=ntra
2687          ELSE
2688             nbtr_tmp=nbtr
2689          ENDIF
2690          !jyg   iflag_con est dans clesphys
2691          !c          CALL concvl (iflag_con,iflag_clos,
2692          CALL concvl (iflag_clos, &
2693               dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
2694               t_w,q_w,wake_s, &
2695               u_seri,v_seri,tr_seri,nbtr_tmp, &
2696               ALE,ALP, &
2697               sig1,w01, &
2698               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2699               rain_con, snow_con, ibas_con, itop_con, sigd, &
2700               ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2701               Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2702               pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2703                                ! RomP >>>
2704                                !!     .        pmflxr,pmflxs,da,phi,mp,
2705                                !!     .        ftd,fqd,lalim_conv,wght_th)
2706               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2707               ftd,fqd,lalim_conv,wght_th, &
2708               ev, ep,epmlmMm,eplaMm, &
2709               wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2710               tau_cld_cv,coefw_cld_cv,epmax_diag)
2711
2712          ! RomP <<<
2713
2714          !IM begin
2715          !       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2716          !    .dnwd0(1,1),ftd(1,1),fqd(1,1)
2717          !IM end
2718          !IM cf. FH
2719          clwcon0=qcondc
2720          pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2721          !
2722          !jyg<
2723          ! If convective tendencies are too large, then call convection
2724          !  every time step
2725          cvpas = cvpas_0
2726          DO k=1,k_upper_cv
2727             DO i=1,klon
2728               IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
2729                   d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
2730                     dtcon_multistep_max = 3.
2731                     dqcon_multistep_max = 0.02
2732               ENDIF
2733             ENDDO
2734          ENDDO
2735!
2736          DO k=1,k_upper_cv
2737             DO i=1,klon
2738!!               IF (abs(d_t_con(i,k)) > 0.24 .OR. &
2739!!                   abs(d_q_con(i,k)) > 2.e-2) THEN
2740               IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
2741                   abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
2742                 cvpas = 1
2743!!                 print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
2744!!                                   i,k,d_t_con(i,k),d_q_con(i,k)
2745               ENDIF
2746             ENDDO
2747          ENDDO
2748!!!   Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
2749!!!          call bcast(cvpas)
2750!!!   ------------------------------------------------------------
2751          !>jyg
2752          !
2753          DO i = 1, klon
2754             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
2755          ENDDO
2756          !
2757          !jyg<
2758          !    Add the tendency due to the dry adjustment of the wake profile
2759          IF (iflag_wake>=1) THEN
2760            IF (iflag_adjwk == 2) THEN
2761              DO k=1,klev
2762                 DO i=1,klon
2763                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
2764                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
2765                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2766                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2767                 ENDDO
2768              ENDDO
2769            ENDIF  ! (iflag_adjwk = 2)
2770          ENDIF   ! (iflag_wake>=1)
2771          !>jyg
2772          !
2773       ELSE ! ok_cvl
2774
2775          ! MAF conema3 ne contient pas les traceurs
2776          CALL conema3 (dtime, &
2777               paprs,pplay,t_seri,q_seri, &
2778               u_seri,v_seri,tr_seri,ntra, &
2779               sig1,w01, &
2780               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2781               rain_con, snow_con, ibas_con, itop_con, &
2782               upwd,dnwd,dnwd0,bas,top, &
2783               Ma,cape,tvp,rflag, &
2784               pbase &
2785               ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2786               ,clwcon0)
2787
2788       ENDIF ! ok_cvl
2789
2790       !
2791       ! Correction precip
2792       rain_con = rain_con * cvl_corr
2793       snow_con = snow_con * cvl_corr
2794       !
2795
2796       IF (.NOT. ok_gust) THEN
2797          do i = 1, klon
2798             wd(i)=0.0
2799          enddo
2800       ENDIF
2801
2802       ! =================================================================== c
2803       ! Calcul des proprietes des nuages convectifs
2804       !
2805
2806       !   calcul des proprietes des nuages convectifs
2807       clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2808       IF (iflag_cld_cv == 0) THEN
2809          CALL clouds_gno &
2810               (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2811       ELSE
2812          CALL clouds_bigauss &
2813               (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2814       ENDIF
2815
2816
2817       ! =================================================================== c
2818
2819       DO i = 1, klon
2820          itop_con(i) = min(max(itop_con(i),1),klev)
2821          ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2822       ENDDO
2823
2824       DO i = 1, klon
2825          ema_pcb(i)  = paprs(i,ibas_con(i))
2826       ENDDO
2827       DO i = 1, klon
2828          ! L'idicage de itop_con peut cacher un pb potentiel
2829          ! FH sous la dictee de JYG, CR
2830          ema_pct(i)  = paprs(i,itop_con(i)+1)
2831
2832          IF (itop_con(i).gt.klev-3) THEN
2833             IF (prt_level >= 9) THEN
2834                write(lunout,*)'La convection monte trop haut '
2835                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2836             ENDIF
2837          ENDIF
2838       ENDDO
2839    ELSE IF (iflag_con.eq.0) THEN
2840       write(lunout,*) 'On n appelle pas la convection'
2841       clwcon0=0.
2842       rnebcon0=0.
2843       d_t_con=0.
2844       d_q_con=0.
2845       d_u_con=0.
2846       d_v_con=0.
2847       rain_con=0.
2848       snow_con=0.
2849       bas=1
2850       top=1
2851    ELSE
2852       WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2853       CALL abort_physic("physiq", "", 1)
2854    ENDIF
2855
2856    !     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2857    !    .              d_u_con, d_v_con)
2858
2859!jyg    Reinitialize proba_notrig and itapcv when convection has been called
2860    proba_notrig(:) = 1.
2861    itapcv = 0
2862    ENDIF !  (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
2863!
2864    itapcv = itapcv+1
2865    !
2866    ! Compter les steps ou cvpas=1
2867    IF (cvpas == 1) THEN
2868      Ncvpaseq1 = Ncvpaseq1+1
2869    ENDIF
2870    IF (mod(itap,1000) == 0) THEN
2871      print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
2872    ENDIF
2873
2874!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de
2875!!!     l'energie dans les courants satures.
2876!!    d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
2877!!    d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
2878!!    dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
2879!!    CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat,   &
2880!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,& 
2881!!                     itap, 1)
2882!!    call prt_enerbil('convection_sat',itap)
2883!!
2884!!
2885    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2886         'convection',abortphy,flag_inhib_tend,itap,0)
2887    call prt_enerbil('convection',itap)
2888
2889    !-------------------------------------------------------------------------
2890
2891    IF (mydebug) THEN
2892       CALL writefield_phy('u_seri',u_seri,nbp_lev)
2893       CALL writefield_phy('v_seri',v_seri,nbp_lev)
2894       CALL writefield_phy('t_seri',t_seri,nbp_lev)
2895       CALL writefield_phy('q_seri',q_seri,nbp_lev)
2896    ENDIF
2897
2898    IF (check) THEN
2899       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2900       WRITE(lunout,*)"aprescon=", za
2901       zx_t = 0.0
2902       za = 0.0
2903       DO i = 1, klon
2904          za = za + cell_area(i)/REAL(klon)
2905          zx_t = zx_t + (rain_con(i)+ &
2906               snow_con(i))*cell_area(i)/REAL(klon)
2907       ENDDO
2908       zx_t = zx_t/za*dtime
2909       WRITE(lunout,*)"Precip=", zx_t
2910    ENDIF
2911    IF (zx_ajustq) THEN
2912       DO i = 1, klon
2913          z_apres(i) = 0.0
2914       ENDDO
2915       DO k = 1, klev
2916          DO i = 1, klon
2917             z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2918                  *(paprs(i,k)-paprs(i,k+1))/RG
2919          ENDDO
2920       ENDDO
2921       DO i = 1, klon
2922          z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2923               /z_apres(i)
2924       ENDDO
2925       DO k = 1, klev
2926          DO i = 1, klon
2927             IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2928                  z_factor(i).LT.(1.0-1.0E-08)) THEN
2929                q_seri(i,k) = q_seri(i,k) * z_factor(i)
2930             ENDIF
2931          ENDDO
2932       ENDDO
2933    ENDIF
2934    zx_ajustq=.FALSE.
2935
2936    !
2937    !==========================================================================
2938    !RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2939    !pour la couche limite diffuse pour l instant
2940    !
2941    !
2942    ! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
2943    ! il faut rajouter cette tendance calcul\'ee hors des poches
2944    ! froides
2945    !
2946    IF (iflag_wake>=1) THEN
2947       !
2948       !
2949       ! Call wakes every "wkpas" step
2950       !
2951       IF (MOD(itapwk,wkpas).EQ.0) THEN
2952          !
2953          DO k=1,klev
2954             DO i=1,klon
2955                dt_dwn(i,k)  = ftd(i,k)
2956                dq_dwn(i,k)  = fqd(i,k)
2957                M_dwn(i,k)   = dnwd0(i,k)
2958                M_up(i,k)    = upwd(i,k)
2959                dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
2960                dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
2961             ENDDO
2962          ENDDO
2963         
2964          IF (iflag_wake==2) THEN
2965             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2966             DO k = 1,klev
2967                dt_dwn(:,k)= dt_dwn(:,k)+ &
2968                     ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2969                dq_dwn(:,k)= dq_dwn(:,k)+ &
2970                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2971             ENDDO
2972          ELSEIF (iflag_wake==3) THEN
2973             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2974             DO k = 1,klev
2975                DO i=1,klon
2976                   IF (rneb(i,k)==0.) THEN
2977                      ! On ne tient compte des tendances qu'en dehors des
2978                      ! nuages (c'est-\`a-dire a priri dans une region ou
2979                      ! l'eau se reevapore).
2980                      dt_dwn(i,k)= dt_dwn(i,k)+ &
2981                           ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2982                      dq_dwn(i,k)= dq_dwn(i,k)+ &
2983                           ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2984                   ENDIF
2985                ENDDO
2986             ENDDO
2987          ENDIF
2988         
2989          !
2990          !calcul caracteristiques de la poche froide
2991          CALL calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
2992               t_seri, q_seri, omega,  &
2993               dt_dwn, dq_dwn, M_dwn, M_up,  &
2994               dt_a, dq_a,  &
2995               sigd,  &
2996               wake_deltat, wake_deltaq, wake_s, wake_dens,  &
2997               wake_dth, wake_h,  &
2998!!               wake_pe, wake_fip, wake_gfl,  &
2999               wake_pe, wake_fip_0, wake_gfl,  &   !! jyg
3000               d_t_wake, d_q_wake,  &
3001               wake_k, t_x, q_x,  &
3002               wake_omgbdth, wake_dp_omgb,  &
3003               wake_dtKE, wake_dqKE,  &
3004               wake_omg, wake_dp_deltomg,  &
3005               wake_spread, wake_Cstar, d_deltat_wk_gw,  &
3006               d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
3007          !
3008          !jyg    Reinitialize itapwk when wakes have been called
3009          itapwk = 0
3010       ENDIF !  (MOD(itapwk,wkpas).EQ.0)
3011       !
3012       itapwk = itapwk+1
3013       !
3014       !-----------------------------------------------------------------------
3015       ! ajout des tendances des poches froides
3016       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
3017            abortphy,flag_inhib_tend,itap,0)
3018       call prt_enerbil('wake',itap)
3019       !------------------------------------------------------------------------
3020
3021       ! Increment Wake state variables
3022       IF (iflag_wake_tend .GT. 0.) THEN
3023
3024         CALL add_wake_tend &
3025            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
3026             'wake', abortphy)
3027          call prt_enerbil('wake',itap)
3028       ENDIF   ! (iflag_wake_tend .GT. 0.)
3029
3030       IF (iflag_alp_wk_cond .GT. 0.) THEN
3031
3032         CALL alpale_wk(dtime, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3033                        wake_fip)
3034       ELSE
3035         wake_fip(:) = wake_fip_0(:)
3036       ENDIF   ! (iflag_alp_wk_cond .GT. 0.)
3037
3038    ENDIF  ! (iflag_wake>=1)
3039    !
3040    !===================================================================
3041    ! Convection seche (thermiques ou ajustement)
3042    !===================================================================
3043    !
3044    CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3045         ,seuil_inversion,weak_inversion,dthmin)
3046
3047
3048
3049    d_t_ajsb(:,:)=0.
3050    d_q_ajsb(:,:)=0.
3051    d_t_ajs(:,:)=0.
3052    d_u_ajs(:,:)=0.
3053    d_v_ajs(:,:)=0.
3054    d_q_ajs(:,:)=0.
3055    clwcon0th(:,:)=0.
3056    !
3057    !      fm_therm(:,:)=0.
3058    !      entr_therm(:,:)=0.
3059    !      detr_therm(:,:)=0.
3060    !
3061    IF (prt_level>9) WRITE(lunout,*) &
3062         'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3063         ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
3064    IF (iflag_thermals<0) THEN
3065       !  Rien
3066       !  ====
3067       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3068
3069
3070    ELSE
3071
3072       !  Thermiques
3073       !  ==========
3074       IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3075            ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
3076
3077
3078       !cc nrlmd le 10/04/2012
3079       DO k=1,klev+1
3080          DO i=1,klon
3081             pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3082             pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3083             pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3084             pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3085          ENDDO
3086       ENDDO
3087       !cc fin nrlmd le 10/04/2012
3088
3089       IF (iflag_thermals>=1) THEN
3090          !jyg<
3091!!       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3092       IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3093             !  Appel des thermiques avec les profils exterieurs aux poches
3094             DO k=1,klev
3095                DO i=1,klon
3096                   t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3097                   q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3098                   u_therm(i,k) = u_seri(i,k)
3099                   v_therm(i,k) = v_seri(i,k)
3100                ENDDO
3101             ENDDO
3102          ELSE
3103             !  Appel des thermiques avec les profils moyens
3104             DO k=1,klev
3105                DO i=1,klon
3106                   t_therm(i,k) = t_seri(i,k)
3107                   q_therm(i,k) = q_seri(i,k)
3108                   u_therm(i,k) = u_seri(i,k)
3109                   v_therm(i,k) = v_seri(i,k)
3110                ENDDO
3111             ENDDO
3112          ENDIF
3113          !>jyg
3114          CALL calltherm(pdtphys &
3115               ,pplay,paprs,pphi,weak_inversion &
3116                        ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3117               ,u_therm,v_therm,t_therm,q_therm,zqsat,debut &  !jyg
3118               ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3119               ,fm_therm,entr_therm,detr_therm &
3120               ,zqasc,clwcon0th,lmax_th,ratqscth &
3121               ,ratqsdiff,zqsatth &
3122                                !on rajoute ale et alp, et les
3123                                !caracteristiques de la couche alim
3124               ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3125               ,ztv,zpspsk,ztla,zthl &
3126                                !cc nrlmd le 10/04/2012
3127               ,pbl_tke_input,pctsrf,omega,cell_area &
3128               ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3129               ,n2,s2,ale_bl_stat &
3130               ,therm_tke_max,env_tke_max &
3131               ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3132               ,alp_bl_conv,alp_bl_stat &
3133                                !cc fin nrlmd le 10/04/2012
3134               ,zqla,ztva )
3135          !
3136          !jyg<
3137!!jyg          IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3138          IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3139             !  Si les thermiques ne sont presents que hors des
3140             !  poches, la tendance moyenne associ\'ee doit etre
3141             !  multipliee par la fraction surfacique qu'ils couvrent.
3142             DO k=1,klev
3143                DO i=1,klon
3144                   !
3145                   d_deltat_the(i,k) = - d_t_ajs(i,k)
3146                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
3147                   !
3148                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3149                   d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3150                   d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3151                   d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3152                   !
3153                ENDDO
3154             ENDDO
3155          !
3156             CALL add_wake_tend &
3157                 (d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
3158             call prt_enerbil('the',itap)
3159          !
3160          ENDIF  ! (mod(iflag_pbl_split/10,10) .GE. 1)
3161          !
3162          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
3163                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3164          call prt_enerbil('thermals',itap)
3165          !
3166!
3167          CALL alpale_th( dtime, lmax_th, t_seri, cell_area,  &
3168                          cin, s2, n2,  &
3169                          ale_bl_trig, ale_bl_stat, ale_bl,  &
3170                          alp_bl, alp_bl_stat, &
3171                          proba_notrig, random_notrig)
3172          !>jyg
3173
3174          ! ------------------------------------------------------------------
3175          ! Transport de la TKE par les panaches thermiques.
3176          ! FH : 2010/02/01
3177          !     if (iflag_pbl.eq.10) then
3178          !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3179          !    s           rg,paprs,pbl_tke)
3180          !     endif
3181          ! -------------------------------------------------------------------
3182
3183          DO i=1,klon
3184             !           zmax_th(i)=pphi(i,lmax_th(i))/rg
3185             !CR:04/05/12:correction calcul zmax
3186             zmax_th(i)=zmax0(i)
3187          ENDDO
3188
3189       ENDIF
3190
3191       !  Ajustement sec
3192       !  ==============
3193
3194       ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3195       ! a partir du sommet des thermiques.
3196       ! Dans le cas contraire, on demarre au niveau 1.
3197
3198       IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3199
3200          IF (iflag_thermals.eq.0) THEN
3201             IF (prt_level>9) WRITE(lunout,*)'ajsec'
3202             limbas(:)=1
3203          ELSE
3204             limbas(:)=lmax_th(:)
3205          ENDIF
3206
3207          ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3208          ! pour des test de convergence numerique.
3209          ! Le nouveau ajsec est a priori mieux, meme pour le cas
3210          ! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3211          ! non nulles numeriquement pour des mailles non concernees.
3212
3213          IF (iflag_thermals==0) THEN
3214             ! Calling adjustment alone (but not the thermal plume model)
3215             CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3216                  , d_t_ajsb, d_q_ajsb)
3217          ELSE IF (iflag_thermals>0) THEN
3218             ! Calling adjustment above the top of thermal plumes
3219             CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3220                  , d_t_ajsb, d_q_ajsb)
3221          ENDIF
3222
3223          !--------------------------------------------------------------------
3224          ! ajout des tendances de l'ajustement sec ou des thermiques
3225          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3226               'ajsb',abortphy,flag_inhib_tend,itap,0)
3227          call prt_enerbil('ajsb',itap)
3228          d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3229          d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3230
3231          !---------------------------------------------------------------------
3232
3233       ENDIF
3234
3235    ENDIF
3236    !
3237    !===================================================================
3238    ! Computation of ratqs, the width (normalized) of the subrid scale
3239    ! water distribution
3240    CALL  calcratqs(klon,klev,prt_level,lunout,        &
3241         iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
3242         ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3243         tau_ratqs,fact_cldcon,   &
3244         ptconv,ptconvth,clwcon0th, rnebcon0th,     &
3245         paprs,pplay,q_seri,zqsat,fm_therm, &
3246         ratqs,ratqsc)
3247
3248
3249    !
3250    ! Appeler le processus de condensation a grande echelle
3251    ! et le processus de precipitation
3252    !-------------------------------------------------------------------------
3253    IF (prt_level .GE.10) THEN
3254       print *,'itap, ->fisrtilp ',itap
3255    ENDIF
3256    !
3257    CALL fisrtilp(dtime,paprs,pplay, &
3258         t_seri, q_seri,ptconv,ratqs, &
3259         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3260         rain_lsc, snow_lsc, &
3261         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3262         frac_impa, frac_nucl, beta_prec_fisrt, &
3263         prfl, psfl, rhcl,  &
3264         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3265         iflag_ice_thermo)
3266    !
3267    WHERE (rain_lsc < 0) rain_lsc = 0.
3268    WHERE (snow_lsc < 0) snow_lsc = 0.
3269
3270!+JLD
3271!    write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3272!        & ,((rcw-rcpd)*rain_lsc + (rcs-rcpd)*snow_lsc)*t_seri(1,1)-rlvtt*rain_lsc+rlstt*snow_lsc &
3273!        & ,rain_lsc,snow_lsc
3274!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3275!-JLD
3276    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3277         'lsc',abortphy,flag_inhib_tend,itap,0)
3278    call prt_enerbil('lsc',itap)
3279    rain_num(:)=0.
3280    DO k = 1, klev
3281       DO i = 1, klon
3282          IF (ql_seri(i,k)>oliqmax) THEN
3283             rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3284             ql_seri(i,k)=oliqmax
3285          ENDIF
3286       ENDDO
3287    ENDDO
3288    IF (nqo==3) THEN
3289    DO k = 1, klev
3290       DO i = 1, klon
3291          IF (qs_seri(i,k)>oicemax) THEN
3292             rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3293             qs_seri(i,k)=oicemax
3294          ENDIF
3295       ENDDO
3296    ENDDO
3297    ENDIF
3298
3299    !---------------------------------------------------------------------------
3300    DO k = 1, klev
3301       DO i = 1, klon
3302          cldfra(i,k) = rneb(i,k)
3303          !CR: a quoi ca sert? Faut-il ajouter qs_seri?
3304          IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3305       ENDDO
3306    ENDDO
3307    IF (check) THEN
3308       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3309       WRITE(lunout,*)"apresilp=", za
3310       zx_t = 0.0
3311       za = 0.0
3312       DO i = 1, klon
3313          za = za + cell_area(i)/REAL(klon)
3314          zx_t = zx_t + (rain_lsc(i) &
3315               + snow_lsc(i))*cell_area(i)/REAL(klon)
3316       ENDDO
3317       zx_t = zx_t/za*dtime
3318       WRITE(lunout,*)"Precip=", zx_t
3319    ENDIF
3320
3321    IF (mydebug) THEN
3322       CALL writefield_phy('u_seri',u_seri,nbp_lev)
3323       CALL writefield_phy('v_seri',v_seri,nbp_lev)
3324       CALL writefield_phy('t_seri',t_seri,nbp_lev)
3325       CALL writefield_phy('q_seri',q_seri,nbp_lev)
3326    ENDIF
3327
3328    !
3329    !-------------------------------------------------------------------
3330    !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3331    !-------------------------------------------------------------------
3332
3333    ! 1. NUAGES CONVECTIFS
3334    !
3335    !IM cf FH
3336    !     IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3337    IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3338       snow_tiedtke=0.
3339       !     print*,'avant calcul de la pseudo precip '
3340       !     print*,'iflag_cld_th',iflag_cld_th
3341       IF (iflag_cld_th.eq.-1) THEN
3342          rain_tiedtke=rain_con
3343       ELSE
3344          !       print*,'calcul de la pseudo precip '
3345          rain_tiedtke=0.
3346          !         print*,'calcul de la pseudo precip 0'
3347          DO k=1,klev
3348             DO i=1,klon
3349                IF (d_q_con(i,k).lt.0.) THEN
3350                   rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3351                        *(paprs(i,k)-paprs(i,k+1))/rg
3352                ENDIF
3353             ENDDO
3354          ENDDO
3355       ENDIF
3356       !
3357       !     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3358       !
3359
3360       ! Nuages diagnostiques pour Tiedtke
3361       CALL diagcld1(paprs,pplay, &
3362                                !IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3363            rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3364            diafra,dialiq)
3365       DO k = 1, klev
3366          DO i = 1, klon
3367             IF (diafra(i,k).GT.cldfra(i,k)) THEN
3368                cldliq(i,k) = dialiq(i,k)
3369                cldfra(i,k) = diafra(i,k)
3370             ENDIF
3371          ENDDO
3372       ENDDO
3373
3374    ELSE IF (iflag_cld_th.ge.3) THEN
3375       !  On prend pour les nuages convectifs le max du calcul de la
3376       !  convection et du calcul du pas de temps precedent diminue d'un facteur
3377       !  facttemps
3378       facteur = pdtphys *facttemps
3379       DO k=1,klev
3380          DO i=1,klon
3381             rnebcon(i,k)=rnebcon(i,k)*facteur
3382             IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
3383                rnebcon(i,k)=rnebcon0(i,k)
3384                clwcon(i,k)=clwcon0(i,k)
3385             ENDIF
3386          ENDDO
3387       ENDDO
3388
3389       !   On prend la somme des fractions nuageuses et des contenus en eau
3390
3391       IF (iflag_cld_th>=5) THEN
3392
3393          DO k=1,klev
3394             ptconvth(:,k)=fm_therm(:,k+1)>0.
3395          ENDDO
3396
3397          IF (iflag_coupl==4) THEN
3398
3399             ! Dans le cas iflag_coupl==4, on prend la somme des convertures
3400             ! convectives et lsc dans la partie des thermiques
3401             ! Le controle par iflag_coupl est peut etre provisoire.
3402             DO k=1,klev
3403                DO i=1,klon
3404                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3405                      cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3406                      cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3407                   ELSE IF (ptconv(i,k)) THEN
3408                      cldfra(i,k)=rnebcon(i,k)
3409                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3410                   ENDIF
3411                ENDDO
3412             ENDDO
3413
3414          ELSE IF (iflag_coupl==5) THEN
3415             DO k=1,klev
3416                DO i=1,klon
3417                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3418                   cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3419                ENDDO
3420             ENDDO
3421
3422          ELSE
3423
3424             ! Si on est sur un point touche par la convection
3425             ! profonde et pas par les thermiques, on prend la
3426             ! couverture nuageuse et l'eau nuageuse de la convection
3427             ! profonde.
3428
3429             !IM/FH: 2011/02/23
3430             ! definition des points sur lesquels ls thermiques sont actifs
3431
3432             DO k=1,klev
3433                DO i=1,klon
3434                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3435                      cldfra(i,k)=rnebcon(i,k)
3436                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3437                   ENDIF
3438                ENDDO
3439             ENDDO
3440
3441          ENDIF
3442
3443       ELSE
3444
3445          ! Ancienne version
3446          cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3447          cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3448       ENDIF
3449
3450    ENDIF
3451
3452    !     plulsc(:)=0.
3453    !     do k=1,klev,-1
3454    !        do i=1,klon
3455    !              zzz=prfl(:,k)+psfl(:,k)
3456    !           if (.not.ptconvth.zzz.gt.0.)
3457    !        enddo prfl, psfl,
3458    !     enddo
3459    !
3460    ! 2. NUAGES STARTIFORMES
3461    !
3462    IF (ok_stratus) THEN
3463       CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3464       DO k = 1, klev
3465          DO i = 1, klon
3466             IF (diafra(i,k).GT.cldfra(i,k)) THEN
3467                cldliq(i,k) = dialiq(i,k)
3468                cldfra(i,k) = diafra(i,k)
3469             ENDIF
3470          ENDDO
3471       ENDDO
3472    ENDIF
3473    !
3474    ! Precipitation totale
3475    !
3476    DO i = 1, klon
3477       rain_fall(i) = rain_con(i) + rain_lsc(i)
3478       snow_fall(i) = snow_con(i) + snow_lsc(i)
3479    ENDDO
3480    !
3481    ! Calculer l'humidite relative pour diagnostique
3482    !
3483    DO k = 1, klev
3484       DO i = 1, klon
3485          zx_t = t_seri(i,k)
3486          IF (thermcep) THEN
3487             !!           if (iflag_ice_thermo.eq.0) then                 !jyg
3488             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3489             !!           else                                            !jyg
3490             !!           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))      !jyg
3491             !!           endif                                           !jyg
3492             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3493             zx_qs  = MIN(0.5,zx_qs)
3494             zcor   = 1./(1.-retv*zx_qs)
3495             zx_qs  = zx_qs*zcor
3496          ELSE
3497             !!           IF (zx_t.LT.t_coup) THEN             !jyg
3498             IF (zx_t.LT.rtt) THEN                  !jyg
3499                zx_qs = qsats(zx_t)/pplay(i,k)
3500             ELSE
3501                zx_qs = qsatl(zx_t)/pplay(i,k)
3502             ENDIF
3503          ENDIF
3504          zx_rh(i,k) = q_seri(i,k)/zx_qs
3505          zqsat(i,k)=zx_qs
3506       ENDDO
3507    ENDDO
3508
3509    !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3510    !   equivalente a 2m (tpote) pour diagnostique
3511    !
3512    DO i = 1, klon
3513       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
3514       IF (thermcep) THEN
3515          IF(zt2m(i).LT.RTT) then
3516             Lheat=RLSTT
3517          ELSE
3518             Lheat=RLVTT
3519          ENDIF
3520       ELSE
3521          IF (zt2m(i).LT.RTT) THEN
3522             Lheat=RLSTT
3523          ELSE
3524             Lheat=RLVTT
3525          ENDIF
3526       ENDIF
3527       tpote(i) = tpot(i)*      &
3528            EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
3529    ENDDO
3530
3531    IF (type_trac == 'inca') THEN
3532#ifdef INCA
3533       CALL VTe(VTphysiq)
3534       CALL VTb(VTinca)
3535       calday = REAL(days_elapsed + 1) + jH_cur
3536
3537       CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
3538       CALL aerosol_meteo_calc( &
3539            calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3540            prfl,psfl,pctsrf,cell_area, &
3541            latitude_deg,longitude_deg,u10m,v10m)
3542
3543       zxsnow_dummy(:) = 0.0
3544
3545       CALL chemhook_begin (calday, &
3546            days_elapsed+1, &
3547            jH_cur, &
3548            pctsrf(1,1), &
3549            latitude_deg, &
3550            longitude_deg, &
3551            cell_area, &
3552            paprs, &
3553            pplay, &
3554            coefh(1:klon,1:klev,is_ave), &
3555            pphi, &
3556            t_seri, &
3557            u, &
3558            v, &
3559            wo(:, :, 1), &
3560            q_seri, &
3561            zxtsol, &
3562            zxsnow_dummy, &
3563            solsw, &
3564            albsol1, &
3565            rain_fall, &
3566            snow_fall, &
3567            itop_con, &
3568            ibas_con, &
3569            cldfra, &
3570            nbp_lon, &
3571            nbp_lat-1, &
3572            tr_seri, &
3573            ftsol, &
3574            paprs, &
3575            cdragh, &
3576            cdragm, &
3577            pctsrf, &
3578            pdtphys, &
3579            itap)
3580
3581       CALL VTe(VTinca)
3582       CALL VTb(VTphysiq)
3583#endif
3584    ENDIF !type_trac = inca
3585
3586
3587    !
3588    ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3589    !
3590    IF (MOD(itaprad,radpas).EQ.0) THEN
3591
3592       !
3593       !jq - introduce the aerosol direct and first indirect radiative forcings
3594       !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3595       IF (flag_aerosol .GT. 0) THEN
3596          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3597             IF (.NOT. aerosol_couple) THEN
3598                !
3599                CALL readaerosol_optic( &
3600                     debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3601                     pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
3602                     mass_solu_aero, mass_solu_aero_pi,  &
3603                     tau_aero, piz_aero, cg_aero,  &
3604                     tausum_aero, tau3d_aero)
3605             ENDIF
3606          ELSE                       ! RRTM radiation
3607             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3608                abort_message='config_inca=aero et rrtm=1 impossible'
3609                CALL abort_physic(modname,abort_message,1)
3610             ELSE
3611                !
3612#ifdef CPP_RRTM
3613                IF (NSW.EQ.6) THEN
3614                   !--new aerosol properties SW and LW
3615                   !
3616#ifdef CPP_Dust
3617                   !--SPL aerosol model
3618                   CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3619                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
3620                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
3621                        tausum_aero, tau3d_aero)
3622#else
3623                   !--climatologies or INCA aerosols
3624                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
3625                        new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3626                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
3627                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
3628                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
3629                        tausum_aero, drytausum_aero, tau3d_aero)
3630#endif
3631
3632                   IF (flag_aerosol .EQ. 7) THEN
3633                      CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
3634                                   tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3635                   ENDIF
3636
3637                   !
3638                ELSE IF (NSW.EQ.2) THEN
3639                   !--for now we use the old aerosol properties
3640                   !
3641                   CALL readaerosol_optic( &
3642                        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3643                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
3644                        mass_solu_aero, mass_solu_aero_pi,  &
3645                        tau_aero, piz_aero, cg_aero,  &
3646                        tausum_aero, tau3d_aero)
3647                   !
3648                   !--natural aerosols
3649                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3650                   piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3651                   cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3652                   !--all aerosols
3653                   tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3654                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3655                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3656                   !
3657                   !--no LW optics
3658                   tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3659                   !
3660                ELSE
3661                   abort_message='Only NSW=2 or 6 are possible with ' &
3662                        // 'aerosols and iflag_rrtm=1'
3663                   CALL abort_physic(modname,abort_message,1)
3664                ENDIF
3665#else
3666                abort_message='You should compile with -rrtm if running ' &
3667                     // 'with iflag_rrtm=1'
3668                CALL abort_physic(modname,abort_message,1)
3669#endif
3670                !
3671             ENDIF
3672          ENDIF
3673       ELSE   !--flag_aerosol = 0
3674          tausum_aero(:,:,:) = 0.
3675          drytausum_aero(:,:) = 0.
3676          mass_solu_aero(:,:) = 0.
3677          mass_solu_aero_pi(:,:) = 0.
3678          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3679             tau_aero(:,:,:,:) = 1.e-15
3680             piz_aero(:,:,:,:) = 1.
3681             cg_aero(:,:,:,:)  = 0.
3682          ELSE
3683             tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3684             tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3685             piz_aero_sw_rrtm(:,:,:,:) = 1.0
3686             cg_aero_sw_rrtm(:,:,:,:)  = 0.0
3687          ENDIF
3688       ENDIF
3689       !
3690       !--WMO criterion to determine tropopause
3691       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3692       !
3693       !--STRAT AEROSOL
3694       !--updates tausum_aero,tau_aero,piz_aero,cg_aero
3695       IF (flag_aerosol_strat.GT.0) THEN
3696          IF (prt_level .GE.10) THEN
3697             PRINT *,'appel a readaerosolstrat', mth_cur
3698          ENDIF
3699          IF (iflag_rrtm.EQ.0) THEN
3700           IF (flag_aerosol_strat.EQ.1) THEN
3701             CALL readaerosolstrato(debut)
3702           ELSE
3703             abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3704             CALL abort_physic(modname,abort_message,1)
3705           ENDIF
3706          ELSE
3707#ifdef CPP_RRTM
3708#ifndef CPP_StratAer
3709          !--prescribed strat aerosols
3710          !--only in the case of non-interactive strat aerosols
3711            IF (flag_aerosol_strat.EQ.1) THEN
3712             CALL readaerosolstrato1_rrtm(debut)
3713            ELSEIF (flag_aerosol_strat.EQ.2) THEN
3714             CALL readaerosolstrato2_rrtm(debut, ok_volcan)
3715            ELSE
3716             abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3717             CALL abort_physic(modname,abort_message,1)
3718            ENDIF
3719#endif
3720#else
3721             abort_message='You should compile with -rrtm if running ' &
3722                  // 'with iflag_rrtm=1'
3723             CALL abort_physic(modname,abort_message,1)
3724#endif
3725          ENDIF
3726       ENDIF
3727!
3728#ifdef CPP_RRTM
3729#ifdef CPP_StratAer
3730       !--compute stratospheric mask
3731       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3732       !--interactive strat aerosols
3733       CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3734#endif
3735#endif
3736       !--fin STRAT AEROSOL
3737       !     
3738
3739       ! Calculer les parametres optiques des nuages et quelques
3740       ! parametres pour diagnostiques:
3741       !
3742       IF (aerosol_couple.AND.config_inca=='aero') THEN
3743          mass_solu_aero(:,:)    = ccm(:,:,1)
3744          mass_solu_aero_pi(:,:) = ccm(:,:,2)
3745       ENDIF
3746
3747       IF (ok_newmicro) then
3748          IF (iflag_rrtm.NE.0) THEN
3749#ifdef CPP_RRTM
3750             IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3751             abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3752                  // 'pour ok_cdnc'
3753             CALL abort_physic(modname,abort_message,1)
3754             ENDIF
3755#else
3756
3757             abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3758             CALL abort_physic(modname,abort_message,1)
3759#endif
3760          ENDIF
3761          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
3762               paprs, pplay, t_seri, cldliq, cldfra, &
3763               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3764               flwp, fiwp, flwc, fiwc, &
3765               mass_solu_aero, mass_solu_aero_pi, &
3766               cldtaupi, re, fl, ref_liq, ref_ice, &
3767               ref_liq_pi, ref_ice_pi)
3768       ELSE
3769          CALL nuage (paprs, pplay, &
3770               t_seri, cldliq, cldfra, cldtau, cldemi, &
3771               cldh, cldl, cldm, cldt, cldq, &
3772               ok_aie, &
3773               mass_solu_aero, mass_solu_aero_pi, &
3774               bl95_b0, bl95_b1, &
3775               cldtaupi, re, fl)
3776       ENDIF
3777       !
3778       !IM betaCRF
3779       !
3780       cldtaurad   = cldtau
3781       cldtaupirad = cldtaupi
3782       cldemirad   = cldemi
3783       cldfrarad   = cldfra
3784
3785       !
3786       IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3787           lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3788          !
3789          ! global
3790          !
3791!IM 251017 begin
3792!                print*,'physiq betaCRF global zdtime=',zdtime
3793!IM 251017 end
3794          DO k=1, klev
3795             DO i=1, klon
3796                IF (pplay(i,k).GE.pfree) THEN
3797                   beta(i,k) = beta_pbl
3798                ELSE
3799                   beta(i,k) = beta_free
3800                ENDIF
3801                IF (mskocean_beta) THEN
3802                   beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3803                ENDIF
3804                cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3805                cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3806                cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3807                cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3808             ENDDO
3809          ENDDO
3810          !
3811       ELSE
3812          !
3813          ! regional
3814          !
3815          DO k=1, klev
3816             DO i=1,klon
3817                !
3818                IF (longitude_deg(i).ge.lon1_beta.AND. &
3819                    longitude_deg(i).le.lon2_beta.AND. &
3820                    latitude_deg(i).le.lat1_beta.AND.  &
3821                    latitude_deg(i).ge.lat2_beta) THEN
3822                   IF (pplay(i,k).GE.pfree) THEN
3823                      beta(i,k) = beta_pbl
3824                   ELSE
3825                      beta(i,k) = beta_free
3826                   ENDIF
3827                   IF (mskocean_beta) THEN
3828                      beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3829                   ENDIF
3830                   cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3831                   cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3832                   cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3833                   cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3834                ENDIF
3835             !
3836             ENDDO
3837          ENDDO
3838       !
3839       ENDIF
3840
3841       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3842       IF (ok_chlorophyll) THEN
3843          print*,"-- reading chlorophyll"
3844          CALL readchlorophyll(debut)
3845       ENDIF
3846
3847!--if ok_suntime_rrtm we use ancillay data for RSUN
3848!--previous values are therefore overwritten
3849!--this is needed for CMIP6 runs
3850!--and only possible for new radiation scheme
3851       IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
3852#ifdef CPP_RRTM
3853         CALL read_rsun_rrtm(debut)
3854#endif
3855       ENDIF
3856
3857       IF (mydebug) THEN
3858          CALL writefield_phy('u_seri',u_seri,nbp_lev)
3859          CALL writefield_phy('v_seri',v_seri,nbp_lev)
3860          CALL writefield_phy('t_seri',t_seri,nbp_lev)
3861          CALL writefield_phy('q_seri',q_seri,nbp_lev)
3862       ENDIF
3863
3864       !
3865       !sonia : If Iflag_radia >=2, pertubation of some variables
3866       !input to radiation (DICE)
3867       !
3868       IF (iflag_radia .ge. 2) THEN
3869          zsav_tsol (:) = zxtsol(:)
3870          CALL perturb_radlwsw(zxtsol,iflag_radia)
3871       ENDIF
3872
3873       IF (aerosol_couple.AND.config_inca=='aero') THEN
3874#ifdef INCA
3875          CALL radlwsw_inca  &
3876               (chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
3877               paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3878               size(wo,3), wo, &
3879               cldfrarad, cldemirad, cldtaurad, &
3880               heat,heat0,cool,cool0,albpla, &
3881               topsw,toplw,solsw,sollw, &
3882               sollwdown, &
3883               topsw0,toplw0,solsw0,sollw0, &
3884               lwdn0, lwdn, lwup0, lwup,  &
3885               swdn0, swdn, swup0, swup, &
3886               ok_ade, ok_aie, &
3887               tau_aero, piz_aero, cg_aero, &
3888               topswad_aero, solswad_aero, &
3889               topswad0_aero, solswad0_aero, &
3890               topsw_aero, topsw0_aero, &
3891               solsw_aero, solsw0_aero, &
3892               cldtaupirad, &
3893               topswai_aero, solswai_aero)
3894#endif
3895       ELSE
3896          !
3897          !IM calcul radiatif pour le cas actuel
3898          !
3899          RCO2 = RCO2_act
3900          RCH4 = RCH4_act
3901          RN2O = RN2O_act
3902          RCFC11 = RCFC11_act
3903          RCFC12 = RCFC12_act
3904          !
3905          IF (prt_level .GE.10) THEN
3906             print *,' ->radlwsw, number 1 '
3907          ENDIF
3908
3909          !
3910          CALL radlwsw &
3911               (dist, rmu0, fract,  &
3912                                !albedo SB >>>
3913                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
3914               paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif,  &
3915                                !albedo SB <<<
3916               t_seri,q_seri,wo, &
3917               cldfrarad, cldemirad, cldtaurad, &
3918               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
3919               flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
3920               tau_aero, piz_aero, cg_aero, &
3921               tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3922               ! Rajoute par OB pour RRTM
3923               tau_aero_lw_rrtm, &
3924               cldtaupirad,new_aod, &
3925!              zqsat, flwcrad, fiwcrad, &
3926               zqsat, flwc, fiwc, &
3927               ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3928               heat,heat0,cool,cool0,albpla, &
3929               heat_volc,cool_volc, &
3930               topsw,toplw,solsw,sollw, &
3931               sollwdown, &
3932               topsw0,toplw0,solsw0,sollw0, &
3933               lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,  &
3934               swdnc0, swdn0, swdn, swupc0, swup0, swup, &
3935               topswad_aero, solswad_aero, &
3936               topswai_aero, solswai_aero, &
3937               topswad0_aero, solswad0_aero, &
3938               topsw_aero, topsw0_aero, &
3939               solsw_aero, solsw0_aero, &
3940               topswcf_aero, solswcf_aero, &
3941                                !-C. Kleinschmitt for LW diagnostics
3942               toplwad_aero, sollwad_aero,&
3943               toplwai_aero, sollwai_aero, &
3944               toplwad0_aero, sollwad0_aero,&
3945                                !-end
3946               ZLWFT0_i, ZFLDN0, ZFLUP0, &
3947               ZSWFT0_i, ZFSDN0, ZFSUP0)
3948
3949          !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
3950          !schemes
3951          toplw = toplw + betalwoff * (toplw0 - toplw)
3952          sollw = sollw + betalwoff * (sollw0 - sollw)
3953          lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
3954          lwup = lwup + betalwoff * (lwup0 - lwup)
3955          sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
3956                        sollwdown(:))
3957          cool = cool + betalwoff * (cool0 - cool)
3958 
3959#ifndef CPP_XIOS
3960          !--OB 30/05/2016 modified 21/10/2016
3961          !--here we return swaero_diag and dryaod_diag to FALSE
3962          !--and histdef will switch it back to TRUE if necessary
3963          !--this is necessary to get the right swaero at first step
3964          !--but only in the case of no XIOS as XIOS is covered elsewhere
3965          IF (debut) swaerofree_diag = .FALSE.
3966          IF (debut) swaero_diag = .FALSE.
3967          IF (debut) dryaod_diag = .FALSE.
3968          !--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
3969          !--as for swaero_diag, see above
3970          IF (debut) ok_4xCO2atm = .FALSE.
3971
3972          !
3973          !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3974          !IM des taux doit etre different du taux actuel
3975          !IM Par defaut on a les taux perturbes egaux aux taux actuels
3976          !
3977          IF (RCO2_per.NE.RCO2_act.OR. &
3978              RCH4_per.NE.RCH4_act.OR. &
3979              RN2O_per.NE.RN2O_act.OR. &
3980              RCFC11_per.NE.RCFC11_act.OR. &
3981              RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
3982#endif
3983   !
3984          IF (ok_4xCO2atm) THEN
3985                !
3986                RCO2 = RCO2_per
3987                RCH4 = RCH4_per
3988                RN2O = RN2O_per
3989                RCFC11 = RCFC11_per
3990                RCFC12 = RCFC12_per
3991                !
3992                IF (prt_level .GE.10) THEN
3993                   print *,' ->radlwsw, number 2 '
3994                ENDIF
3995                !
3996                CALL radlwsw &
3997                     (dist, rmu0, fract,  &
3998                                !albedo SB >>>
3999                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
4000                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4001                                !albedo SB <<<
4002                     t_seri,q_seri,wo, &
4003                     cldfrarad, cldemirad, cldtaurad, &
4004                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4005                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4006                     tau_aero, piz_aero, cg_aero, &
4007                     tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4008                                ! Rajoute par OB pour RRTM
4009                     tau_aero_lw_rrtm, &
4010                     cldtaupi,new_aod, &
4011!                    zqsat, flwcrad, fiwcrad, &
4012                     zqsat, flwc, fiwc, &
4013                     ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4014                     heatp,heat0p,coolp,cool0p,albplap, &
4015                     heat_volc,cool_volc, &
4016                     topswp,toplwp,solswp,sollwp, &
4017                     sollwdownp, &
4018                     topsw0p,toplw0p,solsw0p,sollw0p, &
4019                     lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp,  &
4020                     swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4021                     topswad_aerop, solswad_aerop, &
4022                     topswai_aerop, solswai_aerop, &
4023                     topswad0_aerop, solswad0_aerop, &
4024                     topsw_aerop, topsw0_aerop, &
4025                     solsw_aerop, solsw0_aerop, &
4026                     topswcf_aerop, solswcf_aerop, &
4027                                !-C. Kleinschmitt for LW diagnostics
4028                     toplwad_aerop, sollwad_aerop,&
4029                     toplwai_aerop, sollwai_aerop, &
4030                     toplwad0_aerop, sollwad0_aerop,&
4031                                !-end
4032                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
4033                     ZSWFT0_i, ZFSDN0, ZFSUP0)
4034          endif !ok_4xCO2atm
4035       ENDIF ! aerosol_couple
4036       itaprad = 0
4037       !
4038       !  If Iflag_radia >=2, reset pertubed variables
4039       !
4040       IF (iflag_radia .ge. 2) THEN
4041          zxtsol(:) = zsav_tsol (:)
4042       ENDIF
4043    ENDIF ! MOD(itaprad,radpas)
4044    itaprad = itaprad + 1
4045
4046    IF (iflag_radia.eq.0) THEN
4047       IF (prt_level.ge.9) THEN
4048          PRINT *,'--------------------------------------------------'
4049          PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4050          PRINT *,'>>>>           heat et cool mis a zero '
4051          PRINT *,'--------------------------------------------------'
4052       ENDIF
4053       heat=0.
4054       cool=0.
4055       sollw=0.   ! MPL 01032011
4056       solsw=0.
4057       radsol=0.
4058       swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4059       swup0=0.
4060       lwup=0.
4061       lwup0=0.
4062       lwdn=0.
4063       lwdn0=0.
4064    ENDIF
4065
4066    !
4067    ! Calculer radsol a l'exterieur de radlwsw
4068    ! pour prendre en compte le cycle diurne
4069    ! recode par Olivier Boucher en sept 2015
4070    !
4071    radsol=solsw*swradcorr+sollw
4072
4073    IF (ok_4xCO2atm) THEN
4074       radsolp=solswp*swradcorr+sollwp
4075    ENDIF
4076
4077    !
4078    ! Ajouter la tendance des rayonnements (tous les pas)
4079    ! avec une correction pour le cycle diurne dans le SW
4080    !
4081
4082    DO k=1, klev
4083       d_t_swr(:,k)=swradcorr(:)*heat(:,k)*dtime/RDAY
4084       d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*dtime/RDAY
4085       d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
4086       d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
4087    ENDDO
4088
4089    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4090    call prt_enerbil('SW',itap)
4091    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4092    call prt_enerbil('LW',itap)
4093
4094    !
4095    IF (mydebug) THEN
4096       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4097       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4098       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4099       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4100    ENDIF
4101
4102    ! Calculer l'hydrologie de la surface
4103    !
4104    !      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4105    !     .            agesno, ftsol,fqsurf,fsnow, ruis)
4106    !
4107
4108    !
4109    ! Calculer le bilan du sol et la derive de temperature (couplage)
4110    !
4111    DO i = 1, klon
4112       !         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4113       ! a la demande de JLD
4114       bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4115    ENDDO
4116    !
4117    !moddeblott(jan95)
4118    ! Appeler le programme de parametrisation de l'orographie
4119    ! a l'echelle sous-maille:
4120    !
4121    IF (prt_level .GE.10) THEN
4122       print *,' call orography ? ', ok_orodr
4123    ENDIF
4124    !
4125    IF (ok_orodr) THEN
4126       !
4127       !  selection des points pour lesquels le shema est actif:
4128       igwd=0
4129       DO i=1,klon
4130          itest(i)=0
4131          !        IF ((zstd(i).gt.10.0)) THEN
4132          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4133             itest(i)=1
4134             igwd=igwd+1
4135             idx(igwd)=i
4136          ENDIF
4137       ENDDO
4138       !        igwdim=MAX(1,igwd)
4139       !
4140       IF (ok_strato) THEN
4141
4142          CALL drag_noro_strato(0,klon,klev,dtime,paprs,pplay, &
4143               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4144               igwd,idx,itest, &
4145               t_seri, u_seri, v_seri, &
4146               zulow, zvlow, zustrdr, zvstrdr, &
4147               d_t_oro, d_u_oro, d_v_oro)
4148
4149       ELSE
4150          CALL drag_noro(klon,klev,dtime,paprs,pplay, &
4151               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4152               igwd,idx,itest, &
4153               t_seri, u_seri, v_seri, &
4154               zulow, zvlow, zustrdr, zvstrdr, &
4155               d_t_oro, d_u_oro, d_v_oro)
4156       ENDIF
4157       !
4158       !  ajout des tendances
4159       !-----------------------------------------------------------------------
4160       ! ajout des tendances de la trainee de l'orographie
4161       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4162            abortphy,flag_inhib_tend,itap,0)
4163       call prt_enerbil('oro',itap)
4164       !----------------------------------------------------------------------
4165       !
4166    ENDIF ! fin de test sur ok_orodr
4167    !
4168    IF (mydebug) THEN
4169       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4170       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4171       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4172       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4173    ENDIF
4174
4175    IF (ok_orolf) THEN
4176       !
4177       !  selection des points pour lesquels le shema est actif:
4178       igwd=0
4179       DO i=1,klon
4180          itest(i)=0
4181          IF ((zpic(i)-zmea(i)).GT.100.) THEN
4182             itest(i)=1
4183             igwd=igwd+1
4184             idx(igwd)=i
4185          ENDIF
4186       ENDDO
4187       !        igwdim=MAX(1,igwd)
4188       !
4189       IF (ok_strato) THEN
4190
4191          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
4192               latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4193               igwd,idx,itest, &
4194               t_seri, u_seri, v_seri, &
4195               zulow, zvlow, zustrli, zvstrli, &
4196               d_t_lif, d_u_lif, d_v_lif               )
4197
4198       ELSE
4199          CALL lift_noro(klon,klev,dtime,paprs,pplay, &
4200               latitude_deg,zmea,zstd,zpic, &
4201               itest, &
4202               t_seri, u_seri, v_seri, &
4203               zulow, zvlow, zustrli, zvstrli, &
4204               d_t_lif, d_u_lif, d_v_lif)
4205       ENDIF
4206
4207       ! ajout des tendances de la portance de l'orographie
4208       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4209            'lif', abortphy,flag_inhib_tend,itap,0)
4210       call prt_enerbil('lif',itap)
4211    ENDIF ! fin de test sur ok_orolf
4212
4213    IF (ok_hines) then
4214       !  HINES GWD PARAMETRIZATION
4215       east_gwstress=0.
4216       west_gwstress=0.
4217       du_gwd_hines=0.
4218       dv_gwd_hines=0.
4219       CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
4220            u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4221            du_gwd_hines, dv_gwd_hines)
4222       zustr_gwd_hines=0.
4223       zvstr_gwd_hines=0.
4224       DO k = 1, klev
4225          zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
4226               * (paprs(:, k)-paprs(:, k+1))/rg
4227          zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
4228               * (paprs(:, k)-paprs(:, k+1))/rg
4229       ENDDO
4230
4231       d_t_hin(:, :)=0.
4232       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4233            dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4234       call prt_enerbil('hin',itap)
4235    ENDIF
4236
4237    IF (.not. ok_hines .and. ok_gwd_rando) then
4238       CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
4239            v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4240            dv_gwd_front, east_gwstress, west_gwstress)
4241       zustr_gwd_front=0.
4242       zvstr_gwd_front=0.
4243       DO k = 1, klev
4244          zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
4245               * (paprs(:, k)-paprs(:, k+1))/rg
4246          zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
4247               * (paprs(:, k)-paprs(:, k+1))/rg
4248       ENDDO
4249
4250       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4251            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4252       call prt_enerbil('front_gwd_rando',itap)
4253    ENDIF
4254
4255    IF (ok_gwd_rando) THEN
4256       CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
4257            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4258            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4259       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4260            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4261       call prt_enerbil('flott_gwd_rando',itap)
4262       zustr_gwd_rando=0.
4263       zvstr_gwd_rando=0.
4264       DO k = 1, klev
4265          zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
4266               * (paprs(:, k)-paprs(:, k+1))/rg
4267          zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
4268               * (paprs(:, k)-paprs(:, k+1))/rg
4269       ENDDO
4270    ENDIF
4271
4272    ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4273
4274    IF (mydebug) THEN
4275       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4276       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4277       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4278       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4279    ENDIF
4280
4281    DO i = 1, klon
4282       zustrph(i)=0.
4283       zvstrph(i)=0.
4284    ENDDO
4285    DO k = 1, klev
4286       DO i = 1, klon
4287          zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
4288               (paprs(i,k)-paprs(i,k+1))/rg
4289          zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
4290               (paprs(i,k)-paprs(i,k+1))/rg
4291       ENDDO
4292    ENDDO
4293    !
4294    !IM calcul composantes axiales du moment angulaire et couple des montagnes
4295    !
4296    IF (is_sequential .and. ok_orodr) THEN
4297       CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4298            ra,rg,romega, &
4299            latitude_deg,longitude_deg,pphis, &
4300            zustrdr,zustrli,zustrph, &
4301            zvstrdr,zvstrli,zvstrph, &
4302            paprs,u,v, &
4303            aam, torsfc)
4304    ENDIF
4305    !IM cf. FLott END
4306    !DC Calcul de la tendance due au methane
4307    IF(ok_qch4) THEN
4308       CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4309       ! ajout de la tendance d'humidite due au methane
4310       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
4311       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4312            'q_ch4', abortphy,flag_inhib_tend,itap,0)
4313       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
4314    ENDIF
4315    !
4316    !
4317
4318!===============================================================
4319!            Additional tendency of TKE due to orography
4320!===============================================================
4321!
4322! Inititialization
4323!------------------
4324
4325   
4326
4327       addtkeoro=0   
4328       CALL getin_p('addtkeoro',addtkeoro)
4329     
4330       IF (prt_level.ge.5) &
4331            print*,'addtkeoro', addtkeoro
4332           
4333       alphatkeoro=1.   
4334       CALL getin_p('alphatkeoro',alphatkeoro)
4335       alphatkeoro=min(max(0.,alphatkeoro),1.)
4336
4337       smallscales_tkeoro=.false.   
4338       CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4339
4340
4341        dtadd(:,:)=0.
4342        duadd(:,:)=0.
4343        dvadd(:,:)=0.
4344
4345
4346
4347! Choices for addtkeoro:
4348!      ** 0 no TKE tendency from orography   
4349!      ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4350!      ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4351!
4352
4353       IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4354!      -------------------------------------------
4355
4356
4357       !  selection des points pour lesquels le schema est actif:
4358
4359
4360
4361  IF (addtkeoro .EQ. 1 ) THEN
4362
4363            duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4364            dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4365
4366  ELSE IF (addtkeoro .EQ. 2) THEN
4367
4368
4369
4370       IF (smallscales_tkeoro) THEN
4371       igwd=0
4372       DO i=1,klon
4373          itest(i)=0
4374! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4375! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4376! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4377          IF (zstd(i).GT.1.0) THEN
4378             itest(i)=1
4379             igwd=igwd+1
4380             idx(igwd)=i
4381          ENDIF
4382       ENDDO
4383
4384     ELSE
4385
4386       igwd=0
4387       DO i=1,klon
4388          itest(i)=0
4389        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4390             itest(i)=1
4391             igwd=igwd+1
4392             idx(igwd)=i
4393          ENDIF
4394       ENDDO
4395
4396       END IF
4397
4398
4399
4400
4401       CALL drag_noro_strato(addtkeoro,klon,klev,dtime,paprs,pplay, &
4402               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4403               igwd,idx,itest, &
4404               t_seri, u_seri, v_seri, &
4405               zulow, zvlow, zustrdr, zvstrdr, &
4406               d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4407
4408            zustrdr(:)=0.
4409            zvstrdr(:)=0.
4410            zulow(:)=0.
4411            zvlow(:)=0.
4412
4413            duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4414            dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4415 END IF
4416   
4417
4418
4419   ! TKE update from subgrid temperature and wind tendencies
4420   !----------------------------------------------------------
4421    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4422
4423
4424    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4425
4426
4427
4428       ENDIF
4429!      -----
4430!===============================================================
4431
4432
4433
4434    !====================================================================
4435    ! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4436    !====================================================================
4437    ! Abderrahmane 24.08.09
4438
4439    IF (ok_cosp) THEN
4440       ! adeclarer
4441#ifdef CPP_COSP
4442       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4443
4444          IF (prt_level .GE.10) THEN
4445             print*,'freq_cosp',freq_cosp
4446          ENDIF
4447          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4448          !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4449          !     s        ref_liq,ref_ice
4450          CALL phys_cosp(itap,dtime,freq_cosp, &
4451               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4452               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4453               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4454               JrNt,ref_liq,ref_ice, &
4455               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4456               zu10m,zv10m,pphis, &
4457               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4458               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4459               prfl(:,1:klev),psfl(:,1:klev), &
4460               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4461               mr_ozone,cldtau, cldemi)
4462
4463          !     L         calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4464          !     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
4465          !     M          clMISR,
4466          !     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4467          !     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
4468
4469       ENDIF
4470
4471#endif
4472    ENDIF  !ok_cosp
4473
4474
4475! Marine
4476
4477  IF (ok_airs) then
4478
4479  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
4480     write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4481     CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4482        & map_prop_hc,map_prop_hist,&
4483        & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4484        & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4485        & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4486        & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4487        & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4488        & map_ntot,map_hc,map_hist,&
4489        & map_Cb,map_ThCi,map_Anv,&
4490        & alt_tropo )
4491  ENDIF
4492
4493  ENDIF  ! ok_airs
4494
4495
4496    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4497    !AA
4498    !AA Installation de l'interface online-offline pour traceurs
4499    !AA
4500    !====================================================================
4501    !   Calcul  des tendances traceurs
4502    !====================================================================
4503    !
4504
4505    IF (type_trac=='repr') THEN
4506       sh_in(:,:) = q_seri(:,:)
4507    ELSE
4508       sh_in(:,:) = qx(:,:,ivap)
4509       ch_in(:,:) = qx(:,:,iliq)
4510    ENDIF
4511
4512    IF (iflag_phytrac == 1 ) THEN
4513
4514#ifdef CPP_Dust
4515      CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
4516                      pdtphys,ftsol,                                   &  ! I
4517                      t,q_seri,paprs,pplay,RHcl,                  &  ! I
4518                      pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,          &  ! I
4519                      coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1,                 &  ! I
4520                      u_seri, v_seri, latitude_deg, longitude_deg,  &
4521                      pphis,pctsrf,pmflxr,pmflxs,prfl,psfl,            &  ! I
4522                      da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij,     &  ! I
4523                      epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con,      &  ! I
4524                      ev,wdtrainA,  wdtrainM,wght_cvfd,              &  ! I
4525                      fm_therm, entr_therm, rneb,                      &  ! I
4526                      beta_prec_fisrt,beta_prec, & !I
4527                      zu10m,zv10m,wstar,ale_bl,ale_wake,               &  ! I
4528                      d_tr_dyn,tr_seri)
4529
4530#else
4531
4532    CALL phytrac ( &
4533         itap,     days_elapsed+1,    jH_cur,   debut, &
4534         lafin,    dtime,     u, v,     t, &
4535         paprs,    pplay,     pmfu,     pmfd, &
4536         pen_u,    pde_u,     pen_d,    pde_d, &
4537         cdragh,   coefh(1:klon,1:klev,is_ave),   fm_therm, entr_therm, &
4538         u1,       v1,        ftsol,    pctsrf, &
4539         zustar,   zu10m,     zv10m, &
4540         wstar(:,is_ave),    ale_bl,         ale_wake, &
4541         latitude_deg, longitude_deg, &
4542         frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4543         presnivs, pphis,     pphi,     albsol1, &
4544         sh_in,   ch_in,    rhcl,      cldfra,   rneb, &
4545         diafra,   cldliq,    itop_con, ibas_con, &
4546         pmflxr,   pmflxs,    prfl,     psfl, &
4547         da,       phi,       mp,       upwd, &
4548         phi2,     d1a,       dam,      sij, wght_cvfd, &        !<<RomP+RL
4549         wdtrainA, wdtrainM,  sigd,     clw,elij, &   !<<RomP
4550         ev,       ep,        epmlmMm,  eplaMm, &     !<<RomP
4551         dnwd,     aerosol_couple,      flxmass_w, &
4552         tau_aero, piz_aero,  cg_aero,  ccm, &
4553         rfname, &
4554         d_tr_dyn, &                                 !<<RomP
4555         tr_seri, init_source)
4556#endif
4557    ENDIF    ! (iflag_phytrac=1)
4558
4559    IF (offline) THEN
4560
4561       IF (prt_level.ge.9) &
4562            print*,'Attention on met a 0 les thermiques pour phystoke'
4563       CALL phystokenc ( &
4564            nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4565            t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4566            fm_therm,entr_therm, &
4567            cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4568            frac_impa, frac_nucl, &
4569            pphis,cell_area,dtime,itap, &
4570            qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4571
4572
4573    ENDIF
4574
4575    !
4576    ! Calculer le transport de l'eau et de l'energie (diagnostique)
4577    !
4578    CALL transp (paprs,zxtsol, &
4579         t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4580         ve, vq, ue, uq, vwat, uwat)
4581    !
4582    !IM global posePB BEG
4583    IF(1.EQ.0) THEN
4584       !
4585       CALL transp_lay (paprs,zxtsol, &
4586            t_seri, q_seri, u_seri, v_seri, zphi, &
4587            ve_lay, vq_lay, ue_lay, uq_lay)
4588       !
4589    ENDIF !(1.EQ.0) THEN
4590    !IM global posePB END
4591    ! Accumuler les variables a stocker dans les fichiers histoire:
4592    !
4593
4594    !================================================================
4595    ! Conversion of kinetic and potential energy into heat, for
4596    ! parameterisation of subgrid-scale motions
4597    !================================================================
4598
4599    d_t_ec(:,:)=0.
4600    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4601    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4602         u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4603         zmasse,exner,d_t_ec)
4604    t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4605
4606    !=======================================================================
4607    !   SORTIES
4608    !=======================================================================
4609    !
4610    !IM initialisation + calculs divers diag AMIP2
4611    !
4612    include "calcul_divers.h"
4613    !
4614    !IM Interpolation sur les niveaux de pression du NMC
4615    !   -------------------------------------------------
4616    !
4617!    include "calcul_STDlev.h"
4618    !
4619    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4620    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4621    !
4622    !cc prw  = eau precipitable
4623    !   prlw = colonne eau liquide
4624    !   prlw = colonne eau solide
4625    prw(:) = 0.
4626    prlw(:) = 0.
4627    prsw(:) = 0.
4628    DO k = 1, klev
4629       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
4630       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4631       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4632    ENDDO
4633    !
4634    IF (type_trac == 'inca') THEN
4635#ifdef INCA
4636       CALL VTe(VTphysiq)
4637       CALL VTb(VTinca)
4638
4639       CALL chemhook_end ( &
4640            dtime, &
4641            pplay, &
4642            t_seri, &
4643            tr_seri, &
4644            nbtr, &
4645            paprs, &
4646            q_seri, &
4647            cell_area, &
4648            pphi, &
4649            pphis, &
4650            zx_rh, &
4651            aps, bps)
4652
4653       CALL VTe(VTinca)
4654       CALL VTb(VTphysiq)
4655#endif
4656    ENDIF
4657
4658
4659    !
4660    ! Convertir les incrementations en tendances
4661    !
4662    IF (prt_level .GE.10) THEN
4663       print *,'Convertir les incrementations en tendances '
4664    ENDIF
4665    !
4666    IF (mydebug) THEN
4667       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4668       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4669       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4670       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4671    ENDIF
4672
4673    DO k = 1, klev
4674       DO i = 1, klon
4675          d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4676          d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4677          d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4678          d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4679          d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4680          !CR: on ajoute le contenu en glace
4681          IF (nqo.eq.3) THEN
4682             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
4683          ENDIF
4684       ENDDO
4685    ENDDO
4686    !
4687    !CR: nb de traceurs eau: nqo
4688    !  IF (nqtot.GE.3) THEN
4689    IF (nqtot.GE.(nqo+1)) THEN
4690       !     DO iq = 3, nqtot
4691       DO iq = nqo+1, nqtot
4692          DO  k = 1, klev
4693             DO  i = 1, klon
4694                ! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4695                d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4696             ENDDO
4697          ENDDO
4698       ENDDO
4699    ENDIF
4700    !
4701    !IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4702    !IM global posePB      include "write_bilKP_ins.h"
4703    !IM global posePB      include "write_bilKP_ave.h"
4704    !
4705
4706    !--OB mass fixer
4707    !--profile is corrected to force mass conservation of water
4708    IF (mass_fixer) THEN
4709    qql2(:)=0.0
4710    DO k = 1, klev
4711      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4712    ENDDO
4713    DO i = 1, klon
4714      !--compute ratio of what q+ql should be with conservation to what it is
4715      corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4716      DO k = 1, klev
4717        q_seri(i,k) =q_seri(i,k)*corrqql
4718        ql_seri(i,k)=ql_seri(i,k)*corrqql
4719      ENDDO
4720    ENDDO
4721    ENDIF
4722    !--fin mass fixer
4723
4724    ! Sauvegarder les valeurs de t et q a la fin de la physique:
4725    !
4726    u_ancien(:,:)  = u_seri(:,:)
4727    v_ancien(:,:)  = v_seri(:,:)
4728    t_ancien(:,:)  = t_seri(:,:)
4729    q_ancien(:,:)  = q_seri(:,:)
4730    ql_ancien(:,:) = ql_seri(:,:)
4731    qs_ancien(:,:) = qs_seri(:,:)
4732    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4733    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4734    CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4735    ! !! RomP >>>
4736    !CR: nb de traceurs eau: nqo
4737    IF (nqtot.GT.nqo) THEN
4738       DO iq = nqo+1, nqtot
4739          tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4740       ENDDO
4741    ENDIF
4742    ! !! RomP <<<
4743    !==========================================================================
4744    ! Sorties des tendances pour un point particulier
4745    ! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4746    ! pour le debug
4747    ! La valeur de igout est attribuee plus haut dans le programme
4748    !==========================================================================
4749
4750    IF (prt_level.ge.1) THEN
4751       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4752       write(lunout,*) &
4753            'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4754       write(lunout,*) &
4755            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4756            pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4757            pctsrf(igout,is_sic)
4758       write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4759       DO k=1,klev
4760          write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4761               d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4762               d_t_eva(igout,k)
4763       ENDDO
4764       write(lunout,*) 'cool,heat'
4765       DO k=1,klev
4766          write(lunout,*) cool(igout,k),heat(igout,k)
4767       ENDDO
4768
4769       !jyg<     (En attendant de statuer sur le sort de d_t_oli)
4770       !jyg!     write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4771       !jyg!     do k=1,klev
4772       !jyg!        write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4773       !jyg!             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4774       !jyg!     enddo
4775       write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4776       DO k=1,klev
4777          write(lunout,*) d_t_vdf(igout,k), &
4778               d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4779       ENDDO
4780       !>jyg
4781
4782       write(lunout,*) 'd_ps ',d_ps(igout)
4783       write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4784       DO k=1,klev
4785          write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4786               d_qx(igout,k,1),d_qx(igout,k,2)
4787       ENDDO
4788    ENDIF
4789
4790    !============================================================
4791    !   Calcul de la temperature potentielle
4792    !============================================================
4793    DO k = 1, klev
4794       DO i = 1, klon
4795          !JYG/IM theta en debut du pas de temps
4796          !JYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4797          !JYG/IM theta en fin de pas de temps de physique
4798          theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4799          ! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4800          !     MPL 20130625
4801          ! fth_fonctions.F90 et parkind1.F90
4802          ! sinon thetal=theta
4803          !       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4804          !    :         ql_seri(i,k))
4805          thetal(i,k)=theta(i,k)
4806       ENDDO
4807    ENDDO
4808    !
4809
4810    ! 22.03.04 BEG
4811    !=============================================================
4812    !   Ecriture des sorties
4813    !=============================================================
4814#ifdef CPP_IOIPSL
4815
4816    ! Recupere des varibles calcule dans differents modules
4817    ! pour ecriture dans histxxx.nc
4818
4819    ! Get some variables from module fonte_neige_mod
4820    CALL fonte_neige_get_vars(pctsrf,  &
4821         zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4822
4823
4824    !=============================================================
4825    ! Separation entre thermiques et non thermiques dans les sorties
4826    ! de fisrtilp
4827    !=============================================================
4828
4829    IF (iflag_thermals>=1) THEN
4830       d_t_lscth=0.
4831       d_t_lscst=0.
4832       d_q_lscth=0.
4833       d_q_lscst=0.
4834       DO k=1,klev
4835          DO i=1,klon
4836             IF (ptconvth(i,k)) THEN
4837                d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4838                d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4839             ELSE
4840                d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4841                d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4842             ENDIF
4843          ENDDO
4844       ENDDO
4845
4846       DO i=1,klon
4847          plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4848          plul_th(i)=prfl(i,1)+psfl(i,1)
4849       ENDDO
4850    ENDIF
4851
4852    !On effectue les sorties:
4853
4854#ifdef CPP_Dust
4855  CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,  &
4856       pplay, lmax_th, aerosol_couple,                 &
4857       ok_ade, ok_aie, ivap, new_aod, ok_sync,         &
4858       ptconv, read_climoz, clevSTD,                   &
4859       ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
4860       flag_aerosol, flag_aerosol_strat, ok_cdnc)
4861#else
4862    CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
4863         pplay, lmax_th, aerosol_couple,                 &
4864         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod,      &
4865         ok_sync, ptconv, read_climoz, clevSTD,          &
4866         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
4867         flag_aerosol, flag_aerosol_strat, ok_cdnc)
4868#endif
4869
4870#ifndef CPP_XIOS
4871    CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4872#endif
4873
4874#endif
4875
4876! On remet des variables a .false. apres un premier appel
4877!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4878    if (debut) then
4879#ifdef CPP_XIOS
4880      swaero_diag=.FALSE.
4881      swaerofree_diag=.FALSE.
4882      dryaod_diag=.FALSE.
4883      ok_4xCO2atm= .FALSE.
4884!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4885
4886      IF (is_master) then
4887        !--setting up swaero_diag to TRUE in XIOS case
4888        IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
4889           xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
4890           xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  &
4891             (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
4892                                 xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  &
4893           !!!--for now these fields are not in the XML files so they are omitted
4894           !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
4895           swaero_diag=.TRUE.
4896
4897        !--setting up swaerofree_diag to TRUE in XIOS case
4898        IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
4899           xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
4900           xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
4901           xios_field_is_active("LWupTOAcleanclr")) &
4902           swaerofree_diag=.TRUE.
4903
4904        !--setting up dryaod_diag to TRUE in XIOS case
4905        DO naero = 1, naero_tot-1
4906         IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
4907        ENDDO
4908        !
4909        !--setting up ok_4xCO2atm to TRUE in XIOS case
4910        IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
4911           xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
4912           xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
4913           xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
4914           xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
4915           xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
4916           ok_4xCO2atm=.TRUE.
4917      endif
4918      !$OMP BARRIER
4919      call bcast(swaero_diag)
4920      call bcast(swaerofree_diag)
4921      call bcast(dryaod_diag)
4922      call bcast(ok_4xCO2atm)
4923!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4924#endif
4925    endif
4926!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4927
4928    !====================================================================
4929    ! Arret du modele apres hgardfou en cas de detection d'un
4930    ! plantage par hgardfou
4931    !====================================================================
4932
4933    IF (abortphy==1) THEN
4934       abort_message ='Plantage hgardfou'
4935       CALL abort_physic (modname,abort_message,1)
4936    ENDIF
4937
4938    ! 22.03.04 END
4939    !
4940    !====================================================================
4941    ! Si c'est la fin, il faut conserver l'etat de redemarrage
4942    !====================================================================
4943    !
4944
4945    IF (lafin) THEN
4946       itau_phy = itau_phy + itap
4947       CALL phyredem ("restartphy.nc")
4948       !         open(97,form="unformatted",file="finbin")
4949       !         write(97) u_seri,v_seri,t_seri,q_seri
4950       !         close(97)
4951       !$OMP MASTER
4952       IF (read_climoz >= 1) THEN
4953          IF (is_mpi_root) THEN
4954             CALL nf95_close(ncid_climoz)
4955          ENDIF
4956          DEALLOCATE(press_edg_climoz) ! pointer
4957          DEALLOCATE(press_cen_climoz) ! pointer
4958       ENDIF
4959       !$OMP END MASTER
4960       print *,' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
4961    ENDIF
4962
4963    !      first=.false.
4964
4965
4966  END SUBROUTINE physiq
4967
4968END MODULE physiq_mod
Note: See TracBrowser for help on using the repository browser.