source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/physiq_mod.F90 @ 5306

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