Context Navigation

physiq_mod.F90 @ 3981

Last change on this file since 3981 was 3981, checked in by Laurent Fairhead, 3 years ago

Adds a routine 'prt_alerte' to enable developpers to print out informative messages on the first pass through physics
and from the master process only. This can be used to remind oneself of potential problems or further enhancements.
Messages can be differentiated by a 'priority' code 0/1/2 (corresponding to GREEN/ORANGE/RED alerts)
By default, the messages are output in a file called ALERTES.txt but changing the unit number to 6 in the
print_control_mod.F90 file allows you to print out your messages to the screen.

To use you simply need to:

have this USE statement at the start of your routine:

USE print_control_mod, ONLY: prt_level, lunout, call_alert, prt_alerte

ensure that the modname variable is defined and contains the name of your routine
then you just need to add the following lines in your routine wherever you want to ouput your message (preferably not in a do loop)

message='your informative message'
IF (call_alert) CALL prt_alerte(message,modname,CODE)

where

message is the message to print out
modname is the routine name
CODE is an integer representing your priority code (0/1/2)

and you should get messages of the form

ALERTE ROUGE cva_driver! ym missing init, need to have a look by developpers
ALERTE VERTE orosetup_strato! ym correction en attendant mieux

in your output file.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 186.7 KB

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

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

Download in other formats:

Original Format