Context Navigation

physiq.F90 @ 1912

Last change on this file since 1912 was 1912, checked in by musat, 10 years ago

1) Modifications pour faire des simulations par an avec un calendrier realiste (365 jours ou autre).

Il faut mettre une frequence de sortie de -1 (variable phys_out_filetimesteps dans config.def) pour
que IOIPSL calcule les moyennes mensuels en prenant en compte des longuers variables de chaque
mois. Par exemple, pour le fichier histmth (1er fichier) et histmthNMC (7eme).
phys_out_filetimesteps= -1 1day 6hr 6hr 6hr 1d -1 1day 6hr

2) Corrections titres variables niveaux de pression des fichiers histmth, histday, etc

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

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

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

Download in other formats:

Original Format