Context Navigation

physiq_mod.F90 @ 3153

Last change on this file since 3153 was 3153, checked in by jyg, 6 years ago
bug fix on commit 3150, print added
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 168.1 KB

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

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

Download in other formats:

Original Format