Context Navigation

physiq.F90 @ 2201

Last change on this file since 2201 was 2200, checked in by acozic, 9 years ago
Aco : add some variable in inca routine call to correct a bug on time_axis with xios output
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: 133.3 KB

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

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

Download in other formats:

Original Format