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

Last change on this file since 5484 was 5484, checked in by aclsce, 2 years ago

Updated LMDZ sources modified to be used in IPSLCM7 configuration.

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