Context Navigation

physiq_mod.F90 @ 3180

Last change on this file since 3180 was 3180, checked in by jyg, 6 years ago
Correction to commit 3178: the bug is corrected only if ok_bug_split_th is "false". The default is "true".
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 `Author Date Id Revision`
File size: 169.4 KB

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

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

Download in other formats:

Original Format