Context Navigation

physiq_mod.F90 @ 2731

Last change on this file since 2731 was 2731, checked in by jyg, 8 years ago
Test to guarantee that dtime (in seconds) is an integer.
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: 152.3 KB

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

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

Download in other formats:

Original Format