Context Navigation

physiq_mod.F90 @ 3609

Last change on this file since 3609 was 3609, checked in by jghattas, 5 years ago
Put back changes done in commit [3604] which were removed temporary in commit [3607].
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 3609 2019-11-29 08:25:41Z jghattas $
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
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	wo(:, :, 1), &
3560	q_seri, &
3561	zxtsol, &
3562	zxsnow_dummy, &
3563	solsw, &
3564	albsol1, &
3565	rain_fall, &
3566	snow_fall, &
3567	itop_con, &
3568	ibas_con, &
3569	cldfra, &
3570	nbp_lon, &
3571	nbp_lat-1, &
3572	tr_seri, &
3573	ftsol, &
3574	paprs, &
3575	cdragh, &
3576	cdragm, &
3577	pctsrf, &
3578	pdtphys, &
3579	itap)
3580
3581	CALL VTe(VTinca)
3582	CALL VTb(VTphysiq)
3583	#endif
3584	ENDIF !type_trac = inca
3585
3586
3587	!
3588	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3589	!
3590	IF (MOD(itaprad,radpas).EQ.0) THEN
3591
3592	!
3593	!jq - introduce the aerosol direct and first indirect radiative forcings
3594	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3595	IF (flag_aerosol .GT. 0) THEN
3596	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3597	IF (.NOT. aerosol_couple) THEN
3598	!
3599	CALL readaerosol_optic( &
3600	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3601	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3602	mass_solu_aero, mass_solu_aero_pi, &
3603	tau_aero, piz_aero, cg_aero, &
3604	tausum_aero, tau3d_aero)
3605	ENDIF
3606	ELSE ! RRTM radiation
3607	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3608	abort_message='config_inca=aero et rrtm=1 impossible'
3609	CALL abort_physic(modname,abort_message,1)
3610	ELSE
3611	!
3612	#ifdef CPP_RRTM
3613	IF (NSW.EQ.6) THEN
3614	!--new aerosol properties SW and LW
3615	!
3616	#ifdef CPP_Dust
3617	!--SPL aerosol model
3618	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3619	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3620	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3621	tausum_aero, tau3d_aero)
3622	#else
3623	!--climatologies or INCA aerosols
3624	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
3625	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3626	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3627	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3628	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3629	tausum_aero, drytausum_aero, tau3d_aero)
3630	#endif
3631
3632	IF (flag_aerosol .EQ. 7) THEN
3633	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
3634	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3635	ENDIF
3636
3637	!
3638	ELSE IF (NSW.EQ.2) THEN
3639	!--for now we use the old aerosol properties
3640	!
3641	CALL readaerosol_optic( &
3642	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3643	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3644	mass_solu_aero, mass_solu_aero_pi, &
3645	tau_aero, piz_aero, cg_aero, &
3646	tausum_aero, tau3d_aero)
3647	!
3648	!--natural aerosols
3649	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3650	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3651	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3652	!--all aerosols
3653	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3654	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3655	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3656	!
3657	!--no LW optics
3658	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3659	!
3660	ELSE
3661	abort_message='Only NSW=2 or 6 are possible with ' &
3662	// 'aerosols and iflag_rrtm=1'
3663	CALL abort_physic(modname,abort_message,1)
3664	ENDIF
3665	#else
3666	abort_message='You should compile with -rrtm if running ' &
3667	// 'with iflag_rrtm=1'
3668	CALL abort_physic(modname,abort_message,1)
3669	#endif
3670	!
3671	ENDIF
3672	ENDIF
3673	ELSE !--flag_aerosol = 0
3674	tausum_aero(:,:,:) = 0.
3675	drytausum_aero(:,:) = 0.
3676	mass_solu_aero(:,:) = 0.
3677	mass_solu_aero_pi(:,:) = 0.
3678	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3679	tau_aero(:,:,:,:) = 1.e-15
3680	piz_aero(:,:,:,:) = 1.
3681	cg_aero(:,:,:,:) = 0.
3682	ELSE
3683	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3684	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3685	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3686	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3687	ENDIF
3688	ENDIF
3689	!
3690	!--WMO criterion to determine tropopause
3691	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3692	!
3693	!--STRAT AEROSOL
3694	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3695	IF (flag_aerosol_strat.GT.0) THEN
3696	IF (prt_level .GE.10) THEN
3697	PRINT *,'appel a readaerosolstrat', mth_cur
3698	ENDIF
3699	IF (iflag_rrtm.EQ.0) THEN
3700	IF (flag_aerosol_strat.EQ.1) THEN
3701	CALL readaerosolstrato(debut)
3702	ELSE
3703	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3704	CALL abort_physic(modname,abort_message,1)
3705	ENDIF
3706	ELSE
3707	#ifdef CPP_RRTM
3708	#ifndef CPP_StratAer
3709	!--prescribed strat aerosols
3710	!--only in the case of non-interactive strat aerosols
3711	IF (flag_aerosol_strat.EQ.1) THEN
3712	CALL readaerosolstrato1_rrtm(debut)
3713	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3714	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
3715	ELSE
3716	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3717	CALL abort_physic(modname,abort_message,1)
3718	ENDIF
3719	#endif
3720	#else
3721	abort_message='You should compile with -rrtm if running ' &
3722	// 'with iflag_rrtm=1'
3723	CALL abort_physic(modname,abort_message,1)
3724	#endif
3725	ENDIF
3726	ENDIF
3727	!
3728	#ifdef CPP_RRTM
3729	#ifdef CPP_StratAer
3730	!--compute stratospheric mask
3731	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3732	!--interactive strat aerosols
3733	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3734	#endif
3735	#endif
3736	!--fin STRAT AEROSOL
3737	!
3738
3739	! Calculer les parametres optiques des nuages et quelques
3740	! parametres pour diagnostiques:
3741	!
3742	IF (aerosol_couple.AND.config_inca=='aero') THEN
3743	mass_solu_aero(:,:) = ccm(:,:,1)
3744	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3745	ENDIF
3746
3747	IF (ok_newmicro) then
3748	IF (iflag_rrtm.NE.0) THEN
3749	#ifdef CPP_RRTM
3750	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3751	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3752	// 'pour ok_cdnc'
3753	CALL abort_physic(modname,abort_message,1)
3754	ENDIF
3755	#else
3756
3757	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3758	CALL abort_physic(modname,abort_message,1)
3759	#endif
3760	ENDIF
3761	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
3762	paprs, pplay, t_seri, cldliq, cldfra, &
3763	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3764	flwp, fiwp, flwc, fiwc, &
3765	mass_solu_aero, mass_solu_aero_pi, &
3766	cldtaupi, re, fl, ref_liq, ref_ice, &
3767	ref_liq_pi, ref_ice_pi)
3768	ELSE
3769	CALL nuage (paprs, pplay, &
3770	t_seri, cldliq, cldfra, cldtau, cldemi, &
3771	cldh, cldl, cldm, cldt, cldq, &
3772	ok_aie, &
3773	mass_solu_aero, mass_solu_aero_pi, &
3774	bl95_b0, bl95_b1, &
3775	cldtaupi, re, fl)
3776	ENDIF
3777	!
3778	!IM betaCRF
3779	!
3780	cldtaurad = cldtau
3781	cldtaupirad = cldtaupi
3782	cldemirad = cldemi
3783	cldfrarad = cldfra
3784
3785	!
3786	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3787	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3788	!
3789	! global
3790	!
3791	!IM 251017 begin
3792	! print*,'physiq betaCRF global zdtime=',zdtime
3793	!IM 251017 end
3794	DO k=1, klev
3795	DO i=1, klon
3796	IF (pplay(i,k).GE.pfree) THEN
3797	beta(i,k) = beta_pbl
3798	ELSE
3799	beta(i,k) = beta_free
3800	ENDIF
3801	IF (mskocean_beta) THEN
3802	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3803	ENDIF
3804	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3805	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3806	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3807	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3808	ENDDO
3809	ENDDO
3810	!
3811	ELSE
3812	!
3813	! regional
3814	!
3815	DO k=1, klev
3816	DO i=1,klon
3817	!
3818	IF (longitude_deg(i).ge.lon1_beta.AND. &
3819	longitude_deg(i).le.lon2_beta.AND. &
3820	latitude_deg(i).le.lat1_beta.AND. &
3821	latitude_deg(i).ge.lat2_beta) THEN
3822	IF (pplay(i,k).GE.pfree) THEN
3823	beta(i,k) = beta_pbl
3824	ELSE
3825	beta(i,k) = beta_free
3826	ENDIF
3827	IF (mskocean_beta) THEN
3828	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3829	ENDIF
3830	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3831	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3832	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3833	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3834	ENDIF
3835	!
3836	ENDDO
3837	ENDDO
3838	!
3839	ENDIF
3840
3841	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3842	IF (ok_chlorophyll) THEN
3843	print*,"-- reading chlorophyll"
3844	CALL readchlorophyll(debut)
3845	ENDIF
3846
3847	!--if ok_suntime_rrtm we use ancillay data for RSUN
3848	!--previous values are therefore overwritten
3849	!--this is needed for CMIP6 runs
3850	!--and only possible for new radiation scheme
3851	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
3852	#ifdef CPP_RRTM
3853	CALL read_rsun_rrtm(debut)
3854	#endif
3855	ENDIF
3856
3857	IF (mydebug) THEN
3858	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3859	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3860	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3861	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3862	ENDIF
3863
3864	!
3865	!sonia : If Iflag_radia >=2, pertubation of some variables
3866	!input to radiation (DICE)
3867	!
3868	IF (iflag_radia .ge. 2) THEN
3869	zsav_tsol (:) = zxtsol(:)
3870	CALL perturb_radlwsw(zxtsol,iflag_radia)
3871	ENDIF
3872
3873	IF (aerosol_couple.AND.config_inca=='aero') THEN
3874	#ifdef INCA
3875	CALL radlwsw_inca &
3876	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
3877	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3878	size(wo,3), wo, &
3879	cldfrarad, cldemirad, cldtaurad, &
3880	heat,heat0,cool,cool0,albpla, &
3881	topsw,toplw,solsw,sollw, &
3882	sollwdown, &
3883	topsw0,toplw0,solsw0,sollw0, &
3884	lwdn0, lwdn, lwup0, lwup, &
3885	swdn0, swdn, swup0, swup, &
3886	ok_ade, ok_aie, &
3887	tau_aero, piz_aero, cg_aero, &
3888	topswad_aero, solswad_aero, &
3889	topswad0_aero, solswad0_aero, &
3890	topsw_aero, topsw0_aero, &
3891	solsw_aero, solsw0_aero, &
3892	cldtaupirad, &
3893	topswai_aero, solswai_aero)
3894	#endif
3895	ELSE
3896	!
3897	!IM calcul radiatif pour le cas actuel
3898	!
3899	RCO2 = RCO2_act
3900	RCH4 = RCH4_act
3901	RN2O = RN2O_act
3902	RCFC11 = RCFC11_act
3903	RCFC12 = RCFC12_act
3904	!
3905	IF (prt_level .GE.10) THEN
3906	print *,' ->radlwsw, number 1 '
3907	ENDIF
3908
3909	!
3910	CALL radlwsw &
3911	(dist, rmu0, fract, &
3912	!albedo SB >>>
3913	! paprs, pplay,zxtsol,albsol1, albsol2, &
3914	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3915	!albedo SB <<<
3916	t_seri,q_seri,wo, &
3917	cldfrarad, cldemirad, cldtaurad, &
3918	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
3919	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
3920	tau_aero, piz_aero, cg_aero, &
3921	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3922	! Rajoute par OB pour RRTM
3923	tau_aero_lw_rrtm, &
3924	cldtaupirad,new_aod, &
3925	! zqsat, flwcrad, fiwcrad, &
3926	zqsat, flwc, fiwc, &
3927	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3928	heat,heat0,cool,cool0,albpla, &
3929	heat_volc,cool_volc, &
3930	topsw,toplw,solsw,sollw, &
3931	sollwdown, &
3932	topsw0,toplw0,solsw0,sollw0, &
3933	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
3934	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
3935	topswad_aero, solswad_aero, &
3936	topswai_aero, solswai_aero, &
3937	topswad0_aero, solswad0_aero, &
3938	topsw_aero, topsw0_aero, &
3939	solsw_aero, solsw0_aero, &
3940	topswcf_aero, solswcf_aero, &
3941	!-C. Kleinschmitt for LW diagnostics
3942	toplwad_aero, sollwad_aero,&
3943	toplwai_aero, sollwai_aero, &
3944	toplwad0_aero, sollwad0_aero,&
3945	!-end
3946	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3947	ZSWFT0_i, ZFSDN0, ZFSUP0)
3948
3949	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
3950	!schemes
3951	toplw = toplw + betalwoff * (toplw0 - toplw)
3952	sollw = sollw + betalwoff * (sollw0 - sollw)
3953	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
3954	lwup = lwup + betalwoff * (lwup0 - lwup)
3955	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
3956	sollwdown(:))
3957	cool = cool + betalwoff * (cool0 - cool)
3958
3959	#ifndef CPP_XIOS
3960	!--OB 30/05/2016 modified 21/10/2016
3961	!--here we return swaero_diag and dryaod_diag to FALSE
3962	!--and histdef will switch it back to TRUE if necessary
3963	!--this is necessary to get the right swaero at first step
3964	!--but only in the case of no XIOS as XIOS is covered elsewhere
3965	IF (debut) swaerofree_diag = .FALSE.
3966	IF (debut) swaero_diag = .FALSE.
3967	IF (debut) dryaod_diag = .FALSE.
3968	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
3969	!--as for swaero_diag, see above
3970	IF (debut) ok_4xCO2atm = .FALSE.
3971
3972	!
3973	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3974	!IM des taux doit etre different du taux actuel
3975	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3976	!
3977	IF (RCO2_per.NE.RCO2_act.OR. &
3978	RCH4_per.NE.RCH4_act.OR. &
3979	RN2O_per.NE.RN2O_act.OR. &
3980	RCFC11_per.NE.RCFC11_act.OR. &
3981	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
3982	#endif
3983	!
3984	IF (ok_4xCO2atm) THEN
3985	!
3986	RCO2 = RCO2_per
3987	RCH4 = RCH4_per
3988	RN2O = RN2O_per
3989	RCFC11 = RCFC11_per
3990	RCFC12 = RCFC12_per
3991	!
3992	IF (prt_level .GE.10) THEN
3993	print *,' ->radlwsw, number 2 '
3994	ENDIF
3995	!
3996	CALL radlwsw &
3997	(dist, rmu0, fract, &
3998	!albedo SB >>>
3999	! paprs, pplay,zxtsol,albsol1, albsol2, &
4000	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4001	!albedo SB <<<
4002	t_seri,q_seri,wo, &
4003	cldfrarad, cldemirad, cldtaurad, &
4004	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4005	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4006	tau_aero, piz_aero, cg_aero, &
4007	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4008	! Rajoute par OB pour RRTM
4009	tau_aero_lw_rrtm, &
4010	cldtaupi,new_aod, &
4011	! zqsat, flwcrad, fiwcrad, &
4012	zqsat, flwc, fiwc, &
4013	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4014	heatp,heat0p,coolp,cool0p,albplap, &
4015	heat_volc,cool_volc, &
4016	topswp,toplwp,solswp,sollwp, &
4017	sollwdownp, &
4018	topsw0p,toplw0p,solsw0p,sollw0p, &
4019	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4020	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4021	topswad_aerop, solswad_aerop, &
4022	topswai_aerop, solswai_aerop, &
4023	topswad0_aerop, solswad0_aerop, &
4024	topsw_aerop, topsw0_aerop, &
4025	solsw_aerop, solsw0_aerop, &
4026	topswcf_aerop, solswcf_aerop, &
4027	!-C. Kleinschmitt for LW diagnostics
4028	toplwad_aerop, sollwad_aerop,&
4029	toplwai_aerop, sollwai_aerop, &
4030	toplwad0_aerop, sollwad0_aerop,&
4031	!-end
4032	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4033	ZSWFT0_i, ZFSDN0, ZFSUP0)
4034	endif !ok_4xCO2atm
4035	ENDIF ! aerosol_couple
4036	itaprad = 0
4037	!
4038	! If Iflag_radia >=2, reset pertubed variables
4039	!
4040	IF (iflag_radia .ge. 2) THEN
4041	zxtsol(:) = zsav_tsol (:)
4042	ENDIF
4043	ENDIF ! MOD(itaprad,radpas)
4044	itaprad = itaprad + 1
4045
4046	IF (iflag_radia.eq.0) THEN
4047	IF (prt_level.ge.9) THEN
4048	PRINT *,'--------------------------------------------------'
4049	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4050	PRINT *,'>>>> heat et cool mis a zero '
4051	PRINT *,'--------------------------------------------------'
4052	ENDIF
4053	heat=0.
4054	cool=0.
4055	sollw=0. ! MPL 01032011
4056	solsw=0.
4057	radsol=0.
4058	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4059	swup0=0.
4060	lwup=0.
4061	lwup0=0.
4062	lwdn=0.
4063	lwdn0=0.
4064	ENDIF
4065
4066	!
4067	! Calculer radsol a l'exterieur de radlwsw
4068	! pour prendre en compte le cycle diurne
4069	! recode par Olivier Boucher en sept 2015
4070	!
4071	radsol=solsw*swradcorr+sollw
4072
4073	IF (ok_4xCO2atm) THEN
4074	radsolp=solswp*swradcorr+sollwp
4075	ENDIF
4076
4077	!
4078	! Ajouter la tendance des rayonnements (tous les pas)
4079	! avec une correction pour le cycle diurne dans le SW
4080	!
4081
4082	DO k=1, klev
4083	d_t_swr(:,k)=swradcorr(:)heat(:,k)dtime/RDAY
4084	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)dtime/RDAY
4085	d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
4086	d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
4087	ENDDO
4088
4089	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4090	call prt_enerbil('SW',itap)
4091	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4092	call prt_enerbil('LW',itap)
4093
4094	!
4095	IF (mydebug) THEN
4096	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4097	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4098	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4099	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4100	ENDIF
4101
4102	! Calculer l'hydrologie de la surface
4103	!
4104	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4105	! . agesno, ftsol,fqsurf,fsnow, ruis)
4106	!
4107
4108	!
4109	! Calculer le bilan du sol et la derive de temperature (couplage)
4110	!
4111	DO i = 1, klon
4112	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4113	! a la demande de JLD
4114	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4115	ENDDO
4116	!
4117	!moddeblott(jan95)
4118	! Appeler le programme de parametrisation de l'orographie
4119	! a l'echelle sous-maille:
4120	!
4121	IF (prt_level .GE.10) THEN
4122	print *,' call orography ? ', ok_orodr
4123	ENDIF
4124	!
4125	IF (ok_orodr) THEN
4126	!
4127	! selection des points pour lesquels le shema est actif:
4128	igwd=0
4129	DO i=1,klon
4130	itest(i)=0
4131	! IF ((zstd(i).gt.10.0)) THEN
4132	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4133	itest(i)=1
4134	igwd=igwd+1
4135	idx(igwd)=i
4136	ENDIF
4137	ENDDO
4138	! igwdim=MAX(1,igwd)
4139	!
4140	IF (ok_strato) THEN
4141
4142	CALL drag_noro_strato(0,klon,klev,dtime,paprs,pplay, &
4143	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4144	igwd,idx,itest, &
4145	t_seri, u_seri, v_seri, &
4146	zulow, zvlow, zustrdr, zvstrdr, &
4147	d_t_oro, d_u_oro, d_v_oro)
4148
4149	ELSE
4150	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
4151	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4152	igwd,idx,itest, &
4153	t_seri, u_seri, v_seri, &
4154	zulow, zvlow, zustrdr, zvstrdr, &
4155	d_t_oro, d_u_oro, d_v_oro)
4156	ENDIF
4157	!
4158	! ajout des tendances
4159	!-----------------------------------------------------------------------
4160	! ajout des tendances de la trainee de l'orographie
4161	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4162	abortphy,flag_inhib_tend,itap,0)
4163	call prt_enerbil('oro',itap)
4164	!----------------------------------------------------------------------
4165	!
4166	ENDIF ! fin de test sur ok_orodr
4167	!
4168	IF (mydebug) THEN
4169	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4170	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4171	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4172	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4173	ENDIF
4174
4175	IF (ok_orolf) THEN
4176	!
4177	! selection des points pour lesquels le shema est actif:
4178	igwd=0
4179	DO i=1,klon
4180	itest(i)=0
4181	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4182	itest(i)=1
4183	igwd=igwd+1
4184	idx(igwd)=i
4185	ENDIF
4186	ENDDO
4187	! igwdim=MAX(1,igwd)
4188	!
4189	IF (ok_strato) THEN
4190
4191	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
4192	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4193	igwd,idx,itest, &
4194	t_seri, u_seri, v_seri, &
4195	zulow, zvlow, zustrli, zvstrli, &
4196	d_t_lif, d_u_lif, d_v_lif )
4197
4198	ELSE
4199	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
4200	latitude_deg,zmea,zstd,zpic, &
4201	itest, &
4202	t_seri, u_seri, v_seri, &
4203	zulow, zvlow, zustrli, zvstrli, &
4204	d_t_lif, d_u_lif, d_v_lif)
4205	ENDIF
4206
4207	! ajout des tendances de la portance de l'orographie
4208	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4209	'lif', abortphy,flag_inhib_tend,itap,0)
4210	call prt_enerbil('lif',itap)
4211	ENDIF ! fin de test sur ok_orolf
4212
4213	IF (ok_hines) then
4214	! HINES GWD PARAMETRIZATION
4215	east_gwstress=0.
4216	west_gwstress=0.
4217	du_gwd_hines=0.
4218	dv_gwd_hines=0.
4219	CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
4220	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4221	du_gwd_hines, dv_gwd_hines)
4222	zustr_gwd_hines=0.
4223	zvstr_gwd_hines=0.
4224	DO k = 1, klev
4225	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
4226	* (paprs(:, k)-paprs(:, k+1))/rg
4227	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
4228	* (paprs(:, k)-paprs(:, k+1))/rg
4229	ENDDO
4230
4231	d_t_hin(:, :)=0.
4232	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4233	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4234	call prt_enerbil('hin',itap)
4235	ENDIF
4236
4237	IF (.not. ok_hines .and. ok_gwd_rando) then
4238	CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
4239	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4240	dv_gwd_front, east_gwstress, west_gwstress)
4241	zustr_gwd_front=0.
4242	zvstr_gwd_front=0.
4243	DO k = 1, klev
4244	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
4245	* (paprs(:, k)-paprs(:, k+1))/rg
4246	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
4247	* (paprs(:, k)-paprs(:, k+1))/rg
4248	ENDDO
4249
4250	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4251	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4252	call prt_enerbil('front_gwd_rando',itap)
4253	ENDIF
4254
4255	IF (ok_gwd_rando) THEN
4256	CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
4257	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4258	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4259	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4260	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4261	call prt_enerbil('flott_gwd_rando',itap)
4262	zustr_gwd_rando=0.
4263	zvstr_gwd_rando=0.
4264	DO k = 1, klev
4265	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
4266	* (paprs(:, k)-paprs(:, k+1))/rg
4267	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
4268	* (paprs(:, k)-paprs(:, k+1))/rg
4269	ENDDO
4270	ENDIF
4271
4272	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4273
4274	IF (mydebug) THEN
4275	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4276	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4277	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4278	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4279	ENDIF
4280
4281	DO i = 1, klon
4282	zustrph(i)=0.
4283	zvstrph(i)=0.
4284	ENDDO
4285	DO k = 1, klev
4286	DO i = 1, klon
4287	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
4288	(paprs(i,k)-paprs(i,k+1))/rg
4289	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
4290	(paprs(i,k)-paprs(i,k+1))/rg
4291	ENDDO
4292	ENDDO
4293	!
4294	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4295	!
4296	IF (is_sequential .and. ok_orodr) THEN
4297	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4298	ra,rg,romega, &
4299	latitude_deg,longitude_deg,pphis, &
4300	zustrdr,zustrli,zustrph, &
4301	zvstrdr,zvstrli,zvstrph, &
4302	paprs,u,v, &
4303	aam, torsfc)
4304	ENDIF
4305	!IM cf. FLott END
4306	!DC Calcul de la tendance due au methane
4307	IF(ok_qch4) THEN
4308	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4309	! ajout de la tendance d'humidite due au methane
4310	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
4311	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4312	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4313	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
4314	ENDIF
4315	!
4316	!
4317
4318	!===============================================================
4319	! Additional tendency of TKE due to orography
4320	!===============================================================
4321	!
4322	! Inititialization
4323	!------------------
4324
4325
4326
4327	addtkeoro=0
4328	CALL getin_p('addtkeoro',addtkeoro)
4329
4330	IF (prt_level.ge.5) &
4331	print*,'addtkeoro', addtkeoro
4332
4333	alphatkeoro=1.
4334	CALL getin_p('alphatkeoro',alphatkeoro)
4335	alphatkeoro=min(max(0.,alphatkeoro),1.)
4336
4337	smallscales_tkeoro=.false.
4338	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4339
4340
4341	dtadd(:,:)=0.
4342	duadd(:,:)=0.
4343	dvadd(:,:)=0.
4344
4345
4346
4347	! Choices for addtkeoro:
4348	! ** 0 no TKE tendency from orography
4349	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4350	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4351	!
4352
4353	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4354	! -------------------------------------------
4355
4356
4357	! selection des points pour lesquels le schema est actif:
4358
4359
4360
4361	IF (addtkeoro .EQ. 1 ) THEN
4362
4363	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4364	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4365
4366	ELSE IF (addtkeoro .EQ. 2) THEN
4367
4368
4369
4370	IF (smallscales_tkeoro) THEN
4371	igwd=0
4372	DO i=1,klon
4373	itest(i)=0
4374	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4375	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4376	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4377	IF (zstd(i).GT.1.0) THEN
4378	itest(i)=1
4379	igwd=igwd+1
4380	idx(igwd)=i
4381	ENDIF
4382	ENDDO
4383
4384	ELSE
4385
4386	igwd=0
4387	DO i=1,klon
4388	itest(i)=0
4389	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4390	itest(i)=1
4391	igwd=igwd+1
4392	idx(igwd)=i
4393	ENDIF
4394	ENDDO
4395
4396	END IF
4397
4398
4399
4400
4401	CALL drag_noro_strato(addtkeoro,klon,klev,dtime,paprs,pplay, &
4402	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4403	igwd,idx,itest, &
4404	t_seri, u_seri, v_seri, &
4405	zulow, zvlow, zustrdr, zvstrdr, &
4406	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4407
4408	zustrdr(:)=0.
4409	zvstrdr(:)=0.
4410	zulow(:)=0.
4411	zvlow(:)=0.
4412
4413	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4414	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4415	END IF
4416
4417
4418
4419	! TKE update from subgrid temperature and wind tendencies
4420	!----------------------------------------------------------
4421	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4422
4423
4424	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4425
4426
4427
4428	ENDIF
4429	! -----
4430	!===============================================================
4431
4432
4433
4434	!====================================================================
4435	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4436	!====================================================================
4437	! Abderrahmane 24.08.09
4438
4439	IF (ok_cosp) THEN
4440	! adeclarer
4441	#ifdef CPP_COSP
4442	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4443
4444	IF (prt_level .GE.10) THEN
4445	print*,'freq_cosp',freq_cosp
4446	ENDIF
4447	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4448	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4449	! s ref_liq,ref_ice
4450	CALL phys_cosp(itap,dtime,freq_cosp, &
4451	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4452	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4453	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4454	JrNt,ref_liq,ref_ice, &
4455	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4456	zu10m,zv10m,pphis, &
4457	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4458	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4459	prfl(:,1:klev),psfl(:,1:klev), &
4460	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4461	mr_ozone,cldtau, cldemi)
4462
4463	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4464	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4465	! M clMISR,
4466	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4467	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4468
4469	ENDIF
4470
4471	#endif
4472	ENDIF !ok_cosp
4473
4474
4475	! Marine
4476
4477	IF (ok_airs) then
4478
4479	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
4480	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4481	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4482	& map_prop_hc,map_prop_hist,&
4483	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4484	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4485	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4486	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4487	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4488	& map_ntot,map_hc,map_hist,&
4489	& map_Cb,map_ThCi,map_Anv,&
4490	& alt_tropo )
4491	ENDIF
4492
4493	ENDIF ! ok_airs
4494
4495
4496	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4497	!AA
4498	!AA Installation de l'interface online-offline pour traceurs
4499	!AA
4500	!====================================================================
4501	! Calcul des tendances traceurs
4502	!====================================================================
4503	!
4504
4505	IF (type_trac=='repr') THEN
4506	sh_in(:,:) = q_seri(:,:)
4507	ELSE
4508	sh_in(:,:) = qx(:,:,ivap)
4509	ch_in(:,:) = qx(:,:,iliq)
4510	ENDIF
4511
4512	IF (iflag_phytrac == 1 ) THEN
4513
4514	#ifdef CPP_Dust
4515	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4516	pdtphys,ftsol, & ! I
4517	t,q_seri,paprs,pplay,RHcl, & ! I
4518	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4519	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4520	u_seri, v_seri, latitude_deg, longitude_deg, &
4521	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4522	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4523	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4524	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4525	fm_therm, entr_therm, rneb, & ! I
4526	beta_prec_fisrt,beta_prec, & !I
4527	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4528	d_tr_dyn,tr_seri)
4529
4530	#else
4531
4532	CALL phytrac ( &
4533	itap, days_elapsed+1, jH_cur, debut, &
4534	lafin, dtime, u, v, t, &
4535	paprs, pplay, pmfu, pmfd, &
4536	pen_u, pde_u, pen_d, pde_d, &
4537	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4538	u1, v1, ftsol, pctsrf, &
4539	zustar, zu10m, zv10m, &
4540	wstar(:,is_ave), ale_bl, ale_wake, &
4541	latitude_deg, longitude_deg, &
4542	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4543	presnivs, pphis, pphi, albsol1, &
4544	sh_in, ch_in, rhcl, cldfra, rneb, &
4545	diafra, cldliq, itop_con, ibas_con, &
4546	pmflxr, pmflxs, prfl, psfl, &
4547	da, phi, mp, upwd, &
4548	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4549	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4550	ev, ep, epmlmMm, eplaMm, & !<<RomP
4551	dnwd, aerosol_couple, flxmass_w, &
4552	tau_aero, piz_aero, cg_aero, ccm, &
4553	rfname, &
4554	d_tr_dyn, & !<<RomP
4555	tr_seri, init_source)
4556	#endif
4557	ENDIF ! (iflag_phytrac=1)
4558
4559	IF (offline) THEN
4560
4561	IF (prt_level.ge.9) &
4562	print*,'Attention on met a 0 les thermiques pour phystoke'
4563	CALL phystokenc ( &
4564	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4565	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4566	fm_therm,entr_therm, &
4567	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4568	frac_impa, frac_nucl, &
4569	pphis,cell_area,dtime,itap, &
4570	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4571
4572
4573	ENDIF
4574
4575	!
4576	! Calculer le transport de l'eau et de l'energie (diagnostique)
4577	!
4578	CALL transp (paprs,zxtsol, &
4579	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4580	ve, vq, ue, uq, vwat, uwat)
4581	!
4582	!IM global posePB BEG
4583	IF(1.EQ.0) THEN
4584	!
4585	CALL transp_lay (paprs,zxtsol, &
4586	t_seri, q_seri, u_seri, v_seri, zphi, &
4587	ve_lay, vq_lay, ue_lay, uq_lay)
4588	!
4589	ENDIF !(1.EQ.0) THEN
4590	!IM global posePB END
4591	! Accumuler les variables a stocker dans les fichiers histoire:
4592	!
4593
4594	!================================================================
4595	! Conversion of kinetic and potential energy into heat, for
4596	! parameterisation of subgrid-scale motions
4597	!================================================================
4598
4599	d_t_ec(:,:)=0.
4600	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4601	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4602	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4603	zmasse,exner,d_t_ec)
4604	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4605
4606	!=======================================================================
4607	! SORTIES
4608	!=======================================================================
4609	!
4610	!IM initialisation + calculs divers diag AMIP2
4611	!
4612	include "calcul_divers.h"
4613	!
4614	!IM Interpolation sur les niveaux de pression du NMC
4615	! -------------------------------------------------
4616	!
4617	! include "calcul_STDlev.h"
4618	!
4619	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4620	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4621	!
4622	!cc prw = eau precipitable
4623	! prlw = colonne eau liquide
4624	! prlw = colonne eau solide
4625	prw(:) = 0.
4626	prlw(:) = 0.
4627	prsw(:) = 0.
4628	DO k = 1, klev
4629	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4630	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4631	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4632	ENDDO
4633	!
4634	IF (type_trac == 'inca') THEN
4635	#ifdef INCA
4636	CALL VTe(VTphysiq)
4637	CALL VTb(VTinca)
4638
4639	CALL chemhook_end ( &
4640	dtime, &
4641	pplay, &
4642	t_seri, &
4643	tr_seri, &
4644	nbtr, &
4645	paprs, &
4646	q_seri, &
4647	cell_area, &
4648	pphi, &
4649	pphis, &
4650	zx_rh, &
4651	aps, bps)
4652
4653	CALL VTe(VTinca)
4654	CALL VTb(VTphysiq)
4655	#endif
4656	ENDIF
4657
4658
4659	!
4660	! Convertir les incrementations en tendances
4661	!
4662	IF (prt_level .GE.10) THEN
4663	print *,'Convertir les incrementations en tendances '
4664	ENDIF
4665	!
4666	IF (mydebug) THEN
4667	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4668	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4669	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4670	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4671	ENDIF
4672
4673	DO k = 1, klev
4674	DO i = 1, klon
4675	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4676	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4677	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4678	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4679	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4680	!CR: on ajoute le contenu en glace
4681	IF (nqo.eq.3) THEN
4682	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
4683	ENDIF
4684	ENDDO
4685	ENDDO
4686	!
4687	!CR: nb de traceurs eau: nqo
4688	! IF (nqtot.GE.3) THEN
4689	IF (nqtot.GE.(nqo+1)) THEN
4690	! DO iq = 3, nqtot
4691	DO iq = nqo+1, nqtot
4692	DO k = 1, klev
4693	DO i = 1, klon
4694	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4695	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4696	ENDDO
4697	ENDDO
4698	ENDDO
4699	ENDIF
4700	!
4701	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4702	!IM global posePB include "write_bilKP_ins.h"
4703	!IM global posePB include "write_bilKP_ave.h"
4704	!
4705
4706	!--OB mass fixer
4707	!--profile is corrected to force mass conservation of water
4708	IF (mass_fixer) THEN
4709	qql2(:)=0.0
4710	DO k = 1, klev
4711	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4712	ENDDO
4713	DO i = 1, klon
4714	!--compute ratio of what q+ql should be with conservation to what it is
4715	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4716	DO k = 1, klev
4717	q_seri(i,k) =q_seri(i,k)*corrqql
4718	ql_seri(i,k)=ql_seri(i,k)*corrqql
4719	ENDDO
4720	ENDDO
4721	ENDIF
4722	!--fin mass fixer
4723
4724	! Sauvegarder les valeurs de t et q a la fin de la physique:
4725	!
4726	u_ancien(:,:) = u_seri(:,:)
4727	v_ancien(:,:) = v_seri(:,:)
4728	t_ancien(:,:) = t_seri(:,:)
4729	q_ancien(:,:) = q_seri(:,:)
4730	ql_ancien(:,:) = ql_seri(:,:)
4731	qs_ancien(:,:) = qs_seri(:,:)
4732	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4733	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4734	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4735	! !! RomP >>>
4736	!CR: nb de traceurs eau: nqo
4737	IF (nqtot.GT.nqo) THEN
4738	DO iq = nqo+1, nqtot
4739	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4740	ENDDO
4741	ENDIF
4742	! !! RomP <<<
4743	!==========================================================================
4744	! Sorties des tendances pour un point particulier
4745	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4746	! pour le debug
4747	! La valeur de igout est attribuee plus haut dans le programme
4748	!==========================================================================
4749
4750	IF (prt_level.ge.1) THEN
4751	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4752	write(lunout,*) &
4753	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4754	write(lunout,*) &
4755	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4756	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4757	pctsrf(igout,is_sic)
4758	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4759	DO k=1,klev
4760	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4761	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4762	d_t_eva(igout,k)
4763	ENDDO
4764	write(lunout,*) 'cool,heat'
4765	DO k=1,klev
4766	write(lunout,*) cool(igout,k),heat(igout,k)
4767	ENDDO
4768
4769	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4770	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4771	!jyg! do k=1,klev
4772	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4773	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4774	!jyg! enddo
4775	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4776	DO k=1,klev
4777	write(lunout,*) d_t_vdf(igout,k), &
4778	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4779	ENDDO
4780	!>jyg
4781
4782	write(lunout,*) 'd_ps ',d_ps(igout)
4783	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4784	DO k=1,klev
4785	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4786	d_qx(igout,k,1),d_qx(igout,k,2)
4787	ENDDO
4788	ENDIF
4789
4790	!============================================================
4791	! Calcul de la temperature potentielle
4792	!============================================================
4793	DO k = 1, klev
4794	DO i = 1, klon
4795	!JYG/IM theta en debut du pas de temps
4796	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4797	!JYG/IM theta en fin de pas de temps de physique
4798	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4799	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4800	! MPL 20130625
4801	! fth_fonctions.F90 et parkind1.F90
4802	! sinon thetal=theta
4803	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4804	! : ql_seri(i,k))
4805	thetal(i,k)=theta(i,k)
4806	ENDDO
4807	ENDDO
4808	!
4809
4810	! 22.03.04 BEG
4811	!=============================================================
4812	! Ecriture des sorties
4813	!=============================================================
4814	#ifdef CPP_IOIPSL
4815
4816	! Recupere des varibles calcule dans differents modules
4817	! pour ecriture dans histxxx.nc
4818
4819	! Get some variables from module fonte_neige_mod
4820	CALL fonte_neige_get_vars(pctsrf, &
4821	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4822
4823
4824	!=============================================================
4825	! Separation entre thermiques et non thermiques dans les sorties
4826	! de fisrtilp
4827	!=============================================================
4828
4829	IF (iflag_thermals>=1) THEN
4830	d_t_lscth=0.
4831	d_t_lscst=0.
4832	d_q_lscth=0.
4833	d_q_lscst=0.
4834	DO k=1,klev
4835	DO i=1,klon
4836	IF (ptconvth(i,k)) THEN
4837	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4838	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4839	ELSE
4840	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4841	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4842	ENDIF
4843	ENDDO
4844	ENDDO
4845
4846	DO i=1,klon
4847	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4848	plul_th(i)=prfl(i,1)+psfl(i,1)
4849	ENDDO
4850	ENDIF
4851
4852	!On effectue les sorties:
4853
4854	#ifdef CPP_Dust
4855	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
4856	pplay, lmax_th, aerosol_couple, &
4857	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4858	ptconv, read_climoz, clevSTD, &
4859	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
4860	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4861	#else
4862	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4863	pplay, lmax_th, aerosol_couple, &
4864	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, &
4865	ok_sync, ptconv, read_climoz, clevSTD, &
4866	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
4867	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4868	#endif
4869
4870	#ifndef CPP_XIOS
4871	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4872	#endif
4873
4874	#endif
4875
4876	! On remet des variables a .false. apres un premier appel
4877	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4878	if (debut) then
4879	#ifdef CPP_XIOS
4880	swaero_diag=.FALSE.
4881	swaerofree_diag=.FALSE.
4882	dryaod_diag=.FALSE.
4883	ok_4xCO2atm= .FALSE.
4884	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4885
4886	IF (is_master) then
4887	!--setting up swaero_diag to TRUE in XIOS case
4888	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
4889	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
4890	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
4891	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
4892	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
4893	!!!--for now these fields are not in the XML files so they are omitted
4894	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
4895	swaero_diag=.TRUE.
4896
4897	!--setting up swaerofree_diag to TRUE in XIOS case
4898	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
4899	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
4900	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
4901	xios_field_is_active("LWupTOAcleanclr")) &
4902	swaerofree_diag=.TRUE.
4903
4904	!--setting up dryaod_diag to TRUE in XIOS case
4905	DO naero = 1, naero_tot-1
4906	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
4907	ENDDO
4908	!
4909	!--setting up ok_4xCO2atm to TRUE in XIOS case
4910	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
4911	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
4912	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
4913	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
4914	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
4915	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
4916	ok_4xCO2atm=.TRUE.
4917	endif
4918	!$OMP BARRIER
4919	call bcast(swaero_diag)
4920	call bcast(swaerofree_diag)
4921	call bcast(dryaod_diag)
4922	call bcast(ok_4xCO2atm)
4923	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4924	#endif
4925	endif
4926	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4927
4928	!====================================================================
4929	! Arret du modele apres hgardfou en cas de detection d'un
4930	! plantage par hgardfou
4931	!====================================================================
4932
4933	IF (abortphy==1) THEN
4934	abort_message ='Plantage hgardfou'
4935	CALL abort_physic (modname,abort_message,1)
4936	ENDIF
4937
4938	! 22.03.04 END
4939	!
4940	!====================================================================
4941	! Si c'est la fin, il faut conserver l'etat de redemarrage
4942	!====================================================================
4943	!
4944
4945	IF (lafin) THEN
4946	itau_phy = itau_phy + itap
4947	CALL phyredem ("restartphy.nc")
4948	! open(97,form="unformatted",file="finbin")
4949	! write(97) u_seri,v_seri,t_seri,q_seri
4950	! close(97)
4951	!$OMP MASTER
4952	IF (read_climoz >= 1) THEN
4953	IF (is_mpi_root) THEN
4954	CALL nf95_close(ncid_climoz)
4955	ENDIF
4956	DEALLOCATE(press_edg_climoz) ! pointer
4957	DEALLOCATE(press_cen_climoz) ! pointer
4958	ENDIF
4959	!$OMP END MASTER
4960	print *,' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
4961	ENDIF
4962
4963	! first=.false.
4964
4965
4966	END SUBROUTINE physiq
4967
4968	END MODULE physiq_mod

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

Download in other formats:

Original Format