Context Navigation

physiq.F90 @ 2098

Last change on this file since 2098 was 2098, checked in by lguez, 10 years ago

Replaced 360 in calbeta_clim by length of current year according to
chosen calendar. Length of current year is given by
ioget_year_len. Since we already need this for ozone, moved the call
to ioget_year_len from physiq to phys_cal_mod and created variable
year_len of module phys_cal_mod.

Control the output from minmaxqfi.

Non-ASCII characters in comments are not always rendered properly and
they risk being lost. See revision 1740.

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: 125.9 KB

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

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

Download in other formats:

Original Format