Context Navigation

physiq_mod.F90 @ 3601

Last change on this file since 3601 was 3418, checked in by acozic, 6 years ago
Add modification to take into account value read in INCA restart file
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: 173.3 KB

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

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

Download in other formats:

Original Format