Context Navigation

physiq_mod.F90 @ 3199

Last change on this file since 3199 was 3199, checked in by fhourdin, 6 years ago
Pour la convergence independamment du decoupage de domaine. Pour la convection, le niveau maximum k_upper_cv etait fixe sur un profil d'altitude au milieu du domaine klon/2+1/klon dependant donc du decoupage. Remplace par un test en presnivs -7*log(presnivs(k)/presnivs(1)) > 25. independant du decoupage. Garantit la convergence numerique pour l'ancienne physique. Pas forcement pour la nouvelle si on a des valeurs non nulles des variables convectives vers le sommet.
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.6 KB

Line
1	!
2	! $Id: physiq_mod.F90 3199 2018-02-12 00:40:09Z fhourdin $
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	k_upper_cv = klev
2638	!izero = klon/2+1/klon
2639	!DO k = klev,1,-1
2640	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2641	!ENDDO
2642	! FH : nouveau calcul base sur un profil global sans quoi
2643	! le modele etait sensible au decoupage de domaines
2644	DO k = klev,1,-1
2645	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
2646	ENDDO
2647	IF (prt_level .ge. 5) THEN
2648	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2649	k_upper_cv
2650	ENDIF
2651	!
2652	!>jyg
2653	IF (type_trac == 'repr') THEN
2654	nbtr_tmp=ntra
2655	ELSE
2656	nbtr_tmp=nbtr
2657	ENDIF
2658	!jyg iflag_con est dans clesphys
2659	!c CALL concvl (iflag_con,iflag_clos,
2660	CALL concvl (iflag_clos, &
2661	dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
2662	t_w,q_w,wake_s, &
2663	u_seri,v_seri,tr_seri,nbtr_tmp, &
2664	ALE,ALP, &
2665	sig1,w01, &
2666	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2667	rain_con, snow_con, ibas_con, itop_con, sigd, &
2668	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2669	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2670	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2671	! RomP >>>
2672	!! . pmflxr,pmflxs,da,phi,mp,
2673	!! . ftd,fqd,lalim_conv,wght_th)
2674	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2675	ftd,fqd,lalim_conv,wght_th, &
2676	ev, ep,epmlmMm,eplaMm, &
2677	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2678	tau_cld_cv,coefw_cld_cv,epmax_diag)
2679
2680	! RomP <<<
2681
2682	!IM begin
2683	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2684	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2685	!IM end
2686	!IM cf. FH
2687	clwcon0=qcondc
2688	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2689	!
2690	!jyg<
2691	! If convective tendencies are too large, then call convection
2692	! every time step
2693	cvpas = cvpas_0
2694	DO k=1,k_upper_cv
2695	DO i=1,klon
2696	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
2697	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
2698	dtcon_multistep_max = 3.
2699	dqcon_multistep_max = 0.02
2700	ENDIF
2701	ENDDO
2702	ENDDO
2703	!
2704	DO k=1,k_upper_cv
2705	DO i=1,klon
2706	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
2707	!! abs(d_q_con(i,k)) > 2.e-2) THEN
2708	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
2709	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
2710	cvpas = 1
2711	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
2712	!! i,k,d_t_con(i,k),d_q_con(i,k)
2713	ENDIF
2714	ENDDO
2715	ENDDO
2716	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
2717	!!! call bcast(cvpas)
2718	!!! ------------------------------------------------------------
2719	!>jyg
2720	!
2721	DO i = 1, klon
2722	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
2723	ENDDO
2724	!
2725	!jyg<
2726	! Add the tendency due to the dry adjustment of the wake profile
2727	IF (iflag_wake>=1) THEN
2728	IF (iflag_adjwk == 2) THEN
2729	DO k=1,klev
2730	DO i=1,klon
2731	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
2732	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
2733	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2734	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2735	ENDDO
2736	ENDDO
2737	ENDIF ! (iflag_adjwk = 2)
2738	ENDIF ! (iflag_wake>=1)
2739	!>jyg
2740	!
2741	ELSE ! ok_cvl
2742
2743	! MAF conema3 ne contient pas les traceurs
2744	CALL conema3 (dtime, &
2745	paprs,pplay,t_seri,q_seri, &
2746	u_seri,v_seri,tr_seri,ntra, &
2747	sig1,w01, &
2748	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2749	rain_con, snow_con, ibas_con, itop_con, &
2750	upwd,dnwd,dnwd0,bas,top, &
2751	Ma,cape,tvp,rflag, &
2752	pbase &
2753	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2754	,clwcon0)
2755
2756	ENDIF ! ok_cvl
2757
2758	!
2759	! Correction precip
2760	rain_con = rain_con * cvl_corr
2761	snow_con = snow_con * cvl_corr
2762	!
2763
2764	IF (.NOT. ok_gust) THEN
2765	do i = 1, klon
2766	wd(i)=0.0
2767	enddo
2768	ENDIF
2769
2770	! =================================================================== c
2771	! Calcul des proprietes des nuages convectifs
2772	!
2773
2774	! calcul des proprietes des nuages convectifs
2775	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2776	IF (iflag_cld_cv == 0) THEN
2777	CALL clouds_gno &
2778	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2779	ELSE
2780	CALL clouds_bigauss &
2781	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2782	ENDIF
2783
2784
2785	! =================================================================== c
2786
2787	DO i = 1, klon
2788	itop_con(i) = min(max(itop_con(i),1),klev)
2789	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2790	ENDDO
2791
2792	DO i = 1, klon
2793	ema_pcb(i) = paprs(i,ibas_con(i))
2794	ENDDO
2795	DO i = 1, klon
2796	! L'idicage de itop_con peut cacher un pb potentiel
2797	! FH sous la dictee de JYG, CR
2798	ema_pct(i) = paprs(i,itop_con(i)+1)
2799
2800	IF (itop_con(i).gt.klev-3) THEN
2801	IF (prt_level >= 9) THEN
2802	write(lunout,*)'La convection monte trop haut '
2803	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2804	ENDIF
2805	ENDIF
2806	ENDDO
2807	ELSE IF (iflag_con.eq.0) THEN
2808	write(lunout,*) 'On n appelle pas la convection'
2809	clwcon0=0.
2810	rnebcon0=0.
2811	d_t_con=0.
2812	d_q_con=0.
2813	d_u_con=0.
2814	d_v_con=0.
2815	rain_con=0.
2816	snow_con=0.
2817	bas=1
2818	top=1
2819	ELSE
2820	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2821	CALL abort_physic("physiq", "", 1)
2822	ENDIF
2823
2824	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2825	! . d_u_con, d_v_con)
2826
2827	!jyg Reinitialize proba_notrig and itapcv when convection has been called
2828	proba_notrig(:) = 1.
2829	itapcv = 0
2830	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
2831	!
2832	itapcv = itapcv+1
2833	!
2834	! Compter les steps ou cvpas=1
2835	IF (cvpas == 1) THEN
2836	Ncvpaseq1 = Ncvpaseq1+1
2837	ENDIF
2838	IF (mod(itap,1000) == 0) THEN
2839	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
2840	ENDIF
2841
2842	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
2843	!!! l'energie dans les courants satures.
2844	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
2845	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
2846	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
2847	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
2848	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
2849	!! itap, 1)
2850	!! call prt_enerbil('convection_sat',itap)
2851	!!
2852	!!
2853	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2854	'convection',abortphy,flag_inhib_tend,itap,0)
2855	call prt_enerbil('convection',itap)
2856
2857	!-------------------------------------------------------------------------
2858
2859	IF (mydebug) THEN
2860	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2861	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2862	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2863	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2864	ENDIF
2865
2866	IF (check) THEN
2867	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2868	WRITE(lunout,*)"aprescon=", za
2869	zx_t = 0.0
2870	za = 0.0
2871	DO i = 1, klon
2872	za = za + cell_area(i)/REAL(klon)
2873	zx_t = zx_t + (rain_con(i)+ &
2874	snow_con(i))*cell_area(i)/REAL(klon)
2875	ENDDO
2876	zx_t = zx_t/za*dtime
2877	WRITE(lunout,*)"Precip=", zx_t
2878	ENDIF
2879	IF (zx_ajustq) THEN
2880	DO i = 1, klon
2881	z_apres(i) = 0.0
2882	ENDDO
2883	DO k = 1, klev
2884	DO i = 1, klon
2885	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2886	*(paprs(i,k)-paprs(i,k+1))/RG
2887	ENDDO
2888	ENDDO
2889	DO i = 1, klon
2890	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2891	/z_apres(i)
2892	ENDDO
2893	DO k = 1, klev
2894	DO i = 1, klon
2895	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2896	z_factor(i).LT.(1.0-1.0E-08)) THEN
2897	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2898	ENDIF
2899	ENDDO
2900	ENDDO
2901	ENDIF
2902	zx_ajustq=.FALSE.
2903
2904	!
2905	!==========================================================================
2906	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2907	!pour la couche limite diffuse pour l instant
2908	!
2909	!
2910	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
2911	! il faut rajouter cette tendance calcul\'ee hors des poches
2912	! froides
2913	!
2914	IF (iflag_wake>=1) THEN
2915	!
2916	!
2917	! Call wakes every "wkpas" step
2918	!
2919	IF (MOD(itapwk,wkpas).EQ.0) THEN
2920	!
2921	DO k=1,klev
2922	DO i=1,klon
2923	dt_dwn(i,k) = ftd(i,k)
2924	dq_dwn(i,k) = fqd(i,k)
2925	M_dwn(i,k) = dnwd0(i,k)
2926	M_up(i,k) = upwd(i,k)
2927	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2928	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2929	ENDDO
2930	ENDDO
2931
2932	IF (iflag_wake==2) THEN
2933	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2934	DO k = 1,klev
2935	dt_dwn(:,k)= dt_dwn(:,k)+ &
2936	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2937	dq_dwn(:,k)= dq_dwn(:,k)+ &
2938	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2939	ENDDO
2940	ELSEIF (iflag_wake==3) THEN
2941	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2942	DO k = 1,klev
2943	DO i=1,klon
2944	IF (rneb(i,k)==0.) THEN
2945	! On ne tient compte des tendances qu'en dehors des
2946	! nuages (c'est-\`a-dire a priri dans une region ou
2947	! l'eau se reevapore).
2948	dt_dwn(i,k)= dt_dwn(i,k)+ &
2949	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2950	dq_dwn(i,k)= dq_dwn(i,k)+ &
2951	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2952	ENDIF
2953	ENDDO
2954	ENDDO
2955	ENDIF
2956
2957	!
2958	!calcul caracteristiques de la poche froide
2959	CALL calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
2960	t_seri, q_seri, omega, &
2961	dt_dwn, dq_dwn, M_dwn, M_up, &
2962	dt_a, dq_a, &
2963	sigd, &
2964	wake_deltat, wake_deltaq, wake_s, wake_dens, &
2965	wake_dth, wake_h, &
2966	!! wake_pe, wake_fip, wake_gfl, &
2967	wake_pe, wake_fip_0, wake_gfl, & !! jyg
2968	d_t_wake, d_q_wake, &
2969	wake_k, t_x, q_x, &
2970	wake_omgbdth, wake_dp_omgb, &
2971	wake_dtKE, wake_dqKE, &
2972	wake_omg, wake_dp_deltomg, &
2973	wake_spread, wake_Cstar, d_deltat_wk_gw, &
2974	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
2975	!
2976	!jyg Reinitialize itapwk when wakes have been called
2977	itapwk = 0
2978	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
2979	!
2980	itapwk = itapwk+1
2981	!
2982	!-----------------------------------------------------------------------
2983	! ajout des tendances des poches froides
2984	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
2985	abortphy,flag_inhib_tend,itap,0)
2986	call prt_enerbil('wake',itap)
2987	!------------------------------------------------------------------------
2988
2989	! Increment Wake state variables
2990	IF (iflag_wake_tend .GT. 0.) THEN
2991
2992	CALL add_wake_tend &
2993	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
2994	'wake', abortphy)
2995	call prt_enerbil('wake',itap)
2996	ENDIF ! (iflag_wake_tend .GT. 0.)
2997	!
2998	IF (prt_level .GE. 10) THEN
2999	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3000	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3001	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3002	ENDIF
3003
3004	IF (iflag_alp_wk_cond .GT. 0.) THEN
3005
3006	CALL alpale_wk(dtime, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3007	wake_fip)
3008	ELSE
3009	wake_fip(:) = wake_fip_0(:)
3010	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3011
3012	ENDIF ! (iflag_wake>=1)
3013	!
3014	!===================================================================
3015	! Convection seche (thermiques ou ajustement)
3016	!===================================================================
3017	!
3018	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3019	,seuil_inversion,weak_inversion,dthmin)
3020
3021
3022
3023	d_t_ajsb(:,:)=0.
3024	d_q_ajsb(:,:)=0.
3025	d_t_ajs(:,:)=0.
3026	d_u_ajs(:,:)=0.
3027	d_v_ajs(:,:)=0.
3028	d_q_ajs(:,:)=0.
3029	clwcon0th(:,:)=0.
3030	!
3031	! fm_therm(:,:)=0.
3032	! entr_therm(:,:)=0.
3033	! detr_therm(:,:)=0.
3034	!
3035	IF (prt_level>9) WRITE(lunout,*) &
3036	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3037	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3038	IF (iflag_thermals<0) THEN
3039	! Rien
3040	! ====
3041	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3042
3043
3044	ELSE
3045
3046	! Thermiques
3047	! ==========
3048	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3049	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3050
3051
3052	!cc nrlmd le 10/04/2012
3053	DO k=1,klev+1
3054	DO i=1,klon
3055	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3056	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3057	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3058	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3059	ENDDO
3060	ENDDO
3061	!cc fin nrlmd le 10/04/2012
3062
3063	IF (iflag_thermals>=1) THEN
3064	!jyg<
3065	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3066	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3067	! Appel des thermiques avec les profils exterieurs aux poches
3068	DO k=1,klev
3069	DO i=1,klon
3070	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3071	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3072	u_therm(i,k) = u_seri(i,k)
3073	v_therm(i,k) = v_seri(i,k)
3074	ENDDO
3075	ENDDO
3076	ELSE
3077	! Appel des thermiques avec les profils moyens
3078	DO k=1,klev
3079	DO i=1,klon
3080	t_therm(i,k) = t_seri(i,k)
3081	q_therm(i,k) = q_seri(i,k)
3082	u_therm(i,k) = u_seri(i,k)
3083	v_therm(i,k) = v_seri(i,k)
3084	ENDDO
3085	ENDDO
3086	ENDIF
3087	!>jyg
3088	CALL calltherm(pdtphys &
3089	,pplay,paprs,pphi,weak_inversion &
3090	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3091	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
3092	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3093	,fm_therm,entr_therm,detr_therm &
3094	,zqasc,clwcon0th,lmax_th,ratqscth &
3095	,ratqsdiff,zqsatth &
3096	!on rajoute ale et alp, et les
3097	!caracteristiques de la couche alim
3098	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3099	,ztv,zpspsk,ztla,zthl &
3100	!cc nrlmd le 10/04/2012
3101	,pbl_tke_input,pctsrf,omega,cell_area &
3102	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3103	,n2,s2,ale_bl_stat &
3104	,therm_tke_max,env_tke_max &
3105	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3106	,alp_bl_conv,alp_bl_stat &
3107	!cc fin nrlmd le 10/04/2012
3108	,zqla,ztva )
3109	!
3110	!jyg<
3111	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3112	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3113	! Si les thermiques ne sont presents que hors des
3114	! poches, la tendance moyenne associ\'ee doit etre
3115	! multipliee par la fraction surfacique qu'ils couvrent.
3116	DO k=1,klev
3117	DO i=1,klon
3118	!
3119	d_deltat_the(i,k) = - d_t_ajs(i,k)
3120	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3121	!
3122	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3123	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3124	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3125	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3126	!
3127	ENDDO
3128	ENDDO
3129	!
3130	IF (ok_bug_split_th) THEN
3131	CALL add_wake_tend &
3132	(d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
3133	ELSE
3134	CALL add_wake_tend &
3135	(d_deltat_the, d_deltaq_the, dsig0, ddens0, wake_k, 'the', abortphy)
3136	ENDIF
3137	call prt_enerbil('the',itap)
3138	!
3139	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3140	!
3141	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3142	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3143	call prt_enerbil('thermals',itap)
3144	!
3145	!
3146	CALL alpale_th( dtime, lmax_th, t_seri, cell_area, &
3147	cin, s2, n2, &
3148	ale_bl_trig, ale_bl_stat, ale_bl, &
3149	alp_bl, alp_bl_stat, &
3150	proba_notrig, random_notrig)
3151	!>jyg
3152
3153	! ------------------------------------------------------------------
3154	! Transport de la TKE par les panaches thermiques.
3155	! FH : 2010/02/01
3156	! if (iflag_pbl.eq.10) then
3157	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3158	! s rg,paprs,pbl_tke)
3159	! endif
3160	! -------------------------------------------------------------------
3161
3162	DO i=1,klon
3163	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3164	!CR:04/05/12:correction calcul zmax
3165	zmax_th(i)=zmax0(i)
3166	ENDDO
3167
3168	ENDIF
3169
3170	! Ajustement sec
3171	! ==============
3172
3173	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3174	! a partir du sommet des thermiques.
3175	! Dans le cas contraire, on demarre au niveau 1.
3176
3177	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3178
3179	IF (iflag_thermals.eq.0) THEN
3180	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3181	limbas(:)=1
3182	ELSE
3183	limbas(:)=lmax_th(:)
3184	ENDIF
3185
3186	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3187	! pour des test de convergence numerique.
3188	! Le nouveau ajsec est a priori mieux, meme pour le cas
3189	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3190	! non nulles numeriquement pour des mailles non concernees.
3191
3192	IF (iflag_thermals==0) THEN
3193	! Calling adjustment alone (but not the thermal plume model)
3194	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3195	, d_t_ajsb, d_q_ajsb)
3196	ELSE IF (iflag_thermals>0) THEN
3197	! Calling adjustment above the top of thermal plumes
3198	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3199	, d_t_ajsb, d_q_ajsb)
3200	ENDIF
3201
3202	!--------------------------------------------------------------------
3203	! ajout des tendances de l'ajustement sec ou des thermiques
3204	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3205	'ajsb',abortphy,flag_inhib_tend,itap,0)
3206	call prt_enerbil('ajsb',itap)
3207	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3208	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3209
3210	!---------------------------------------------------------------------
3211
3212	ENDIF
3213
3214	ENDIF
3215	!
3216	!===================================================================
3217	! Computation of ratqs, the width (normalized) of the subrid scale
3218	! water distribution
3219	CALL calcratqs(klon,klev,prt_level,lunout, &
3220	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3221	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3222	tau_ratqs,fact_cldcon, &
3223	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3224	paprs,pplay,q_seri,zqsat,fm_therm, &
3225	ratqs,ratqsc)
3226
3227
3228	!
3229	! Appeler le processus de condensation a grande echelle
3230	! et le processus de precipitation
3231	!-------------------------------------------------------------------------
3232	IF (prt_level .GE.10) THEN
3233	print *,'itap, ->fisrtilp ',itap
3234	ENDIF
3235	!
3236	CALL fisrtilp(dtime,paprs,pplay, &
3237	t_seri, q_seri,ptconv,ratqs, &
3238	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3239	rain_lsc, snow_lsc, &
3240	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3241	frac_impa, frac_nucl, beta_prec_fisrt, &
3242	prfl, psfl, rhcl, &
3243	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3244	iflag_ice_thermo)
3245	!
3246	WHERE (rain_lsc < 0) rain_lsc = 0.
3247	WHERE (snow_lsc < 0) snow_lsc = 0.
3248
3249	!+JLD
3250	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3251	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3252	! & ,rain_lsc,snow_lsc
3253	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3254	!-JLD
3255	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3256	'lsc',abortphy,flag_inhib_tend,itap,0)
3257	call prt_enerbil('lsc',itap)
3258	rain_num(:)=0.
3259	DO k = 1, klev
3260	DO i = 1, klon
3261	IF (ql_seri(i,k)>oliqmax) THEN
3262	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3263	ql_seri(i,k)=oliqmax
3264	ENDIF
3265	ENDDO
3266	ENDDO
3267	IF (nqo==3) THEN
3268	DO k = 1, klev
3269	DO i = 1, klon
3270	IF (qs_seri(i,k)>oicemax) THEN
3271	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3272	qs_seri(i,k)=oicemax
3273	ENDIF
3274	ENDDO
3275	ENDDO
3276	ENDIF
3277
3278	!---------------------------------------------------------------------------
3279	DO k = 1, klev
3280	DO i = 1, klon
3281	cldfra(i,k) = rneb(i,k)
3282	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3283	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3284	ENDDO
3285	ENDDO
3286	IF (check) THEN
3287	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3288	WRITE(lunout,*)"apresilp=", za
3289	zx_t = 0.0
3290	za = 0.0
3291	DO i = 1, klon
3292	za = za + cell_area(i)/REAL(klon)
3293	zx_t = zx_t + (rain_lsc(i) &
3294	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3295	ENDDO
3296	zx_t = zx_t/za*dtime
3297	WRITE(lunout,*)"Precip=", zx_t
3298	ENDIF
3299
3300	IF (mydebug) THEN
3301	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3302	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3303	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3304	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3305	ENDIF
3306
3307	!
3308	!-------------------------------------------------------------------
3309	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3310	!-------------------------------------------------------------------
3311
3312	! 1. NUAGES CONVECTIFS
3313	!
3314	!IM cf FH
3315	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3316	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3317	snow_tiedtke=0.
3318	! print*,'avant calcul de la pseudo precip '
3319	! print*,'iflag_cld_th',iflag_cld_th
3320	IF (iflag_cld_th.eq.-1) THEN
3321	rain_tiedtke=rain_con
3322	ELSE
3323	! print*,'calcul de la pseudo precip '
3324	rain_tiedtke=0.
3325	! print*,'calcul de la pseudo precip 0'
3326	DO k=1,klev
3327	DO i=1,klon
3328	IF (d_q_con(i,k).lt.0.) THEN
3329	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3330	*(paprs(i,k)-paprs(i,k+1))/rg
3331	ENDIF
3332	ENDDO
3333	ENDDO
3334	ENDIF
3335	!
3336	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3337	!
3338
3339	! Nuages diagnostiques pour Tiedtke
3340	CALL diagcld1(paprs,pplay, &
3341	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3342	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3343	diafra,dialiq)
3344	DO k = 1, klev
3345	DO i = 1, klon
3346	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3347	cldliq(i,k) = dialiq(i,k)
3348	cldfra(i,k) = diafra(i,k)
3349	ENDIF
3350	ENDDO
3351	ENDDO
3352
3353	ELSE IF (iflag_cld_th.ge.3) THEN
3354	! On prend pour les nuages convectifs le max du calcul de la
3355	! convection et du calcul du pas de temps precedent diminue d'un facteur
3356	! facttemps
3357	facteur = pdtphys *facttemps
3358	DO k=1,klev
3359	DO i=1,klon
3360	rnebcon(i,k)=rnebcon(i,k)*facteur
3361	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3362	rnebcon(i,k)=rnebcon0(i,k)
3363	clwcon(i,k)=clwcon0(i,k)
3364	ENDIF
3365	ENDDO
3366	ENDDO
3367
3368	! On prend la somme des fractions nuageuses et des contenus en eau
3369
3370	IF (iflag_cld_th>=5) THEN
3371
3372	DO k=1,klev
3373	ptconvth(:,k)=fm_therm(:,k+1)>0.
3374	ENDDO
3375
3376	IF (iflag_coupl==4) THEN
3377
3378	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3379	! convectives et lsc dans la partie des thermiques
3380	! Le controle par iflag_coupl est peut etre provisoire.
3381	DO k=1,klev
3382	DO i=1,klon
3383	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3384	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3385	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3386	ELSE IF (ptconv(i,k)) THEN
3387	cldfra(i,k)=rnebcon(i,k)
3388	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3389	ENDIF
3390	ENDDO
3391	ENDDO
3392
3393	ELSE IF (iflag_coupl==5) THEN
3394	DO k=1,klev
3395	DO i=1,klon
3396	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3397	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3398	ENDDO
3399	ENDDO
3400
3401	ELSE
3402
3403	! Si on est sur un point touche par la convection
3404	! profonde et pas par les thermiques, on prend la
3405	! couverture nuageuse et l'eau nuageuse de la convection
3406	! profonde.
3407
3408	!IM/FH: 2011/02/23
3409	! definition des points sur lesquels ls thermiques sont actifs
3410
3411	DO k=1,klev
3412	DO i=1,klon
3413	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3414	cldfra(i,k)=rnebcon(i,k)
3415	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3416	ENDIF
3417	ENDDO
3418	ENDDO
3419
3420	ENDIF
3421
3422	ELSE
3423
3424	! Ancienne version
3425	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3426	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3427	ENDIF
3428
3429	ENDIF
3430
3431	! plulsc(:)=0.
3432	! do k=1,klev,-1
3433	! do i=1,klon
3434	! zzz=prfl(:,k)+psfl(:,k)
3435	! if (.not.ptconvth.zzz.gt.0.)
3436	! enddo prfl, psfl,
3437	! enddo
3438	!
3439	! 2. NUAGES STARTIFORMES
3440	!
3441	IF (ok_stratus) THEN
3442	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3443	DO k = 1, klev
3444	DO i = 1, klon
3445	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3446	cldliq(i,k) = dialiq(i,k)
3447	cldfra(i,k) = diafra(i,k)
3448	ENDIF
3449	ENDDO
3450	ENDDO
3451	ENDIF
3452	!
3453	! Precipitation totale
3454	!
3455	DO i = 1, klon
3456	rain_fall(i) = rain_con(i) + rain_lsc(i)
3457	snow_fall(i) = snow_con(i) + snow_lsc(i)
3458	ENDDO
3459	!
3460	! Calculer l'humidite relative pour diagnostique
3461	!
3462	DO k = 1, klev
3463	DO i = 1, klon
3464	zx_t = t_seri(i,k)
3465	IF (thermcep) THEN
3466	!! if (iflag_ice_thermo.eq.0) then !jyg
3467	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3468	!! else !jyg
3469	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3470	!! endif !jyg
3471	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3472	zx_qs = MIN(0.5,zx_qs)
3473	zcor = 1./(1.-retv*zx_qs)
3474	zx_qs = zx_qs*zcor
3475	ELSE
3476	!! IF (zx_t.LT.t_coup) THEN !jyg
3477	IF (zx_t.LT.rtt) THEN !jyg
3478	zx_qs = qsats(zx_t)/pplay(i,k)
3479	ELSE
3480	zx_qs = qsatl(zx_t)/pplay(i,k)
3481	ENDIF
3482	ENDIF
3483	zx_rh(i,k) = q_seri(i,k)/zx_qs
3484	zqsat(i,k)=zx_qs
3485	ENDDO
3486	ENDDO
3487
3488	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3489	! equivalente a 2m (tpote) pour diagnostique
3490	!
3491	DO i = 1, klon
3492	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3493	IF (thermcep) THEN
3494	IF(zt2m(i).LT.RTT) then
3495	Lheat=RLSTT
3496	ELSE
3497	Lheat=RLVTT
3498	ENDIF
3499	ELSE
3500	IF (zt2m(i).LT.RTT) THEN
3501	Lheat=RLSTT
3502	ELSE
3503	Lheat=RLVTT
3504	ENDIF
3505	ENDIF
3506	tpote(i) = tpot(i)* &
3507	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3508	ENDDO
3509
3510	IF (type_trac == 'inca') THEN
3511	#ifdef INCA
3512	CALL VTe(VTphysiq)
3513	CALL VTb(VTinca)
3514	calday = REAL(days_elapsed + 1) + jH_cur
3515
3516	CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
3517	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3518	CALL AEROSOL_METEO_CALC( &
3519	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3520	prfl,psfl,pctsrf,cell_area, &
3521	latitude_deg,longitude_deg,u10m,v10m)
3522	ENDIF
3523
3524	zxsnow_dummy(:) = 0.0
3525
3526	CALL chemhook_begin (calday, &
3527	days_elapsed+1, &
3528	jH_cur, &
3529	pctsrf(1,1), &
3530	latitude_deg, &
3531	longitude_deg, &
3532	cell_area, &
3533	paprs, &
3534	pplay, &
3535	coefh(1:klon,1:klev,is_ave), &
3536	pphi, &
3537	t_seri, &
3538	u, &
3539	v, &
3540	wo(:, :, 1), &
3541	q_seri, &
3542	zxtsol, &
3543	zxsnow_dummy, &
3544	solsw, &
3545	albsol1, &
3546	rain_fall, &
3547	snow_fall, &
3548	itop_con, &
3549	ibas_con, &
3550	cldfra, &
3551	nbp_lon, &
3552	nbp_lat-1, &
3553	tr_seri, &
3554	ftsol, &
3555	paprs, &
3556	cdragh, &
3557	cdragm, &
3558	pctsrf, &
3559	pdtphys, &
3560	itap)
3561
3562	CALL VTe(VTinca)
3563	CALL VTb(VTphysiq)
3564	#endif
3565	ENDIF !type_trac = inca
3566
3567
3568	!
3569	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3570	!
3571	IF (MOD(itaprad,radpas).EQ.0) THEN
3572
3573	!
3574	!jq - introduce the aerosol direct and first indirect radiative forcings
3575	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3576	IF (flag_aerosol .GT. 0) THEN
3577	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3578	IF (.NOT. aerosol_couple) THEN
3579	!
3580	CALL readaerosol_optic( &
3581	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3582	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3583	mass_solu_aero, mass_solu_aero_pi, &
3584	tau_aero, piz_aero, cg_aero, &
3585	tausum_aero, tau3d_aero)
3586	ENDIF
3587	ELSE ! RRTM radiation
3588	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3589	abort_message='config_inca=aero et rrtm=1 impossible'
3590	CALL abort_physic(modname,abort_message,1)
3591	ELSE
3592	!
3593	#ifdef CPP_RRTM
3594	IF (NSW.EQ.6) THEN
3595	!--new aerosol properties SW and LW
3596	!
3597	#ifdef CPP_Dust
3598	!--SPL aerosol model
3599	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3600	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3601	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3602	tausum_aero, tau3d_aero)
3603	#else
3604	!--climatologies or INCA aerosols
3605	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
3606	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3607	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3608	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3609	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3610	tausum_aero, drytausum_aero, tau3d_aero)
3611	#endif
3612	!
3613	ELSE IF (NSW.EQ.2) THEN
3614	!--for now we use the old aerosol properties
3615	!
3616	CALL readaerosol_optic( &
3617	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3618	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3619	mass_solu_aero, mass_solu_aero_pi, &
3620	tau_aero, piz_aero, cg_aero, &
3621	tausum_aero, tau3d_aero)
3622	!
3623	!--natural aerosols
3624	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3625	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3626	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3627	!--all aerosols
3628	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3629	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3630	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3631	!
3632	!--no LW optics
3633	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3634	!
3635	ELSE
3636	abort_message='Only NSW=2 or 6 are possible with ' &
3637	// 'aerosols and iflag_rrtm=1'
3638	CALL abort_physic(modname,abort_message,1)
3639	ENDIF
3640	#else
3641	abort_message='You should compile with -rrtm if running ' &
3642	// 'with iflag_rrtm=1'
3643	CALL abort_physic(modname,abort_message,1)
3644	#endif
3645	!
3646	ENDIF
3647	ENDIF
3648	ELSE !--flag_aerosol = 0
3649	tausum_aero(:,:,:) = 0.
3650	drytausum_aero(:,:) = 0.
3651	mass_solu_aero(:,:) = 0.
3652	mass_solu_aero_pi(:,:) = 0.
3653	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3654	tau_aero(:,:,:,:) = 1.e-15
3655	piz_aero(:,:,:,:) = 1.
3656	cg_aero(:,:,:,:) = 0.
3657	ELSE
3658	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3659	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3660	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3661	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3662	ENDIF
3663	ENDIF
3664	!
3665	!--WMO criterion to determine tropopause
3666	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3667	!
3668	!--STRAT AEROSOL
3669	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3670	IF (flag_aerosol_strat.GT.0) THEN
3671	IF (prt_level .GE.10) THEN
3672	PRINT *,'appel a readaerosolstrat', mth_cur
3673	ENDIF
3674	IF (iflag_rrtm.EQ.0) THEN
3675	IF (flag_aerosol_strat.EQ.1) THEN
3676	CALL readaerosolstrato(debut)
3677	ELSE
3678	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3679	CALL abort_physic(modname,abort_message,1)
3680	ENDIF
3681	ELSE
3682	#ifdef CPP_RRTM
3683	#ifndef CPP_StratAer
3684	!--prescribed strat aerosols
3685	!--only in the case of non-interactive strat aerosols
3686	IF (flag_aerosol_strat.EQ.1) THEN
3687	CALL readaerosolstrato1_rrtm(debut)
3688	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3689	CALL readaerosolstrato2_rrtm(debut)
3690	ELSE
3691	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3692	CALL abort_physic(modname,abort_message,1)
3693	ENDIF
3694	#endif
3695	#else
3696	abort_message='You should compile with -rrtm if running ' &
3697	// 'with iflag_rrtm=1'
3698	CALL abort_physic(modname,abort_message,1)
3699	#endif
3700	ENDIF
3701	ENDIF
3702	!
3703	#ifdef CPP_RRTM
3704	#ifdef CPP_StratAer
3705	!--compute stratospheric mask
3706	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3707	!--interactive strat aerosols
3708	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3709	#endif
3710	#endif
3711	!--fin STRAT AEROSOL
3712	!
3713
3714	! Calculer les parametres optiques des nuages et quelques
3715	! parametres pour diagnostiques:
3716	!
3717	IF (aerosol_couple.AND.config_inca=='aero') THEN
3718	mass_solu_aero(:,:) = ccm(:,:,1)
3719	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3720	ENDIF
3721
3722	IF (ok_newmicro) then
3723	IF (iflag_rrtm.NE.0) THEN
3724	#ifdef CPP_RRTM
3725	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3726	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3727	// 'pour ok_cdnc'
3728	CALL abort_physic(modname,abort_message,1)
3729	ENDIF
3730	#else
3731
3732	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3733	CALL abort_physic(modname,abort_message,1)
3734	#endif
3735	ENDIF
3736	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3737	paprs, pplay, t_seri, cldliq, cldfra, &
3738	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3739	flwp, fiwp, flwc, fiwc, &
3740	mass_solu_aero, mass_solu_aero_pi, &
3741	cldtaupi, re, fl, ref_liq, ref_ice, &
3742	ref_liq_pi, ref_ice_pi)
3743	ELSE
3744	CALL nuage (paprs, pplay, &
3745	t_seri, cldliq, cldfra, cldtau, cldemi, &
3746	cldh, cldl, cldm, cldt, cldq, &
3747	ok_aie, &
3748	mass_solu_aero, mass_solu_aero_pi, &
3749	bl95_b0, bl95_b1, &
3750	cldtaupi, re, fl)
3751	ENDIF
3752	!
3753	!IM betaCRF
3754	!
3755	cldtaurad = cldtau
3756	cldtaupirad = cldtaupi
3757	cldemirad = cldemi
3758	cldfrarad = cldfra
3759
3760	!
3761	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3762	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3763	!
3764	! global
3765	!
3766	!IM 251017 begin
3767	print*,'physiq betaCRF global zdtime=',zdtime
3768	!IM 251017 end
3769	DO k=1, klev
3770	DO i=1, klon
3771	IF (pplay(i,k).GE.pfree) THEN
3772	beta(i,k) = beta_pbl
3773	ELSE
3774	beta(i,k) = beta_free
3775	ENDIF
3776	IF (mskocean_beta) THEN
3777	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3778	ENDIF
3779	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3780	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3781	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3782	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3783	ENDDO
3784	ENDDO
3785	!
3786	ELSE
3787	!
3788	! regional
3789	!
3790	DO k=1, klev
3791	DO i=1,klon
3792	!
3793	IF (longitude_deg(i).ge.lon1_beta.AND. &
3794	longitude_deg(i).le.lon2_beta.AND. &
3795	latitude_deg(i).le.lat1_beta.AND. &
3796	latitude_deg(i).ge.lat2_beta) THEN
3797	IF (pplay(i,k).GE.pfree) THEN
3798	beta(i,k) = beta_pbl
3799	ELSE
3800	beta(i,k) = beta_free
3801	ENDIF
3802	IF (mskocean_beta) THEN
3803	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3804	ENDIF
3805	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3806	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3807	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3808	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3809	ENDIF
3810	!
3811	ENDDO
3812	ENDDO
3813	!
3814	ENDIF
3815
3816	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3817	IF (ok_chlorophyll) THEN
3818	print*,"-- reading chlorophyll"
3819	CALL readchlorophyll(debut)
3820	ENDIF
3821
3822	!--if ok_suntime_rrtm we use ancillay data for RSUN
3823	!--previous values are therefore overwritten
3824	!--this is needed for CMIP6 runs
3825	!--and only possible for new radiation scheme
3826	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
3827	#ifdef CPP_RRTM
3828	CALL read_rsun_rrtm(debut)
3829	#endif
3830	ENDIF
3831
3832	IF (mydebug) THEN
3833	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3834	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3835	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3836	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3837	ENDIF
3838
3839	!
3840	!sonia : If Iflag_radia >=2, pertubation of some variables
3841	!input to radiation (DICE)
3842	!
3843	IF (iflag_radia .ge. 2) THEN
3844	zsav_tsol (:) = zxtsol(:)
3845	CALL perturb_radlwsw(zxtsol,iflag_radia)
3846	ENDIF
3847
3848	IF (aerosol_couple.AND.config_inca=='aero') THEN
3849	#ifdef INCA
3850	CALL radlwsw_inca &
3851	(kdlon,kflev,dist, rmu0, fract, solaire, &
3852	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3853	size(wo,3), wo, &
3854	cldfrarad, cldemirad, cldtaurad, &
3855	heat,heat0,cool,cool0,albpla, &
3856	topsw,toplw,solsw,sollw, &
3857	sollwdown, &
3858	topsw0,toplw0,solsw0,sollw0, &
3859	lwdn0, lwdn, lwup0, lwup, &
3860	swdn0, swdn, swup0, swup, &
3861	ok_ade, ok_aie, &
3862	tau_aero, piz_aero, cg_aero, &
3863	topswad_aero, solswad_aero, &
3864	topswad0_aero, solswad0_aero, &
3865	topsw_aero, topsw0_aero, &
3866	solsw_aero, solsw0_aero, &
3867	cldtaupirad, &
3868	topswai_aero, solswai_aero)
3869	#endif
3870	ELSE
3871	!
3872	!IM calcul radiatif pour le cas actuel
3873	!
3874	RCO2 = RCO2_act
3875	RCH4 = RCH4_act
3876	RN2O = RN2O_act
3877	RCFC11 = RCFC11_act
3878	RCFC12 = RCFC12_act
3879	!
3880	IF (prt_level .GE.10) THEN
3881	print *,' ->radlwsw, number 1 '
3882	ENDIF
3883
3884	!
3885	CALL radlwsw &
3886	(dist, rmu0, fract, &
3887	!albedo SB >>>
3888	! paprs, pplay,zxtsol,albsol1, albsol2, &
3889	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3890	!albedo SB <<<
3891	t_seri,q_seri,wo, &
3892	cldfrarad, cldemirad, cldtaurad, &
3893	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3894	flag_aerosol_strat, &
3895	tau_aero, piz_aero, cg_aero, &
3896	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3897	! Rajoute par OB pour RRTM
3898	tau_aero_lw_rrtm, &
3899	cldtaupirad,new_aod, &
3900	! zqsat, flwcrad, fiwcrad, &
3901	zqsat, flwc, fiwc, &
3902	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3903	heat,heat0,cool,cool0,albpla, &
3904	topsw,toplw,solsw,sollw, &
3905	sollwdown, &
3906	topsw0,toplw0,solsw0,sollw0, &
3907	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
3908	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
3909	topswad_aero, solswad_aero, &
3910	topswai_aero, solswai_aero, &
3911	topswad0_aero, solswad0_aero, &
3912	topsw_aero, topsw0_aero, &
3913	solsw_aero, solsw0_aero, &
3914	topswcf_aero, solswcf_aero, &
3915	!-C. Kleinschmitt for LW diagnostics
3916	toplwad_aero, sollwad_aero,&
3917	toplwai_aero, sollwai_aero, &
3918	toplwad0_aero, sollwad0_aero,&
3919	!-end
3920	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3921	ZSWFT0_i, ZFSDN0, ZFSUP0)
3922
3923	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
3924	!schemes
3925	toplw = toplw + betalwoff * (toplw0 - toplw)
3926	sollw = sollw + betalwoff * (sollw0 - sollw)
3927	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
3928	lwup = lwup + betalwoff * (lwup0 - lwup)
3929	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
3930	sollwdown(:))
3931	cool = cool + betalwoff * (cool0 - cool)
3932
3933	#ifndef CPP_XIOS
3934	!--OB 30/05/2016 modified 21/10/2016
3935	!--here we return swaero_diag and dryaod_diag to FALSE
3936	!--and histdef will switch it back to TRUE if necessary
3937	!--this is necessary to get the right swaero at first step
3938	!--but only in the case of no XIOS as XIOS is covered elsewhere
3939	IF (debut) swaerofree_diag = .FALSE.
3940	IF (debut) swaero_diag = .FALSE.
3941	IF (debut) dryaod_diag = .FALSE.
3942	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
3943	!--as for swaero_diag, see above
3944	IF (debut) ok_4xCO2atm = .FALSE.
3945
3946	!
3947	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3948	!IM des taux doit etre different du taux actuel
3949	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3950	!
3951	IF (RCO2_per.NE.RCO2_act.OR. &
3952	RCH4_per.NE.RCH4_act.OR. &
3953	RN2O_per.NE.RN2O_act.OR. &
3954	RCFC11_per.NE.RCFC11_act.OR. &
3955	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
3956	#endif
3957	!
3958	IF (ok_4xCO2atm) THEN
3959	!
3960	RCO2 = RCO2_per
3961	RCH4 = RCH4_per
3962	RN2O = RN2O_per
3963	RCFC11 = RCFC11_per
3964	RCFC12 = RCFC12_per
3965	!
3966	IF (prt_level .GE.10) THEN
3967	print *,' ->radlwsw, number 2 '
3968	ENDIF
3969	!
3970	CALL radlwsw &
3971	(dist, rmu0, fract, &
3972	!albedo SB >>>
3973	! paprs, pplay,zxtsol,albsol1, albsol2, &
3974	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3975	!albedo SB <<<
3976	t_seri,q_seri,wo, &
3977	cldfrarad, cldemirad, cldtaurad, &
3978	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3979	flag_aerosol_strat, &
3980	tau_aero, piz_aero, cg_aero, &
3981	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3982	! Rajoute par OB pour RRTM
3983	tau_aero_lw_rrtm, &
3984	cldtaupi,new_aod, &
3985	! zqsat, flwcrad, fiwcrad, &
3986	zqsat, flwc, fiwc, &
3987	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3988	heatp,heat0p,coolp,cool0p,albplap, &
3989	topswp,toplwp,solswp,sollwp, &
3990	sollwdownp, &
3991	topsw0p,toplw0p,solsw0p,sollw0p, &
3992	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
3993	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
3994	topswad_aerop, solswad_aerop, &
3995	topswai_aerop, solswai_aerop, &
3996	topswad0_aerop, solswad0_aerop, &
3997	topsw_aerop, topsw0_aerop, &
3998	solsw_aerop, solsw0_aerop, &
3999	topswcf_aerop, solswcf_aerop, &
4000	!-C. Kleinschmitt for LW diagnostics
4001	toplwad_aerop, sollwad_aerop,&
4002	toplwai_aerop, sollwai_aerop, &
4003	toplwad0_aerop, sollwad0_aerop,&
4004	!-end
4005	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4006	ZSWFT0_i, ZFSDN0, ZFSUP0)
4007	endif !ok_4xCO2atm
4008	ENDIF ! aerosol_couple
4009	itaprad = 0
4010	!
4011	! If Iflag_radia >=2, reset pertubed variables
4012	!
4013	IF (iflag_radia .ge. 2) THEN
4014	zxtsol(:) = zsav_tsol (:)
4015	ENDIF
4016	ENDIF ! MOD(itaprad,radpas)
4017	itaprad = itaprad + 1
4018
4019	IF (iflag_radia.eq.0) THEN
4020	IF (prt_level.ge.9) THEN
4021	PRINT *,'--------------------------------------------------'
4022	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4023	PRINT *,'>>>> heat et cool mis a zero '
4024	PRINT *,'--------------------------------------------------'
4025	ENDIF
4026	heat=0.
4027	cool=0.
4028	sollw=0. ! MPL 01032011
4029	solsw=0.
4030	radsol=0.
4031	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4032	swup0=0.
4033	lwup=0.
4034	lwup0=0.
4035	lwdn=0.
4036	lwdn0=0.
4037	ENDIF
4038
4039	!
4040	! Calculer radsol a l'exterieur de radlwsw
4041	! pour prendre en compte le cycle diurne
4042	! recode par Olivier Boucher en sept 2015
4043	!
4044	radsol=solsw*swradcorr+sollw
4045
4046	IF (ok_4xCO2atm) THEN
4047	radsolp=solswp*swradcorr+sollwp
4048	ENDIF
4049
4050	!
4051	! Ajouter la tendance des rayonnements (tous les pas)
4052	! avec une correction pour le cycle diurne dans le SW
4053	!
4054
4055	DO k=1, klev
4056	d_t_swr(:,k)=swradcorr(:)heat(:,k)dtime/RDAY
4057	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)dtime/RDAY
4058	d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
4059	d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
4060	ENDDO
4061
4062	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4063	call prt_enerbil('SW',itap)
4064	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4065	call prt_enerbil('LW',itap)
4066
4067	!
4068	IF (mydebug) THEN
4069	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4070	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4071	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4072	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4073	ENDIF
4074
4075	! Calculer l'hydrologie de la surface
4076	!
4077	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4078	! . agesno, ftsol,fqsurf,fsnow, ruis)
4079	!
4080
4081	!
4082	! Calculer le bilan du sol et la derive de temperature (couplage)
4083	!
4084	DO i = 1, klon
4085	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4086	! a la demande de JLD
4087	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4088	ENDDO
4089	!
4090	!moddeblott(jan95)
4091	! Appeler le programme de parametrisation de l'orographie
4092	! a l'echelle sous-maille:
4093	!
4094	IF (prt_level .GE.10) THEN
4095	print *,' call orography ? ', ok_orodr
4096	ENDIF
4097	!
4098	IF (ok_orodr) THEN
4099	!
4100	! selection des points pour lesquels le shema est actif:
4101	igwd=0
4102	DO i=1,klon
4103	itest(i)=0
4104	! IF ((zstd(i).gt.10.0)) THEN
4105	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4106	itest(i)=1
4107	igwd=igwd+1
4108	idx(igwd)=i
4109	ENDIF
4110	ENDDO
4111	! igwdim=MAX(1,igwd)
4112	!
4113	IF (ok_strato) THEN
4114
4115	CALL drag_noro_strato(0,klon,klev,dtime,paprs,pplay, &
4116	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4117	igwd,idx,itest, &
4118	t_seri, u_seri, v_seri, &
4119	zulow, zvlow, zustrdr, zvstrdr, &
4120	d_t_oro, d_u_oro, d_v_oro)
4121
4122	ELSE
4123	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
4124	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4125	igwd,idx,itest, &
4126	t_seri, u_seri, v_seri, &
4127	zulow, zvlow, zustrdr, zvstrdr, &
4128	d_t_oro, d_u_oro, d_v_oro)
4129	ENDIF
4130	!
4131	! ajout des tendances
4132	!-----------------------------------------------------------------------
4133	! ajout des tendances de la trainee de l'orographie
4134	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4135	abortphy,flag_inhib_tend,itap,0)
4136	call prt_enerbil('oro',itap)
4137	!----------------------------------------------------------------------
4138	!
4139	ENDIF ! fin de test sur ok_orodr
4140	!
4141	IF (mydebug) THEN
4142	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4143	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4144	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4145	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4146	ENDIF
4147
4148	IF (ok_orolf) THEN
4149	!
4150	! selection des points pour lesquels le shema est actif:
4151	igwd=0
4152	DO i=1,klon
4153	itest(i)=0
4154	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4155	itest(i)=1
4156	igwd=igwd+1
4157	idx(igwd)=i
4158	ENDIF
4159	ENDDO
4160	! igwdim=MAX(1,igwd)
4161	!
4162	IF (ok_strato) THEN
4163
4164	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
4165	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4166	igwd,idx,itest, &
4167	t_seri, u_seri, v_seri, &
4168	zulow, zvlow, zustrli, zvstrli, &
4169	d_t_lif, d_u_lif, d_v_lif )
4170
4171	ELSE
4172	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
4173	latitude_deg,zmea,zstd,zpic, &
4174	itest, &
4175	t_seri, u_seri, v_seri, &
4176	zulow, zvlow, zustrli, zvstrli, &
4177	d_t_lif, d_u_lif, d_v_lif)
4178	ENDIF
4179
4180	! ajout des tendances de la portance de l'orographie
4181	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4182	'lif', abortphy,flag_inhib_tend,itap,0)
4183	call prt_enerbil('lif',itap)
4184	ENDIF ! fin de test sur ok_orolf
4185
4186	IF (ok_hines) then
4187	! HINES GWD PARAMETRIZATION
4188	east_gwstress=0.
4189	west_gwstress=0.
4190	du_gwd_hines=0.
4191	dv_gwd_hines=0.
4192	CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
4193	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4194	du_gwd_hines, dv_gwd_hines)
4195	zustr_gwd_hines=0.
4196	zvstr_gwd_hines=0.
4197	DO k = 1, klev
4198	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
4199	* (paprs(:, k)-paprs(:, k+1))/rg
4200	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
4201	* (paprs(:, k)-paprs(:, k+1))/rg
4202	ENDDO
4203
4204	d_t_hin(:, :)=0.
4205	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4206	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4207	call prt_enerbil('hin',itap)
4208	ENDIF
4209
4210	IF (.not. ok_hines .and. ok_gwd_rando) then
4211	CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
4212	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4213	dv_gwd_front, east_gwstress, west_gwstress)
4214	zustr_gwd_front=0.
4215	zvstr_gwd_front=0.
4216	DO k = 1, klev
4217	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
4218	* (paprs(:, k)-paprs(:, k+1))/rg
4219	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
4220	* (paprs(:, k)-paprs(:, k+1))/rg
4221	ENDDO
4222
4223	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4224	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4225	call prt_enerbil('front_gwd_rando',itap)
4226	ENDIF
4227
4228	IF (ok_gwd_rando) THEN
4229	CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
4230	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4231	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4232	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4233	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4234	call prt_enerbil('flott_gwd_rando',itap)
4235	zustr_gwd_rando=0.
4236	zvstr_gwd_rando=0.
4237	DO k = 1, klev
4238	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
4239	* (paprs(:, k)-paprs(:, k+1))/rg
4240	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
4241	* (paprs(:, k)-paprs(:, k+1))/rg
4242	ENDDO
4243	ENDIF
4244
4245	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4246
4247	IF (mydebug) THEN
4248	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4249	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4250	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4251	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4252	ENDIF
4253
4254	DO i = 1, klon
4255	zustrph(i)=0.
4256	zvstrph(i)=0.
4257	ENDDO
4258	DO k = 1, klev
4259	DO i = 1, klon
4260	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
4261	(paprs(i,k)-paprs(i,k+1))/rg
4262	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
4263	(paprs(i,k)-paprs(i,k+1))/rg
4264	ENDDO
4265	ENDDO
4266	!
4267	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4268	!
4269	IF (is_sequential .and. ok_orodr) THEN
4270	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4271	ra,rg,romega, &
4272	latitude_deg,longitude_deg,pphis, &
4273	zustrdr,zustrli,zustrph, &
4274	zvstrdr,zvstrli,zvstrph, &
4275	paprs,u,v, &
4276	aam, torsfc)
4277	ENDIF
4278	!IM cf. FLott END
4279	!DC Calcul de la tendance due au methane
4280	IF(ok_qch4) THEN
4281	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4282	! ajout de la tendance d'humidite due au methane
4283	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
4284	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4285	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4286	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
4287	ENDIF
4288	!
4289	!
4290
4291	!===============================================================
4292	! Additional tendency of TKE due to orography
4293	!===============================================================
4294	!
4295	! Inititialization
4296	!------------------
4297
4298
4299
4300	addtkeoro=0
4301	CALL getin_p('addtkeoro',addtkeoro)
4302
4303	IF (prt_level.ge.5) &
4304	print*,'addtkeoro', addtkeoro
4305
4306	alphatkeoro=1.
4307	CALL getin_p('alphatkeoro',alphatkeoro)
4308	alphatkeoro=min(max(0.,alphatkeoro),1.)
4309
4310	smallscales_tkeoro=.false.
4311	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4312
4313
4314	dtadd(:,:)=0.
4315	duadd(:,:)=0.
4316	dvadd(:,:)=0.
4317
4318
4319
4320	! Choices for addtkeoro:
4321	! ** 0 no TKE tendency from orography
4322	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4323	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4324	!
4325
4326	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4327	! -------------------------------------------
4328
4329
4330	! selection des points pour lesquels le schema est actif:
4331
4332
4333
4334	IF (addtkeoro .EQ. 1 ) THEN
4335
4336	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4337	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4338
4339	ELSE IF (addtkeoro .EQ. 2) THEN
4340
4341
4342
4343	IF (smallscales_tkeoro) THEN
4344	igwd=0
4345	DO i=1,klon
4346	itest(i)=0
4347	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4348	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4349	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4350	IF (zstd(i).GT.1.0) THEN
4351	itest(i)=1
4352	igwd=igwd+1
4353	idx(igwd)=i
4354	ENDIF
4355	ENDDO
4356
4357	ELSE
4358
4359	igwd=0
4360	DO i=1,klon
4361	itest(i)=0
4362	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4363	itest(i)=1
4364	igwd=igwd+1
4365	idx(igwd)=i
4366	ENDIF
4367	ENDDO
4368
4369	END IF
4370
4371
4372
4373
4374	CALL drag_noro_strato(addtkeoro,klon,klev,dtime,paprs,pplay, &
4375	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4376	igwd,idx,itest, &
4377	t_seri, u_seri, v_seri, &
4378	zulow, zvlow, zustrdr, zvstrdr, &
4379	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4380
4381	zustrdr(:)=0.
4382	zvstrdr(:)=0.
4383	zulow(:)=0.
4384	zvlow(:)=0.
4385
4386	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4387	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4388	END IF
4389
4390
4391
4392	! TKE update from subgrid temperature and wind tendencies
4393	!----------------------------------------------------------
4394	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4395
4396
4397	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4398
4399
4400
4401	ENDIF
4402	! -----
4403	!===============================================================
4404
4405
4406
4407	!====================================================================
4408	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4409	!====================================================================
4410	! Abderrahmane 24.08.09
4411
4412	IF (ok_cosp) THEN
4413	! adeclarer
4414	#ifdef CPP_COSP
4415	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4416
4417	IF (prt_level .GE.10) THEN
4418	print*,'freq_cosp',freq_cosp
4419	ENDIF
4420	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4421	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4422	! s ref_liq,ref_ice
4423	CALL phys_cosp(itap,dtime,freq_cosp, &
4424	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4425	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4426	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4427	JrNt,ref_liq,ref_ice, &
4428	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4429	zu10m,zv10m,pphis, &
4430	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4431	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4432	prfl(:,1:klev),psfl(:,1:klev), &
4433	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4434	mr_ozone,cldtau, cldemi)
4435
4436	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4437	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4438	! M clMISR,
4439	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4440	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4441
4442	ENDIF
4443
4444	#endif
4445	ENDIF !ok_cosp
4446
4447
4448	! Marine
4449
4450	IF (ok_airs) then
4451
4452	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
4453	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4454	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4455	& map_prop_hc,map_prop_hist,&
4456	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4457	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4458	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4459	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4460	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4461	& map_ntot,map_hc,map_hist,&
4462	& map_Cb,map_ThCi,map_Anv,&
4463	& alt_tropo )
4464	ENDIF
4465
4466	ENDIF ! ok_airs
4467
4468
4469	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4470	!AA
4471	!AA Installation de l'interface online-offline pour traceurs
4472	!AA
4473	!====================================================================
4474	! Calcul des tendances traceurs
4475	!====================================================================
4476	!
4477
4478	IF (type_trac=='repr') THEN
4479	sh_in(:,:) = q_seri(:,:)
4480	ELSE
4481	sh_in(:,:) = qx(:,:,ivap)
4482	ch_in(:,:) = qx(:,:,iliq)
4483	ENDIF
4484
4485	IF (iflag_phytrac == 1 ) THEN
4486
4487	#ifdef CPP_Dust
4488	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4489	pdtphys,ftsol, & ! I
4490	t,q_seri,paprs,pplay,RHcl, & ! I
4491	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4492	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4493	u_seri, v_seri, latitude_deg, longitude_deg, &
4494	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4495	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4496	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4497	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4498	fm_therm, entr_therm, rneb, & ! I
4499	beta_prec_fisrt,beta_prec, & !I
4500	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4501	d_tr_dyn,tr_seri)
4502
4503	#else
4504
4505	CALL phytrac ( &
4506	itap, days_elapsed+1, jH_cur, debut, &
4507	lafin, dtime, u, v, t, &
4508	paprs, pplay, pmfu, pmfd, &
4509	pen_u, pde_u, pen_d, pde_d, &
4510	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4511	u1, v1, ftsol, pctsrf, &
4512	zustar, zu10m, zv10m, &
4513	wstar(:,is_ave), ale_bl, ale_wake, &
4514	latitude_deg, longitude_deg, &
4515	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4516	presnivs, pphis, pphi, albsol1, &
4517	sh_in, ch_in, rhcl, cldfra, rneb, &
4518	diafra, cldliq, itop_con, ibas_con, &
4519	pmflxr, pmflxs, prfl, psfl, &
4520	da, phi, mp, upwd, &
4521	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4522	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4523	ev, ep, epmlmMm, eplaMm, & !<<RomP
4524	dnwd, aerosol_couple, flxmass_w, &
4525	tau_aero, piz_aero, cg_aero, ccm, &
4526	rfname, &
4527	d_tr_dyn, & !<<RomP
4528	tr_seri)
4529	#endif
4530	ENDIF ! (iflag_phytrac=1)
4531
4532	IF (offline) THEN
4533
4534	IF (prt_level.ge.9) &
4535	print*,'Attention on met a 0 les thermiques pour phystoke'
4536	CALL phystokenc ( &
4537	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4538	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4539	fm_therm,entr_therm, &
4540	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4541	frac_impa, frac_nucl, &
4542	pphis,cell_area,dtime,itap, &
4543	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4544
4545
4546	ENDIF
4547
4548	!
4549	! Calculer le transport de l'eau et de l'energie (diagnostique)
4550	!
4551	CALL transp (paprs,zxtsol, &
4552	t_seri, q_seri, u_seri, v_seri, zphi, &
4553	ve, vq, ue, uq)
4554	!
4555	!IM global posePB BEG
4556	IF(1.EQ.0) THEN
4557	!
4558	CALL transp_lay (paprs,zxtsol, &
4559	t_seri, q_seri, u_seri, v_seri, zphi, &
4560	ve_lay, vq_lay, ue_lay, uq_lay)
4561	!
4562	ENDIF !(1.EQ.0) THEN
4563	!IM global posePB END
4564	! Accumuler les variables a stocker dans les fichiers histoire:
4565	!
4566
4567	!================================================================
4568	! Conversion of kinetic and potential energy into heat, for
4569	! parameterisation of subgrid-scale motions
4570	!================================================================
4571
4572	d_t_ec(:,:)=0.
4573	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4574	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4575	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4576	zmasse,exner,d_t_ec)
4577	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4578
4579	!=======================================================================
4580	! SORTIES
4581	!=======================================================================
4582	!
4583	!IM initialisation + calculs divers diag AMIP2
4584	!
4585	include "calcul_divers.h"
4586	!
4587	!IM Interpolation sur les niveaux de pression du NMC
4588	! -------------------------------------------------
4589	!
4590	include "calcul_STDlev.h"
4591	!
4592	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4593	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4594	!
4595	!cc prw = eau precipitable
4596	! prlw = colonne eau liquide
4597	! prlw = colonne eau solide
4598	prw(:) = 0.
4599	prlw(:) = 0.
4600	prsw(:) = 0.
4601	DO k = 1, klev
4602	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4603	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4604	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4605	ENDDO
4606	!
4607	IF (type_trac == 'inca') THEN
4608	#ifdef INCA
4609	CALL VTe(VTphysiq)
4610	CALL VTb(VTinca)
4611
4612	CALL chemhook_end ( &
4613	dtime, &
4614	pplay, &
4615	t_seri, &
4616	tr_seri, &
4617	nbtr, &
4618	paprs, &
4619	q_seri, &
4620	cell_area, &
4621	pphi, &
4622	pphis, &
4623	zx_rh, &
4624	aps, bps)
4625
4626	CALL VTe(VTinca)
4627	CALL VTb(VTphysiq)
4628	#endif
4629	ENDIF
4630
4631
4632	!
4633	! Convertir les incrementations en tendances
4634	!
4635	IF (prt_level .GE.10) THEN
4636	print *,'Convertir les incrementations en tendances '
4637	ENDIF
4638	!
4639	IF (mydebug) THEN
4640	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4641	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4642	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4643	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4644	ENDIF
4645
4646	DO k = 1, klev
4647	DO i = 1, klon
4648	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4649	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4650	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4651	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4652	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4653	!CR: on ajoute le contenu en glace
4654	IF (nqo.eq.3) THEN
4655	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
4656	ENDIF
4657	ENDDO
4658	ENDDO
4659	!
4660	!CR: nb de traceurs eau: nqo
4661	! IF (nqtot.GE.3) THEN
4662	IF (nqtot.GE.(nqo+1)) THEN
4663	! DO iq = 3, nqtot
4664	DO iq = nqo+1, nqtot
4665	DO k = 1, klev
4666	DO i = 1, klon
4667	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4668	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4669	ENDDO
4670	ENDDO
4671	ENDDO
4672	ENDIF
4673	!
4674	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4675	!IM global posePB include "write_bilKP_ins.h"
4676	!IM global posePB include "write_bilKP_ave.h"
4677	!
4678
4679	!--OB mass fixer
4680	!--profile is corrected to force mass conservation of water
4681	IF (mass_fixer) THEN
4682	qql2(:)=0.0
4683	DO k = 1, klev
4684	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4685	ENDDO
4686	DO i = 1, klon
4687	!--compute ratio of what q+ql should be with conservation to what it is
4688	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4689	DO k = 1, klev
4690	q_seri(i,k) =q_seri(i,k)*corrqql
4691	ql_seri(i,k)=ql_seri(i,k)*corrqql
4692	ENDDO
4693	ENDDO
4694	ENDIF
4695	!--fin mass fixer
4696
4697	! Sauvegarder les valeurs de t et q a la fin de la physique:
4698	!
4699	u_ancien(:,:) = u_seri(:,:)
4700	v_ancien(:,:) = v_seri(:,:)
4701	t_ancien(:,:) = t_seri(:,:)
4702	q_ancien(:,:) = q_seri(:,:)
4703	ql_ancien(:,:) = ql_seri(:,:)
4704	qs_ancien(:,:) = qs_seri(:,:)
4705	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4706	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4707	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4708	! !! RomP >>>
4709	!CR: nb de traceurs eau: nqo
4710	IF (nqtot.GT.nqo) THEN
4711	DO iq = nqo+1, nqtot
4712	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4713	ENDDO
4714	ENDIF
4715	! !! RomP <<<
4716	!==========================================================================
4717	! Sorties des tendances pour un point particulier
4718	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4719	! pour le debug
4720	! La valeur de igout est attribuee plus haut dans le programme
4721	!==========================================================================
4722
4723	IF (prt_level.ge.1) THEN
4724	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4725	write(lunout,*) &
4726	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4727	write(lunout,*) &
4728	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4729	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4730	pctsrf(igout,is_sic)
4731	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4732	DO k=1,klev
4733	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4734	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4735	d_t_eva(igout,k)
4736	ENDDO
4737	write(lunout,*) 'cool,heat'
4738	DO k=1,klev
4739	write(lunout,*) cool(igout,k),heat(igout,k)
4740	ENDDO
4741
4742	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4743	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4744	!jyg! do k=1,klev
4745	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4746	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4747	!jyg! enddo
4748	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4749	DO k=1,klev
4750	write(lunout,*) d_t_vdf(igout,k), &
4751	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4752	ENDDO
4753	!>jyg
4754
4755	write(lunout,*) 'd_ps ',d_ps(igout)
4756	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4757	DO k=1,klev
4758	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4759	d_qx(igout,k,1),d_qx(igout,k,2)
4760	ENDDO
4761	ENDIF
4762
4763	!============================================================
4764	! Calcul de la temperature potentielle
4765	!============================================================
4766	DO k = 1, klev
4767	DO i = 1, klon
4768	!JYG/IM theta en debut du pas de temps
4769	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4770	!JYG/IM theta en fin de pas de temps de physique
4771	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4772	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4773	! MPL 20130625
4774	! fth_fonctions.F90 et parkind1.F90
4775	! sinon thetal=theta
4776	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4777	! : ql_seri(i,k))
4778	thetal(i,k)=theta(i,k)
4779	ENDDO
4780	ENDDO
4781	!
4782
4783	! 22.03.04 BEG
4784	!=============================================================
4785	! Ecriture des sorties
4786	!=============================================================
4787	#ifdef CPP_IOIPSL
4788
4789	! Recupere des varibles calcule dans differents modules
4790	! pour ecriture dans histxxx.nc
4791
4792	! Get some variables from module fonte_neige_mod
4793	CALL fonte_neige_get_vars(pctsrf, &
4794	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4795
4796
4797	!=============================================================
4798	! Separation entre thermiques et non thermiques dans les sorties
4799	! de fisrtilp
4800	!=============================================================
4801
4802	IF (iflag_thermals>=1) THEN
4803	d_t_lscth=0.
4804	d_t_lscst=0.
4805	d_q_lscth=0.
4806	d_q_lscst=0.
4807	DO k=1,klev
4808	DO i=1,klon
4809	IF (ptconvth(i,k)) THEN
4810	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4811	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4812	ELSE
4813	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4814	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4815	ENDIF
4816	ENDDO
4817	ENDDO
4818
4819	DO i=1,klon
4820	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4821	plul_th(i)=prfl(i,1)+psfl(i,1)
4822	ENDDO
4823	ENDIF
4824
4825	!On effectue les sorties:
4826
4827	#ifdef CPP_Dust
4828	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
4829	pplay, lmax_th, aerosol_couple, &
4830	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4831	ptconv, read_climoz, clevSTD, &
4832	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
4833	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4834	#else
4835	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4836	pplay, lmax_th, aerosol_couple, &
4837	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
4838	ok_sync, ptconv, read_climoz, clevSTD, &
4839	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
4840	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4841	#endif
4842
4843	#ifndef CPP_XIOS
4844	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4845	#endif
4846
4847	#endif
4848
4849
4850	!====================================================================
4851	! Arret du modele apres hgardfou en cas de detection d'un
4852	! plantage par hgardfou
4853	!====================================================================
4854
4855	IF (abortphy==1) THEN
4856	abort_message ='Plantage hgardfou'
4857	CALL abort_physic (modname,abort_message,1)
4858	ENDIF
4859
4860	! 22.03.04 END
4861	!
4862	!====================================================================
4863	! Si c'est la fin, il faut conserver l'etat de redemarrage
4864	!====================================================================
4865	!
4866
4867	IF (lafin) THEN
4868	itau_phy = itau_phy + itap
4869	CALL phyredem ("restartphy.nc")
4870	! open(97,form="unformatted",file="finbin")
4871	! write(97) u_seri,v_seri,t_seri,q_seri
4872	! close(97)
4873	!$OMP MASTER
4874	IF (read_climoz >= 1) THEN
4875	IF (is_mpi_root) THEN
4876	CALL nf95_close(ncid_climoz)
4877	ENDIF
4878	DEALLOCATE(press_edg_climoz) ! pointer
4879	DEALLOCATE(press_cen_climoz) ! pointer
4880	ENDIF
4881	!$OMP END MASTER
4882	print *,' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
4883	ENDIF
4884
4885	! first=.false.
4886
4887
4888	END SUBROUTINE physiq
4889
4890	END MODULE physiq_mod

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

Download in other formats:

Original Format