Context Navigation

physiq.F90 @ 2344

Last change on this file since 2344 was 2344, checked in by Ehouarn Millour, 9 years ago
Physics/dynamics separation: get rid of all the 'include "temps.h"' in the physics; variables in module time_phylmdz_mod must be used instead. Also added JD_cur, JH_cur and JD_ref in module phys_cal_mod, in preparation for having physics handle its calendar internally. EM
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: 142.6 KB

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

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

Download in other formats:

Original Format