Context Navigation

physiq_mod.F90 @ 3719

Last change on this file since 3719 was 3633, checked in by acozic, 5 years ago
merge with rev 3632 on LMDZ6/trunk
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: 172.7 KB

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

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

Download in other formats:

Original Format