Context Navigation

physiq.F90 @ 2160

Last change on this file since 2160 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
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: 128.2 KB

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

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

Download in other formats:

Original Format