Context Navigation

physiq_mod.F90 @ 2586

Last change on this file since 2586 was 2580, checked in by idelkadi, 8 years ago
Implementation du simulateur AIRS: Le but du simulateur est de permettre la comparaison de proprietes macro-et microphysiques des nuages hauts de LMDZ avec celles restituees par les observations du satellite AIRS (Atmospheric IR Sounder). La methode est decrite dans Hendricks et al. Meteorol. Z., 2010
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: 150.0 KB

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

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

Download in other formats:

Original Format