Context Navigation

physiq.F @ 1149

Last change on this file since 1149 was 1149, checked in by Laurent Fairhead, 15 years ago
Erreur de langage LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 114.9 KB

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

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

Download in other formats:

Original Format