Context Navigation

physiq_mod.F90 @ 3134

Last change on this file since 3134 was 3134, checked in by jyg, 6 years ago

Bug fix: some communication variables between the
convective scheme and phytrac were set to zero at
the beginning of physiq (and others were not
saved) which led to an erroneous behaviour when
convection was not called at every time step. Now
all theses variables are in phys_local_var_mod and
are no longer reset to zero at the beginning of
physiq provided ok_bug_cv_trac=n .

If ok_bug_cv_trac=y, then the reset to zero is
still performed.

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: 166.4 KB

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

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

Download in other formats:

Original Format