Context Navigation

physiq.F90 @ 1915

Last change on this file since 1915 was 1915, checked in by musat, 11 years ago
D'autres oublis pour les sorties annuelles IM
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.3 KB

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

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

Download in other formats:

Original Format