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

Last change on this file since 3400 was 3400, checked in by Laurent Fairhead, 6 years ago

Due to a difficult to correct bug, near surface fields can take abnormal values in very specific circumstances.
While waiting to find a correction, these modifications will force the values of these fields within reasonable bounds.
To get back to the previous behaviour of the model, use iflag_bug_t2m_ipslcm61 = 1 (it is 0 by default)

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