Context Navigation

physiq_mod.F90 @ 2448

Last change on this file since 2448 was 2448, checked in by jyg, 8 years ago
Correction in the computation of qsat in physiq_mod.F90: the transition of humidity from saturation over ice to saturation over liquid water occurs now at rtt (=273.15 K) instead of t_glace_min or t_coup.
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: 144.7 KB

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

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

Download in other formats:

Original Format