Context Navigation

physiq.F @ 955

Last change on this file since 955 was 955, checked in by lsce, 16 years ago

ACo+JG

Ajout du flag aerosol_couple :
false=lecture des sulfates dans un fichier(par defaut) - configuration existante
true=calcul des aerosol par INCA

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

File size: 109.7 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c#define IO_DEBUG
6
7	SUBROUTINE physiq (nlon,nlev,nqmax,
8	. debut,lafin,rjourvrai,gmtime,pdtphys,
9	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
10	. u,v,t,qx,
11	. omega,
12	#ifdef INCA
13	. flxmass_w,
14	#endif
15	. d_u, d_v, d_t, d_qx, d_ps
16	. , dudyn
17	. , PVteta)
18
19	USE ioipsl
20	USE comgeomphy
21	USE write_field_phy
22	USE dimphy
23	USE mod_grid_phy_lmdz
24	USE mod_phys_lmdz_para
25	USE iophy
26	USE misc_mod, mydebug=>debug
27	USE vampir
28	USE pbl_surface_mod, ONLY : pbl_surface
29	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
30	USE phys_state_var_mod ! Variables sauvegardees de la physique
31
32
33	USE ocean_slab_mod, ONLY : ocean_slab_get_vars
34	USE ocean_cpl_mod, ONLY : ocean_cpl_get_vars
35	USE ocean_forced_mod, ONLY : ocean_forced_get_vars
36	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
37	USE phys_output_mod
38
39	IMPLICIT none
40	c======================================================================
41	c
42	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
43	c
44	c Objet: Moniteur general de la physique du modele
45	cAA Modifications quant aux traceurs :
46	cAA - uniformisation des parametrisations ds phytrac
47	cAA - stockage des moyennes des champs necessaires
48	cAA en mode traceur off-line
49	c======================================================================
50	c CLEFS CPP POUR LES IO
51	c =====================
52	c#define histhf
53	#define histday
54	#define histmth
55	c#define histmthNMC
56	c#define histins
57	c#define histISCCP
58	c======================================================================
59	c modif ( P. Le Van , 12/10/98 )
60	c
61	c Arguments:
62	c
63	c nlon----input-I-nombre de points horizontaux
64	c nlev----input-I-nombre de couches verticales
65	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
66	c debut---input-L-variable logique indiquant le premier passage
67	c lafin---input-L-variable logique indiquant le dernier passage
68	c rjour---input-R-numero du jour de l'experience
69	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
70	c pdtphys-input-R-pas d'integration pour la physique (seconde)
71	c paprs---input-R-pression pour chaque inter-couche (en Pa)
72	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
73	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
74	c pphis---input-R-geopotentiel du sol
75	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
76	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
77	c v-------input-R-vitesse Y (de S a N) en m/s
78	c t-------input-R-temperature (K)
79	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
80	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
81	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
82	c omega---input-R-vitesse verticale en Pa/s
83	c d_u-----output-R-tendance physique de "u" (m/s/s)
84	c d_v-----output-R-tendance physique de "v" (m/s/s)
85	c d_t-----output-R-tendance physique de "t" (K/s)
86	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
87	c d_ps----output-R-tendance physique de la pression au sol
88	cIM
89	c PVteta--output-R-vorticite potentielle a des thetas constantes
90	c======================================================================
91	#include "dimensions.h"
92	integer jjmp1
93	parameter (jjmp1=jjm+1-1/jjm)
94	integer iip1
95	parameter (iip1=iim+1)
96
97	#include "regdim.h"
98	#include "indicesol.h"
99	#include "dimsoil.h"
100	#include "clesphys.h"
101	#include "control.h"
102	#include "logic.h"
103	#include "temps.h"
104	cym#include "comgeomphy.h"
105	#include "advtrac.h"
106	#include "iniprint.h"
107	#include "thermcell.h"
108	c======================================================================
109	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
110	PARAMETER (ok_cvl=.TRUE.)
111	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
112	PARAMETER (ok_gust=.FALSE.)
113	integer iflag_radia ! active ou non le rayonnement (MPL)
114	save iflag_radia
115	c======================================================================
116	LOGICAL check ! Verifier la conservation du modele en eau
117	PARAMETER (check=.FALSE.)
118	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
119	PARAMETER (ok_stratus=.FALSE.)
120	c======================================================================
121	logical rnpb
122	#ifdef INCA
123	parameter(rnpb=.false.)
124	#else
125	parameter(rnpb=.true.)
126	#endif
127
128	cIM "slab" ocean
129	REAL tslab(klon) !Temperature du slab-ocean
130	REAL seaice(klon) !glace de mer (kg/m2)
131	REAL fluxo(klon) !flux turbulents ocean-glace de mer
132	REAL fluxg(klon) !flux turbulents ocean-atmosphere
133	REAL amn, amx
134	INTEGER igout
135	c======================================================================
136	c Clef controlant l'activation du cycle diurne:
137	ccc LOGICAL cycle_diurne
138	ccc PARAMETER (cycle_diurne=.FALSE.)
139	c======================================================================
140	c Modele thermique du sol, a activer pour le cycle diurne:
141	ccc LOGICAL soil_model
142	ccc PARAMETER (soil_model=.FALSE.)
143	c======================================================================
144	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
145	c le calcul du rayonnement est celle apres la precipitation des nuages.
146	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
147	c la condensation et la precipitation. Cette cle augmente les impacts
148	c radiatifs des nuages.
149	ccc LOGICAL new_oliq
150	ccc PARAMETER (new_oliq=.FALSE.)
151	c======================================================================
152	c Clefs controlant deux parametrisations de l'orographie:
153	cc LOGICAL ok_orodr
154	ccc PARAMETER (ok_orodr=.FALSE.)
155	ccc LOGICAL ok_orolf
156	ccc PARAMETER (ok_orolf=.FALSE.)
157	c======================================================================
158	LOGICAL ok_journe ! sortir le fichier journalier
159	save ok_journe
160	c$OMP THREADPRIVATE(ok_journe)
161	c
162	LOGICAL ok_mensuel ! sortir le fichier mensuel
163	save ok_mensuel
164	c$OMP THREADPRIVATE(ok_mensuel)
165	c
166	LOGICAL ok_instan ! sortir le fichier instantane
167	save ok_instan
168	c$OMP THREADPRIVATE(ok_instan)
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
181	integer lmax_th(klon)
182	integer limbas(klon)
183	real ratqscth(klon,klev)
184	real ratqsdiff(klon,klev)
185	real zqsatth(klon,klev)
186
187	c======================================================================
188	c
189	INTEGER ivap ! indice de traceurs pour vapeur d'eau
190	PARAMETER (ivap=1)
191	INTEGER iliq ! indice de traceurs pour eau liquide
192	PARAMETER (iliq=2)
193
194	c
195	c
196	c Variables argument:
197	c
198	INTEGER nlon
199	INTEGER nlev
200	INTEGER nqmax
201	REAL rjourvrai
202	REAL gmtime
203	REAL pdtphys
204	LOGICAL debut, lafin
205	REAL paprs(klon,klev+1)
206	REAL pplay(klon,klev)
207	REAL pphi(klon,klev)
208	REAL pphis(klon)
209	REAL presnivs(klev)
210	REAL znivsig(klev)
211	real pir
212
213	REAL u(klon,klev)
214	REAL v(klon,klev)
215	REAL t(klon,klev),theta(klon,klev)
216	REAL qx(klon,klev,nqmax)
217
218	REAL omega(klon,klev)
219
220	#ifdef INCA
221	REAL flxmass_w(klon,klev)
222	#endif
223	REAL d_u(klon,klev)
224	REAL d_v(klon,klev)
225	REAL d_t(klon,klev)
226	REAL d_qx(klon,klev,nqmax)
227	REAL d_ps(klon)
228	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
229	c
230	cIM Amip2 PV a theta constante
231	c
232	INTEGER nbteta
233	PARAMETER(nbteta=3)
234	CHARACTER*3 ctetaSTD(nbteta)
235	DATA ctetaSTD/'350','380','405'/
236	SAVE ctetaSTD
237	c$OMP THREADPRIVATE(ctetaSTD)
238	REAL rtetaSTD(nbteta)
239	DATA rtetaSTD/350., 380., 405./
240	SAVE rtetaSTD
241	c$OMP THREADPRIVATE(rtetaSTD)
242	c
243	REAL PVteta(klon,nbteta)
244	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
245	c
246	cMI Amip2 PV a theta constante
247
248	cym INTEGER klevp1, klevm1
249	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
250	cym#include "raddim.h"
251	c
252	c
253	cIM Amip2
254	c variables a une pression donnee
255	c
256	real rlevSTD(nlevSTD)
257	DATA rlevSTD/100000., 92500., 85000., 70000.,
258	.60000., 50000., 40000., 30000., 25000., 20000.,
259	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
260	SAVE rlevstd
261	CHARACTER*4 clevSTD(nlevSTD)
262	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
263	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
264	.'70 ','50 ','30 ','20 ','10 '/
265	SAVE clevSTD
266	c
267	CHARACTER*4 bb2
268	CHARACTER*2 bb3
269	c
270	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
271	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
272	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
273	real wlevSTD(klon,nlevSTD)
274	c
275	c nout : niveau de output des variables a une pression donnee
276	logical oknondef(klon,nlevSTD,nout)
277	c
278	c les produits uvSTD, vqSTD, .., T2STD sont calcules
279	c a partir des valeurs instantannees toutes les 6 h
280	c qui sont moyennees sur le mois
281	c
282	real uvSTD(klon,nlevSTD)
283	real vqSTD(klon,nlevSTD)
284	real vTSTD(klon,nlevSTD)
285	real wqSTD(klon,nlevSTD)
286	c
287	real vphiSTD(klon,nlevSTD)
288	real wTSTD(klon,nlevSTD)
289	real u2STD(klon,nlevSTD)
290	real v2STD(klon,nlevSTD)
291	real T2STD(klon,nlevSTD)
292	c
293	#include "radepsi.h"
294	#include "radopt.h"
295	c
296	c
297	c prw: precipitable water
298	real prw(klon)
299
300	REAL convliq(klon,klev) ! eau liquide nuageuse convective
301	REAL convfra(klon,klev) ! fraction nuageuse convective
302
303	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
304	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
305	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
306	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
307
308	INTEGER linv, kp1
309	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
310	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
311	REAL flwp(klon), fiwp(klon)
312	REAL flwc(klon,klev), fiwc(klon,klev)
313	REAL flwp_c(klon), fiwp_c(klon)
314	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
315	REAL flwp_s(klon), fiwp_s(klon)
316	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
317
318	cIM ISCCP simulator v3.4
319	c dans clesphys.h top_height, overlap
320	cv3.4
321	INTEGER debug, debugcol
322	cym INTEGER npoints
323	cym PARAMETER(npoints=klon)
324	c
325	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
326	INTEGER nregISCtot
327	PARAMETER(nregISCtot=1)
328	c
329	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
330	c y compris pour 1 point
331	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
332	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
333	INTEGER imin_debut, nbpti
334	INTEGER jmin_debut, nbptj
335	cIM parametres ISCCP BEG
336	INTEGER nbapp_isccp
337	! INTEGER nbapp_isccp,isccppas
338	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
339	! !i.e. toutes les 3 heures
340	INTEGER n
341	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
342	DATA ifreq_isccp/3/
343	SAVE ifreq_isccp
344	c$OMP THREADPRIVATE(ifreq_isccp)
345	CHARACTER*5 typinout(napisccp)
346	DATA typinout/'i3od'/
347	SAVE typinout
348	c$OMP THREADPRIVATE(typinout)
349	cIM verif boxptop BEG
350	CHARACTER*1 verticaxe(napisccp)
351	DATA verticaxe/'1'/
352	SAVE verticaxe
353	c$OMP THREADPRIVATE(verticaxe)
354	cIM verif boxptop END
355	INTEGER nvlev(napisccp)
356	c INTEGER nvlev
357	REAL t1, aa
358	REAL seed_re(klon,napisccp)
359	cym !!!! A voir plus tard
360	cym INTEGER iphy(iim,jjmp1)
361	cIM parametres ISCCP END
362	c
363	c ncol = nb. de sous-colonnes pour chaque maille du GCM
364	c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
365	c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
366	INTEGER,SAVE :: ncol(napisccp)
367	INTEGER ncolmx, seed(klon,napisccp)
368	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
369	c PARAMETER(ncolmx=1500)
370	PARAMETER(ncolmx=300)
371	c
372	cIM verif boxptop BEG
373	REAL vertlev(ncolmx,napisccp)
374	cIM verif boxptop END
375	c
376	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
377	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
378	c$OMP THREADPRIVATE(tautab,invtau)
379	REAL emsfc_lw
380	PARAMETER(emsfc_lw=0.99)
381	c REAL ran0 ! type for random number fuction
382	c
383	REAL cldtot(klon,klev)
384	c variables de haut en bas pour le simulateur ISCCP
385	REAL dtau_s(klon,klev) !tau nuages startiformes
386	REAL dtau_c(klon,klev) !tau nuages convectifs
387	REAL dem_s(klon,klev) !emissivite nuages startiformes
388	REAL dem_c(klon,klev) !emissivite nuages convectifs
389	c
390	c variables de haut en bas pour le simulateur ISCCP
391	REAL pfull(klon,klev)
392	REAL phalf(klon,klev+1)
393	REAL qv(klon,klev)
394	REAL cc(klon,klev)
395	REAL conv(klon,klev)
396	REAL dtau_sH2B(klon,klev)
397	REAL dtau_cH2B(klon,klev)
398	REAL at(klon,klev)
399	REAL dem_sH2B(klon,klev)
400	REAL dem_cH2B(klon,klev)
401	c
402	INTEGER kmax, lmax, lmax3
403	PARAMETER(kmax=8, lmax=8, lmax3=3)
404	INTEGER kmaxm1, lmaxm1
405	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
406	INTEGER iimx7, jjmx7, jjmp1x7
407	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
408	.jjmp1x7=jjmp1*lmaxm1)
409	c
410	c output from ISCCP simulator
411	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
412	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
413	REAL totalcldarea(klon,napisccp)
414	REAL meanptop(klon,napisccp)
415	REAL meantaucld(klon,napisccp)
416	REAL boxtau(klon,ncolmx,napisccp)
417	REAL boxptop(klon,ncolmx,napisccp)
418	REAL zx_tmp_fi3d_bx(klon,ncolmx)
419	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
420	c
421	REAL cld_fi3d(klon,lmax3)
422	REAL cld_3d(iim,jjmp1,lmax3)
423	c
424	INTEGER iw, iwmax
425	REAL wmin, pas_w
426	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
427	cIM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
428	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
429	REAL o500(klon)
430	c
431
432	c sorties ISCCP
433
434	integer nid_isccp
435	save nid_isccp
436	c$OMP THREADPRIVATE(nid_isccp)
437
438	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
439	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
440	SAVE zx_tau
441	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
442	SAVE zx_pc
443	c$OMP THREADPRIVATE(zx_tau,zx_pc)
444	c cldtopres pression au sommet des nuages
445	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
446	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
447	DATA cldtopres3/440., 680., 1000./
448	SAVE cldtopres,cldtopres3
449	c$OMP THREADPRIVATE(cldtopres,cldtopres3)
450	cIM 051005 BEG
451	INTEGER komega, nhoriRD
452
453	c taulev: numero du niveau de tau dans les sorties ISCCP
454	CHARACTER *4 taulev(kmaxm1)
455	c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
456	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
457	CHARACTER *3 pclev(lmaxm1)
458	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
459	SAVE taulev,pclev
460	c$OMP THREADPRIVATE(taulev,pclev)
461	c
462	c cnameisccp
463	CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
464	cIM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
465	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
466	. 'pc= 180-310hPa, tau< 0.3',
467	. 'pc= 310-440hPa, tau< 0.3',
468	. 'pc= 440-560hPa, tau< 0.3',
469	. 'pc= 560-680hPa, tau< 0.3',
470	. 'pc= 680-800hPa, tau< 0.3',
471	. 'pc= 800-1000hPa, tau< 0.3',
472	. 'pc= 50-180hPa, tau= 0.3-1.3',
473	. 'pc= 180-310hPa, tau= 0.3-1.3',
474	. 'pc= 310-440hPa, tau= 0.3-1.3',
475	. 'pc= 440-560hPa, tau= 0.3-1.3',
476	. 'pc= 560-680hPa, tau= 0.3-1.3',
477	. 'pc= 680-800hPa, tau= 0.3-1.3',
478	. 'pc= 800-1000hPa, tau= 0.3-1.3',
479	. 'pc= 50-180hPa, tau= 1.3-3.6',
480	. 'pc= 180-310hPa, tau= 1.3-3.6',
481	. 'pc= 310-440hPa, tau= 1.3-3.6',
482	. 'pc= 440-560hPa, tau= 1.3-3.6',
483	. 'pc= 560-680hPa, tau= 1.3-3.6',
484	. 'pc= 680-800hPa, tau= 1.3-3.6',
485	. 'pc= 800-1000hPa, tau= 1.3-3.6',
486	. 'pc= 50-180hPa, tau= 3.6-9.4',
487	. 'pc= 180-310hPa, tau= 3.6-9.4',
488	. 'pc= 310-440hPa, tau= 3.6-9.4',
489	. 'pc= 440-560hPa, tau= 3.6-9.4',
490	. 'pc= 560-680hPa, tau= 3.6-9.4',
491	. 'pc= 680-800hPa, tau= 3.6-9.4',
492	. 'pc= 800-1000hPa, tau= 3.6-9.4',
493	. 'pc= 50-180hPa, tau= 9.4-23',
494	. 'pc= 180-310hPa, tau= 9.4-23',
495	. 'pc= 310-440hPa, tau= 9.4-23',
496	. 'pc= 440-560hPa, tau= 9.4-23',
497	. 'pc= 560-680hPa, tau= 9.4-23',
498	. 'pc= 680-800hPa, tau= 9.4-23',
499	. 'pc= 800-1000hPa, tau= 9.4-23',
500	. 'pc= 50-180hPa, tau= 23-60',
501	. 'pc= 180-310hPa, tau= 23-60',
502	. 'pc= 310-440hPa, tau= 23-60',
503	. 'pc= 440-560hPa, tau= 23-60',
504	. 'pc= 560-680hPa, tau= 23-60',
505	. 'pc= 680-800hPa, tau= 23-60',
506	. 'pc= 800-1000hPa, tau= 23-60',
507	. 'pc= 50-180hPa, tau> 60.',
508	. 'pc= 180-310hPa, tau> 60.',
509	. 'pc= 310-440hPa, tau> 60.',
510	. 'pc= 440-560hPa, tau> 60.',
511	. 'pc= 560-680hPa, tau> 60.',
512	. 'pc= 680-800hPa, tau> 60.',
513	. 'pc= 800-1000hPa, tau> 60.'/
514	SAVE cnameisccp
515	c$OMP THREADPRIVATE(cnameisccp)
516	c
517	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
518	c INTEGER nhorix7
519	cIM: region='3d' <==> sorties en global
520	CHARACTER*3 region
521	PARAMETER(region='3d')
522	c
523	cIM ISCCP simulator v3.4
524	c
525	logical ok_hf
526	c
527	integer nid_hf, nid_hf3d
528	save ok_hf, nid_hf, nid_hf3d
529	c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
530	c QUESTION : noms de variables ?
531
532	c#ifdef histhf
533	c data ok_hf/.true./
534	c#else
535	c data ok_hf/.false./
536	c#endif
537	INTEGER longcles
538	PARAMETER ( longcles = 20 )
539	REAL clesphy0( longcles )
540	c
541	c Variables quasi-arguments
542	c
543	REAL xjour
544	SAVE xjour
545	c$OMP THREADPRIVATE(xjour)
546	c
547	c
548	c Variables propres a la physique
549	c
550	INTEGER radpas
551	SAVE radpas ! frequence d'appel rayonnement
552	c$OMP THREADPRIVATE(radpas)
553	c
554	cc INTEGER iflag_con
555	cc SAVE iflag_con ! indicateur de la convection
556	c
557	INTEGER itap
558	SAVE itap ! compteur pour la physique
559	c$OMP THREADPRIVATE(itap)
560	c
561	real slp(klon) ! sea level pressure
562	c
563	REAL fevap(klon,nbsrf)
564	REAL fluxlat(klon,nbsrf)
565	c
566	REAL qsol(klon)
567	REAL,save :: solarlong0
568	c
569	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
570	c
571	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
572	REAL zulow(klon),zvlow(klon)
573	c
574	INTEGER igwd,idx(klon),itest(klon)
575	c
576	REAL agesno(klon,nbsrf)
577	c
578	c REAL,allocatable,save :: run_off_lic_0(:)
579	cc$OMP THREADPRIVATE(run_off_lic_0)
580	cym SAVE run_off_lic_0
581	cKE43
582	c Variables liees a la convection de K. Emanuel (sb):
583	c
584	REAL bas, top ! cloud base and top levels
585	SAVE bas
586	SAVE top
587	c$OMP THREADPRIVATE(bas, top)
588
589	REAL wdn(klon), tdn(klon), qdn(klon)
590	c
591	c=================================================================================================
592	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
593	c Variables liées à la poche froide (jyg)
594
595	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
596	REAL Vprecip(klon,klev) ! precipitation vertical profile
597	c
598	REAL wape_prescr, fip_prescr
599	INTEGER it_wape_prescr
600	SAVE wape_prescr, fip_prescr, it_wape_prescr
601	c
602	c variables supplementaires de concvl
603	REAL Tconv(klon,klev)
604	REAL ment(klon,klev,klev),sij(klon,klev,klev)
605	REAL dd_t(klon,klev),dd_q(klon,klev)
606	real alp_bl_prescr
607	save alp_bl_prescr
608	real ale_bl_prescr
609	save ale_bl_prescr
610	real ale_wake(klon)
611	real alp_wake(klon)
612	cRC
613	c Variables liées à la poche froide (jyg et rr)
614	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
615
616	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
617
618	REAL wake_dth(klon,klev) ! wake : temp pot difference
619
620	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
621	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
622	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
623	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
624	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
625	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
626	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
627	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
628	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
629	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
630	REAL wake_spread(klon,klev) ! spreading term in wake_delt
631	c
632	cpourquoi y'a pas de save??
633	REAL wake_h(klon) ! Wake : hauteur de la poche froide
634	c
635	INTEGER wake_k(klon) ! Wake sommet
636	c
637	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
638	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
639	c
640	REAL wake_pe(klon) ! Wake potential energy - WAPE
641
642	REAL wake_gfl(klon) ! Gust Front Length
643	REAL wake_dens(klon)
644	c
645	c
646	REAL dt_dwn(klon,klev)
647	REAL dq_dwn(klon,klev)
648	REAL wdt_PBL(klon,klev)
649	REAL udt_PBL(klon,klev)
650	REAL wdq_PBL(klon,klev)
651	REAL udq_PBL(klon,klev)
652	REAL M_dwn(klon,klev)
653	REAL M_up(klon,klev)
654	REAL dt_a(klon,klev)
655	REAL dq_a(klon,klev)
656	c
657	cRR:fin declarations poches froides
658	c=======================================================================================================
659
660	c Variables locales pour la couche limite (al1):
661	c
662	cAl1 REAL pblh(klon) ! Hauteur de couche limite
663	cAl1 SAVE pblh
664	c34EK
665	c
666	c Variables locales:
667	c
668	REAL cdragh(klon) ! drag coefficient pour T and Q
669	REAL cdragm(klon) ! drag coefficient pour vent
670	cAA
671	cAA Pour phytrac
672	cAA
673	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
674	REAL yu1(klon) ! vents dans la premiere couche U
675	REAL yv1(klon) ! vents dans la premiere couche V
676
677	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
678
679	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
680	c@$$ PARAMETER (offline=.false.)
681	c@$$ INTEGER physid
682	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
683	REAL frac_nucl(klon,klev) ! idem (nucleation)
684	#ifdef INCA
685	INTEGER :: iii
686	REAL :: calday
687	#endif
688
689	cIM cf FH pour Tiedtke 080604
690	REAL rain_tiedtke(klon),snow_tiedtke(klon)
691	c
692	cIM 050204 END
693	REAL evap(klon), devap(klon) ! evaporation et sa derivee
694	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
695
696	REAL bils(klon) ! bilan de chaleur au sol
697	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
698	C ! type de sous-surface et pondere par la fraction
699	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
700	C ! type de sous-surface et pondere par la fraction
701	REAL fder(klon)
702	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
703	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
704	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
705	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
706	c
707	REAL frugs(klon,nbsrf)
708	REAL zxrugs(klon) ! longueur de rugosite
709	c
710	c Conditions aux limites
711	c
712	INTEGER julien
713	c
714	INTEGER lmt_pas
715	SAVE lmt_pas ! frequence de mise a jour
716	c$OMP THREADPRIVATE(lmt_pas)
717	cIM
718	REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
719
720	cym SAVE pctsrf ! sous-fraction du sol
721
722	cIM sorties
723	REAL un_jour
724	PARAMETER(un_jour=86400.)
725	c======================================================================
726	c
727	c Declaration des procedures appelees
728	c
729	EXTERNAL angle ! calculer angle zenithal du soleil
730	EXTERNAL alboc ! calculer l'albedo sur ocean
731	EXTERNAL ajsec ! ajustement sec
732	EXTERNAL conlmd ! convection (schema LMD)
733	cKE43
734	EXTERNAL conema3 ! convect4.3
735	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
736	cAA
737	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
738	c ! stockage des coefficients necessaires au
739	c ! lessivage OFF-LINE et ON-LINE
740	EXTERNAL hgardfou ! verifier les temperatures
741	EXTERNAL nuage ! calculer les proprietes radiatives
742	EXTERNAL o3cm ! initialiser l'ozone
743	EXTERNAL orbite ! calculer l'orbite terrestre
744	EXTERNAL ozonecm ! prescrire l'ozone
745	EXTERNAL phyetat0 ! lire l'etat initial de la physique
746	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
747	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
748	EXTERNAL suphel ! initialiser certaines constantes
749	EXTERNAL transp ! transport total de l'eau et de l'energie
750	EXTERNAL ecribina ! ecrire le fichier binaire global
751	EXTERNAL ecribins ! ecrire le fichier binaire global
752	EXTERNAL ecrirega ! ecrire le fichier binaire regional
753	EXTERNAL ecriregs ! ecrire le fichier binaire regional
754	cIM
755	EXTERNAL haut2bas !variables de haut en bas
756	INTEGER lnblnk1
757	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
758	!caracter
759	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
760	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
761	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
762	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
763	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
764	c
765	c Variables locales
766	c
767	REAL rhcl(klon,klev) ! humiditi relative ciel clair
768	REAL dialiq(klon,klev) ! eau liquide nuageuse
769	REAL diafra(klon,klev) ! fraction nuageuse
770	REAL cldliq(klon,klev) ! eau liquide nuageuse
771	REAL cldfra(klon,klev) ! fraction nuageuse
772	REAL cldtau(klon,klev) ! epaisseur optique
773	REAL cldemi(klon,klev) ! emissivite infrarouge
774	c
775	CXXX PB
776	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
777	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
778	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
779	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
780	c
781	REAL zxfluxt(klon, klev)
782	REAL zxfluxq(klon, klev)
783	REAL zxfluxu(klon, klev)
784	REAL zxfluxv(klon, klev)
785	CXXX
786	c
787	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
788	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
789	c Le rayonnement n'est pas calcule tous les pas, il faut donc
790	c sauvegarder les sorties du rayonnement
791	cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
792	cym SAVE sollwdownclr, toplwdown, toplwdownclr
793	cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
794	c
795	INTEGER itaprad
796	SAVE itaprad
797	c$OMP THREADPRIVATE(itaprad)
798	c
799	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
800	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
801	c
802	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
803	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
804	c
805	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
806	REAL zxsnow_dummy(klon)
807	c
808	REAL dist, rmu0(klon), fract(klon)
809	REAL zdtime, zlongi
810	c
811	CHARACTER*2 str2
812	CHARACTER*2 iqn
813	c
814	REAL qcheck
815	REAL z_avant(klon), z_apres(klon), z_factor(klon)
816	LOGICAL zx_ajustq
817	c
818	REAL za, zb
819	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
820	real zqsat(klon,klev)
821	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
822	REAL t_coup
823	PARAMETER (t_coup=234.0)
824	c
825	REAL zphi(klon,klev)
826	cym A voir plus tard !!
827	cym REAL zx_relief(iim,jjmp1)
828	cym REAL zx_aire(iim,jjmp1)
829	c
830	c Grandeurs de sorties
831	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
832	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
833	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
834	REAL s_trmb3(klon)
835	cKE43
836	c Variables locales pour la convection de K. Emanuel (sb):
837	c
838	REAL upwd(klon,klev) ! saturated updraft mass flux
839	REAL dnwd(klon,klev) ! saturated downdraft mass flux
840	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
841	REAL tvp(klon,klev) ! virtual temp of lifted parcel
842	CHARACTER*40 capemaxcels !max(CAPE)
843
844	REAL rflag(klon) ! flag fonctionnement de convect
845	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
846	c -- convect43:
847	INTEGER ntra ! nb traceurs pour convect4.3
848	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
849	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
850	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
851	REAL dplcldt(klon), dplcldr(klon)
852	c? . condm_con(klon,klev),conda_con(klon,klev),
853	c? . mr_con(klon,klev),ep_con(klon,klev)
854	c? . ,sadiab(klon,klev),wadiab(klon,klev)
855	c --
856	c34EK
857	c
858	c Variables du changement
859	c
860	c con: convection
861	c lsc: condensation a grande echelle (Large-Scale-Condensation)
862	c ajs: ajustement sec
863	c eva: evaporation de l'eau liquide nuageuse
864	c vdf: couche limite (Vertical DiFfusion)
865	REAL rneb(klon,klev)
866
867	! tendance nulles
868	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
869
870	c
871	*********************************************************
872	* declarations
873
874	*********************************************************
875	cIM 081204 END
876	c
877	REAL pmfu(klon,klev), pmfd(klon,klev)
878	REAL pen_u(klon,klev), pen_d(klon,klev)
879	REAL pde_u(klon,klev), pde_d(klon,klev)
880	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
881	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
882	REAL prfl(klon,klev+1), psfl(klon,klev+1)
883	c
884	REAL rain_lsc(klon)
885	REAL snow_lsc(klon)
886	c
887	REAL ratqss(klon,klev),ratqsc(klon,klev)
888	real ratqsbas,ratqshaut
889	save ratqsbas,ratqshaut
890	real zpt_conv(klon,klev)
891
892	c Parametres lies au nouveau schema de nuages (SB, PDF)
893	real fact_cldcon
894	real facttemps
895	logical ok_newmicro
896	save ok_newmicro
897	c$OMP THREADPRIVATE(ok_newmicro)
898	save fact_cldcon,facttemps
899	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
900	real facteur
901
902	integer iflag_cldcon
903	save iflag_cldcon
904	c$OMP THREADPRIVATE(iflag_cldcon)
905	logical ptconv(klon,klev)
906	cIM cf. AM 081204 BEG
907	logical ptconvth(klon,klev)
908	cIM cf. AM 081204 END
909	c
910	c Variables liees a l'ecriture de la bande histoire physique
911	c
912	c======================================================================
913	c
914	cIM cf. AM 081204 BEG
915	c declarations pour sortir sur une sous-region
916	integer imin_ins,imax_ins,jmin_ins,jmax_ins
917	save imin_ins,imax_ins,jmin_ins,jmax_ins
918	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
919	c real lonmin_ins,lonmax_ins,latmin_ins
920	c s ,latmax_ins
921	c data lonmin_ins,lonmax_ins,latmin_ins
922	c s ,latmax_ins/
923	c valeurs initiales s -5.,20.,41.,55./
924	c s 100.,130.,-20.,20./
925	c s -180.,180.,-90.,90./
926	c======================================================================
927	cIM cf. AM 081204 END
928
929	c
930	integer itau_w ! pas de temps ecriture = itap + itau_phy
931	c
932	c
933	c Variables locales pour effectuer les appels en serie
934	c
935	REAL zx_rh(klon,klev)
936	cIM RH a 2m (la surface)
937	REAL rh2m(klon), qsat2m(klon)
938	REAL tpot(klon), tpote(klon)
939	REAL Lheat
940
941	INTEGER length
942	PARAMETER ( length = 100 )
943	REAL tabcntr0( length )
944	c
945	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
946	cIM
947	INTEGER ndex2d1(iwmax)
948	c
949	cIM AMIP2 BEG
950	REAL moyglo, mountor
951	cIM 141004 BEG
952	REAL zustrdr(klon), zvstrdr(klon)
953	REAL zustrli(klon), zvstrli(klon)
954	REAL zustrph(klon), zvstrph(klon)
955	REAL aam, torsfc
956	cIM 141004 END
957	cIM 190504 BEG
958	INTEGER ij, imp1jmp1
959	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
960	cym A voir plus tard
961	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
962	REAL padyn(iim+1,jjmp1,klev+1)
963	REAL dudyn(iim+1,jjmp1,klev)
964	REAL rlatdyn(iim+1,jjmp1)
965	cIM 190504 END
966	LOGICAL ok_msk
967	REAL msk(klon)
968	cIM
969	REAL airetot, pi
970	cym A voir plus tard
971	cym REAL zm_wo(jjmp1, klev)
972	cIM AMIP2 END
973	c
974	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
975	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
976	c#ifdef histmthNMC
977	cym A voir plus tard !!!!
978	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
979	REAL zx_tmp_fiNC(klon,nlevSTD)
980	c#endif
981	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
982	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
983	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
984	c
985	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
986	INTEGER nid_ctesGCM
987	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
988	SAVE nid_ctesGCM
989	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
990	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
991	c
992	cIM 280405 BEG
993	INTEGER nid_bilKPins, nid_bilKPave
994	SAVE nid_bilKPins, nid_bilKPave
995	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
996	c
997	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
998	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
999	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1000	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1001	c
1002	cIM 280405 END
1003	c
1004	INTEGER nhori, nvert, nvert1, nvert3
1005	REAL zsto, zsto1, zsto2
1006	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1007	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1008	REAL zout_isccp(napisccp)
1009	SAVE zcals, zcalh, zoutj, zout_isccp
1010	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1011
1012	real zjulian
1013	save zjulian
1014	c$OMP THREADPRIVATE(zjulian)
1015
1016	character*20 modname
1017	character*80 abort_message
1018	logical ok_sync
1019	real date0
1020	integer idayref
1021
1022	C essai writephys
1023	integer fid_day, fid_mth, fid_ins
1024	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1025	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1026	parameter (prof2d_on = 1, prof3d_on = 2,
1027	. prof2d_av = 3, prof3d_av = 4)
1028	character*30 nom_fichier
1029	character*10 varname
1030	character*40 vartitle
1031	character*20 varunits
1032	C Variables liees au bilan d'energie et d'enthalpi
1033	REAL ztsol(klon)
1034	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1035	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1036	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1037	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1038	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1039	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1040	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1041	REAL d_h_vcol_phy
1042	REAL fs_bound, fq_bound
1043	SAVE d_h_vcol_phy
1044	c$OMP THREADPRIVATE(d_h_vcol_phy)
1045	REAL zero_v(klon)
1046	CHARACTER*15 ztit
1047	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1048	SAVE ip_ebil
1049	DATA ip_ebil/0/
1050	c$OMP THREADPRIVATE(ip_ebil)
1051	INTEGER if_ebil ! level for energy conserv. dignostics
1052	SAVE if_ebil
1053	c$OMP THREADPRIVATE(if_ebil)
1054	c+jld ec_conser
1055	REAL ZRCPD
1056	c-jld ec_conser
1057	REAL t2m(klon,nbsrf) ! temperature a 2m
1058	REAL q2m(klon,nbsrf) ! humidite a 2m
1059
1060	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1061	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1062	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1063	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1064	CHARACTER*40 tinst, tave, typeval
1065	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
1066	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
1067
1068	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1069
1070	REAL re(klon, klev) ! Cloud droplet effective radius
1071	REAL fl(klon, klev) ! denominator of re
1072
1073	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1074
1075	! Aerosol optical properties
1076	#ifdef INCA
1077	! Aerosol optical properties by INCA model
1078	CHARACTER*4 :: rfname(9)
1079	#endif
1080	REAL aerindex(klon) ! POLDER aerosol index
1081
1082	! Parameters
1083	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1084	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1085	cym
1086	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1087	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1088	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1089	! false : lecture des aerosol dans un fichier
1090	c$OMP THREADPRIVATE(aerosol_couple)
1091
1092	c
1093	c Declaration des constantes et des fonctions thermodynamiques
1094	c
1095	LOGICAL,SAVE :: first=.true.
1096	c$OMP THREADPRIVATE(first)
1097	#include "YOMCST.h"
1098	#include "YOETHF.h"
1099	#include "FCTTRE.h"
1100	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1101	#include "conema3.h"
1102	#include "fisrtilp.h"
1103	#include "nuage.h"
1104	#include "compbl.h"
1105	cIM 100106 END : pouvoir sortir les ctes de la physique
1106	c
1107	c======================================================================
1108	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1109	! Utilise notamment en 1D mais peut etre active egalement en 3D
1110	! en imposant la valeur de igout.
1111	c======================================================================
1112
1113	if (prt_level.ge.1) then
1114	igout=klon/2+1/klon
1115	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1116	write(lunout,*)
1117	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
1118	write(lunout,*)
1119	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
1120
1121	write(lunout,*) 'papers, play, phi, u, v, t, omega'
1122	do k=1,nlev
1123	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1124	s u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
1125	enddo
1126	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1127	do k=1,nlev
1128	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1129	enddo
1130	endif
1131
1132	c======================================================================
1133
1134	cym => necessaire pour iflag_con != 2
1135	pmfd(:,:) = 0.
1136	pen_u(:,:) = 0.
1137	pen_d(:,:) = 0.
1138	pde_d(:,:) = 0.
1139	pde_u(:,:) = 0.
1140	aam=0.
1141	torsfc=0.
1142
1143	if (first) then
1144
1145	cCR:nvelles variables convection/poches froides
1146
1147	print*, '================================================='
1148	print*, 'Allocation des variables locales et sauvegardees'
1149	call phys_local_var_init
1150	call phys_state_var_init
1151	print*, '================================================='
1152
1153	paire_ter(:)=0.
1154	clwcon(:,:)=0.
1155	rnebcon(:,:)=0.
1156	ratqs(:,:)=0.
1157	sollw(:)=0.
1158	cym Attention pbase pas initialise dans concvl !!!!
1159	pbase(:)=0
1160
1161	first=.false.
1162
1163	endif ! fisrt
1164
1165	modname = 'physiq'
1166	cIM
1167	IF (ip_ebil_phy.ge.1) THEN
1168	DO i=1,klon
1169	zero_v(i)=0.
1170	END DO
1171	END IF
1172	ok_sync=.TRUE.
1173	IF (nqmax .LT. 2) THEN
1174	abort_message = 'eaux vapeur et liquide sont indispensables'
1175	CALL abort_gcm (modname,abort_message,1)
1176	ENDIF
1177	IF (debut) THEN
1178	CALL suphel ! initialiser constantes et parametres phys.
1179	ENDIF
1180
1181	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1182
1183
1184	c======================================================================
1185	xjour = rjourvrai
1186	c
1187	c Si c'est le debut, il faut initialiser plusieurs choses
1188	c ********
1189	c
1190	IF (debut) THEN
1191	C
1192	!rv
1193	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1194	cde la convection a partir des caracteristiques du thermique
1195	wght_th(:,:)=1.
1196	lalim_conv(:)=1
1197	cRC
1198	u10m(:,:)=0.
1199	v10m(:,:)=0.
1200	piz_ae(:,:,:)=0.
1201	tau_ae(:,:,:)=0.
1202	cg_ae(:,:,:)=0.
1203	rain_con(:)=0.
1204	snow_con(:)=0.
1205	bl95_b0=0.
1206	bl95_b1=0.
1207	topswai(:)=0.
1208	topswad(:)=0.
1209	solswai(:)=0.
1210	solswad(:)=0.
1211	#ifdef INCA
1212	tau_inca(:,:,:,:) = 0.
1213	piz_inca(:,:,:,:) = 0.
1214	cg_inca(:,:,:,:) = 0.
1215	ccm(:,:,:) = 0.
1216	topswai_inca(:) = 0.
1217	topswad_inca(:) = 0.
1218	topswad0_inca(:) = 0.
1219	topsw_inca(:,:) = 0.
1220	topsw0_inca(:,:) = 0.
1221	solswai_inca(:) = 0.
1222	solswad_inca(:) = 0.
1223	solswad0_inca(:) = 0.
1224	solsw_inca(:,:) = 0.
1225	solsw0_inca(:,:) = 0.
1226	#endif
1227	!rv
1228	!ACo
1229	rnebcon0(:,:) = 0.0
1230	clwcon0(:,:) = 0.0
1231	rnebcon(:,:) = 0.0
1232	clwcon(:,:) = 0.0
1233	! fin anne
1234
1235	cIM
1236	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1237	c
1238	c appel a la lecture du run.def physique
1239	c
1240	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1241	. ok_instan, ok_hf,
1242	. solarlong0,seuil_inversion,
1243	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1244	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
1245	. ok_ade, ok_aie, aerosol_couple,
1246	. bl95_b0, bl95_b1,
1247	. iflag_thermals,nsplit_thermals,
1248	cnv flags pour la convection et les poches froides
1249	. iflag_coupl,iflag_clos,iflag_wake)
1250
1251	print*,'iflag_coupl,iflag_clos,iflag_wake',
1252	. iflag_coupl,iflag_clos,iflag_wake
1253
1254	c
1255	c
1256	c Initialiser les compteurs:
1257	c
1258	itap = 0
1259	itaprad = 0
1260
1261	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1262	!! Un petit travail à faire ici.
1263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1264
1265	if (iflag_pbl>1) then
1266	PRINT*, "Using method MELLOR&YAMADA"
1267	endif
1268
1269	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1270	. rlat,rlon,pctsrf, ftsol,
1271	. ocean, ok_veget,
1272	. falb1, falb2, rain_fall,snow_fall,
1273	. solsw, sollw,
1274	. radsol,clesphy0,
1275	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1276	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1277	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
1278
1279
1280
1281
1282	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1283
1284
1285	DO i=1,klon
1286	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1287	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1288	$ THEN
1289	WRITE(,)
1290	$ 'physiq apres lecture de restart: pb sous surface au point ',
1291	$ i, pctsrf(i, 1 : nbsrf)
1292	ENDIF
1293	ENDDO
1294
1295	radpas = NINT( 86400./dtime/nbapp_rad)
1296	c
1297	C on remet le calendrier a zero
1298	c
1299	IF (raz_date .eq. 1) THEN
1300	itau_phy = 0
1301	ENDIF
1302
1303	cIM cf. AM 081204 BEG
1304	PRINT*,'cycle_diurne3 =',cycle_diurne
1305	cIM cf. AM 081204 END
1306	c
1307	CALL printflag( tabcntr0,radpas,ok_journe,
1308	, ok_instan, ok_region )
1309	c
1310	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1311	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1312	. pdtphys
1313	abort_message='Pas physique n est pas correct '
1314	! call abort_gcm(modname,abort_message,1)
1315	dtime=pdtphys
1316	ENDIF
1317	IF (nlon .NE. klon) THEN
1318	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1319	. klon
1320	abort_message='nlon et klon ne sont pas coherents'
1321	call abort_gcm(modname,abort_message,1)
1322	ENDIF
1323	IF (nlev .NE. klev) THEN
1324	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1325	. klev
1326	abort_message='nlev et klev ne sont pas coherents'
1327	call abort_gcm(modname,abort_message,1)
1328	ENDIF
1329	c
1330	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1331	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1332	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1333	abort_message='Nbre d appels au rayonnement insuffisant'
1334	call abort_gcm(modname,abort_message,1)
1335	ENDIF
1336	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1337	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1338	. ok_cvl
1339	c
1340	cKE43
1341	c Initialisation pour la convection de K.E. (sb):
1342	IF (iflag_con.GE.3) THEN
1343
1344	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1345	WRITE(lunout,*)
1346	. "On va utiliser le melange convectif des traceurs qui"
1347	WRITE(lunout,*)"est calcule dans convect4.3"
1348	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1349
1350	DO i = 1, klon
1351	ema_cbmf(i) = 0.
1352	ema_pcb(i) = 0.
1353	ema_pct(i) = 0.
1354	ema_workcbmf(i) = 0.
1355	ENDDO
1356	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1357	DO i = 1, klon
1358	ibas_con(i) = 1
1359	itop_con(i) = 1
1360	ENDDO
1361	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1362	c===============================================================================
1363	cCR:04.12.07: initialisations poches froides
1364	c Controle de ALE et ALP pour la fermeture convective (jyg)
1365	if (iflag_wake.eq.1) then
1366	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1367	s ,alp_bl_prescr, ale_bl_prescr)
1368	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1369	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1370	endif
1371
1372	do i = 1,klon
1373	wake_s(i) = 0.
1374	wake_fip(i) = 0.
1375	wake_cstar(i) = 0.
1376	DO k=1,klev
1377	wake_deltat(i,k)=0.
1378	wake_deltaq(i,k)=0.
1379	ENDDO
1380	enddo
1381	c================================================================================
1382
1383	ENDIF
1384
1385	rugoro=0.
1386	c34EK
1387	IF (ok_orodr) THEN
1388
1389	rugoro=0.
1390
1391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1393	! justement quand ok_orodr = false.
1394	! ce rugoro est utilise par la couche limite et fait double emploi
1395	! avec les paramétrisations spécifiques de Francois Lott.
1396	! DO i=1,klon
1397	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1398	! ENDDO
1399	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1400
1401	CALL SUGWD(klon,klev,paprs,pplay)
1402	DO i=1,klon
1403	zuthe(i)=0.
1404	zvthe(i)=0.
1405	if(zstd(i).gt.10.)then
1406	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1407	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1408	endif
1409	ENDDO
1410	ENDIF
1411	c
1412	c
1413	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1414	WRITE(lunout,*)'La frequence de lecture surface est de ',
1415	. lmt_pas
1416	c
1417	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1418	c IF (ok_mensuel) THEN
1419	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1420	c . ecrit_mth
1421	c ENDIF
1422	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1423	c IF (ok_journe) THEN
1424	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1425	c . ecrit_day
1426	c ENDIF
1427	cIM 130904 BEG
1428	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1429	cIM 170305
1430	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1431	cIM 230305
1432	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1433	c
1434	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1435	c
1436	cIM200505 IF (ok_journe) THEN
1437	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1438	cIM200505 . ecrit_hf
1439	cIM200505 ENDIF
1440	cIM 130904 END
1441	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1442	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1443	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1444	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1445	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1446	cIM200505 IF (ok_instan) THEN
1447	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1448	cIM200505 . ecrit_ins
1449	cIM200505 ENDIF
1450	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1451	cIM200505 IF (ok_region) THEN
1452	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1453	cIM200505 . ecrit_reg
1454	cIM200505 ENDIF
1455	cIM 030306 BEG
1456	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1457	cIM : ne pas modifier ecrit_hf2mth
1458	c
1459	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1460	ecrit_hf2mth = ecrit_mth/ecrit_hf
1461	c ecrit_ins en secondes, chaque pas de temps de la physique
1462	ecrit_ins = dtime
1463	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1464	ecrit_hf = ecrit_hf * un_jour
1465	!IM
1466	IF(ecrit_day.LE.1.) THEN
1467	ecrit_day = ecrit_day * un_jour !en secondes
1468	ENDIF
1469	!IM
1470	ecrit_mth = ecrit_mth * un_jour
1471	ecrit_reg = ecrit_reg * un_jour
1472	ecrit_tra = ecrit_tra * un_jour
1473	ecrit_ISCCP = ecrit_ISCCP * un_jour
1474	c
1475	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1476	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1477	. ecrit_hf2mth
1478	cIM 030306 END
1479
1480	capemaxcels = 't_max(X)'
1481	t2mincels = 't_min(X)'
1482	t2maxcels = 't_max(X)'
1483	tinst = 'inst(X)'
1484	tave = 'ave(X)'
1485	cIM cf. AM 081204 BEG
1486	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1487	cIM cf. AM 081204 END
1488	c
1489	c=============================================================
1490	c Initialisation des sorties
1491	c=============================================================
1492
1493	#ifdef CPP_IOIPSL
1494
1495	c Commente par abderrahmane 11 2 08
1496	c#ifdef histhf
1497	c#include "ini_histhf.h"
1498	c#endif
1499
1500	c#ifdef histday
1501	c#include "ini_histday.h"
1502	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1503	c#include "ini_bilKP_ins.h"
1504	c#include "ini_bilKP_ave.h"
1505	c#endif
1506
1507	c#ifdef histmth
1508	c#include "ini_histmth.h"
1509	c#endif
1510
1511	c#ifdef histins
1512	c#include "ini_histins.h"
1513	c#endif
1514
1515	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
1516	& ctetaSTD,dtime,presnivs,ok_veget,
1517	& ocean,iflag_pbl,ok_mensuel,ok_journe,
1518	& ok_hf,ok_instan,nid_files)
1519
1520	#ifdef histISCCP
1521	#include "ini_histISCCP.h"
1522	#endif
1523
1524	#ifdef histmthNMC
1525	#include "ini_histmthNMC.h"
1526	#endif
1527
1528	#include "ini_histday_seri.h"
1529
1530	#include "ini_paramLMDZ_phy.h"
1531
1532	#endif
1533
1534	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1535	date0 = zjulian
1536	C date0 = day_ini
1537	WRITE(,) 'physiq date0 : ',date0
1538	c
1539	c
1540	c
1541	c Prescrire l'ozone dans l'atmosphere
1542	c
1543	c
1544	cc DO i = 1, klon
1545	cc DO k = 1, klev
1546	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1547	cc ENDDO
1548	cc ENDDO
1549	c
1550	#ifdef INCA
1551	call VTe(VTphysiq)
1552	call VTb(VTinca)
1553	iii = MOD(NINT(xjour),360)
1554	calday = FLOAT(iii) + gmtime
1555	WRITE(lunout,*) 'initial time ', xjour, calday
1556	#ifdef INCAINFO
1557	WRITE(lunout,*) 'Appel CHEMINI ...'
1558	#endif
1559	CALL chemini(
1560	$ rg,
1561	$ ra,
1562	$ airephy,
1563	$ rlat,
1564	$ rlon,
1565	$ presnivs,
1566	$ calday,
1567	$ klon,
1568	$ nqmax,
1569	$ pdtphys,
1570	$ annee_ref,
1571	$ day_ini)
1572	#ifdef INCAINFO
1573	WRITE(lunout,*) 'OK.'
1574	#endif
1575	call VTe(VTinca)
1576	call VTb(VTphysiq)
1577	#endif
1578	c
1579	ENDIF
1580	c
1581	c ************** Fin de IF ( debut ) *************
1582	c
1583
1584	! Tendances bidons pour les processus qui n'affectent pas certaines
1585	! variables.
1586	du0(:,:)=0.
1587	dv0(:,:)=0.
1588	dq0(:,:)=0.
1589	dql0(:,:)=0.
1590	c
1591	c Mettre a zero des variables de sortie (pour securite)
1592	c
1593	DO i = 1, klon
1594	d_ps(i) = 0.0
1595	ENDDO
1596	DO k = 1, klev
1597	DO i = 1, klon
1598	d_t(i,k) = 0.0
1599	d_u(i,k) = 0.0
1600	d_v(i,k) = 0.0
1601	ENDDO
1602	ENDDO
1603	DO iq = 1, nqmax
1604	DO k = 1, klev
1605	DO i = 1, klon
1606	d_qx(i,k,iq) = 0.0
1607	ENDDO
1608	ENDDO
1609	ENDDO
1610	da(:,:)=0.
1611	mp(:,:)=0.
1612	phi(:,:,:)=0.
1613	c
1614	c Ne pas affecter les valeurs entrees de u, v, h, et q
1615	c
1616	DO k = 1, klev
1617	DO i = 1, klon
1618	t_seri(i,k) = t(i,k)
1619	u_seri(i,k) = u(i,k)
1620	v_seri(i,k) = v(i,k)
1621	q_seri(i,k) = qx(i,k,ivap)
1622	ql_seri(i,k) = qx(i,k,iliq)
1623	qs_seri(i,k) = 0.
1624	ENDDO
1625	ENDDO
1626	IF (nqmax.GE.3) THEN
1627	DO iq = 3, nqmax
1628	DO k = 1, klev
1629	DO i = 1, klon
1630	tr_seri(i,k,iq-2) = qx(i,k,iq)
1631	ENDDO
1632	ENDDO
1633	ENDDO
1634	ELSE
1635	DO k = 1, klev
1636	DO i = 1, klon
1637	tr_seri(i,k,1) = 0.0
1638	ENDDO
1639	ENDDO
1640	ENDIF
1641	C
1642	DO i = 1, klon
1643	ztsol(i) = 0.
1644	ENDDO
1645	DO nsrf = 1, nbsrf
1646	DO i = 1, klon
1647	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1648	ENDDO
1649	ENDDO
1650	cIM
1651	IF (ip_ebil_phy.ge.1) THEN
1652	ztit='after dynamic'
1653	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1654	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1655	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1656	C Comme les tendances de la physique sont ajoute dans la dynamique,
1657	C on devrait avoir que la variation d'entalpie par la dynamique
1658	C est egale a la variation de la physique au pas de temps precedent.
1659	C Donc la somme de ces 2 variations devrait etre nulle.
1660	call diagphy(airephy,ztit,ip_ebil_phy
1661	e , zero_v, zero_v, zero_v, zero_v, zero_v
1662	e , zero_v, zero_v, zero_v, ztsol
1663	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1664	s , fs_bound, fq_bound )
1665	END IF
1666
1667	c Diagnostiquer la tendance dynamique
1668	c
1669	IF (ancien_ok) THEN
1670	DO k = 1, klev
1671	DO i = 1, klon
1672	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1673	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1674	ENDDO
1675	ENDDO
1676	ELSE
1677	DO k = 1, klev
1678	DO i = 1, klon
1679	d_t_dyn(i,k) = 0.0
1680	d_q_dyn(i,k) = 0.0
1681	ENDDO
1682	ENDDO
1683	ancien_ok = .TRUE.
1684	ENDIF
1685	c
1686	c Ajouter le geopotentiel du sol:
1687	c
1688	DO k = 1, klev
1689	DO i = 1, klon
1690	zphi(i,k) = pphi(i,k) + pphis(i)
1691	ENDDO
1692	ENDDO
1693	c
1694	c Verifier les temperatures
1695	c
1696	cIM BEG
1697	IF (check) THEN
1698	amn=MIN(ftsol(1,is_ter),1000.)
1699	amx=MAX(ftsol(1,is_ter),-1000.)
1700	DO i=2, klon
1701	amn=MIN(ftsol(i,is_ter),amn)
1702	amx=MAX(ftsol(i,is_ter),amx)
1703	ENDDO
1704	c
1705	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1706	ENDIF !(check) THEN
1707	cIM END
1708	c
1709	CALL hgardfou(t_seri,ftsol,'debutphy')
1710	c
1711	cIM BEG
1712	IF (check) THEN
1713	amn=MIN(ftsol(1,is_ter),1000.)
1714	amx=MAX(ftsol(1,is_ter),-1000.)
1715	DO i=2, klon
1716	amn=MIN(ftsol(i,is_ter),amn)
1717	amx=MAX(ftsol(i,is_ter),amx)
1718	ENDDO
1719	c
1720	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1721	ENDIF !(check) THEN
1722	cIM END
1723	c
1724	c Incrementer le compteur de la physique
1725	c
1726	itap = itap + 1
1727	julien = MOD(NINT(xjour),360)
1728	if (julien .eq. 0) julien = 360
1729	c
1730	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1731	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1732	c
1733	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1734	WRITE(lunout,*)' PHYS cond julien ',julien
1735	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1736	ENDIF
1737	c
1738	c Re-evaporer l'eau liquide nuageuse
1739	c
1740	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1741	DO i = 1, klon
1742	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1743	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1744	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1745	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1746	zb = MAX(0.0,ql_seri(i,k))
1747	za = - MAX(0.0,ql_seri(i,k))
1748	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1749	t_seri(i,k) = t_seri(i,k) + za
1750	q_seri(i,k) = q_seri(i,k) + zb
1751	ql_seri(i,k) = 0.0
1752	d_t_eva(i,k) = za
1753	d_q_eva(i,k) = zb
1754	ENDDO
1755	ENDDO
1756	cIM
1757	IF (ip_ebil_phy.ge.2) THEN
1758	ztit='after reevap'
1759	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1760	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1761	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1762	call diagphy(airephy,ztit,ip_ebil_phy
1763	e , zero_v, zero_v, zero_v, zero_v, zero_v
1764	e , zero_v, zero_v, zero_v, ztsol
1765	e , d_h_vcol, d_qt, d_ec
1766	s , fs_bound, fq_bound )
1767	C
1768	END IF
1769
1770	c
1771	c=========================================================================
1772	! Calculs de l'orbite.
1773	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1774	! doit donc etre placé avant radlwsw et pbl_surface
1775
1776	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1777	! solarlong0
1778
1779	if (solarlong0<-999.) then
1780	CALL orbite(FLOAT(julien),zlongi,dist)
1781	else
1782	zlongi=solarlong0 ! longitude solaire vraie
1783	dist=1. ! distance au soleil / moyenne
1784	endif
1785
1786	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
1787
1788	! Avec ou sans cycle diurne
1789	IF (cycle_diurne) THEN
1790	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1791	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1792	ELSE
1793	rmu0 = -999.999
1794	ENDIF
1795
1796	if (mydebug) then
1797	call writefield_phy('u_seri',u_seri,llm)
1798	call writefield_phy('v_seri',v_seri,llm)
1799	call writefield_phy('t_seri',t_seri,llm)
1800	call writefield_phy('q_seri',q_seri,llm)
1801	endif
1802
1803	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1804	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1805	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1806	c turbulent du couche limit.
1807	c
1808	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1809	c ecriture des fihiers hist_XXXX.nc, ces sont :
1810	c qsol, zq2m, s_pblh, s_lcl,
1811	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1812	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1813	c zxrugs, zu10m, zv10m, fder,
1814	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1815	c frugs, agesno, fsollw, fsolsw,
1816	c d_ts, fevap, fluxlat, t2m,
1817	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1818	c
1819	c Certains ne sont pas utiliser du tout :
1820	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1821	c
1822
1823	CALL pbl_surface(
1824	e dtime, date0, itap, julien,
1825	e debut, lafin,
1826	e rlon, rlat, rugoro, rmu0,
1827	e rain_fall, snow_fall, solsw, sollw,
1828	e t_seri, q_seri, u_seri, v_seri,
1829	e pplay, paprs, pctsrf,
1830	+ ftsol, falb1, falb2, u10m, v10m,
1831	s sollwdown, cdragh, cdragm, yu1, yv1,
1832	s albsol1, albsol2, sens, evap,
1833	s zxtsol, zxfluxlat, zt2m, qsat2m,
1834	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1835	s ycoefh, pctsrf_new,
1836	d qsol, zq2m, s_pblh, s_lcl,
1837	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1838	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1839	d zxrugs, zu10m, zv10m, fder,
1840	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1841	d frugs, agesno, fsollw, fsolsw,
1842	d d_ts, fevap, fluxlat, t2m,
1843	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1844	- dsens, devap, zxsnow,
1845	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1846	c
1847	c
1848	pctsrf(:,:) = pctsrf_new(:,:)
1849
1850	!-----------------------------------------------------------------------------------------
1851	! ajout des tendances de la diffusion turbulente
1852	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1853	!-----------------------------------------------------------------------------------------
1854
1855	if (mydebug) then
1856	call writefield_phy('u_seri',u_seri,llm)
1857	call writefield_phy('v_seri',v_seri,llm)
1858	call writefield_phy('t_seri',t_seri,llm)
1859	call writefield_phy('q_seri',q_seri,llm)
1860	endif
1861
1862
1863	IF (ip_ebil_phy.ge.2) THEN
1864	ztit='after surface_main'
1865	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1866	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1867	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1868	call diagphy(airephy,ztit,ip_ebil_phy
1869	e , zero_v, zero_v, zero_v, zero_v, sens
1870	e , evap , zero_v, zero_v, ztsol
1871	e , d_h_vcol, d_qt, d_ec
1872	s , fs_bound, fq_bound )
1873	END IF
1874
1875	c =================================================================== c
1876	c Calcul de Qsat
1877
1878	DO k = 1, klev
1879	DO i = 1, klon
1880	zx_t = t_seri(i,k)
1881	IF (thermcep) THEN
1882	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1883	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1884	zx_qs = MIN(0.5,zx_qs)
1885	zcor = 1./(1.-retv*zx_qs)
1886	zx_qs = zx_qs*zcor
1887	ELSE
1888	IF (zx_t.LT.t_coup) THEN
1889	zx_qs = qsats(zx_t)/pplay(i,k)
1890	ELSE
1891	zx_qs = qsatl(zx_t)/pplay(i,k)
1892	ENDIF
1893	ENDIF
1894	zqsat(i,k)=zx_qs
1895	ENDDO
1896	ENDDO
1897
1898	if (prt_level.ge.1) then
1899	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1900	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1901	endif
1902	c
1903	c Appeler la convection (au choix)
1904	c
1905	DO k = 1, klev
1906	DO i = 1, klon
1907	conv_q(i,k) = d_q_dyn(i,k)
1908	. + d_q_vdf(i,k)/dtime
1909	conv_t(i,k) = d_t_dyn(i,k)
1910	. + d_t_vdf(i,k)/dtime
1911	ENDDO
1912	ENDDO
1913	IF (check) THEN
1914	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1915	WRITE(lunout,*) "avantcon=", za
1916	ENDIF
1917	zx_ajustq = .FALSE.
1918	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1919	IF (zx_ajustq) THEN
1920	DO i = 1, klon
1921	z_avant(i) = 0.0
1922	ENDDO
1923	DO k = 1, klev
1924	DO i = 1, klon
1925	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1926	. *(paprs(i,k)-paprs(i,k+1))/RG
1927	ENDDO
1928	ENDDO
1929	ENDIF
1930	IF (iflag_con.EQ.1) THEN
1931	stop'reactiver le call conlmd dans physiq.F'
1932	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1933	c . d_t_con, d_q_con,
1934	c . rain_con, snow_con, ibas_con, itop_con)
1935	ELSE IF (iflag_con.EQ.2) THEN
1936	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1937	e conv_t, conv_q, -evap, omega,
1938	s d_t_con, d_q_con, rain_con, snow_con,
1939	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1940	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1941	WHERE (rain_con < 0.) rain_con = 0.
1942	WHERE (snow_con < 0.) snow_con = 0.
1943	DO i = 1, klon
1944	ibas_con(i) = klev+1 - kcbot(i)
1945	itop_con(i) = klev+1 - kctop(i)
1946	ENDDO
1947	ELSE IF (iflag_con.GE.3) THEN
1948	c nb of tracers for the KE convection:
1949	c MAF la partie traceurs est faite dans phytrac
1950	c on met ntra=1 pour limiter les appels mais on peut
1951	c supprimer les calculs / ftra.
1952	ntra = 1
1953
1954	c=====================================================================================
1955	cajout pour la parametrisation des poches froides:
1956	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1957	do k=1,klev
1958	do i=1,klon
1959	if (iflag_wake.eq.1) then
1960	t_wake(i,k) = t_seri(i,k)
1961	. +(1-wake_s(i))*wake_deltat(i,k)
1962	q_wake(i,k) = q_seri(i,k)
1963	. +(1-wake_s(i))*wake_deltaq(i,k)
1964	t_undi(i,k) = t_seri(i,k)
1965	. -wake_s(i)*wake_deltat(i,k)
1966	q_undi(i,k) = q_seri(i,k)
1967	. -wake_s(i)*wake_deltaq(i,k)
1968	else
1969	t_wake(i,k) = t_seri(i,k)
1970	q_wake(i,k) = q_seri(i,k)
1971	t_undi(i,k) = t_seri(i,k)
1972	q_undi(i,k) = q_seri(i,k)
1973	endif
1974	enddo
1975	enddo
1976
1977	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1978	cc-- pour le soulevement des particules dans le modele convectif
1979	c
1980	do i = 1,klon
1981	ALE(i) = 0.
1982	ALP(i) = 0.
1983	enddo
1984	c
1985	ccalcul de ale_wake et alp_wake
1986	do i = 1,klon
1987	if (iflag_wake.eq.1) then
1988	ale_wake(i) = 0.5wake_cstar(i)*2
1989	alp_wake(i) = wake_fip(i)
1990	else
1991	ale_wake(i) = 0.
1992	alp_wake(i) = 0.
1993	endif
1994	enddo
1995	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
1996	cdans le thermique sinon
1997	if (iflag_coupl.eq.0) then
1998	if (debut) print*,'ALE et ALP imposes'
1999	do i = 1,klon
2000	con ne couple que ale
2001	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2002	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2003	con ne couple que alp
2004	c ALP(i) = alp_wake(i) + Alp_bl(i)
2005	ALP(i) = alp_wake(i) + alp_bl_prescr
2006	enddo
2007	else
2008	print*,'ALE et ALP couples au thermique'
2009	do i = 1,klon
2010	ALE(i) = max(ale_wake(i),Ale_bl(i))
2011	ALP(i) = alp_wake(i) + Alp_bl(i)
2012	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2013	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2014	enddo
2015	endif
2016
2017	cfin calcul ale et alp
2018	c=================================================================================================
2019
2020
2021	c sb, oct02:
2022	c Schema de convection modularise et vectorise:
2023	c (driver commun aux versions 3 et 4)
2024	c
2025	IF (ok_cvl) THEN ! new driver for convectL
2026
2027	CALL concvl (iflag_con,iflag_clos,
2028	. dtime,paprs,pplay,t_undi,q_undi,
2029	. t_wake,q_wake,
2030	. u_seri,v_seri,tr_seri,nbtr,
2031	. ALE,ALP,
2032	. ema_work1,ema_work2,
2033	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2034	. rain_con, snow_con, ibas_con, itop_con, sigd,
2035	. upwd,dnwd,dnwd0,
2036	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2037	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2038	. pmflxr,pmflxs,da,phi,mp,
2039	. ftd,fqd,lalim_conv,wght_th)
2040
2041	cIM cf. FH
2042	clwcon0=qcondc
2043	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2044
2045	ELSE ! ok_cvl
2046	c MAF conema3 ne contient pas les traceurs
2047	CALL conema3 (dtime,
2048	. paprs,pplay,t_seri,q_seri,
2049	. u_seri,v_seri,tr_seri,ntra,
2050	. ema_work1,ema_work2,
2051	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2052	. rain_con, snow_con, ibas_con, itop_con,
2053	. upwd,dnwd,dnwd0,bas,top,
2054	. Ma,cape,tvp,rflag,
2055	. pbase
2056	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2057	. ,clwcon0)
2058
2059	ENDIF ! ok_cvl
2060
2061	c
2062	c Correction precip
2063	rain_con = rain_con * cvl_corr
2064	snow_con = snow_con * cvl_corr
2065	c
2066
2067	IF (.NOT. ok_gust) THEN
2068	do i = 1, klon
2069	wd(i)=0.0
2070	enddo
2071	ENDIF
2072
2073	c =================================================================== c
2074	c Calcul des proprietes des nuages convectifs
2075	c
2076
2077	c calcul des proprietes des nuages convectifs
2078	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2079	call clouds_gno
2080	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2081
2082	c =================================================================== c
2083
2084	DO i = 1, klon
2085	ema_pcb(i) = pbase(i)
2086	ENDDO
2087	DO i = 1, klon
2088
2089	! L'idicage de itop_con peut cacher un pb potentiel
2090	! FH sous la dictee de JYG, CR
2091	ema_pct(i) = paprs(i,itop_con(i)+1)
2092
2093	if (itop_con(i).gt.klev-3) then
2094	print*,'La convection monte trop haut '
2095	print*,'itop_con(,',i,',)=',itop_con(i)
2096	endif
2097	ENDDO
2098	DO i = 1, klon
2099	ema_cbmf(i) = ema_workcbmf(i)
2100	ENDDO
2101	ELSE IF (iflag_con.eq.0) THEN
2102	write(lunout,*) 'On n appelle pas la convection'
2103	clwcon0=0.
2104	rnebcon0=0.
2105	d_t_con=0.
2106	d_q_con=0.
2107	d_u_con=0.
2108	d_v_con=0.
2109	rain_con=0.
2110	snow_con=0.
2111	bas=1
2112	top=1
2113	ELSE
2114	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2115	CALL abort
2116	ENDIF
2117
2118	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2119	c . d_u_con, d_v_con)
2120
2121	!-----------------------------------------------------------------------------------------
2122	! ajout des tendances de la diffusion turbulente
2123	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2124	!-----------------------------------------------------------------------------------------
2125
2126	if (mydebug) then
2127	call writefield_phy('u_seri',u_seri,llm)
2128	call writefield_phy('v_seri',v_seri,llm)
2129	call writefield_phy('t_seri',t_seri,llm)
2130	call writefield_phy('q_seri',q_seri,llm)
2131	endif
2132
2133	cIM
2134	IF (ip_ebil_phy.ge.2) THEN
2135	ztit='after convect'
2136	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2137	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2138	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2139	call diagphy(airephy,ztit,ip_ebil_phy
2140	e , zero_v, zero_v, zero_v, zero_v, zero_v
2141	e , zero_v, rain_con, snow_con, ztsol
2142	e , d_h_vcol, d_qt, d_ec
2143	s , fs_bound, fq_bound )
2144	END IF
2145	C
2146	IF (check) THEN
2147	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2148	WRITE(lunout,*)"aprescon=", za
2149	zx_t = 0.0
2150	za = 0.0
2151	DO i = 1, klon
2152	za = za + airephy(i)/FLOAT(klon)
2153	zx_t = zx_t + (rain_con(i)+
2154	. snow_con(i))*airephy(i)/FLOAT(klon)
2155	ENDDO
2156	zx_t = zx_t/za*dtime
2157	WRITE(lunout,*)"Precip=", zx_t
2158	ENDIF
2159	IF (zx_ajustq) THEN
2160	DO i = 1, klon
2161	z_apres(i) = 0.0
2162	ENDDO
2163	DO k = 1, klev
2164	DO i = 1, klon
2165	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2166	. *(paprs(i,k)-paprs(i,k+1))/RG
2167	ENDDO
2168	ENDDO
2169	DO i = 1, klon
2170	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2171	. /z_apres(i)
2172	ENDDO
2173	DO k = 1, klev
2174	DO i = 1, klon
2175	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2176	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2177	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2178	ENDIF
2179	ENDDO
2180	ENDDO
2181	ENDIF
2182	zx_ajustq=.FALSE.
2183
2184	c
2185	c=============================================================================
2186	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2187	cpour la couche limite diffuse pour l instant
2188	c
2189	if (iflag_wake.eq.1) then
2190	DO k=1,klev
2191	DO i=1,klon
2192	dt_dwn(i,k) = ftd(i,k)
2193	wdt_PBL(i,k) = 0.
2194	dq_dwn(i,k) = fqd(i,k)
2195	wdq_PBL(i,k) = 0.
2196	M_dwn(i,k) = dnwd0(i,k)
2197	M_up(i,k) = upwd(i,k)
2198	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2199	udt_PBL(i,k) = 0.
2200	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2201	udq_PBL(i,k) = 0.
2202	ENDDO
2203	ENDDO
2204	c
2205	ccalcul caracteristiques de la poche froide
2206	call calWAKE (paprs,pplay,dtime
2207	: ,t_seri,q_seri,omega
2208	: ,dt_dwn,dq_dwn,M_dwn,M_up
2209	: ,dt_a,dq_a,sigd
2210	: ,wdt_PBL,wdq_PBL
2211	: ,udt_PBL,udq_PBL
2212	o ,wake_deltat,wake_deltaq,wake_dth
2213	o ,wake_h,wake_s,wake_dens
2214	o ,wake_pe,wake_fip,wake_gfl
2215	o ,dt_wake,dq_wake
2216	o ,wake_k, t_undi,q_undi
2217	o ,wake_omgbdth,wake_dp_omgb
2218	o ,wake_dtKE,wake_dqKE
2219	o ,wake_dtPBL,wake_dqPBL
2220	o ,wake_omg,wake_dp_deltomg
2221	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2222	o ,wake_ddeltat,wake_ddeltaq)
2223	c
2224	!-----------------------------------------------------------------------------------------
2225	! ajout des tendances des poches froides
2226	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2227	! coherent avec les autres d_t_...
2228	d_t_wake(:,:)=dt_wake(:,:)*dtime
2229	d_q_wake(:,:)=dq_wake(:,:)*dtime
2230	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2231	!-----------------------------------------------------------------------------------------
2232
2233	endif
2234	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2235	c
2236	c===================================================================
2237	c Convection seche (thermiques ou ajustement)
2238	c===================================================================
2239	c
2240	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2241	s ,seuil_inversion,weak_inversion,dthmin)
2242
2243
2244
2245	d_t_ajsb(:,:)=0.
2246	d_q_ajsb(:,:)=0.
2247	d_t_ajs(:,:)=0.
2248	d_u_ajs(:,:)=0.
2249	d_v_ajs(:,:)=0.
2250	d_q_ajs(:,:)=0.
2251	clwcon0th(:,:)=0.
2252	c
2253	IF(prt_level>9)WRITE(lunout,*)
2254	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2255	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2256	if(iflag_thermals.lt.0) then
2257	c Rien
2258	c ====
2259	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2260
2261
2262	else
2263
2264	c Thermiques
2265	c ==========
2266	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2267	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2268
2269
2270	if (iflag_thermals.gt.1) then
2271	call calltherm(pdtphys
2272	s ,pplay,paprs,pphi,weak_inversion
2273	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2274	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2275	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
2276	s ,ratqsdiff,zqsatth
2277	con rajoute ale et alp, et les caracteristiques de la couche alim
2278	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
2279	endif
2280
2281	! call calltherm(pdtphys
2282	! s ,pplay,paprs,pphi
2283	! s ,u_seri,v_seri,t_seri,q_seri
2284	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2285	! s ,fm_therm,entr_therm)
2286
2287	c Ajustement sec
2288	c ==============
2289
2290	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2291	! a partir du sommet des thermiques.
2292	! Dans le cas contraire, on demarre au niveau 1.
2293
2294	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2295
2296	if(iflag_thermals.eq.0) then
2297	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2298	limbas(:)=1
2299	else
2300	limbas(:)=lmax_th(:)
2301	endif
2302
2303	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2304	! pour des test de convergence numerique.
2305	! Le nouveau ajsec est a priori mieux, meme pour le cas
2306	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2307	! non nulles numeriquement pour des mailles non concernees.
2308
2309	if (iflag_thermals.eq.0) then
2310	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2311	s , d_t_ajsb, d_q_ajsb)
2312	else
2313	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2314	s , d_t_ajsb, d_q_ajsb)
2315	endif
2316
2317	!-----------------------------------------------------------------------------------------
2318	! ajout des tendances de l'ajustement sec ou des thermiques
2319	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2320	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2321	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2322
2323	!-----------------------------------------------------------------------------------------
2324
2325	endif
2326
2327	endif
2328	c
2329	c===================================================================
2330	cIM
2331	IF (ip_ebil_phy.ge.2) THEN
2332	ztit='after dry_adjust'
2333	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2334	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2335	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2336	END IF
2337
2338
2339	c-------------------------------------------------------------------------
2340	c Caclul des ratqs
2341	c-------------------------------------------------------------------------
2342
2343	c print*,'calcul des ratqs'
2344	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2345	c ----------------
2346	c on ecrase le tableau ratqsc calcule par clouds_gno
2347	if (iflag_cldcon.eq.1) then
2348	do k=1,klev
2349	do i=1,klon
2350	if(ptconv(i,k)) then
2351	ratqsc(i,k)=ratqsbas
2352	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2353	else
2354	ratqsc(i,k)=0.
2355	endif
2356	enddo
2357	enddo
2358	else if (iflag_cldcon.eq.4) then
2359	ptconvth(:,:)=.false.
2360	ratqsc(:,:)=0.
2361	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2362	call clouds_gno
2363	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2364	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2365	endif
2366
2367	c ratqs stables
2368	c -------------
2369
2370	if (iflag_ratqs.eq.0) then
2371
2372	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2373	do k=1,klev
2374	do i=1, klon
2375	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2376	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2377	enddo
2378	enddo
2379
2380	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2381	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2382	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2383	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2384	! Il s'agit de differents tests dans la phase de reglage du modele
2385	! avec thermiques.
2386
2387	else if (iflag_ratqs.eq.1) then
2388
2389	do k=1,klev
2390	do i=1, klon
2391	if (pplay(i,k).ge.60000.) then
2392	ratqss(i,k)=ratqsbas
2393	else if ((pplay(i,k).ge.30000.).and.
2394	s (pplay(i,k).lt.60000.)) then
2395	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2396	s (60000.-pplay(i,k))/(60000.-30000.)
2397	else
2398	ratqss(i,k)=ratqshaut
2399	endif
2400	enddo
2401	enddo
2402
2403	else
2404
2405	do k=1,klev
2406	do i=1, klon
2407	if (pplay(i,k).ge.60000.) then
2408	ratqss(i,k)=ratqsbas
2409	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2410	else if ((pplay(i,k).ge.30000.).and.
2411	s (pplay(i,k).lt.60000.)) then
2412	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2413	s (60000.-pplay(i,k))/(60000.-30000.)
2414	else
2415	ratqss(i,k)=ratqshaut
2416	endif
2417	enddo
2418	enddo
2419	endif
2420
2421
2422
2423
2424	c ratqs final
2425	c -----------
2426
2427	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2428	s .or.iflag_cldcon.eq.4) then
2429
2430	! On ajoute une constante au ratqsc*2 pour tenir compte de
2431	! fluctuations turbulentes de petite echelle
2432
2433	do k=1,klev
2434	do i=1,klon
2435	if ((fm_therm(i,k).gt.1.e-10)) then
2436	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2437	endif
2438	enddo
2439	enddo
2440
2441	! les ratqs sont une conbinaison de ratqss et ratqsc
2442	! 1800s, en dur pour le moment, est le temps de
2443	! relaxation des ratqs
2444
2445	facteur=exp(-pdtphys/1800.)
2446
2447	print*,'WARNING ratqs a revoir '
2448	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2449	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2450
2451	else
2452	! on ne prend que le ratqs stable pour fisrtilp
2453	ratqs(:,:)=ratqss(:,:)
2454	endif
2455
2456
2457	c
2458	c Appeler le processus de condensation a grande echelle
2459	c et le processus de precipitation
2460	c-------------------------------------------------------------------------
2461	CALL fisrtilp(dtime,paprs,pplay,
2462	. t_seri, q_seri,ptconv,ratqs,
2463	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2464	. rain_lsc, snow_lsc,
2465	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2466	. frac_impa, frac_nucl,
2467	. prfl, psfl, rhcl)
2468
2469	WHERE (rain_lsc < 0) rain_lsc = 0.
2470	WHERE (snow_lsc < 0) snow_lsc = 0.
2471	!-----------------------------------------------------------------------------------------
2472	! ajout des tendances de la diffusion turbulente
2473	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2474	!-----------------------------------------------------------------------------------------
2475	DO k = 1, klev
2476	DO i = 1, klon
2477	cldfra(i,k) = rneb(i,k)
2478	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2479	ENDDO
2480	ENDDO
2481	IF (check) THEN
2482	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2483	WRITE(lunout,*)"apresilp=", za
2484	zx_t = 0.0
2485	za = 0.0
2486	DO i = 1, klon
2487	za = za + airephy(i)/FLOAT(klon)
2488	zx_t = zx_t + (rain_lsc(i)
2489	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2490	ENDDO
2491	zx_t = zx_t/za*dtime
2492	WRITE(lunout,*)"Precip=", zx_t
2493	ENDIF
2494	cIM
2495	IF (ip_ebil_phy.ge.2) THEN
2496	ztit='after fisrt'
2497	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2498	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2499	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2500	call diagphy(airephy,ztit,ip_ebil_phy
2501	e , zero_v, zero_v, zero_v, zero_v, zero_v
2502	e , zero_v, rain_lsc, snow_lsc, ztsol
2503	e , d_h_vcol, d_qt, d_ec
2504	s , fs_bound, fq_bound )
2505	END IF
2506
2507	if (mydebug) then
2508	call writefield_phy('u_seri',u_seri,llm)
2509	call writefield_phy('v_seri',v_seri,llm)
2510	call writefield_phy('t_seri',t_seri,llm)
2511	call writefield_phy('q_seri',q_seri,llm)
2512	endif
2513
2514	c
2515	c-------------------------------------------------------------------
2516	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2517	c-------------------------------------------------------------------
2518
2519	c 1. NUAGES CONVECTIFS
2520	c
2521	cIM cf FH
2522	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2523	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2524	snow_tiedtke=0.
2525	c print*,'avant calcul de la pseudo precip '
2526	c print*,'iflag_cldcon',iflag_cldcon
2527	if (iflag_cldcon.eq.-1) then
2528	rain_tiedtke=rain_con
2529	else
2530	c print*,'calcul de la pseudo precip '
2531	rain_tiedtke=0.
2532	c print*,'calcul de la pseudo precip 0'
2533	do k=1,klev
2534	do i=1,klon
2535	if (d_q_con(i,k).lt.0.) then
2536	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2537	s *(paprs(i,k)-paprs(i,k+1))/rg
2538	endif
2539	enddo
2540	enddo
2541	endif
2542	c
2543	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2544	c
2545
2546	c Nuages diagnostiques pour Tiedtke
2547	CALL diagcld1(paprs,pplay,
2548	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2549	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2550	. diafra,dialiq)
2551	DO k = 1, klev
2552	DO i = 1, klon
2553	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2554	cldliq(i,k) = dialiq(i,k)
2555	cldfra(i,k) = diafra(i,k)
2556	ENDIF
2557	ENDDO
2558	ENDDO
2559
2560	ELSE IF (iflag_cldcon.ge.3) THEN
2561	c On prend pour les nuages convectifs le max du calcul de la
2562	c convection et du calcul du pas de temps precedent diminue d'un facteur
2563	c facttemps
2564	facteur = pdtphys *facttemps
2565	do k=1,klev
2566	do i=1,klon
2567	rnebcon(i,k)=rnebcon(i,k)*facteur
2568	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2569	s then
2570	rnebcon(i,k)=rnebcon0(i,k)
2571	clwcon(i,k)=clwcon0(i,k)
2572	endif
2573	enddo
2574	enddo
2575
2576	c
2577	cjq - introduce the aerosol direct and first indirect radiative forcings
2578	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2579	IF (ok_ade.OR.ok_aie) THEN
2580	IF ( .NOT. aerosol_couple ) THEN
2581	! Get sulfate aerosol distribution
2582	CALL readsulfate(rjourvrai, debut, sulfate)
2583	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2584
2585	! Calculate aerosol optical properties (Olivier Boucher)
2586	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2587	. tau_ae, piz_ae, cg_ae, aerindex)
2588	ENDIF
2589	ELSE
2590	tau_ae(:,:,:)=0.0
2591	piz_ae(:,:,:)=0.0
2592	cg_ae(:,:,:)=0.0
2593	ENDIF
2594
2595	c
2596	cIM calcul nuages par le simulateur ISCCP
2597	c
2598	#ifdef histISCCP
2599	IF (ok_isccp) THEN
2600	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2601	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2602	#include "calcul_simulISCCP.h"
2603	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2604	ENDIF !ok_isccp
2605	#endif
2606
2607	c On prend la somme des fractions nuageuses et des contenus en eau
2608	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2609	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2610
2611	ENDIF
2612
2613	c
2614	c 2. NUAGES STARTIFORMES
2615	c
2616	IF (ok_stratus) THEN
2617	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2618	DO k = 1, klev
2619	DO i = 1, klon
2620	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2621	cldliq(i,k) = dialiq(i,k)
2622	cldfra(i,k) = diafra(i,k)
2623	ENDIF
2624	ENDDO
2625	ENDDO
2626	ENDIF
2627	c
2628	c Precipitation totale
2629	c
2630	DO i = 1, klon
2631	rain_fall(i) = rain_con(i) + rain_lsc(i)
2632	snow_fall(i) = snow_con(i) + snow_lsc(i)
2633	ENDDO
2634	cIM
2635	IF (ip_ebil_phy.ge.2) THEN
2636	ztit="after diagcld"
2637	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2638	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2639	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2640	END IF
2641	c
2642	c Calculer l'humidite relative pour diagnostique
2643	c
2644	DO k = 1, klev
2645	DO i = 1, klon
2646	zx_t = t_seri(i,k)
2647	IF (thermcep) THEN
2648	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2649	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2650	zx_qs = MIN(0.5,zx_qs)
2651	zcor = 1./(1.-retv*zx_qs)
2652	zx_qs = zx_qs*zcor
2653	ELSE
2654	IF (zx_t.LT.t_coup) THEN
2655	zx_qs = qsats(zx_t)/pplay(i,k)
2656	ELSE
2657	zx_qs = qsatl(zx_t)/pplay(i,k)
2658	ENDIF
2659	ENDIF
2660	zx_rh(i,k) = q_seri(i,k)/zx_qs
2661	zqsat(i,k)=zx_qs
2662	ENDDO
2663	ENDDO
2664
2665	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2666	c equivalente a 2m (tpote) pour diagnostique
2667	c
2668	DO i = 1, klon
2669	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2670	IF (thermcep) THEN
2671	IF(zt2m(i).LT.RTT) then
2672	Lheat=RLSTT
2673	ELSE
2674	Lheat=RLVTT
2675	ENDIF
2676	ELSE
2677	IF (zt2m(i).LT.RTT) THEN
2678	Lheat=RLSTT
2679	ELSE
2680	Lheat=RLVTT
2681	ENDIF
2682	ENDIF
2683	tpote(i) = tpot(i)*
2684	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2685	ENDDO
2686
2687
2688	#ifdef INCA
2689	call VTe(VTphysiq)
2690	call VTb(VTinca)
2691	calday = FLOAT(julien) + gmtime
2692
2693	#ifdef INCA_AER
2694	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
2695	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2696	#endif
2697
2698	#ifdef INCAINFO
2699	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2700	#endif
2701
2702	zxsnow_dummy(:) = 0.0
2703
2704	CALL chemhook_begin (calday,
2705	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
2706	$ julien,
2707	$ gmtime,
2708	#endif
2709	$ pctsrf(1,1),
2710	$ rlat,
2711	$ rlon,
2712	$ airephy,
2713	$ paprs,
2714	$ pplay,
2715	$ ycoefh,
2716	$ pphi,
2717	$ t_seri,
2718	$ u,
2719	$ v,
2720	$ wo,
2721	$ q_seri,
2722	$ zxtsol,
2723	$ zxsnow_dummy,
2724	$ solsw,
2725	$ albsol1,
2726	$ rain_fall,
2727	$ snow_fall,
2728	$ itop_con,
2729	$ ibas_con,
2730	$ cldfra,
2731	$ iim,
2732	$ jjm,
2733	#ifdef INCA_AER
2734	$ tr_seri,
2735	$ ftsol,
2736	$ paprs,
2737	$ cdragh,
2738	$ cdragm,
2739	$ pctsrf,
2740	$ pdtphys,
2741	$ itap)
2742	#else
2743	$ tr_seri)
2744	#endif
2745
2746
2747	#ifdef INCAINFO
2748	WRITE(lunout,*)'OK.'
2749	#endif
2750	call VTe(VTinca)
2751	call VTb(VTphysiq)
2752	#endif
2753	c
2754	c Calculer les parametres optiques des nuages et quelques
2755	c parametres pour diagnostiques:
2756	c
2757	IF (aerosol_couple ) THEN
2758	#ifdef INCA
2759	sulfate(:,:) = ccm(:,:,1)
2760	sulfate_pi(:,:) = ccm(:,:,2)
2761	#endif
2762	ENDIF
2763
2764	if (ok_newmicro) then
2765	CALL newmicro (paprs, pplay,ok_newmicro,
2766	. t_seri, cldliq, cldfra, cldtau, cldemi,
2767	. cldh, cldl, cldm, cldt, cldq,
2768	. flwp, fiwp, flwc, fiwc,
2769	e ok_aie,
2770	e sulfate, sulfate_pi,
2771	e bl95_b0, bl95_b1,
2772	s cldtaupi, re, fl)
2773	else
2774	CALL nuage (paprs, pplay,
2775	. t_seri, cldliq, cldfra, cldtau, cldemi,
2776	. cldh, cldl, cldm, cldt, cldq,
2777	e ok_aie,
2778	e sulfate, sulfate_pi,
2779	e bl95_b0, bl95_b1,
2780	s cldtaupi, re, fl)
2781
2782	endif
2783	c
2784	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2785	c
2786	IF (MOD(itaprad,radpas).EQ.0) THEN
2787
2788	DO i = 1, klon
2789	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2790	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2791	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2792	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2793	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2794	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2795	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2796	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2797	ENDDO
2798
2799	if (mydebug) then
2800	call writefield_phy('u_seri',u_seri,llm)
2801	call writefield_phy('v_seri',v_seri,llm)
2802	call writefield_phy('t_seri',t_seri,llm)
2803	call writefield_phy('q_seri',q_seri,llm)
2804	endif
2805
2806	IF (aerosol_couple) THEN
2807	#ifdef INCA_AER
2808	CALL radlwsw_inca
2809	e (kdlon,kflev,dist, rmu0, fract, solaire,
2810	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2811	e wo,
2812	e cldfra, cldemi, cldtau,
2813	s heat,heat0,cool,cool0,radsol,albpla,
2814	s topsw,toplw,solsw,sollw,
2815	s sollwdown,
2816	s topsw0,toplw0,solsw0,sollw0,
2817	s lwdn0, lwdn, lwup0, lwup,
2818	s swdn0, swdn, swup0, swup,
2819	e ok_ade, ok_aie,
2820	e tau_inca, piz_inca, cg_inca,
2821	s topswad_inca, solswad_inca,
2822	s topswad0_inca, solswad0_inca,
2823	s topsw_inca, topsw0_inca,
2824	s solsw_inca, solsw0_inca,
2825	e cldtaupi,
2826	s topswai_inca, solswai_inca)
2827	#endif
2828	ELSE
2829	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2830	e (dist, rmu0, fract,
2831	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2832	e wo,
2833	e cldfra, cldemi, cldtau,
2834	s heat,heat0,cool,cool0,radsol,albpla,
2835	s topsw,toplw,solsw,sollw,
2836	s sollwdown,
2837	s topsw0,toplw0,solsw0,sollw0,
2838	s lwdn0, lwdn, lwup0, lwup,
2839	s swdn0, swdn, swup0, swup,
2840	e ok_ade, ok_aie, ! new for aerosol radiative effects
2841	e tau_ae, piz_ae, cg_ae, ! ="=
2842	s topswad, solswad, ! ="=
2843	e cldtaupi, ! ="=
2844	s topswai, solswai) ! ="=
2845	ENDIF
2846
2847	itaprad = 0
2848	ENDIF
2849	itaprad = itaprad + 1
2850
2851	if (iflag_radia.eq.0) then
2852	print *,'--------------------------------------------------'
2853	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2854	print *,'>>>> heat et cool mis a zero '
2855	print *,'--------------------------------------------------'
2856	heat=0.
2857	cool=0.
2858	endif
2859
2860
2861	c
2862	c Ajouter la tendance des rayonnements (tous les pas)
2863	c
2864	DO k = 1, klev
2865	DO i = 1, klon
2866	t_seri(i,k) = t_seri(i,k)
2867	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2868	ENDDO
2869	ENDDO
2870	c
2871	if (mydebug) then
2872	call writefield_phy('u_seri',u_seri,llm)
2873	call writefield_phy('v_seri',v_seri,llm)
2874	call writefield_phy('t_seri',t_seri,llm)
2875	call writefield_phy('q_seri',q_seri,llm)
2876	endif
2877
2878	cIM
2879	IF (ip_ebil_phy.ge.2) THEN
2880	ztit='after rad'
2881	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2882	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2883	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2884	call diagphy(airephy,ztit,ip_ebil_phy
2885	e , topsw, toplw, solsw, sollw, zero_v
2886	e , zero_v, zero_v, zero_v, ztsol
2887	e , d_h_vcol, d_qt, d_ec
2888	s , fs_bound, fq_bound )
2889	END IF
2890	c
2891	c
2892	c Calculer l'hydrologie de la surface
2893	c
2894	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2895	c . agesno, ftsol,fqsurf,fsnow, ruis)
2896	c
2897
2898	c
2899	c Calculer le bilan du sol et la derive de temperature (couplage)
2900	c
2901	DO i = 1, klon
2902	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2903	c a la demande de JLD
2904	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2905	ENDDO
2906	c
2907	cmoddeblott(jan95)
2908	c Appeler le programme de parametrisation de l'orographie
2909	c a l'echelle sous-maille:
2910	c
2911	IF (ok_orodr) THEN
2912	c
2913	c selection des points pour lesquels le shema est actif:
2914	igwd=0
2915	DO i=1,klon
2916	itest(i)=0
2917	c IF ((zstd(i).gt.10.0)) THEN
2918	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2919	itest(i)=1
2920	igwd=igwd+1
2921	idx(igwd)=i
2922	ENDIF
2923	ENDDO
2924	c igwdim=MAX(1,igwd)
2925	c
2926	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2927	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2928	e igwd,idx,itest,
2929	e t_seri, u_seri, v_seri,
2930	cIM 141004 s zulow, zvlow, zustr, zvstr,
2931	s zulow, zvlow, zustrdr, zvstrdr,
2932	s d_t_oro, d_u_oro, d_v_oro)
2933	c
2934	c ajout des tendances
2935	!-----------------------------------------------------------------------------------------
2936	! ajout des tendances de la trainee de l'orographie
2937	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
2938	!-----------------------------------------------------------------------------------------
2939	c
2940	ENDIF ! fin de test sur ok_orodr
2941	c
2942	if (mydebug) then
2943	call writefield_phy('u_seri',u_seri,llm)
2944	call writefield_phy('v_seri',v_seri,llm)
2945	call writefield_phy('t_seri',t_seri,llm)
2946	call writefield_phy('q_seri',q_seri,llm)
2947	endif
2948
2949	IF (ok_orolf) THEN
2950	c
2951	c selection des points pour lesquels le shema est actif:
2952	igwd=0
2953	DO i=1,klon
2954	itest(i)=0
2955	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2956	itest(i)=1
2957	igwd=igwd+1
2958	idx(igwd)=i
2959	ENDIF
2960	ENDDO
2961	c igwdim=MAX(1,igwd)
2962	c
2963	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2964	e rlat,zmea,zstd,zpic,
2965	e itest,
2966	e t_seri, u_seri, v_seri,
2967	s zulow, zvlow, zustrli, zvstrli,
2968	s d_t_lif, d_u_lif, d_v_lif)
2969	c
2970	!-----------------------------------------------------------------------------------------
2971	! ajout des tendances de la portance de l'orographie
2972	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
2973	!-----------------------------------------------------------------------------------------
2974	c
2975	ENDIF ! fin de test sur ok_orolf
2976	c
2977	cIM cf. FLott BEG
2978	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
2979
2980	if (mydebug) then
2981	call writefield_phy('u_seri',u_seri,llm)
2982	call writefield_phy('v_seri',v_seri,llm)
2983	call writefield_phy('t_seri',t_seri,llm)
2984	call writefield_phy('q_seri',q_seri,llm)
2985	endif
2986
2987	DO i = 1, klon
2988	zustrph(i)=0.
2989	zvstrph(i)=0.
2990	ENDDO
2991	DO k = 1, klev
2992	DO i = 1, klon
2993	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
2994	c (paprs(i,k)-paprs(i,k+1))/rg
2995	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
2996	c (paprs(i,k)-paprs(i,k+1))/rg
2997	ENDDO
2998	ENDDO
2999	c
3000	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3001	c
3002	IF (is_sequential) THEN
3003
3004	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3005	C ra,rg,romega,
3006	C rlat,rlon,pphis,
3007	C zustrdr,zustrli,zustrph,
3008	C zvstrdr,zvstrli,zvstrph,
3009	C paprs,u,v,
3010	C aam, torsfc)
3011	ENDIF
3012	cIM cf. FLott END
3013	cIM
3014	IF (ip_ebil_phy.ge.2) THEN
3015	ztit='after orography'
3016	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3017	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3018	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3019	END IF
3020	c
3021	c
3022	cAA
3023	cAA Installation de l'interface online-offline pour traceurs
3024	cAA
3025	c====================================================================
3026	c Calcul des tendances traceurs
3027	c====================================================================
3028	C
3029	call phytrac ( rnpb,
3030	I itap,
3031	I julien,
3032	I gmtime,
3033	I debut,
3034	I lafin,
3035	I nqmax-2,
3036	I nlon,
3037	I nlev,
3038	I dtime,
3039	I u,
3040	I v,
3041	I t,
3042	I paprs,
3043	I pplay,
3044	I pmfu,
3045	I pmfd,
3046	I pen_u,
3047	I pde_u,
3048	I pen_d,
3049	I pde_d,
3050	I ycoefh,
3051	I fm_therm,
3052	I entr_therm,
3053	I yu1,
3054	I yv1,
3055	I ftsol,
3056	I pctsrf,
3057	I rlat,
3058	I frac_impa,
3059	I frac_nucl,
3060	I rlon,
3061	I presnivs,
3062	I pphis,
3063	I pphi,
3064	I albsol1,
3065	I qx(1,1,1),
3066	I rhcl,
3067	I cldfra,
3068	I rneb,
3069	I diafra,
3070	I cldliq,
3071	I itop_con,
3072	I ibas_con,
3073	I pmflxr,
3074	I pmflxs,
3075	I prfl,
3076	I psfl,
3077	I da,
3078	I phi,
3079	I mp,
3080	I upwd,
3081	I dnwd,
3082	I aerosol_couple,
3083	#ifdef INCA
3084	I flxmass_w,
3085	I tau_inca,
3086	I piz_inca,
3087	I cg_inca,
3088	I ccm,
3089	I rfname,
3090	#endif
3091	O tr_seri)
3092
3093	IF (offline) THEN
3094
3095	print*,'Attention on met a 0 les thermiques pour phystoke'
3096	call phystokenc (
3097	I nlon,nlev,pdtphys,rlon,rlat,
3098	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3099	I fm_therm,entr_therm,
3100	I ycoefh,yu1,yv1,ftsol,pctsrf,
3101	I frac_impa, frac_nucl,
3102	I pphis,airephy,dtime,itap)
3103
3104
3105	ENDIF
3106
3107	c
3108	c Calculer le transport de l'eau et de l'energie (diagnostique)
3109	c
3110	CALL transp (paprs,zxtsol,
3111	e t_seri, q_seri, u_seri, v_seri, zphi,
3112	s ve, vq, ue, uq)
3113	c
3114	cIM global posePB BEG
3115	IF(1.EQ.0) THEN
3116	c
3117	CALL transp_lay (paprs,zxtsol,
3118	e t_seri, q_seri, u_seri, v_seri, zphi,
3119	s ve_lay, vq_lay, ue_lay, uq_lay)
3120	c
3121	ENDIF !(1.EQ.0) THEN
3122	cIM global posePB END
3123	c Accumuler les variables a stocker dans les fichiers histoire:
3124	c
3125	c+jld ec_conser
3126	DO k = 1, klev
3127	DO i = 1, klon
3128	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3129	d_t_ec(i,k)=0.5/ZRCPD
3130	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3131	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3132	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3133	END DO
3134	END DO
3135	c-jld ec_conser
3136	cIM
3137	IF (ip_ebil_phy.ge.1) THEN
3138	ztit='after physic'
3139	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3140	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3141	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3142	C Comme les tendances de la physique sont ajoute dans la dynamique,
3143	C on devrait avoir que la variation d'entalpie par la dynamique
3144	C est egale a la variation de la physique au pas de temps precedent.
3145	C Donc la somme de ces 2 variations devrait etre nulle.
3146	call diagphy(airephy,ztit,ip_ebil_phy
3147	e , topsw, toplw, solsw, sollw, sens
3148	e , evap, rain_fall, snow_fall, ztsol
3149	e , d_h_vcol, d_qt, d_ec
3150	s , fs_bound, fq_bound )
3151	C
3152	d_h_vcol_phy=d_h_vcol
3153	C
3154	END IF
3155	C
3156	c=======================================================================
3157	c SORTIES
3158	c=======================================================================
3159
3160	cIM Interpolation sur les niveaux de pression du NMC
3161	c -------------------------------------------------
3162	c
3163	#include "calcul_STDlev.h"
3164	c
3165	c slp sea level pressure
3166	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3167	c
3168	ccc prw = eau precipitable
3169	DO i = 1, klon
3170	prw(i) = 0.
3171	DO k = 1, klev
3172	prw(i) = prw(i) +
3173	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3174	ENDDO
3175	ENDDO
3176	c
3177	cIM initialisation + calculs divers diag AMIP2
3178	c
3179	#include "calcul_divers.h"
3180	c
3181	#ifdef INCA
3182	call VTe(VTphysiq)
3183	call VTb(VTinca)
3184	#ifdef INCAINFO
3185	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3186	#endif
3187	CALL chemhook_end (calday,
3188	$ dtime,
3189	$ pplay,
3190	$ t_seri,
3191	$ tr_seri,
3192	$ nbtr,
3193	$ paprs,
3194	$ q_seri,
3195	$ annee_ref,
3196	$ day_ini,
3197	$ airephy,
3198	#ifdef INCA_AER
3199	$ xjour,
3200	$ pphi,
3201	$ pphis,
3202	$ zx_rh)
3203	#else
3204	$ xjour)
3205	#endif
3206	#ifdef INCAINFO
3207	WRITE(lunout,*)'OK.'
3208	#endif
3209	call VTe(VTinca)
3210	call VTb(VTphysiq)
3211	#endif
3212
3213	c=============================================================
3214	c
3215	c Convertir les incrementations en tendances
3216	c
3217	if (mydebug) then
3218	call writefield_phy('u_seri',u_seri,llm)
3219	call writefield_phy('v_seri',v_seri,llm)
3220	call writefield_phy('t_seri',t_seri,llm)
3221	call writefield_phy('q_seri',q_seri,llm)
3222	endif
3223
3224	DO k = 1, klev
3225	DO i = 1, klon
3226	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3227	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3228	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3229	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3230	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3231	ENDDO
3232	ENDDO
3233	c
3234	IF (nqmax.GE.3) THEN
3235	DO iq = 3, nqmax
3236	DO k = 1, klev
3237	DO i = 1, klon
3238	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3239	ENDDO
3240	ENDDO
3241	ENDDO
3242	ENDIF
3243	c
3244	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3245	cIM global posePB#include "write_bilKP_ins.h"
3246	cIM global posePB#include "write_bilKP_ave.h"
3247	c
3248	c Sauvegarder les valeurs de t et q a la fin de la physique:
3249	c
3250	DO k = 1, klev
3251	DO i = 1, klon
3252	t_ancien(i,k) = t_seri(i,k)
3253	q_ancien(i,k) = q_seri(i,k)
3254	ENDDO
3255	ENDDO
3256	c
3257	!==========================================================================
3258	! Sorties des tendances pour un point particulier
3259	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3260	! pour le debug
3261	! La valeur de igout est attribuee plus haut dans le programme
3262	!==========================================================================
3263
3264	if (prt_level.ge.1) then
3265	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3266	write(lunout,*)
3267	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
3268	write(lunout,*)
3269	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
3270	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3271	s pctsrf(igout,is_sic)
3272	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3273	do k=1,nlev
3274	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3275	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3276	s d_t_eva(igout,k)
3277	enddo
3278	write(lunout,*) 'cool,heat'
3279	do k=1,nlev
3280	write(lunout,*) cool(igout,k),heat(igout,k)
3281	enddo
3282
3283	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3284	do k=1,nlev
3285	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3286	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3287	enddo
3288
3289	write(lunout,*) 'd_ps ',d_ps(igout)
3290	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3291	do k=1,nlev
3292	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3293	s d_qx(igout,k,1),d_qx(igout,k,2)
3294	enddo
3295	endif
3296
3297	!==========================================================================
3298
3299	c============================================================
3300	c Calcul de la temperature potentielle
3301	c============================================================
3302	DO k = 1, klev
3303	DO i = 1, klon
3304	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3305	ENDDO
3306	ENDDO
3307	c
3308
3309	c 22.03.04 BEG
3310	c=============================================================
3311	c Ecriture des sorties
3312	c=============================================================
3313	#ifdef CPP_IOIPSL
3314
3315	c Recupere des varibles calcule dans differents modules
3316	c pour ecriture dans histxxx.nc
3317
3318	! Get some variables from module fonte_neige_mod
3319	CALL fonte_neige_get_vars(pctsrf,
3320	. zxfqcalving, zxfqfonte, zxffonte)
3321
3322	IF (ocean == 'slab') THEN
3323	! Get some variables from module ocean_slab_mod
3324	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
3325	ELSEIF (ocean == 'couple') THEN
3326	! Get some variables from module ocean_cpl_mod
3327	CALL ocean_cpl_get_vars(fluxo, fluxg)
3328	ELSE
3329	! Get some variables from module ocean_forced_mod
3330	CALL ocean_forced_get_vars(fluxo, fluxg)
3331	ENDIF
3332
3333
3334	c Commente par abderrahmane le 11 2 08
3335	c#ifdef histhf
3336	c#include "write_histhf.h"
3337	c#endif
3338
3339	c#ifdef histday
3340	c#include "write_histday.h"
3341	c#endif
3342
3343	c#ifdef histmth
3344	c#include "write_histmth.h"
3345	c#endif
3346
3347	c#ifdef histins
3348	c#include "write_histins.h"
3349	c#endif
3350
3351	#include "phys_output_write.h"
3352
3353	#ifdef histISCCP
3354	#include "write_histISCCP.h"
3355	#endif
3356
3357	#ifdef histmthNMC
3358	#include "write_histmthNMC.h"
3359	#endif
3360
3361	#include "write_histday_seri.h"
3362
3363	#include "write_paramLMDZ_phy.h"
3364
3365	#endif
3366
3367	c 22.03.04 END
3368	c
3369	c====================================================================
3370	c Si c'est la fin, il faut conserver l'etat de redemarrage
3371	c====================================================================
3372	c
3373
3374
3375	IF (lafin) THEN
3376	itau_phy = itau_phy + itap
3377	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
3378	. rlat, rlon, pctsrf, ftsol,
3379	. falb1, falb2, rain_fall,
3380	. snow_fall,
3381	. solsw, sollw,
3382	. radsol,
3383	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3384	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
3385	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
3386	open(97,form="unformatted",file="finbin")
3387	write(97) u_seri,v_seri,t_seri,q_seri
3388	close(97)
3389	ENDIF
3390
3391
3392	RETURN
3393	END
3394	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3395	IMPLICIT none
3396	c
3397	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3398	c la conservation de l'eau
3399	c
3400	#include "YOMCST.h"
3401	INTEGER klon,klev
3402	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3403	REAL aire(klon)
3404	REAL qtotal, zx, qcheck
3405	INTEGER i, k
3406	c
3407	zx = 0.0
3408	DO i = 1, klon
3409	zx = zx + aire(i)
3410	ENDDO
3411	qtotal = 0.0
3412	DO k = 1, klev
3413	DO i = 1, klon
3414	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3415	. *(paprs(i,k)-paprs(i,k+1))/RG
3416	ENDDO
3417	ENDDO
3418	c
3419	qcheck = qtotal/zx
3420	c
3421	RETURN
3422	END
3423	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3424	IMPLICIT none
3425	c
3426	c Tranformer une variable de la grille physique a
3427	c la grille d'ecriture
3428	c
3429	INTEGER nfield,nlon,iim,jjmp1, jjm
3430	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3431	c
3432	INTEGER i, n, ig
3433	c
3434	jjm = jjmp1 - 1
3435	DO n = 1, nfield
3436	DO i=1,iim
3437	ecrit(i,n) = fi(1,n)
3438	ecrit(i+jjm*iim,n) = fi(nlon,n)
3439	ENDDO
3440	DO ig = 1, nlon - 2
3441	ecrit(iim+ig,n) = fi(1+ig,n)
3442	ENDDO
3443	ENDDO
3444	RETURN
3445	END

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

Download in other formats:

Original Format