Context Navigation

physiq_mod.F90 @ 2496

Last change on this file since 2496 was 2496, checked in by oboucher, 8 years ago
Adding more diagnostics for water conservation dqlphys, dqldyn, dqsphys, dqsdyn, prlw, prsw
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: 157.5 KB

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

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

Download in other formats:

Original Format