Context Navigation

physiq_mod.F90 @ 3201

Last change on this file since 3201 was 3201, checked in by Laurent Fairhead, 6 years ago

Inclusion of r3199 from trunk

Pour la convergence independamment du decoupage de domaine.
Pour la convection, le niveau maximum k_upper_cv etait fixe
sur un profil d'altitude au milieu du domaine klon/2+1/klon
dependant donc du decoupage.
Remplace par un test en presnivs
-7*log(presnivs(k)/presnivs(1)) > 25.
independant du decoupage.
FH

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: 168.8 KB

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

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

Download in other formats:

Original Format