Context Navigation

physiq_mod.F90 @ 5361

Last change on this file since 5361 was 2914, checked in by acozic, 8 years ago
Add an argument to chemini call ( = rev 2906 on trunk LMDZ5)
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: 160.8 KB

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

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

Download in other formats:

Original Format