Context Navigation

physiq_mod.F90 @ 2513

Last change on this file since 2513 was 2513, checked in by jyg, 8 years ago
Creation of two new subroutines containing all the Ale and Alp stuff previously present in physiq_mod.F90: files alpale.F90 and alpale_th.F90.
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: 147.3 KB

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

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

Download in other formats:

Original Format