Context Navigation

physiq.F90 @ 2339

Last change on this file since 2339 was 2333, checked in by lguez, 9 years ago

New parameterization of gravity wave drag due to front/jet systems, by

de la Camara and F. Lott. The new Camara-Lott parameterization

replaces the Hines parameterization so it is activated if not ok_hines
and ok_gwd_rando.

Also changed distribution of phase speeds in FLOTT_GWD_rando, from
uniform to Gaussian. Bug fix in sugwd_strato. Bug fix in the arguments
of the call to add_phys_tend for methane oxydation.

For the new Camara-Lott parameterization, we need to compute relative
vorticity in calfis and pass it as a new argument "rot" to
physiq. Interpolation of relative vorticity to the physics grid is not
optimal for now: it is not weighted by cell areas.

Alvaro de la Camara, Fran\c{}cois Lott

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

File size: 141.8 KB

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

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

Download in other formats:

Original Format