source: CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/LMDZ/physiq_mod.F90 @ 5479

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