Context Navigation

physiq_mod.F90 @ 2492

Last change on this file since 2492 was 2489, checked in by oboucher, 8 years ago
moved the water vapour mass fixer for the physics before q_seri is copied in q_ancien
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: 156.8 KB

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

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

Download in other formats:

Original Format