Context Navigation

physiq_mod.F90 @ 5017

Last change on this file since 5017 was 3413, checked in by Laurent Fairhead, 6 years ago
Inclusion of Yann's latest (summer/fall 2018) modifications for convergence of DYNAMICO/LMDZ physics YM/LF
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 `Id`
File size: 166.2 KB

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

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

Download in other formats:

Original Format