Context Navigation

physiq.F @ 1785

Last change on this file since 1785 was 1785, checked in by Ehouarn Millour, 11 years ago

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG

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

File size: 136.8 KB

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

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

Download in other formats:

Original Format