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

Last change on this file since 6554 was 6554, checked in by acosce, 11 months ago

fixe a bug when we are using aerosol_couple with dynamico (in this case we don't need to read aerosol input files)
commit in LMDZ with rev 4627

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