Context Navigation

physiq.F90 @ 2311

Last change on this file since 2311 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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

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

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

Download in other formats:

Original Format