Context Navigation

physiq_mod.F90 @ 3208

Last change on this file since 3208 was 3208, checked in by jyg, 6 years ago
Implementation of a first crude model of the dynamic of wake population.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 169.7 KB

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

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

Download in other formats:

Original Format