Context Navigation

physiq.F @ 1707

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

Version testing basée sur la r1706

Testing release based on r1706

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

File size: 132.4 KB

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

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

Download in other formats:

Original Format