Context Navigation

physiq.F @ 1086

Last change on this file since 1086 was 1068, checked in by Laurent Fairhead, 16 years ago
Correction sans cycle diurne FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.8 KB

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

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

Download in other formats:

Original Format