Context Navigation

physiq.F90 @ 2259

Last change on this file since 2259 was 2259, checked in by jyg, 9 years ago

Changes in Emanuel's deep convection scheme: the
upper bound of deep convection loops is set at
the first level above 22 km.

Modified files:

physiq.F90,
concvl.F90,
cva_driver.F90,
cv3a_compress.F90,
cv3a_uncompress.F90,
cv3_routines.F90

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

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

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

Download in other formats:

Original Format