Context Navigation

physiq.F @ 1740

Last change on this file since 1740 was 1740, checked in by lguez, 11 years ago

Fixing bug from revision 1539 (see also revision 1604). Allocatable
arrays "tabijgcm", "longcm", "latgcm", "igcm", "jgcm" cannot be
arguments of "phys_output_open" if they have not been
allocated. Allocate them with zero size when iflag_con < 3.

Replaced non-ASCII characters in comments by ASCII sequences. There
are problems with non-ASCII characters, see for example line 2238:
some information is lost here.

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

File size: 133.1 KB

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

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

Download in other formats:

Original Format