Context Navigation

physiq_mod.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by oboucher, 8 years ago
I moved ini_paramLMDZ_phy and write_paramLMDZ_phy into a module called paramlmdz_phy_mod.F90. Output frequency is increased to daily. RSUN is now an output if iflag_rrtm = 1
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: 149.7 KB

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

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

Download in other formats:

Original Format