Context Navigation

physiq_mod.F90 @ 3400

Last change on this file since 3400 was 3400, checked in by Laurent Fairhead, 6 years ago
Due to a difficult to correct bug, near surface fields can take abnormal values in very specific circumstances. While waiting to find a correction, these modifications will force the values of these fields within reasonable bounds. To get back to the previous behaviour of the model, use iflag_bug_t2m_ipslcm61 = 1 (it is 0 by default)
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: 171.0 KB

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

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

Download in other formats:

Original Format