Context Navigation

physiq_mod.F90 @ 3425

Last change on this file since 3425 was 3425, checked in by Laurent Fairhead, 6 years ago
Followup to r3408 modifications for VolMIP NL
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: 172.0 KB

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

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

Download in other formats:

Original Format