Context Navigation

physiq_mod.F90 @ 2476

Last change on this file since 2476 was 2476, checked in by oboucher, 8 years ago
Introduction of a mass fixer for water in physiq int q+ql dz at the end = int q+ql dz at the beginning + evappdtphys - precippdtphys
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.9 KB

Line
1	!
2	! $Id: physiq_mod.F90 2476 2016-03-23 21:48:19Z oboucher $
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),ep(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
940	logical, parameter :: mass_fixer=.true.
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	zdz=(paprs(i,k)-paprs(i,k+1))/RG
1602	qql1(i)=qql1(i)+(q_seri(i,k)+ql_seri(i,k))*zdz
1603	ENDDO
1604	ENDDO
1605	ENDIF
1606	!--fin mass fixer
1607
1608	tke0(:,:)=pbl_tke(:,:,is_ave)
1609	!CR:Nombre de traceurs de l'eau: nqo
1610	! IF (nqtot.GE.3) THEN
1611	IF (nqtot.GE.(nqo+1)) THEN
1612	! DO iq = 3, nqtot
1613	DO iq = nqo+1, nqtot
1614	DO k = 1, klev
1615	DO i = 1, klon
1616	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1617	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1618	ENDDO
1619	ENDDO
1620	ENDDO
1621	ELSE
1622	DO k = 1, klev
1623	DO i = 1, klon
1624	tr_seri(i,k,1) = 0.0
1625	ENDDO
1626	ENDDO
1627	ENDIF
1628	!
1629	DO i = 1, klon
1630	ztsol(i) = 0.
1631	ENDDO
1632	DO nsrf = 1, nbsrf
1633	DO i = 1, klon
1634	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1635	ENDDO
1636	ENDDO
1637	!IM
1638	IF (ip_ebil_phy.ge.1) THEN
1639	ztit='after dynamic'
1640	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
1641	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1642	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1643	! Comme les tendances de la physique sont ajoute dans la dynamique,
1644	! on devrait avoir que la variation d'entalpie par la dynamique
1645	! est egale a la variation de la physique au pas de temps precedent.
1646	! Donc la somme de ces 2 variations devrait etre nulle.
1647	call diagphy(cell_area,ztit,ip_ebil_phy &
1648	, zero_v, zero_v, zero_v, zero_v, zero_v &
1649	, zero_v, zero_v, zero_v, ztsol &
1650	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1651	, fs_bound, fq_bound )
1652	END IF
1653
1654	! Diagnostiquer la tendance dynamique
1655	!
1656	IF (ancien_ok) THEN
1657	DO k = 1, klev
1658	DO i = 1, klon
1659	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1660	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1661	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1662	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1663	ENDDO
1664	ENDDO
1665	! !! RomP >>> td dyn traceur
1666	!! IF (nqtot.GE.3) THEN ! jyg
1667	!! DO iq = 3, nqtot ! jyg
1668	IF (nqtot.GE.nqo+1) THEN ! jyg
1669	DO iq = nqo+1, nqtot ! jyg
1670	DO k = 1, klev
1671	DO i = 1, klon
1672	!! d_tr_dyn(i,k,iq-2)= & ! jyg
1673	!! (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime ! jyg
1674	d_tr_dyn(i,k,iq-nqo)= & ! jyg
1675	(tr_seri(i,k,iq-nqo)-tr_ancien(i,k,iq-nqo))/dtime ! jyg
1676	! iiq=niadv(iq)
1677	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-nqo),"tra:",iq,&
1678	! tname(iiq)
1679	ENDDO
1680	ENDDO
1681	ENDDO
1682	ENDIF
1683	! !! RomP <<<
1684	ELSE
1685	DO k = 1, klev
1686	DO i = 1, klon
1687	d_u_dyn(i,k) = 0.0
1688	d_v_dyn(i,k) = 0.0
1689	d_t_dyn(i,k) = 0.0
1690	d_q_dyn(i,k) = 0.0
1691	ENDDO
1692	ENDDO
1693	! !! RomP >>> td dyn traceur
1694	!! IF (nqtot.GE.3) THEN ! jyg
1695	!! DO iq = 3, nqtot ! jyg
1696	IF (nqtot.GE.nqo+1) THEN ! jyg
1697	DO iq = nqo+1, nqtot ! jyg
1698	DO k = 1, klev
1699	DO i = 1, klon
1700	!! d_tr_dyn(i,k,iq-2)= 0.0 ! jyg
1701	d_tr_dyn(i,k,iq-nqo)= 0.0 ! jyg
1702	ENDDO
1703	ENDDO
1704	ENDDO
1705	ENDIF
1706	! !! RomP <<<
1707	ancien_ok = .TRUE.
1708	ENDIF
1709	!
1710	! Ajouter le geopotentiel du sol:
1711	!
1712	DO k = 1, klev
1713	DO i = 1, klon
1714	zphi(i,k) = pphi(i,k) + pphis(i)
1715	ENDDO
1716	ENDDO
1717	!
1718	! Verifier les temperatures
1719	!
1720	!IM BEG
1721	IF (check) THEN
1722	amn=MIN(ftsol(1,is_ter),1000.)
1723	amx=MAX(ftsol(1,is_ter),-1000.)
1724	DO i=2, klon
1725	amn=MIN(ftsol(i,is_ter),amn)
1726	amx=MAX(ftsol(i,is_ter),amx)
1727	ENDDO
1728	!
1729	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1730	ENDIF !(check) THEN
1731	!IM END
1732	!
1733	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
1734	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
1735
1736	!
1737	!IM BEG
1738	IF (check) THEN
1739	amn=MIN(ftsol(1,is_ter),1000.)
1740	amx=MAX(ftsol(1,is_ter),-1000.)
1741	DO i=2, klon
1742	amn=MIN(ftsol(i,is_ter),amn)
1743	amx=MAX(ftsol(i,is_ter),amx)
1744	ENDDO
1745	!
1746	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1747	ENDIF !(check) THEN
1748	!IM END
1749	!
1750	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1751	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1752	!
1753	if (read_climoz >= 1) then
1754	! Ozone from a file
1755	! Update required ozone index:
1756	ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
1757	if (ro3i == 361) ro3i = 360
1758	! (This should never occur, except perhaps because of roundup
1759	! error. See documentation.)
1760	if (ro3i /= co3i) then
1761	! Update ozone field:
1762	if (read_climoz == 1) then
1763	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1764	press_in_edg=press_climoz, paprs=paprs, v3=wo)
1765	else
1766	! read_climoz == 2
1767	call regr_pr_av(ncid_climoz, (/"tro3 ", &
1768	"tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
1769	paprs=paprs, v3=wo)
1770	end if
1771	! Convert from mole fraction of ozone to column density of ozone in a
1772	! cell, in kDU:
1773	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1774	* zmasse / dobson_u / 1e3
1775	! (By regridding ozone values for LMDZ only once every 360th of
1776	! year, we have already neglected the variation of pressure in one
1777	! 360th of year. So do not recompute "wo" at each time step even if
1778	! "zmasse" changes a little.)
1779	co3i = ro3i
1780	end if
1781	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1782	! Once per day, update ozone from Royer:
1783
1784	IF (solarlong0<-999.) then
1785	! Generic case with evolvoing season
1786	zzz=real(days_elapsed+1)
1787	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1788	! Particular case with annual mean insolation
1789	zzz=real(90) ! could be revisited
1790	IF (read_climoz/=-1) THEN
1791	abort_message ='read_climoz=-1 is recommended when ' &
1792	// 'solarlong0=1000.'
1793	CALL abort_physic (modname,abort_message,1)
1794	ENDIF
1795	ELSE
1796	! Case where the season is imposed with solarlong0
1797	zzz=real(90) ! could be revisited
1798	ENDIF
1799	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
1800	ENDIF
1801	!
1802	! Re-evaporer l'eau liquide nuageuse
1803	!
1804	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1805	DO i = 1, klon
1806	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1807	!jyg<
1808	! Attention : Arnaud a propose des formules completement differentes
1809	! A verifier !!!
1810	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1811	IF (iflag_ice_thermo .EQ. 0) THEN
1812	zlsdcp=zlvdcp
1813	ENDIF
1814	!>jyg
1815
1816	if (iflag_ice_thermo.eq.0) then
1817	!pas necessaire a priori
1818
1819	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1820	zb = MAX(0.0,ql_seri(i,k))
1821	za = - MAX(0.0,ql_seri(i,k)) &
1822	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1823	t_seri(i,k) = t_seri(i,k) + za
1824	q_seri(i,k) = q_seri(i,k) + zb
1825	ql_seri(i,k) = 0.0
1826	d_t_eva(i,k) = za
1827	d_q_eva(i,k) = zb
1828
1829	else
1830
1831	!CR: on r\'e-\'evapore eau liquide et glace
1832
1833	! zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1834	! zb = MAX(0.0,ql_seri(i,k))
1835	! za = - MAX(0.0,ql_seri(i,k)) &
1836	! * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1837	zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
1838	za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
1839	- MAX(0.0,qs_seri(i,k))*zlsdcp
1840	t_seri(i,k) = t_seri(i,k) + za
1841	q_seri(i,k) = q_seri(i,k) + zb
1842	ql_seri(i,k) = 0.0
1843	!on \'evapore la glace
1844	qs_seri(i,k) = 0.0
1845	d_t_eva(i,k) = za
1846	d_q_eva(i,k) = zb
1847	endif
1848
1849	ENDDO
1850	ENDDO
1851	!IM
1852	IF (ip_ebil_phy.ge.2) THEN
1853	ztit='after reevap'
1854	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,1,dtime &
1855	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1856	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1857	call diagphy(cell_area,ztit,ip_ebil_phy &
1858	, zero_v, zero_v, zero_v, zero_v, zero_v &
1859	, zero_v, zero_v, zero_v, ztsol &
1860	, d_h_vcol, d_qt, d_ec &
1861	, fs_bound, fq_bound )
1862	!
1863	END IF
1864
1865	!
1866	!=========================================================================
1867	! Calculs de l'orbite.
1868	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1869	! doit donc etre plac\'e avant radlwsw et pbl_surface
1870
1871	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1872	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1873	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1874	!
1875	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1876	! solarlong0
1877	if (solarlong0<-999.) then
1878	if (new_orbit) then
1879	! calcul selon la routine utilisee pour les planetes
1880	call solarlong(day_since_equinox, zlongi, dist)
1881	else
1882	! calcul selon la routine utilisee pour l'AR4
1883	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1884	endif
1885	else
1886	zlongi=solarlong0 ! longitude solaire vraie
1887	dist=1. ! distance au soleil / moyenne
1888	endif
1889	if(prt_level.ge.1) &
1890	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1891
1892
1893	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1894	! Calcul de l'ensoleillement :
1895	! ============================
1896	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1897	! l'annee a partir d'une formule analytique.
1898	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1899	! non nul aux poles.
1900	IF (abs(solarlong0-1000.)<1.e-4) then
1901	call zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
1902	latitude_deg,longitude_deg,rmu0,fract)
1903	JrNt = 1.0
1904	ELSE
1905	! recode par Olivier Boucher en sept 2015
1906	SELECT CASE (iflag_cycle_diurne)
1907	CASE(0)
1908	! Sans cycle diurne
1909	CALL angle(zlongi, latitude_deg, fract, rmu0)
1910	swradcorr = 1.0
1911	JrNt = 1.0
1912	zrmu0 = rmu0
1913	CASE(1)
1914	! Avec cycle diurne sans application des poids
1915	! bit comparable a l ancienne formulation cycle_diurne=true
1916	! on integre entre gmtime et gmtime+radpas
1917	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1918	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
1919	latitude_deg,longitude_deg,rmu0,fract)
1920	zrmu0 = rmu0
1921	swradcorr = 1.0
1922	! Calcul du flag jour-nuit
1923	JrNt = 0.0
1924	WHERE (fract.GT.0.0) JrNt = 1.0
1925	CASE(2)
1926	! Avec cycle diurne sans application des poids
1927	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
1928	! Comme cette routine est appele a tous les pas de temps de
1929	! la physique meme si le rayonnement n'est pas appele je
1930	! remonte en arriere les radpas-1 pas de temps
1931	! suivant. Petite ruse avec MOD pour prendre en compte le
1932	! premier pas de temps de la physique pendant lequel
1933	! itaprad=0
1934	zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)
1935	zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
1936	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
1937	latitude_deg,longitude_deg,rmu0,fract)
1938	!
1939	! Calcul des poids
1940	!
1941	zdtime1=-dtime !--on corrige le rayonnement pour representer le
1942	zdtime2=0.0 !--pas de temps de la physique qui se termine
1943	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
1944	latitude_deg,longitude_deg,zrmu0,zfract)
1945	swradcorr = 0.0
1946	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
1947	swradcorr=zfract/fract*zrmu0/rmu0
1948	! Calcul du flag jour-nuit
1949	JrNt = 0.0
1950	WHERE (zfract.GT.0.0) JrNt = 1.0
1951	END SELECT
1952	ENDIF
1953
1954	if (mydebug) then
1955	call writefield_phy('u_seri',u_seri,nbp_lev)
1956	call writefield_phy('v_seri',v_seri,nbp_lev)
1957	call writefield_phy('t_seri',t_seri,nbp_lev)
1958	call writefield_phy('q_seri',q_seri,nbp_lev)
1959	endif
1960
1961	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1962	! Appel au pbl_surface : Planetary Boudary Layer et Surface
1963	! Cela implique tous les interactions des sous-surfaces et la
1964	! partie diffusion turbulent du couche limit.
1965	!
1966	! Certains varibales de sorties de pbl_surface sont utiliser que pour
1967	! ecriture des fihiers hist_XXXX.nc, ces sont :
1968	! qsol, zq2m, s_pblh, s_lcl,
1969	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1970	! s_therm, s_trmb1, s_trmb2, s_trmb3,
1971	! zu10m, zv10m, fder,
1972	! zxqsurf, rh2m, zxfluxu, zxfluxv,
1973	! frugs, agesno, fsollw, fsolsw,
1974	! d_ts, fevap, fluxlat, t2m,
1975	! wfbils, wfbilo, fluxt, fluxu, fluxv,
1976	!
1977	! Certains ne sont pas utiliser du tout :
1978	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1979	!
1980
1981	! Calcul de l'humidite de saturation au niveau du sol
1982
1983
1984
1985	if (iflag_pbl/=0) then
1986
1987	!jyg+nrlmd<
1988	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
1989	print *,'debut du splitting de la PBL'
1990	ENDIF
1991	! !!
1992	!=================================================================
1993	! PROVISOIRE : DECOUPLAGE PBL/WAKE
1994	! --------------------------------
1995	!
1996	!! wake_deltat_sav(:,:)=wake_deltat(:,:)
1997	!! wake_deltaq_sav(:,:)=wake_deltaq(:,:)
1998	!! wake_deltat(:,:)=0.
1999	!! wake_deltaq(:,:)=0.
2000	!=================================================================
2001	!>jyg+nrlmd
2002	!
2003	!-------gustiness calculation-------!
2004	IF (iflag_gusts==0) THEN
2005	gustiness(1:klon)=0
2006	ELSE IF (iflag_gusts==1) THEN
2007	do i = 1, klon
2008	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
2009	enddo
2010	! ELSE IF (iflag_gusts==2) THEN
2011	! do i = 1, klon
2012	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2013	! *ale_wake(i) !! need to make sigma_wk accessible here
2014	! enddo
2015	! ELSE IF (iflag_gusts==3) THEN
2016	! do i = 1, klon
2017	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2018	! enddo
2019	ENDIF
2020
2021
2022
2023	CALL pbl_surface( &
2024	dtime, date0, itap, days_elapsed+1, &
2025	debut, lafin, &
2026	longitude_deg, latitude_deg, rugoro, zrmu0, &
2027	zsig, sollwdown, pphi, cldt, &
2028	rain_fall, snow_fall, solsw, sollw, &
2029	gustiness, &
2030	t_seri, q_seri, u_seri, v_seri, &
2031	!nrlmd+jyg<
2032	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2033	!>nrlmd+jyg
2034	pplay, paprs, pctsrf, &
2035	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2036	!albedo SB <<<
2037	cdragh, cdragm, u1, v1, &
2038	!albedo SB >>>
2039	! albsol1, albsol2, sens, evap, &
2040	albsol_dir, albsol_dif, sens, evap, &
2041	!albedo SB <<<
2042	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2043	zxtsol, zxfluxlat, zt2m, qsat2m, &
2044	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2045	!nrlmd<
2046	!jyg<
2047	d_t_vdf_w, d_q_vdf_w, &
2048	d_t_vdf_x, d_q_vdf_x, &
2049	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2050	!>jyg
2051	delta_tsurf,wake_dens, &
2052	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2053	kh,kh_x,kh_w, &
2054	!>nrlmd
2055	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2056	slab_wfbils, &
2057	qsol, zq2m, s_pblh, s_lcl, &
2058	!jyg<
2059	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2060	!>jyg
2061	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2062	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2063	zustar, zu10m, zv10m, fder, &
2064	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2065	z0m, z0h, agesno, fsollw, fsolsw, &
2066	d_ts, fevap, fluxlat, t2m, &
2067	wfbils, wfbilo, fluxt, fluxu, fluxv, &
2068	dsens, devap, zxsnow, &
2069	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2070	!nrlmd+jyg<
2071	wake_delta_pbl_TKE &
2072	!>nrlmd+jyg
2073	)
2074	!
2075	!=================================================================
2076	! PROVISOIRE : DECOUPLAGE PBL/WAKE
2077	! --------------------------------
2078	!
2079	!! wake_deltat(:,:)=wake_deltat_sav(:,:)
2080	!! wake_deltaq(:,:)=wake_deltaq_sav(:,:)
2081	!=================================================================
2082	!
2083	! Add turbulent diffusion tendency to the wake difference variables
2084	IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2085	wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
2086	wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
2087	ENDIF
2088
2089
2090	!---------------------------------------------------------------------
2091	! ajout des tendances de la diffusion turbulente
2092	IF (klon_glo==1) THEN
2093	CALL add_pbl_tend &
2094	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2095	'vdf',abortphy)
2096	ELSE
2097	CALL add_phys_tend &
2098	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2099	'vdf',abortphy)
2100	ENDIF
2101	!--------------------------------------------------------------------
2102
2103	if (mydebug) then
2104	call writefield_phy('u_seri',u_seri,nbp_lev)
2105	call writefield_phy('v_seri',v_seri,nbp_lev)
2106	call writefield_phy('t_seri',t_seri,nbp_lev)
2107	call writefield_phy('q_seri',q_seri,nbp_lev)
2108	endif
2109
2110
2111	!albedo SB >>>
2112	albsol1=0.
2113	albsol2=0.
2114	falb1=0.
2115	falb2=0.
2116	select case(nsw)
2117	case(2)
2118	albsol1=albsol_dir(:,1)
2119	albsol2=albsol_dir(:,2)
2120	falb1=falb_dir(:,1,:)
2121	falb2=falb_dir(:,2,:)
2122	case(4)
2123	albsol1=albsol_dir(:,1)
2124	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2125	+albsol_dir(:,4)*SFRWL(4)
2126	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2127	falb1=falb_dir(:,1,:)
2128	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2129	+falb_dir(:,4,:)*SFRWL(4)
2130	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2131	case(6)
2132	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2133	+albsol_dir(:,3)*SFRWL(3)
2134	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2135	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2136	+albsol_dir(:,6)*SFRWL(6)
2137	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2138	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2139	+falb_dir(:,3,:)*SFRWL(3)
2140	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2141	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2142	+falb_dir(:,6,:)*SFRWL(6)
2143	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2144	end select
2145	!albedo SB <<<
2146
2147
2148	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2149	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2150
2151
2152	IF (ip_ebil_phy.ge.2) THEN
2153	ztit='after surface_main'
2154	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2155	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2156	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2157	call diagphy(cell_area,ztit,ip_ebil_phy &
2158	, zero_v, zero_v, zero_v, zero_v, sens &
2159	, evap , zero_v, zero_v, ztsol &
2160	, d_h_vcol, d_qt, d_ec &
2161	, fs_bound, fq_bound )
2162	END IF
2163
2164	ENDIF
2165	! =================================================================== c
2166	! Calcul de Qsat
2167
2168	DO k = 1, klev
2169	DO i = 1, klon
2170	zx_t = t_seri(i,k)
2171	IF (thermcep) THEN
2172	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2173	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2174	zx_qs = MIN(0.5,zx_qs)
2175	zcor = 1./(1.-retv*zx_qs)
2176	zx_qs = zx_qs*zcor
2177	ELSE
2178	!! IF (zx_t.LT.t_coup) THEN !jyg
2179	IF (zx_t.LT.rtt) THEN !jyg
2180	zx_qs = qsats(zx_t)/pplay(i,k)
2181	ELSE
2182	zx_qs = qsatl(zx_t)/pplay(i,k)
2183	ENDIF
2184	ENDIF
2185	zqsat(i,k)=zx_qs
2186	ENDDO
2187	ENDDO
2188
2189	if (prt_level.ge.1) then
2190	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2191	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2192	endif
2193	!
2194	! Appeler la convection (au choix)
2195	!
2196	DO k = 1, klev
2197	DO i = 1, klon
2198	conv_q(i,k) = d_q_dyn(i,k) &
2199	+ d_q_vdf(i,k)/dtime
2200	conv_t(i,k) = d_t_dyn(i,k) &
2201	+ d_t_vdf(i,k)/dtime
2202	ENDDO
2203	ENDDO
2204	IF (check) THEN
2205	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2206	WRITE(lunout,*) "avantcon=", za
2207	ENDIF
2208	zx_ajustq = .FALSE.
2209	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2210	IF (zx_ajustq) THEN
2211	DO i = 1, klon
2212	z_avant(i) = 0.0
2213	ENDDO
2214	DO k = 1, klev
2215	DO i = 1, klon
2216	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2217	*(paprs(i,k)-paprs(i,k+1))/RG
2218	ENDDO
2219	ENDDO
2220	ENDIF
2221
2222	! Calcule de vitesse verticale a partir de flux de masse verticale
2223	DO k = 1, klev
2224	DO i = 1, klon
2225	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2226	END DO
2227	END DO
2228	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2229	omega(igout, :)
2230
2231	IF (iflag_con.EQ.1) THEN
2232	abort_message ='reactiver le call conlmd dans physiq.F'
2233	CALL abort_physic (modname,abort_message,1)
2234	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2235	! . d_t_con, d_q_con,
2236	! . rain_con, snow_con, ibas_con, itop_con)
2237	ELSE IF (iflag_con.EQ.2) THEN
2238	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2239	conv_t, conv_q, -evap, omega, &
2240	d_t_con, d_q_con, rain_con, snow_con, &
2241	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2242	kcbot, kctop, kdtop, pmflxr, pmflxs)
2243	d_u_con = 0.
2244	d_v_con = 0.
2245
2246	WHERE (rain_con < 0.) rain_con = 0.
2247	WHERE (snow_con < 0.) snow_con = 0.
2248	DO i = 1, klon
2249	ibas_con(i) = klev+1 - kcbot(i)
2250	itop_con(i) = klev+1 - kctop(i)
2251	ENDDO
2252	ELSE IF (iflag_con.GE.3) THEN
2253	! nb of tracers for the KE convection:
2254	! MAF la partie traceurs est faite dans phytrac
2255	! on met ntra=1 pour limiter les appels mais on peut
2256	! supprimer les calculs / ftra.
2257	ntra = 1
2258
2259	!=======================================================================
2260	!ajout pour la parametrisation des poches froides: calcul de
2261	!t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2262	do k=1,klev
2263	do i=1,klon
2264	if (iflag_wake>=1) then
2265	t_wake(i,k) = t_seri(i,k) &
2266	+(1-wake_s(i))*wake_deltat(i,k)
2267	q_wake(i,k) = q_seri(i,k) &
2268	+(1-wake_s(i))*wake_deltaq(i,k)
2269	t_undi(i,k) = t_seri(i,k) &
2270	-wake_s(i)*wake_deltat(i,k)
2271	q_undi(i,k) = q_seri(i,k) &
2272	-wake_s(i)*wake_deltaq(i,k)
2273	else
2274	t_wake(i,k) = t_seri(i,k)
2275	q_wake(i,k) = q_seri(i,k)
2276	t_undi(i,k) = t_seri(i,k)
2277	q_undi(i,k) = q_seri(i,k)
2278	endif
2279	enddo
2280	enddo
2281	!
2282	!jyg<
2283	! Perform dry adiabatic adjustment on wake profile
2284	! The corresponding tendencies are added to the convective tendencies
2285	! after the call to the convective scheme.
2286	IF (iflag_wake>=1) then
2287	IF (ok_adjwk) THEN
2288	limbas(:) = 1
2289	CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
2290	d_t_adjwk, d_q_adjwk)
2291	ENDIF
2292	!
2293	DO k=1,klev
2294	DO i=1,klon
2295	IF (wake_s(i) .GT. 1.e-3) THEN
2296	t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
2297	q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
2298	wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
2299	wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
2300	ENDIF
2301	ENDDO
2302	ENDDO
2303	ENDIF ! (iflag_wake>=1)
2304	!>jyg
2305	!
2306
2307	! Calcul de l'energie disponible ALE (J/kg) et de la puissance
2308	! disponible ALP (W/m2) pour le soulevement des particules dans
2309	! le modele convectif
2310	!
2311	do i = 1,klon
2312	ALE(i) = 0.
2313	ALP(i) = 0.
2314	enddo
2315	!
2316	!calcul de ale_wake et alp_wake
2317	if (iflag_wake>=1) then
2318	if (itap .le. it_wape_prescr) then
2319	do i = 1,klon
2320	ale_wake(i) = wape_prescr
2321	alp_wake(i) = fip_prescr
2322	enddo
2323	else
2324	do i = 1,klon
2325	!jyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2326	!cc ale_wake(i) = 0.5wake_cstar(i)*2
2327	ale_wake(i) = wake_pe(i)
2328	alp_wake(i) = wake_fip(i)
2329	enddo
2330	endif
2331	else
2332	do i = 1,klon
2333	ale_wake(i) = 0.
2334	alp_wake(i) = 0.
2335	enddo
2336	endif
2337	!combinaison avec ale et alp de couche limite: constantes si pas
2338	!de couplage, valeurs calculees dans le thermique sinon
2339	if (iflag_coupl.eq.0) then
2340	if (debut.and.prt_level.gt.9) &
2341	WRITE(lunout,*)'ALE et ALP imposes'
2342	do i = 1,klon
2343	!on ne couple que ale
2344	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2345	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2346	!on ne couple que alp
2347	! ALP(i) = alp_wake(i) + Alp_bl(i)
2348	ALP(i) = alp_wake(i) + alp_bl_prescr
2349	enddo
2350	else
2351	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2352	! do i = 1,klon
2353	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2354	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2355	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2356	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2357	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2358	! enddo
2359
2360	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2361	! Modif FH 2010/04/27. Sans doute temporaire.
2362	! Deux options pour le alp_offset : constant si >?? 0 ou
2363	! proportionnel ??a w si <0
2364	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2365	! Estimation d'une vitesse verticale effective pour ALP
2366	if (1==0) THEN
2367	www(1:klon)=0.
2368	do k=2,klev-1
2369	do i=1,klon
2370	www(i)=max(www(i),-omega(i,k)RDt_seri(i,k) &
2371	/(RGpaprs(i,k)) zw2(i,k)*zw2(i,k))
2372	! if (paprs(i,k)>pbase(i)) then
2373	! calcul approche de la vitesse verticale en m/s
2374	! www(i)=max(www(i),-omega(i,k)RDtemp(i,k)/(RG*paprs(i,k))
2375	! endif
2376	! Le 0.1 est en gros H / ps = 1e5 / 1e4
2377	enddo
2378	enddo
2379	do i=1,klon
2380	if (www(i)>0. .and. ale_bl(i)>0. ) www(i)=www(i)/ale_bl(i)
2381	enddo
2382	ENDIF
2383
2384
2385	do i = 1,klon
2386	ALE(i) = max(ale_wake(i),Ale_bl(i))
2387	!cc nrlmd le 10/04/2012----------Stochastic triggering------------
2388	if (iflag_trig_bl.ge.1) then
2389	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2390	endif
2391	!cc fin nrlmd le 10/04/2012
2392	if (alp_offset>=0.) then
2393	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2394	else
2395	abort_message ='Ne pas passer la car www non calcule'
2396	CALL abort_physic (modname,abort_message,1)
2397
2398	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2399	! _ _
2400	! Ajout d'une composante 3 * A * w w'2 a w'3 avec
2401	! w=www : w max sous pbase ou A est la fraction
2402	! couverte par les ascendances w' on utilise le fait
2403	! que A * w'3 = ALP et donc A * w'2 ~ ALP / sqrt(ALE)
2404	! (on ajoute 0.1 pour les singularites)
2405	ALP(i)=alp_wake(i)(1.+3.www(i)/( sqrt(ale_wake(i))+0.1) ) &
2406	+alp_bl(i) (1.+3.www(i)/( sqrt(ale_bl(i)) +0.1) )
2407	! ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2408	! if (alp(i)<0.) then
2409	! print*,'ALP ',alp(i),alp_wake(i) &
2410	! ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2411	! endif
2412	endif
2413	enddo
2414	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2415
2416	endif
2417	do i=1,klon
2418	if (alp(i)>alp_max) then
2419	IF(prt_level>9)WRITE(lunout,*) &
2420	'WARNING SUPER ALP (seuil=',alp_max, &
2421	'): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2422	alp(i)=alp_max
2423	endif
2424	if (ale(i)>ale_max) then
2425	IF(prt_level>9)WRITE(lunout,*) &
2426	'WARNING SUPER ALE (seuil=',ale_max, &
2427	'): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2428	ale(i)=ale_max
2429	endif
2430	enddo
2431
2432	!fin calcul ale et alp
2433	!=======================================================================
2434
2435
2436	! sb, oct02:
2437	! Schema de convection modularise et vectorise:
2438	! (driver commun aux versions 3 et 4)
2439	!
2440	IF (ok_cvl) THEN ! new driver for convectL
2441	!
2442	!jyg<
2443	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2444	! Calculate the upmost level of deep convection loops: k_upper_cv
2445	! (near 22 km)
2446	izero = klon/2+1/klon
2447	k_upper_cv = klev
2448	DO k = klev,1,-1
2449	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2450	ENDDO
2451	IF (prt_level .ge. 5) THEN
2452	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2453	k_upper_cv
2454	ENDIF
2455	!
2456	!>jyg
2457	IF (type_trac == 'repr') THEN
2458	nbtr_tmp=ntra
2459	ELSE
2460	nbtr_tmp=nbtr
2461	END IF
2462	!jyg iflag_con est dans clesphys
2463	!c CALL concvl (iflag_con,iflag_clos,
2464	CALL concvl (iflag_clos, &
2465	dtime, paprs, pplay, k_upper_cv, t_undi,q_undi, &
2466	t_wake,q_wake,wake_s, &
2467	u_seri,v_seri,tr_seri,nbtr_tmp, &
2468	ALE,ALP, &
2469	sig1,w01, &
2470	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2471	rain_con, snow_con, ibas_con, itop_con, sigd, &
2472	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2473	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2474	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2475	! RomP >>>
2476	!! . pmflxr,pmflxs,da,phi,mp,
2477	!! . ftd,fqd,lalim_conv,wght_th)
2478	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2479	ftd,fqd,lalim_conv,wght_th, &
2480	ev, ep,epmlmMm,eplaMm, &
2481	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2482	tau_cld_cv,coefw_cld_cv)
2483	! RomP <<<
2484
2485	!IM begin
2486	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2487	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2488	!IM end
2489	!IM cf. FH
2490	clwcon0=qcondc
2491	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2492
2493	do i = 1, klon
2494	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2495	enddo
2496	!
2497	!jyg<
2498	! Add the tendency due to the dry adjustment of the wake profile
2499	IF (iflag_wake>=1) THEN
2500	DO k=1,klev
2501	DO i=1,klon
2502	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
2503	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
2504	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2505	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2506	ENDDO
2507	ENDDO
2508	ENDIF
2509	!>jyg
2510	!
2511	ELSE ! ok_cvl
2512
2513	! MAF conema3 ne contient pas les traceurs
2514	CALL conema3 (dtime, &
2515	paprs,pplay,t_seri,q_seri, &
2516	u_seri,v_seri,tr_seri,ntra, &
2517	sig1,w01, &
2518	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2519	rain_con, snow_con, ibas_con, itop_con, &
2520	upwd,dnwd,dnwd0,bas,top, &
2521	Ma,cape,tvp,rflag, &
2522	pbase &
2523	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2524	,clwcon0)
2525
2526	ENDIF ! ok_cvl
2527
2528	!
2529	! Correction precip
2530	rain_con = rain_con * cvl_corr
2531	snow_con = snow_con * cvl_corr
2532	!
2533
2534	IF (.NOT. ok_gust) THEN
2535	do i = 1, klon
2536	wd(i)=0.0
2537	enddo
2538	ENDIF
2539
2540	! =================================================================== c
2541	! Calcul des proprietes des nuages convectifs
2542	!
2543
2544	! calcul des proprietes des nuages convectifs
2545	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2546	IF (iflag_cld_cv == 0) THEN
2547	call clouds_gno &
2548	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2549	ELSE
2550	call clouds_bigauss &
2551	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2552	ENDIF
2553
2554
2555	! =================================================================== c
2556
2557	DO i = 1, klon
2558	itop_con(i) = min(max(itop_con(i),1),klev)
2559	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2560	ENDDO
2561
2562	DO i = 1, klon
2563	ema_pcb(i) = paprs(i,ibas_con(i))
2564	ENDDO
2565	DO i = 1, klon
2566	! L'idicage de itop_con peut cacher un pb potentiel
2567	! FH sous la dictee de JYG, CR
2568	ema_pct(i) = paprs(i,itop_con(i)+1)
2569
2570	if (itop_con(i).gt.klev-3) then
2571	if(prt_level >= 9) then
2572	write(lunout,*)'La convection monte trop haut '
2573	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2574	endif
2575	endif
2576	ENDDO
2577	ELSE IF (iflag_con.eq.0) THEN
2578	write(lunout,*) 'On n appelle pas la convection'
2579	clwcon0=0.
2580	rnebcon0=0.
2581	d_t_con=0.
2582	d_q_con=0.
2583	d_u_con=0.
2584	d_v_con=0.
2585	rain_con=0.
2586	snow_con=0.
2587	bas=1
2588	top=1
2589	ELSE
2590	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2591	call abort_physic("physiq", "", 1)
2592	ENDIF
2593
2594	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2595	! . d_u_con, d_v_con)
2596
2597	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2598	'convection',abortphy)
2599
2600	!-------------------------------------------------------------------------
2601
2602	if (mydebug) then
2603	call writefield_phy('u_seri',u_seri,nbp_lev)
2604	call writefield_phy('v_seri',v_seri,nbp_lev)
2605	call writefield_phy('t_seri',t_seri,nbp_lev)
2606	call writefield_phy('q_seri',q_seri,nbp_lev)
2607	endif
2608
2609	!IM
2610	IF (ip_ebil_phy.ge.2) THEN
2611	ztit='after convect'
2612	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2613	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2614	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2615	call diagphy(cell_area,ztit,ip_ebil_phy &
2616	, zero_v, zero_v, zero_v, zero_v, zero_v &
2617	, zero_v, rain_con, snow_con, ztsol &
2618	, d_h_vcol, d_qt, d_ec &
2619	, fs_bound, fq_bound )
2620	END IF
2621	!
2622	IF (check) THEN
2623	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2624	WRITE(lunout,*)"aprescon=", za
2625	zx_t = 0.0
2626	za = 0.0
2627	DO i = 1, klon
2628	za = za + cell_area(i)/REAL(klon)
2629	zx_t = zx_t + (rain_con(i)+ &
2630	snow_con(i))*cell_area(i)/REAL(klon)
2631	ENDDO
2632	zx_t = zx_t/za*dtime
2633	WRITE(lunout,*)"Precip=", zx_t
2634	ENDIF
2635	IF (zx_ajustq) THEN
2636	DO i = 1, klon
2637	z_apres(i) = 0.0
2638	ENDDO
2639	DO k = 1, klev
2640	DO i = 1, klon
2641	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2642	*(paprs(i,k)-paprs(i,k+1))/RG
2643	ENDDO
2644	ENDDO
2645	DO i = 1, klon
2646	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2647	/z_apres(i)
2648	ENDDO
2649	DO k = 1, klev
2650	DO i = 1, klon
2651	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2652	z_factor(i).LT.(1.0-1.0E-08)) THEN
2653	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2654	ENDIF
2655	ENDDO
2656	ENDDO
2657	ENDIF
2658	zx_ajustq=.FALSE.
2659
2660	!
2661	!==========================================================================
2662	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2663	!pour la couche limite diffuse pour l instant
2664	!
2665	!
2666	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
2667	! il faut rajouter cette tendance calcul\'ee hors des poches
2668	! froides
2669	!
2670	if (iflag_wake>=1) then
2671	DO k=1,klev
2672	DO i=1,klon
2673	dt_dwn(i,k) = ftd(i,k)
2674	dq_dwn(i,k) = fqd(i,k)
2675	M_dwn(i,k) = dnwd0(i,k)
2676	M_up(i,k) = upwd(i,k)
2677	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2678	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2679	ENDDO
2680	ENDDO
2681	!nrlmd+jyg<
2682	DO k=1,klev
2683	DO i=1,klon
2684	wdt_PBL(i,k) = 0.
2685	wdq_PBL(i,k) = 0.
2686	udt_PBL(i,k) = 0.
2687	udq_PBL(i,k) = 0.
2688	ENDDO
2689	ENDDO
2690	!
2691	IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
2692	DO k=1,klev
2693	DO i=1,klon
2694	wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
2695	wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
2696	udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
2697	udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
2698	!! dt_dwn(i,k) = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
2699	!! dq_dwn(i,k) = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
2700	!! dt_a (i,k) = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
2701	!! dq_a (i,k) = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
2702	ENDDO
2703	ENDDO
2704	ENDIF
2705	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2706	DO k=1,klev
2707	DO i=1,klon
2708	!! dt_dwn(i,k) = dt_dwn(i,k) + 0.
2709	!! dq_dwn(i,k) = dq_dwn(i,k) + 0.
2710	!! dt_a(i,k) = dt_a(i,k) + d_t_ajs(i,k)/dtime
2711	!! dq_a(i,k) = dq_a(i,k) + d_q_ajs(i,k)/dtime
2712	udt_PBL(i,k) = udt_PBL(i,k) + d_t_ajs(i,k)/dtime
2713	udq_PBL(i,k) = udq_PBL(i,k) + d_q_ajs(i,k)/dtime
2714	ENDDO
2715	ENDDO
2716	ENDIF
2717	!>nrlmd+jyg
2718
2719	IF (iflag_wake==2) THEN
2720	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2721	DO k = 1,klev
2722	dt_dwn(:,k)= dt_dwn(:,k)+ &
2723	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2724	dq_dwn(:,k)= dq_dwn(:,k)+ &
2725	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2726	ENDDO
2727	ELSEIF (iflag_wake==3) THEN
2728	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2729	DO k = 1,klev
2730	DO i=1,klon
2731	IF (rneb(i,k)==0.) THEN
2732	! On ne tient compte des tendances qu'en dehors des
2733	! nuages (c'est-\`a-dire a priri dans une region ou
2734	! l'eau se reevapore).
2735	dt_dwn(i,k)= dt_dwn(i,k)+ &
2736	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2737	dq_dwn(i,k)= dq_dwn(i,k)+ &
2738	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2739	ENDIF
2740	ENDDO
2741	ENDDO
2742	ENDIF
2743
2744	!
2745	!calcul caracteristiques de la poche froide
2746	call calWAKE (paprs,pplay,dtime &
2747	,t_seri,q_seri,omega &
2748	,dt_dwn,dq_dwn,M_dwn,M_up &
2749	,dt_a,dq_a,sigd &
2750	,wdt_PBL,wdq_PBL &
2751	,udt_PBL,udq_PBL &
2752	,wake_deltat,wake_deltaq,wake_dth &
2753	,wake_h,wake_s,wake_dens &
2754	,wake_pe,wake_fip,wake_gfl &
2755	,dt_wake,dq_wake &
2756	,wake_k, t_undi,q_undi &
2757	,wake_omgbdth,wake_dp_omgb &
2758	,wake_dtKE,wake_dqKE &
2759	,wake_dtPBL,wake_dqPBL &
2760	,wake_omg,wake_dp_deltomg &
2761	,wake_spread,wake_Cstar,wake_d_deltat_gw &
2762	,wake_ddeltat,wake_ddeltaq)
2763	!
2764	!-----------------------------------------------------------------------
2765	! ajout des tendances des poches froides
2766	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2767	! coherent avec les autres d_t_...
2768	d_t_wake(:,:)=dt_wake(:,:)*dtime
2769	d_q_wake(:,:)=dq_wake(:,:)*dtime
2770	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
2771	abortphy)
2772	!------------------------------------------------------------------------
2773
2774	endif ! (iflag_wake>=1)
2775	!
2776	!===================================================================
2777	!JYG
2778	IF (ip_ebil_phy.ge.2) THEN
2779	ztit='after wake'
2780	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2781	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2782	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2783	call diagphy(cell_area,ztit,ip_ebil_phy &
2784	, zero_v, zero_v, zero_v, zero_v, zero_v &
2785	, zero_v, zero_v, zero_v, ztsol &
2786	, d_h_vcol, d_qt, d_ec &
2787	, fs_bound, fq_bound )
2788	END IF
2789
2790	! print*,'apres callwake iflag_cld_th=', iflag_cld_th
2791	!
2792	!===================================================================
2793	! Convection seche (thermiques ou ajustement)
2794	!===================================================================
2795	!
2796	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2797	,seuil_inversion,weak_inversion,dthmin)
2798
2799
2800
2801	d_t_ajsb(:,:)=0.
2802	d_q_ajsb(:,:)=0.
2803	d_t_ajs(:,:)=0.
2804	d_u_ajs(:,:)=0.
2805	d_v_ajs(:,:)=0.
2806	d_q_ajs(:,:)=0.
2807	clwcon0th(:,:)=0.
2808	!
2809	! fm_therm(:,:)=0.
2810	! entr_therm(:,:)=0.
2811	! detr_therm(:,:)=0.
2812	!
2813	IF(prt_level>9)WRITE(lunout,*) &
2814	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2815	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2816	if(iflag_thermals<0) then
2817	! Rien
2818	! ====
2819	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2820
2821
2822	else
2823
2824	! Thermiques
2825	! ==========
2826	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2827	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2828
2829
2830	!cc nrlmd le 10/04/2012
2831	DO k=1,klev+1
2832	DO i=1,klon
2833	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2834	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2835	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2836	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2837	ENDDO
2838	ENDDO
2839	!cc fin nrlmd le 10/04/2012
2840
2841	if (iflag_thermals>=1) then
2842	!jyg<
2843	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2844	! Appel des thermiques avec les profils exterieurs aux poches
2845	DO k=1,klev
2846	DO i=1,klon
2847	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2848	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2849	ENDDO
2850	ENDDO
2851	ELSE
2852	! Appel des thermiques avec les profils moyens
2853	DO k=1,klev
2854	DO i=1,klon
2855	t_therm(i,k) = t_seri(i,k)
2856	q_therm(i,k) = q_seri(i,k)
2857	ENDDO
2858	ENDDO
2859	ENDIF
2860	!>jyg
2861	call calltherm(pdtphys &
2862	,pplay,paprs,pphi,weak_inversion &
2863	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
2864	!jyg
2865	,u_seri,v_seri,t_therm,q_therm,zqsat,debut & !jyg
2866	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2867	,fm_therm,entr_therm,detr_therm &
2868	,zqasc,clwcon0th,lmax_th,ratqscth &
2869	,ratqsdiff,zqsatth &
2870	!on rajoute ale et alp, et les
2871	!caracteristiques de la couche alim
2872	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2873	,ztv,zpspsk,ztla,zthl &
2874	!cc nrlmd le 10/04/2012
2875	,pbl_tke_input,pctsrf,omega,cell_area &
2876	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2877	,n2,s2,ale_bl_stat &
2878	,therm_tke_max,env_tke_max &
2879	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2880	,alp_bl_conv,alp_bl_stat &
2881	!cc fin nrlmd le 10/04/2012
2882	,zqla,ztva )
2883	!
2884	!jyg<
2885	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2886	! Si les thermiques ne sont presents que hors des
2887	! poches, la tendance moyenne associ\'ee doit etre
2888	! multipliee par la fraction surfacique qu'ils couvrent.
2889	DO k=1,klev
2890	DO i=1,klon
2891	!
2892	wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
2893	wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
2894	t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
2895	q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
2896	!
2897	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2898	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2899	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2900	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2901	!
2902	ENDDO
2903	ENDDO
2904	ELSE
2905	DO k=1,klev
2906	DO i=1,klon
2907	t_seri(i,k) = t_therm(i,k)
2908	q_seri(i,k) = q_therm(i,k)
2909	ENDDO
2910	ENDDO
2911	ENDIF
2912	!>jyg
2913
2914	!cc nrlmd le 10/04/2012
2915	!-----------Stochastic triggering-----------
2916	if (iflag_trig_bl.ge.1) then
2917	!
2918	IF (prt_level .GE. 10) THEN
2919	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2920	cin, ale_bl_stat, alp_bl_stat
2921	ENDIF
2922
2923
2924	!----Initialisations
2925	do i=1,klon
2926	proba_notrig(i)=1.
2927	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2928	if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
2929	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2930	tau_trig(i)=tau_trig_shallow
2931	else
2932	tau_trig(i)=tau_trig_deep
2933	endif
2934	enddo
2935	!
2936	IF (prt_level .GE. 10) THEN
2937	print *,'random_notrig, tau_trig ', &
2938	random_notrig, tau_trig
2939	print *,'s_trig,s2,n2 ', &
2940	s_trig,s2,n2
2941	ENDIF
2942
2943	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2944	IF (iflag_trig_bl.eq.1) then
2945
2946	!----Tirage al\'eatoire et calcul de ale_bl_trig
2947	do i=1,klon
2948	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2949	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2950	(n2(i)*dtime/tau_trig(i))
2951	! print *, 'proba_notrig(i) ',proba_notrig(i)
2952	if (random_notrig(i) .ge. proba_notrig(i)) then
2953	ale_bl_trig(i)=ale_bl_stat(i)
2954	else
2955	ale_bl_trig(i)=0.
2956	endif
2957	else
2958	proba_notrig(i)=1.
2959	random_notrig(i)=0.
2960	ale_bl_trig(i)=0.
2961	endif
2962	enddo
2963
2964	ELSE IF (iflag_trig_bl.ge.2) then
2965
2966	do i=1,klon
2967	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
2968	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2969	(n2(i)*dtime/tau_trig(i))
2970	! print *, 'proba_notrig(i) ',proba_notrig(i)
2971	if (random_notrig(i) .ge. proba_notrig(i)) then
2972	ale_bl_trig(i)=Ale_bl(i)
2973	else
2974	ale_bl_trig(i)=0.
2975	endif
2976	else
2977	proba_notrig(i)=1.
2978	random_notrig(i)=0.
2979	ale_bl_trig(i)=0.
2980	endif
2981	enddo
2982
2983	ENDIF
2984
2985	!
2986	IF (prt_level .GE. 10) THEN
2987	print *,'proba_notrig, ale_bl_trig ', &
2988	proba_notrig, ale_bl_trig
2989	ENDIF
2990
2991	endif !(iflag_trig_bl)
2992
2993	!-----------Statistical closure-----------
2994	if (iflag_clos_bl.eq.1) then
2995
2996	do i=1,klon
2997	!CR: alp probabiliste
2998	if (ale_bl_trig(i).gt.0.) then
2999	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
3000	endif
3001	enddo
3002
3003	else if (iflag_clos_bl.eq.2) then
3004
3005	!CR: alp calculee dans thermcell_main
3006	do i=1,klon
3007	alp_bl(i)=alp_bl_stat(i)
3008	enddo
3009
3010	else
3011
3012	alp_bl_stat(:)=0.
3013
3014	endif !(iflag_clos_bl)
3015
3016	IF (prt_level .GE. 10) THEN
3017	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
3018	ENDIF
3019
3020	!cc fin nrlmd le 10/04/2012
3021
3022	! ------------------------------------------------------------------
3023	! Transport de la TKE par les panaches thermiques.
3024	! FH : 2010/02/01
3025	! if (iflag_pbl.eq.10) then
3026	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3027	! s rg,paprs,pbl_tke)
3028	! endif
3029	! -------------------------------------------------------------------
3030	!IM/FH: 2011/02/23
3031	! Couplage Thermiques/Emanuel seulement si T<0
3032	if (iflag_coupl==2) then
3033	IF (prt_level .GE. 10) THEN
3034	print*,'Couplage Thermiques/Emanuel seulement si T<0'
3035	ENDIF
3036	do i=1,klon
3037	if (t_seri(i,lmax_th(i))>273.) then
3038	Ale_bl(i)=0.
3039	endif
3040	enddo
3041	endif
3042
3043	do i=1,klon
3044	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3045	!CR:04/05/12:correction calcul zmax
3046	zmax_th(i)=zmax0(i)
3047	enddo
3048
3049	endif
3050
3051
3052	! Ajustement sec
3053	! ==============
3054
3055	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3056	! a partir du sommet des thermiques.
3057	! Dans le cas contraire, on demarre au niveau 1.
3058
3059	if (iflag_thermals>=13.or.iflag_thermals<=0) then
3060
3061	if(iflag_thermals.eq.0) then
3062	IF(prt_level>9)WRITE(lunout,*)'ajsec'
3063	limbas(:)=1
3064	else
3065	limbas(:)=lmax_th(:)
3066	endif
3067
3068	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3069	! pour des test de convergence numerique.
3070	! Le nouveau ajsec est a priori mieux, meme pour le cas
3071	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3072	! non nulles numeriquement pour des mailles non concernees.
3073
3074	if (iflag_thermals==0) then
3075	! Calling adjustment alone (but not the thermal plume model)
3076	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3077	, d_t_ajsb, d_q_ajsb)
3078	else if (iflag_thermals>0) then
3079	! Calling adjustment above the top of thermal plumes
3080	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3081	, d_t_ajsb, d_q_ajsb)
3082	endif
3083
3084	!--------------------------------------------------------------------
3085	! ajout des tendances de l'ajustement sec ou des thermiques
3086	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3087	'ajsb',abortphy)
3088	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3089	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3090
3091	!---------------------------------------------------------------------
3092
3093	endif
3094
3095	endif
3096	!
3097	!===================================================================
3098	!IM
3099	IF (ip_ebil_phy.ge.2) THEN
3100	ztit='after dry_adjust'
3101	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3102	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3103	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3104	call diagphy(cell_area,ztit,ip_ebil_phy &
3105	, zero_v, zero_v, zero_v, zero_v, zero_v &
3106	, zero_v, zero_v, zero_v, ztsol &
3107	, d_h_vcol, d_qt, d_ec &
3108	, fs_bound, fq_bound )
3109	END IF
3110
3111
3112	!-------------------------------------------------------------------------
3113	! Computation of ratqs, the width (normalized) of the subrid scale
3114	! water distribution
3115	CALL calcratqs(klon,klev,prt_level,lunout, &
3116	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3117	ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
3118	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3119	paprs,pplay,q_seri,zqsat,fm_therm, &
3120	ratqs,ratqsc)
3121
3122
3123	!
3124	! Appeler le processus de condensation a grande echelle
3125	! et le processus de precipitation
3126	!-------------------------------------------------------------------------
3127	IF (prt_level .GE.10) THEN
3128	print *,'itap, ->fisrtilp ',itap
3129	ENDIF
3130	!
3131	CALL fisrtilp(dtime,paprs,pplay, &
3132	t_seri, q_seri,ptconv,ratqs, &
3133	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3134	rain_lsc, snow_lsc, &
3135	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3136	frac_impa, frac_nucl, beta_prec_fisrt, &
3137	prfl, psfl, rhcl, &
3138	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3139	iflag_ice_thermo)
3140	!
3141	WHERE (rain_lsc < 0) rain_lsc = 0.
3142	WHERE (snow_lsc < 0) snow_lsc = 0.
3143
3144	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3145	'lsc',abortphy)
3146	!---------------------------------------------------------------------------
3147	DO k = 1, klev
3148	DO i = 1, klon
3149	cldfra(i,k) = rneb(i,k)
3150	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3151	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3152	ENDDO
3153	ENDDO
3154	IF (check) THEN
3155	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3156	WRITE(lunout,*)"apresilp=", za
3157	zx_t = 0.0
3158	za = 0.0
3159	DO i = 1, klon
3160	za = za + cell_area(i)/REAL(klon)
3161	zx_t = zx_t + (rain_lsc(i) &
3162	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3163	ENDDO
3164	zx_t = zx_t/za*dtime
3165	WRITE(lunout,*)"Precip=", zx_t
3166	ENDIF
3167	!IM
3168	IF (ip_ebil_phy.ge.2) THEN
3169	ztit='after fisrt'
3170	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3171	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3172	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3173	call diagphy(cell_area,ztit,ip_ebil_phy &
3174	, zero_v, zero_v, zero_v, zero_v, zero_v &
3175	, zero_v, rain_lsc, snow_lsc, ztsol &
3176	, d_h_vcol, d_qt, d_ec &
3177	, fs_bound, fq_bound )
3178	END IF
3179
3180	if (mydebug) then
3181	call writefield_phy('u_seri',u_seri,nbp_lev)
3182	call writefield_phy('v_seri',v_seri,nbp_lev)
3183	call writefield_phy('t_seri',t_seri,nbp_lev)
3184	call writefield_phy('q_seri',q_seri,nbp_lev)
3185	endif
3186
3187	!
3188	!-------------------------------------------------------------------
3189	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3190	!-------------------------------------------------------------------
3191
3192	! 1. NUAGES CONVECTIFS
3193	!
3194	!IM cf FH
3195	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3196	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3197	snow_tiedtke=0.
3198	! print*,'avant calcul de la pseudo precip '
3199	! print*,'iflag_cld_th',iflag_cld_th
3200	if (iflag_cld_th.eq.-1) then
3201	rain_tiedtke=rain_con
3202	else
3203	! print*,'calcul de la pseudo precip '
3204	rain_tiedtke=0.
3205	! print*,'calcul de la pseudo precip 0'
3206	do k=1,klev
3207	do i=1,klon
3208	if (d_q_con(i,k).lt.0.) then
3209	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3210	*(paprs(i,k)-paprs(i,k+1))/rg
3211	endif
3212	enddo
3213	enddo
3214	endif
3215	!
3216	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3217	!
3218
3219	! Nuages diagnostiques pour Tiedtke
3220	CALL diagcld1(paprs,pplay, &
3221	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3222	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3223	diafra,dialiq)
3224	DO k = 1, klev
3225	DO i = 1, klon
3226	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3227	cldliq(i,k) = dialiq(i,k)
3228	cldfra(i,k) = diafra(i,k)
3229	ENDIF
3230	ENDDO
3231	ENDDO
3232
3233	ELSE IF (iflag_cld_th.ge.3) THEN
3234	! On prend pour les nuages convectifs le max du calcul de la
3235	! convection et du calcul du pas de temps precedent diminue d'un facteur
3236	! facttemps
3237	facteur = pdtphys *facttemps
3238	do k=1,klev
3239	do i=1,klon
3240	rnebcon(i,k)=rnebcon(i,k)*facteur
3241	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
3242	then
3243	rnebcon(i,k)=rnebcon0(i,k)
3244	clwcon(i,k)=clwcon0(i,k)
3245	endif
3246	enddo
3247	enddo
3248
3249	!
3250	!jq - introduce the aerosol direct and first indirect radiative forcings
3251	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3252	IF (flag_aerosol .gt. 0) THEN
3253	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3254	IF (.NOT. aerosol_couple) THEN
3255	!
3256	CALL readaerosol_optic( &
3257	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3258	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3259	mass_solu_aero, mass_solu_aero_pi, &
3260	tau_aero, piz_aero, cg_aero, &
3261	tausum_aero, tau3d_aero)
3262	ENDIF
3263	ELSE ! RRTM radiation
3264	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3265	abort_message='config_inca=aero et rrtm=1 impossible'
3266	call abort_physic(modname,abort_message,1)
3267	ELSE
3268	!
3269	#ifdef CPP_RRTM
3270	IF (NSW.EQ.6) THEN
3271	!--new aerosol properties
3272	!
3273	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
3274	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3275	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3276	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3277	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3278	tausum_aero, tau3d_aero)
3279
3280	ELSE IF (NSW.EQ.2) THEN
3281	!--for now we use the old aerosol properties
3282	!
3283	CALL readaerosol_optic( &
3284	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3285	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3286	mass_solu_aero, mass_solu_aero_pi, &
3287	tau_aero, piz_aero, cg_aero, &
3288	tausum_aero, tau3d_aero)
3289	!
3290	!--natural aerosols
3291	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3292	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3293	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3294	!--all aerosols
3295	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3296	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3297	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3298	ELSE
3299	abort_message='Only NSW=2 or 6 are possible with ' &
3300	// 'aerosols and iflag_rrtm=1'
3301	call abort_physic(modname,abort_message,1)
3302	ENDIF
3303
3304	CALL aeropt_lw_rrtm
3305	!
3306	#else
3307	abort_message='You should compile with -rrtm if running ' &
3308	// 'with iflag_rrtm=1'
3309	call abort_physic(modname,abort_message,1)
3310	#endif
3311	!
3312	ENDIF
3313	ENDIF
3314	ELSE
3315	tausum_aero(:,:,:) = 0.
3316	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3317	tau_aero(:,:,:,:) = 1.e-15
3318	piz_aero(:,:,:,:) = 1.
3319	cg_aero(:,:,:,:) = 0.
3320	ELSE
3321	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3322	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3323	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3324	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3325	ENDIF
3326	ENDIF
3327	!
3328	!--STRAT AEROSOL
3329	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3330	IF (flag_aerosol_strat) THEN
3331	IF (prt_level .GE.10) THEN
3332	PRINT *,'appel a readaerosolstrat', mth_cur
3333	ENDIF
3334	IF (iflag_rrtm.EQ.0) THEN
3335	CALL readaerosolstrato(debut)
3336	ELSE
3337	#ifdef CPP_RRTM
3338	CALL readaerosolstrato_rrtm(debut)
3339	#else
3340
3341	abort_message='You should compile with -rrtm if running ' &
3342	// 'with iflag_rrtm=1'
3343	call abort_physic(modname,abort_message,1)
3344	#endif
3345	ENDIF
3346	ENDIF
3347	!--fin STRAT AEROSOL
3348
3349
3350	! On prend la somme des fractions nuageuses et des contenus en eau
3351
3352	if (iflag_cld_th>=5) then
3353
3354	do k=1,klev
3355	ptconvth(:,k)=fm_therm(:,k+1)>0.
3356	enddo
3357
3358	if (iflag_coupl==4) then
3359
3360	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3361	! convectives et lsc dans la partie des thermiques
3362	! Le controle par iflag_coupl est peut etre provisoire.
3363	do k=1,klev
3364	do i=1,klon
3365	if (ptconv(i,k).and.ptconvth(i,k)) then
3366	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3367	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3368	else if (ptconv(i,k)) then
3369	cldfra(i,k)=rnebcon(i,k)
3370	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3371	endif
3372	enddo
3373	enddo
3374
3375	else if (iflag_coupl==5) then
3376	do k=1,klev
3377	do i=1,klon
3378	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3379	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3380	enddo
3381	enddo
3382
3383	else
3384
3385	! Si on est sur un point touche par la convection
3386	! profonde et pas par les thermiques, on prend la
3387	! couverture nuageuse et l'eau nuageuse de la convection
3388	! profonde.
3389
3390	!IM/FH: 2011/02/23
3391	! definition des points sur lesquels ls thermiques sont actifs
3392
3393	do k=1,klev
3394	do i=1,klon
3395	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3396	cldfra(i,k)=rnebcon(i,k)
3397	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3398	endif
3399	enddo
3400	enddo
3401
3402	endif
3403
3404	else
3405
3406	! Ancienne version
3407	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3408	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3409	endif
3410
3411	ENDIF
3412
3413	! plulsc(:)=0.
3414	! do k=1,klev,-1
3415	! do i=1,klon
3416	! zzz=prfl(:,k)+psfl(:,k)
3417	! if (.not.ptconvth.zzz.gt.0.)
3418	! enddo prfl, psfl,
3419	! enddo
3420	!
3421	! 2. NUAGES STARTIFORMES
3422	!
3423	IF (ok_stratus) THEN
3424	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3425	DO k = 1, klev
3426	DO i = 1, klon
3427	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3428	cldliq(i,k) = dialiq(i,k)
3429	cldfra(i,k) = diafra(i,k)
3430	ENDIF
3431	ENDDO
3432	ENDDO
3433	ENDIF
3434	!
3435	! Precipitation totale
3436	!
3437	DO i = 1, klon
3438	rain_fall(i) = rain_con(i) + rain_lsc(i)
3439	snow_fall(i) = snow_con(i) + snow_lsc(i)
3440	ENDDO
3441	!IM
3442	IF (ip_ebil_phy.ge.2) THEN
3443	ztit="after diagcld"
3444	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3445	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3446	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3447	call diagphy(cell_area,ztit,ip_ebil_phy &
3448	, zero_v, zero_v, zero_v, zero_v, zero_v &
3449	, zero_v, zero_v, zero_v, ztsol &
3450	, d_h_vcol, d_qt, d_ec &
3451	, fs_bound, fq_bound )
3452	END IF
3453	!
3454	! Calculer l'humidite relative pour diagnostique
3455	!
3456	DO k = 1, klev
3457	DO i = 1, klon
3458	zx_t = t_seri(i,k)
3459	IF (thermcep) THEN
3460	!! if (iflag_ice_thermo.eq.0) then !jyg
3461	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3462	!! else !jyg
3463	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3464	!! endif !jyg
3465	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3466	zx_qs = MIN(0.5,zx_qs)
3467	zcor = 1./(1.-retv*zx_qs)
3468	zx_qs = zx_qs*zcor
3469	ELSE
3470	!! IF (zx_t.LT.t_coup) THEN !jyg
3471	IF (zx_t.LT.rtt) THEN !jyg
3472	zx_qs = qsats(zx_t)/pplay(i,k)
3473	ELSE
3474	zx_qs = qsatl(zx_t)/pplay(i,k)
3475	ENDIF
3476	ENDIF
3477	zx_rh(i,k) = q_seri(i,k)/zx_qs
3478	zqsat(i,k)=zx_qs
3479	ENDDO
3480	ENDDO
3481
3482	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3483	! equivalente a 2m (tpote) pour diagnostique
3484	!
3485	DO i = 1, klon
3486	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3487	IF (thermcep) THEN
3488	IF(zt2m(i).LT.RTT) then
3489	Lheat=RLSTT
3490	ELSE
3491	Lheat=RLVTT
3492	ENDIF
3493	ELSE
3494	IF (zt2m(i).LT.RTT) THEN
3495	Lheat=RLSTT
3496	ELSE
3497	Lheat=RLVTT
3498	ENDIF
3499	ENDIF
3500	tpote(i) = tpot(i)* &
3501	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3502	ENDDO
3503
3504	IF (type_trac == 'inca') THEN
3505	#ifdef INCA
3506	CALL VTe(VTphysiq)
3507	CALL VTb(VTinca)
3508	calday = REAL(days_elapsed + 1) + jH_cur
3509
3510	call chemtime(itap+itau_phy-1, date0, dtime, itap)
3511	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3512	CALL AEROSOL_METEO_CALC( &
3513	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3514	prfl,psfl,pctsrf,cell_area, &
3515	latitude_deg,longitude_deg,u10m,v10m)
3516	END IF
3517
3518	zxsnow_dummy(:) = 0.0
3519
3520	CALL chemhook_begin (calday, &
3521	days_elapsed+1, &
3522	jH_cur, &
3523	pctsrf(1,1), &
3524	latitude_deg, &
3525	longitude_deg, &
3526	cell_area, &
3527	paprs, &
3528	pplay, &
3529	coefh(1:klon,1:klev,is_ave), &
3530	pphi, &
3531	t_seri, &
3532	u, &
3533	v, &
3534	wo(:, :, 1), &
3535	q_seri, &
3536	zxtsol, &
3537	zxsnow_dummy, &
3538	solsw, &
3539	albsol1, &
3540	rain_fall, &
3541	snow_fall, &
3542	itop_con, &
3543	ibas_con, &
3544	cldfra, &
3545	nbp_lon, &
3546	nbp_lat-1, &
3547	tr_seri, &
3548	ftsol, &
3549	paprs, &
3550	cdragh, &
3551	cdragm, &
3552	pctsrf, &
3553	pdtphys, &
3554	itap)
3555
3556	CALL VTe(VTinca)
3557	CALL VTb(VTphysiq)
3558	#endif
3559	END IF !type_trac = inca
3560	!
3561	! Calculer les parametres optiques des nuages et quelques
3562	! parametres pour diagnostiques:
3563	!
3564
3565	IF (aerosol_couple.AND.config_inca=='aero') THEN
3566	mass_solu_aero(:,:) = ccm(:,:,1)
3567	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3568	END IF
3569
3570	if (ok_newmicro) then
3571	IF (iflag_rrtm.NE.0) THEN
3572	#ifdef CPP_RRTM
3573	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3574	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3575	// 'pour ok_cdnc'
3576	call abort_physic(modname,abort_message,1)
3577	endif
3578	#else
3579
3580	abort_message='You should compile with -rrtm if running with ' &
3581	// 'iflag_rrtm=1'
3582	call abort_physic(modname,abort_message,1)
3583	#endif
3584	ENDIF
3585	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3586	paprs, pplay, t_seri, cldliq, cldfra, &
3587	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3588	flwp, fiwp, flwc, fiwc, &
3589	mass_solu_aero, mass_solu_aero_pi, &
3590	cldtaupi, re, fl, ref_liq, ref_ice, &
3591	ref_liq_pi, ref_ice_pi)
3592	else
3593	CALL nuage (paprs, pplay, &
3594	t_seri, cldliq, cldfra, cldtau, cldemi, &
3595	cldh, cldl, cldm, cldt, cldq, &
3596	ok_aie, &
3597	mass_solu_aero, mass_solu_aero_pi, &
3598	bl95_b0, bl95_b1, &
3599	cldtaupi, re, fl)
3600	endif
3601	!
3602	!IM betaCRF
3603	!
3604	cldtaurad = cldtau
3605	cldtaupirad = cldtaupi
3606	cldemirad = cldemi
3607	cldfrarad = cldfra
3608
3609	!
3610	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3611	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3612	!
3613	! global
3614	!
3615	DO k=1, klev
3616	DO i=1, klon
3617	if (pplay(i,k).GE.pfree) THEN
3618	beta(i,k) = beta_pbl
3619	else
3620	beta(i,k) = beta_free
3621	endif
3622	if (mskocean_beta) THEN
3623	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3624	endif
3625	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3626	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3627	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3628	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3629	ENDDO
3630	ENDDO
3631	!
3632	else
3633	!
3634	! regional
3635	!
3636	DO k=1, klev
3637	DO i=1,klon
3638	!
3639	if (longitude_deg(i).ge.lon1_beta.AND. &
3640	longitude_deg(i).le.lon2_beta.AND. &
3641	latitude_deg(i).le.lat1_beta.AND. &
3642	latitude_deg(i).ge.lat2_beta) THEN
3643	if (pplay(i,k).GE.pfree) THEN
3644	beta(i,k) = beta_pbl
3645	else
3646	beta(i,k) = beta_free
3647	endif
3648	if (mskocean_beta) THEN
3649	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3650	endif
3651	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3652	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3653	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3654	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3655	endif
3656	!
3657	ENDDO
3658	ENDDO
3659	!
3660	endif
3661	!
3662	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3663	!
3664	IF (MOD(itaprad,radpas).EQ.0) THEN
3665
3666	!albedo SB >>>
3667	if(ok_chlorophyll)then
3668	print*,"-- reading chlorophyll"
3669	call readchlorophyll(debut)
3670	endif
3671	!do i=1,klon
3672	!if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter)
3673	!enddo
3674	!albedo SB <<<
3675
3676
3677	if (mydebug) then
3678	call writefield_phy('u_seri',u_seri,nbp_lev)
3679	call writefield_phy('v_seri',v_seri,nbp_lev)
3680	call writefield_phy('t_seri',t_seri,nbp_lev)
3681	call writefield_phy('q_seri',q_seri,nbp_lev)
3682	endif
3683
3684	!
3685	!sonia : If Iflag_radia >=2, pertubation of some variables
3686	!input to radiation (DICE)
3687	!
3688	IF (iflag_radia .ge. 2) THEN
3689	zsav_tsol (:) = zxtsol(:)
3690	call perturb_radlwsw(zxtsol,iflag_radia)
3691	ENDIF
3692
3693	IF (aerosol_couple.AND.config_inca=='aero') THEN
3694	#ifdef INCA
3695	CALL radlwsw_inca &
3696	(kdlon,kflev,dist, rmu0, fract, solaire, &
3697	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3698	wo(:, :, 1), &
3699	cldfrarad, cldemirad, cldtaurad, &
3700	heat,heat0,cool,cool0,albpla, &
3701	topsw,toplw,solsw,sollw, &
3702	sollwdown, &
3703	topsw0,toplw0,solsw0,sollw0, &
3704	lwdn0, lwdn, lwup0, lwup, &
3705	swdn0, swdn, swup0, swup, &
3706	ok_ade, ok_aie, &
3707	tau_aero, piz_aero, cg_aero, &
3708	topswad_aero, solswad_aero, &
3709	topswad0_aero, solswad0_aero, &
3710	topsw_aero, topsw0_aero, &
3711	solsw_aero, solsw0_aero, &
3712	cldtaupirad, &
3713	topswai_aero, solswai_aero)
3714
3715	#endif
3716	ELSE
3717	!
3718	!IM calcul radiatif pour le cas actuel
3719	!
3720	RCO2 = RCO2_act
3721	RCH4 = RCH4_act
3722	RN2O = RN2O_act
3723	RCFC11 = RCFC11_act
3724	RCFC12 = RCFC12_act
3725	!
3726	IF (prt_level .GE.10) THEN
3727	print *,' ->radlwsw, number 1 '
3728	ENDIF
3729	!
3730	CALL radlwsw &
3731	(dist, rmu0, fract, &
3732	!albedo SB >>>
3733	! paprs, pplay,zxtsol,albsol1, albsol2, &
3734	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3735	!albedo SB <<<
3736	t_seri,q_seri,wo, &
3737	cldfrarad, cldemirad, cldtaurad, &
3738	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3739	flag_aerosol_strat, &
3740	tau_aero, piz_aero, cg_aero, &
3741	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3742	! Rajoute par OB pour RRTM
3743	tau_aero_lw_rrtm, &
3744	cldtaupirad,new_aod, &
3745	zqsat, flwc, fiwc, &
3746	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3747	heat,heat0,cool,cool0,albpla, &
3748	topsw,toplw,solsw,sollw, &
3749	sollwdown, &
3750	topsw0,toplw0,solsw0,sollw0, &
3751	lwdn0, lwdn, lwup0, lwup, &
3752	swdn0, swdn, swup0, swup, &
3753	topswad_aero, solswad_aero, &
3754	topswai_aero, solswai_aero, &
3755	topswad0_aero, solswad0_aero, &
3756	topsw_aero, topsw0_aero, &
3757	solsw_aero, solsw0_aero, &
3758	topswcf_aero, solswcf_aero, &
3759	!-C. Kleinschmitt for LW diagnostics
3760	toplwad_aero, sollwad_aero,&
3761	toplwai_aero, sollwai_aero, &
3762	toplwad0_aero, sollwad0_aero,&
3763	!-end
3764	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3765	ZSWFT0_i, ZFSDN0, ZFSUP0)
3766
3767	!
3768	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3769	!IM des taux doit etre different du taux actuel
3770	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3771	!
3772	if (ok_4xCO2atm) then
3773	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3774	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3775	RCFC12_per.NE.RCFC12_act) THEN
3776	!
3777	RCO2 = RCO2_per
3778	RCH4 = RCH4_per
3779	RN2O = RN2O_per
3780	RCFC11 = RCFC11_per
3781	RCFC12 = RCFC12_per
3782	!
3783	IF (prt_level .GE.10) THEN
3784	print *,' ->radlwsw, number 2 '
3785	ENDIF
3786	!
3787	CALL radlwsw &
3788	(dist, rmu0, fract, &
3789	!albedo SB >>>
3790	! paprs, pplay,zxtsol,albsol1, albsol2, &
3791	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3792	!albedo SB <<<
3793	t_seri,q_seri,wo, &
3794	cldfra, cldemi, cldtau, &
3795	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3796	flag_aerosol_strat, &
3797	tau_aero, piz_aero, cg_aero, &
3798	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3799	! Rajoute par OB pour RRTM
3800	tau_aero_lw_rrtm, &
3801	cldtaupi,new_aod, &
3802	zqsat, flwc, fiwc, &
3803	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3804	heatp,heat0p,coolp,cool0p,albplap, &
3805	topswp,toplwp,solswp,sollwp, &
3806	sollwdownp, &
3807	topsw0p,toplw0p,solsw0p,sollw0p, &
3808	lwdn0p, lwdnp, lwup0p, lwupp, &
3809	swdn0p, swdnp, swup0p, swupp, &
3810	topswad_aerop, solswad_aerop, &
3811	topswai_aerop, solswai_aerop, &
3812	topswad0_aerop, solswad0_aerop, &
3813	topsw_aerop, topsw0_aerop, &
3814	solsw_aerop, solsw0_aerop, &
3815	topswcf_aerop, solswcf_aerop, &
3816	!-C. Kleinschmitt for LW diagnostics
3817	toplwad_aerop, sollwad_aerop,&
3818	toplwai_aerop, sollwai_aerop, &
3819	toplwad0_aerop, sollwad0_aerop,&
3820	!-end
3821	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3822	ZSWFT0_i, ZFSDN0, ZFSUP0)
3823	endif
3824	endif
3825	!
3826	ENDIF ! aerosol_couple
3827	itaprad = 0
3828	!
3829	! If Iflag_radia >=2, reset pertubed variables
3830	!
3831	IF (iflag_radia .ge. 2) THEN
3832	zxtsol(:) = zsav_tsol (:)
3833	ENDIF
3834	ENDIF ! MOD(itaprad,radpas)
3835	itaprad = itaprad + 1
3836
3837	IF (iflag_radia.eq.0) THEN
3838	IF (prt_level.ge.9) THEN
3839	PRINT *,'--------------------------------------------------'
3840	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3841	PRINT *,'>>>> heat et cool mis a zero '
3842	PRINT *,'--------------------------------------------------'
3843	END IF
3844	heat=0.
3845	cool=0.
3846	sollw=0. ! MPL 01032011
3847	solsw=0.
3848	radsol=0.
3849	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3850	swup0=0.
3851	lwup=0.
3852	lwup0=0.
3853	lwdn=0.
3854	lwdn0=0.
3855	END IF
3856
3857	!
3858	! Calculer radsol a l'exterieur de radlwsw
3859	! pour prendre en compte le cycle diurne
3860	! recode par Olivier Boucher en sept 2015
3861	!
3862	radsol=solsw*swradcorr+sollw
3863	if (ok_4xCO2atm) then
3864	radsolp=solswp*swradcorr+sollwp
3865	endif
3866
3867	!
3868	! Ajouter la tendance des rayonnements (tous les pas)
3869	! avec une correction pour le cycle diurne dans le SW
3870	!
3871
3872	DO k=1, klev
3873	d_t_swr(:,k)=swradcorr(:)heat(:,k)dtime/RDAY
3874	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)dtime/RDAY
3875	d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
3876	d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
3877	ENDDO
3878
3879	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy)
3880	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy)
3881
3882	!
3883	if (mydebug) then
3884	call writefield_phy('u_seri',u_seri,nbp_lev)
3885	call writefield_phy('v_seri',v_seri,nbp_lev)
3886	call writefield_phy('t_seri',t_seri,nbp_lev)
3887	call writefield_phy('q_seri',q_seri,nbp_lev)
3888	endif
3889
3890	!IM
3891	IF (ip_ebil_phy.ge.2) THEN
3892	ztit='after rad'
3893	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3894	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3895	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3896	call diagphy(cell_area,ztit,ip_ebil_phy &
3897	, topsw, toplw, solsw, sollw, zero_v &
3898	, zero_v, zero_v, zero_v, ztsol &
3899	, d_h_vcol, d_qt, d_ec &
3900	, fs_bound, fq_bound )
3901	END IF
3902	!
3903	!
3904	! Calculer l'hydrologie de la surface
3905	!
3906	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3907	! . agesno, ftsol,fqsurf,fsnow, ruis)
3908	!
3909
3910	!
3911	! Calculer le bilan du sol et la derive de temperature (couplage)
3912	!
3913	DO i = 1, klon
3914	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3915	! a la demande de JLD
3916	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3917	ENDDO
3918	!
3919	!moddeblott(jan95)
3920	! Appeler le programme de parametrisation de l'orographie
3921	! a l'echelle sous-maille:
3922	!
3923	IF (prt_level .GE.10) THEN
3924	print *,' call orography ? ', ok_orodr
3925	ENDIF
3926	!
3927	IF (ok_orodr) THEN
3928	!
3929	! selection des points pour lesquels le shema est actif:
3930	igwd=0
3931	DO i=1,klon
3932	itest(i)=0
3933	! IF ((zstd(i).gt.10.0)) THEN
3934	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3935	itest(i)=1
3936	igwd=igwd+1
3937	idx(igwd)=i
3938	ENDIF
3939	ENDDO
3940	! igwdim=MAX(1,igwd)
3941	!
3942	IF (ok_strato) THEN
3943
3944	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3945	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3946	igwd,idx,itest, &
3947	t_seri, u_seri, v_seri, &
3948	zulow, zvlow, zustrdr, zvstrdr, &
3949	d_t_oro, d_u_oro, d_v_oro)
3950
3951	ELSE
3952	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3953	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3954	igwd,idx,itest, &
3955	t_seri, u_seri, v_seri, &
3956	zulow, zvlow, zustrdr, zvstrdr, &
3957	d_t_oro, d_u_oro, d_v_oro)
3958	ENDIF
3959	!
3960	! ajout des tendances
3961	!-----------------------------------------------------------------------
3962	! ajout des tendances de la trainee de l'orographie
3963	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
3964	abortphy)
3965	!----------------------------------------------------------------------
3966	!
3967	ENDIF ! fin de test sur ok_orodr
3968	!
3969	if (mydebug) then
3970	call writefield_phy('u_seri',u_seri,nbp_lev)
3971	call writefield_phy('v_seri',v_seri,nbp_lev)
3972	call writefield_phy('t_seri',t_seri,nbp_lev)
3973	call writefield_phy('q_seri',q_seri,nbp_lev)
3974	endif
3975
3976	IF (ok_orolf) THEN
3977	!
3978	! selection des points pour lesquels le shema est actif:
3979	igwd=0
3980	DO i=1,klon
3981	itest(i)=0
3982	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3983	itest(i)=1
3984	igwd=igwd+1
3985	idx(igwd)=i
3986	ENDIF
3987	ENDDO
3988	! igwdim=MAX(1,igwd)
3989	!
3990	IF (ok_strato) THEN
3991
3992	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3993	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3994	igwd,idx,itest, &
3995	t_seri, u_seri, v_seri, &
3996	zulow, zvlow, zustrli, zvstrli, &
3997	d_t_lif, d_u_lif, d_v_lif )
3998
3999	ELSE
4000	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
4001	latitude_deg,zmea,zstd,zpic, &
4002	itest, &
4003	t_seri, u_seri, v_seri, &
4004	zulow, zvlow, zustrli, zvstrli, &
4005	d_t_lif, d_u_lif, d_v_lif)
4006	ENDIF
4007
4008	! ajout des tendances de la portance de l'orographie
4009	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4010	'lif', abortphy)
4011	ENDIF ! fin de test sur ok_orolf
4012
4013	IF (ok_hines) then
4014	! HINES GWD PARAMETRIZATION
4015	east_gwstress=0.
4016	west_gwstress=0.
4017	du_gwd_hines=0.
4018	dv_gwd_hines=0.
4019	CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
4020	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4021	du_gwd_hines, dv_gwd_hines)
4022	zustr_gwd_hines=0.
4023	zvstr_gwd_hines=0.
4024	DO k = 1, klev
4025	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
4026	* (paprs(:, k)-paprs(:, k+1))/rg
4027	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
4028	* (paprs(:, k)-paprs(:, k+1))/rg
4029	ENDDO
4030
4031	d_t_hin(:, :)=0.
4032	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4033	dqi0, paprs, 'hin', abortphy)
4034	ENDIF
4035
4036	IF (.not. ok_hines .and. ok_gwd_rando) then
4037	CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
4038	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4039	dv_gwd_front, east_gwstress, west_gwstress)
4040	zustr_gwd_front=0.
4041	zvstr_gwd_front=0.
4042	DO k = 1, klev
4043	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
4044	* (paprs(:, k)-paprs(:, k+1))/rg
4045	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
4046	* (paprs(:, k)-paprs(:, k+1))/rg
4047	ENDDO
4048
4049	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4050	paprs, 'front_gwd_rando', abortphy)
4051	ENDIF
4052
4053	if (ok_gwd_rando) then
4054	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
4055	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4056	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4057	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4058	paprs, 'flott_gwd_rando', abortphy)
4059	zustr_gwd_rando=0.
4060	zvstr_gwd_rando=0.
4061	DO k = 1, klev
4062	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
4063	* (paprs(:, k)-paprs(:, k+1))/rg
4064	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
4065	* (paprs(:, k)-paprs(:, k+1))/rg
4066	ENDDO
4067	end if
4068
4069	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4070
4071	if (mydebug) then
4072	call writefield_phy('u_seri',u_seri,nbp_lev)
4073	call writefield_phy('v_seri',v_seri,nbp_lev)
4074	call writefield_phy('t_seri',t_seri,nbp_lev)
4075	call writefield_phy('q_seri',q_seri,nbp_lev)
4076	endif
4077
4078	DO i = 1, klon
4079	zustrph(i)=0.
4080	zvstrph(i)=0.
4081	ENDDO
4082	DO k = 1, klev
4083	DO i = 1, klon
4084	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
4085	(paprs(i,k)-paprs(i,k+1))/rg
4086	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
4087	(paprs(i,k)-paprs(i,k+1))/rg
4088	ENDDO
4089	ENDDO
4090	!
4091	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4092	!
4093	IF (is_sequential .and. ok_orodr) THEN
4094	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4095	ra,rg,romega, &
4096	latitude_deg,longitude_deg,pphis, &
4097	zustrdr,zustrli,zustrph, &
4098	zvstrdr,zvstrli,zvstrph, &
4099	paprs,u,v, &
4100	aam, torsfc)
4101	ENDIF
4102	!IM cf. FLott END
4103	!IM
4104	IF (ip_ebil_phy.ge.2) THEN
4105	ztit='after orography'
4106	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
4107	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
4108	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4109	call diagphy(cell_area,ztit,ip_ebil_phy &
4110	, zero_v, zero_v, zero_v, zero_v, zero_v &
4111	, zero_v, zero_v, zero_v, ztsol &
4112	, d_h_vcol, d_qt, d_ec &
4113	, fs_bound, fq_bound )
4114	END IF
4115
4116	!DC Calcul de la tendance due au methane
4117	IF(ok_qch4) THEN
4118	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4119	! ajout de la tendance d'humidite due au methane
4120	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4*dtime, dql0, dqi0, paprs, &
4121	'q_ch4', abortphy)
4122	END IF
4123	!
4124	!
4125	!====================================================================
4126	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4127	!====================================================================
4128	! Abderrahmane 24.08.09
4129
4130	IF (ok_cosp) THEN
4131	! adeclarer
4132	#ifdef CPP_COSP
4133	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4134
4135	IF (prt_level .GE.10) THEN
4136	print*,'freq_cosp',freq_cosp
4137	ENDIF
4138	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4139	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4140	! s ref_liq,ref_ice
4141	call phys_cosp(itap,dtime,freq_cosp, &
4142	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4143	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
4144	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4145	JrNt,ref_liq,ref_ice, &
4146	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4147	zu10m,zv10m,pphis, &
4148	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4149	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4150	prfl(:,1:klev),psfl(:,1:klev), &
4151	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4152	mr_ozone,cldtau, cldemi)
4153
4154	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4155	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4156	! M clMISR,
4157	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4158	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4159
4160	ENDIF
4161
4162	#endif
4163	ENDIF !ok_cosp
4164	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4165	!AA
4166	!AA Installation de l'interface online-offline pour traceurs
4167	!AA
4168	!====================================================================
4169	! Calcul des tendances traceurs
4170	!====================================================================
4171	!
4172
4173	IF (type_trac=='repr') THEN
4174	sh_in(:,:) = q_seri(:,:)
4175	ELSE
4176	sh_in(:,:) = qx(:,:,ivap)
4177	END IF
4178
4179	call phytrac ( &
4180	itap, days_elapsed+1, jH_cur, debut, &
4181	lafin, dtime, u, v, t, &
4182	paprs, pplay, pmfu, pmfd, &
4183	pen_u, pde_u, pen_d, pde_d, &
4184	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4185	u1, v1, ftsol, pctsrf, &
4186	zustar, zu10m, zv10m, &
4187	wstar(:,is_ave), ale_bl, ale_wake, &
4188	latitude_deg, longitude_deg, &
4189	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4190	presnivs, pphis, pphi, albsol1, &
4191	sh_in, rhcl, cldfra, rneb, &
4192	diafra, cldliq, itop_con, ibas_con, &
4193	pmflxr, pmflxs, prfl, psfl, &
4194	da, phi, mp, upwd, &
4195	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4196	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4197	ev, ep, epmlmMm, eplaMm, & !<<RomP
4198	dnwd, aerosol_couple, flxmass_w, &
4199	tau_aero, piz_aero, cg_aero, ccm, &
4200	rfname, &
4201	d_tr_dyn, & !<<RomP
4202	tr_seri)
4203
4204	IF (offline) THEN
4205
4206	IF (prt_level.ge.9) &
4207	print*,'Attention on met a 0 les thermiques pour phystoke'
4208	call phystokenc ( &
4209	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4210	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4211	fm_therm,entr_therm, &
4212	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4213	frac_impa, frac_nucl, &
4214	pphis,cell_area,dtime,itap, &
4215	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4216
4217
4218	ENDIF
4219
4220	!
4221	! Calculer le transport de l'eau et de l'energie (diagnostique)
4222	!
4223	CALL transp (paprs,zxtsol, &
4224	t_seri, q_seri, u_seri, v_seri, zphi, &
4225	ve, vq, ue, uq)
4226	!
4227	!IM global posePB BEG
4228	IF(1.EQ.0) THEN
4229	!
4230	CALL transp_lay (paprs,zxtsol, &
4231	t_seri, q_seri, u_seri, v_seri, zphi, &
4232	ve_lay, vq_lay, ue_lay, uq_lay)
4233	!
4234	ENDIF !(1.EQ.0) THEN
4235	!IM global posePB END
4236	! Accumuler les variables a stocker dans les fichiers histoire:
4237	!
4238
4239	!================================================================
4240	! Conversion of kinetic and potential energy into heat, for
4241	! parameterisation of subgrid-scale motions
4242	!================================================================
4243
4244	d_t_ec(:,:)=0.
4245	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4246	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
4247	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4248	zmasse,exner,d_t_ec)
4249	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4250
4251	!IM
4252	IF (ip_ebil_phy.ge.1) THEN
4253	ztit='after physic'
4254	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
4255	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
4256	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
4257	! Comme les tendances de la physique sont ajoute dans la dynamique,
4258	! on devrait avoir que la variation d'entalpie par la dynamique
4259	! est egale a la variation de la physique au pas de temps precedent.
4260	! Donc la somme de ces 2 variations devrait etre nulle.
4261
4262	call diagphy(cell_area,ztit,ip_ebil_phy &
4263	, topsw, toplw, solsw, sollw, sens &
4264	, evap, rain_fall, snow_fall, ztsol &
4265	, d_h_vcol, d_qt, d_ec &
4266	, fs_bound, fq_bound )
4267	!
4268	d_h_vcol_phy=d_h_vcol
4269	!
4270	END IF
4271	!
4272	!=======================================================================
4273	! SORTIES
4274	!=======================================================================
4275	!
4276	!IM initialisation + calculs divers diag AMIP2
4277	!
4278	include "calcul_divers.h"
4279	!
4280	!IM Interpolation sur les niveaux de pression du NMC
4281	! -------------------------------------------------
4282	#ifdef CPP_XIOS
4283	!$OMP MASTER
4284	!On recupere la valeur de la missing value donnee dans le xml
4285	CALL xios_get_field_attr("t850",default_value=missing_val_omp)
4286	! PRINT *,"ARNAUD value missing ",missing_val_omp
4287	!$OMP END MASTER
4288	!$OMP BARRIER
4289	missing_val=missing_val_omp
4290	#endif
4291	#ifndef CPP_XIOS
4292	missing_val=missing_val_nf90
4293	#endif
4294	!
4295	include "calcul_STDlev.h"
4296	!
4297	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4298	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4299	!
4300	!cc prw = eau precipitable
4301	DO i = 1, klon
4302	prw(i) = 0.
4303	DO k = 1, klev
4304	prw(i) = prw(i) + &
4305	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
4306	ENDDO
4307	ENDDO
4308	!
4309	IF (type_trac == 'inca') THEN
4310	#ifdef INCA
4311	CALL VTe(VTphysiq)
4312	CALL VTb(VTinca)
4313
4314	CALL chemhook_end ( &
4315	dtime, &
4316	pplay, &
4317	t_seri, &
4318	tr_seri, &
4319	nbtr, &
4320	paprs, &
4321	q_seri, &
4322	cell_area, &
4323	pphi, &
4324	pphis, &
4325	zx_rh)
4326
4327	CALL VTe(VTinca)
4328	CALL VTb(VTphysiq)
4329	#endif
4330	END IF
4331
4332
4333	!
4334	! Convertir les incrementations en tendances
4335	!
4336	IF (prt_level .GE.10) THEN
4337	print *,'Convertir les incrementations en tendances '
4338	ENDIF
4339	!
4340	if (mydebug) then
4341	call writefield_phy('u_seri',u_seri,nbp_lev)
4342	call writefield_phy('v_seri',v_seri,nbp_lev)
4343	call writefield_phy('t_seri',t_seri,nbp_lev)
4344	call writefield_phy('q_seri',q_seri,nbp_lev)
4345	endif
4346
4347	DO k = 1, klev
4348	DO i = 1, klon
4349	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4350	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4351	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4352	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4353	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4354	!CR: on ajoute le contenu en glace
4355	if (nqo.eq.3) then
4356	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
4357	endif
4358	ENDDO
4359	ENDDO
4360	!
4361	!CR: nb de traceurs eau: nqo
4362	! IF (nqtot.GE.3) THEN
4363	IF (nqtot.GE.(nqo+1)) THEN
4364	! DO iq = 3, nqtot
4365	DO iq = nqo+1, nqtot
4366	DO k = 1, klev
4367	DO i = 1, klon
4368	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4369	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4370	ENDDO
4371	ENDDO
4372	ENDDO
4373	ENDIF
4374	!
4375	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4376	!IM global posePB include "write_bilKP_ins.h"
4377	!IM global posePB include "write_bilKP_ave.h"
4378	!
4379
4380	! Sauvegarder les valeurs de t et q a la fin de la physique:
4381	!
4382	DO k = 1, klev
4383	DO i = 1, klon
4384	u_ancien(i,k) = u_seri(i,k)
4385	v_ancien(i,k) = v_seri(i,k)
4386	t_ancien(i,k) = t_seri(i,k)
4387	q_ancien(i,k) = q_seri(i,k)
4388	ENDDO
4389	ENDDO
4390
4391	! !! RomP >>>
4392	!CR: nb de traceurs eau: nqo
4393	! IF (nqtot.GE.3) THEN
4394	IF (nqtot.GE.(nqo+1)) THEN
4395	! DO iq = 3, nqtot
4396	DO iq = nqo+1, nqtot
4397	DO k = 1, klev
4398	DO i = 1, klon
4399	! tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4400	tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
4401	ENDDO
4402	ENDDO
4403	ENDDO
4404	ENDIF
4405	! !! RomP <<<
4406	!==========================================================================
4407	! Sorties des tendances pour un point particulier
4408	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4409	! pour le debug
4410	! La valeur de igout est attribuee plus haut dans le programme
4411	!==========================================================================
4412
4413	if (prt_level.ge.1) then
4414	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4415	write(lunout,*) &
4416	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4417	write(lunout,*) &
4418	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4419	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4420	pctsrf(igout,is_sic)
4421	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4422	do k=1,klev
4423	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4424	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4425	d_t_eva(igout,k)
4426	enddo
4427	write(lunout,*) 'cool,heat'
4428	do k=1,klev
4429	write(lunout,*) cool(igout,k),heat(igout,k)
4430	enddo
4431
4432	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4433	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4434	!jyg! do k=1,klev
4435	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4436	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4437	!jyg! enddo
4438	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4439	do k=1,klev
4440	write(lunout,*) d_t_vdf(igout,k), &
4441	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4442	enddo
4443	!>jyg
4444
4445	write(lunout,*) 'd_ps ',d_ps(igout)
4446	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4447	do k=1,klev
4448	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4449	d_qx(igout,k,1),d_qx(igout,k,2)
4450	enddo
4451	endif
4452
4453	!==========================================================================
4454
4455	!============================================================
4456	! Calcul de la temperature potentielle
4457	!============================================================
4458	DO k = 1, klev
4459	DO i = 1, klon
4460	!JYG/IM theta en debut du pas de temps
4461	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4462	!JYG/IM theta en fin de pas de temps de physique
4463	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4464	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4465	! MPL 20130625
4466	! fth_fonctions.F90 et parkind1.F90
4467	! sinon thetal=theta
4468	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4469	! : ql_seri(i,k))
4470	thetal(i,k)=theta(i,k)
4471	ENDDO
4472	ENDDO
4473	!
4474
4475	! 22.03.04 BEG
4476	!=============================================================
4477	! Ecriture des sorties
4478	!=============================================================
4479	#ifdef CPP_IOIPSL
4480
4481	! Recupere des varibles calcule dans differents modules
4482	! pour ecriture dans histxxx.nc
4483
4484	! Get some variables from module fonte_neige_mod
4485	CALL fonte_neige_get_vars(pctsrf, &
4486	zxfqcalving, zxfqfonte, zxffonte)
4487
4488
4489
4490
4491	!=============================================================
4492	! Separation entre thermiques et non thermiques dans les sorties
4493	! de fisrtilp
4494	!=============================================================
4495
4496	if (iflag_thermals>=1) then
4497	d_t_lscth=0.
4498	d_t_lscst=0.
4499	d_q_lscth=0.
4500	d_q_lscst=0.
4501	do k=1,klev
4502	do i=1,klon
4503	if (ptconvth(i,k)) then
4504	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4505	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4506	else
4507	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4508	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4509	endif
4510	enddo
4511	enddo
4512
4513	do i=1,klon
4514	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4515	plul_th(i)=prfl(i,1)+psfl(i,1)
4516	enddo
4517	endif
4518
4519
4520	!On effectue les sorties:
4521
4522	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4523	pplay, lmax_th, aerosol_couple, &
4524	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4525	ptconv, read_climoz, clevSTD, &
4526	ptconvth, d_t, qx, d_qx, zmasse, &
4527	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4528
4529
4530
4531	include "write_histday_seri.h"
4532
4533	include "write_paramLMDZ_phy.h"
4534
4535	#endif
4536
4537
4538	!====================================================================
4539	! Arret du modele apres hgardfou en cas de detection d'un
4540	! plantage par hgardfou
4541	!====================================================================
4542
4543	IF (abortphy==1) THEN
4544	abort_message ='Plantage hgardfou'
4545	CALL abort_physic (modname,abort_message,1)
4546	ENDIF
4547
4548	!--OB mass fixer
4549	!--profile is corrected to force mass conservation of water
4550	IF (mass_fixer) THEN
4551	qql2(:)=0.0
4552	DO i = 1, klon
4553	DO k = 1, klev
4554	zdz=(paprs(i,k)-paprs(i,k+1))/RG
4555	qql2(i)=qql2(i)+(q_seri(i,k)+ql_seri(i,k))*zdz
4556	ENDDO
4557	ENDDO
4558	DO i = 1, klon
4559	!--compute ratio of what q+ql should be with conservation to what it is
4560	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4561	DO k = 1, klev
4562	q_seri(i,k) =q_seri(i,k)*corrqql
4563	ql_seri(i,k)=ql_seri(i,k)*corrqql
4564	ENDDO
4565	ENDDO
4566	ENDIF
4567	!--fin mass fixer
4568
4569	! 22.03.04 END
4570	!
4571	!====================================================================
4572	! Si c'est la fin, il faut conserver l'etat de redemarrage
4573	!====================================================================
4574	!
4575
4576	IF (lafin) THEN
4577	itau_phy = itau_phy + itap
4578	CALL phyredem ("restartphy.nc")
4579	! open(97,form="unformatted",file="finbin")
4580	! write(97) u_seri,v_seri,t_seri,q_seri
4581	! close(97)
4582	!$OMP MASTER
4583	if (read_climoz >= 1) then
4584	if (is_mpi_root) then
4585	call nf95_close(ncid_climoz)
4586	end if
4587	deallocate(press_climoz) ! pointer
4588	end if
4589	!$OMP END MASTER
4590	ENDIF
4591
4592	! first=.false.
4593
4594
4595	END SUBROUTINE physiq
4596
4597	END MODULE physiq_mod

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

Download in other formats:

Original Format