Context Navigation

physiq_mod.F90 @ 5446

Last change on this file since 5446 was 3823, checked in by musat, 4 years ago
Nouveaux calculs a 2m et 10m
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: 183.6 KB

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

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

Download in other formats:

Original Format