Context Navigation

physiq.F90 @ 1955

Last change on this file since 1955 was 1955, checked in by fhourdin, 11 years ago
Options for idealized ozone (read_climoz=0) : Imposing ozone be constant in case were solarlong0> -999 (defaut value is -999.999 when seasonal cycle is compiuted) and to be symetric with respect to equator when read_climoz=-1 FH
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: 120.9 KB

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

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

Download in other formats:

Original Format