Context Navigation

physiq_mod.F90 @ 2525

Last change on this file since 2525 was 2525, checked in by oboucher, 8 years ago
Adding a compilation directive in case the model is not compiled with -rrtm true This fixes the previous commit
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.6 KB

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

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

Download in other formats:

Original Format