Context Navigation

physiq_mod.F90 @ 2421

Last change on this file since 2421 was 2421, checked in by Ehouarn Millour, 8 years ago
Stop sending pducov from dynamics to physics; it is not used in the physics anymore (and what was done was only valid when runnning in serial). EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 144.3 KB

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

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

Download in other formats:

Original Format