Context Navigation

physiq_mod.F90 @ 3148

Last change on this file since 3148 was 3148, checked in by musat, 6 years ago
Bug fix for ratqsc (add save/OMP THREADPRIVATE) and itau_con counter JYG/IM
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: 166.3 KB

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

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

Download in other formats:

Original Format