Context Navigation

physiq.F @ 998

Last change on this file since 998 was 998, checked in by Laurent Fairhead, 16 years ago
Modifications necessaires a la preparation au passage au nouveau rayonnement RRTM MPL LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.7 KB

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

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

Download in other formats:

Original Format