Context Navigation

physiq_mod.F90 @ 3511

Last change on this file since 3511 was 3511, checked in by idelkadi, 5 years ago
Implementation de la version cospv2 dans LMDZ (suite)
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: 180.7 KB

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

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

Download in other formats:

Original Format