Context Navigation

physiq_mod.F90 @ 2502

Last change on this file since 2502 was 2501, checked in by oboucher, 9 years ago
Adding a template for tropospheric aerosol LW optical properties Only works for aerosol_couple = TRUE
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 156.7 KB

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

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

Download in other formats:

Original Format