Context Navigation

physiq.F90 @ 1865

Last change on this file since 1865 was 1865, checked in by Laurent Fairhead, 11 years ago

Inclusion de la bibliothèque SISVAT/MAR à LMDZ pour le traitement des surfaces
"land ice"

Menegoz

Integration of the SISVAT/MAR library to LMDZ to model the land ice surfaces

Menegoz

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 130.0 KB

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

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

Download in other formats:

Original Format