Context Navigation

physiq_mod.F90 @ 4140

Last change on this file since 4140 was 4140, checked in by lfalletti, 2 years ago

Fixed bug for coupling with reprobus

rm of a call Init_chem_rep_phys that was set twice
add of hybrid coordinates in reprobus outputs (for CMIP6 dr2xml)
bug fix to take into account H2Orep tracer name

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 190.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format