Context Navigation

physiq.F90 @ 2294

Last change on this file since 2294 was 2294, checked in by fhourdin, 9 years ago

Introduction d'un seuil sur la probabilité de non déclenchement
stochastique random_notrig_max
Si le tirage aleatoire uniforme entre 0 et 1 est > random_notrig_max
on ne déclenche pas.
random_notrig_max=1-epsilon avec epsilon petit.

Corrections dans grid_atob_m.F90 pour la compilation Linux.
A vérifier.

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

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

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

Download in other formats:

Original Format